#### Abstract

An efficient empirical approach was developed to accurately represent the blast shock wave loading resulting from the launch of a missile from a military aircraft to be used in numerical analyses. Based on experimental test series of missile launches in laboratory environment and from a helicopter, equations were derived to predict the time- and position-dependent overpressure. The method was finally applied and validated in a structural analysis of a helicopter tail boom under missile launch shock wave loading.

#### 1. Introduction

The launch of a space rocket or military missile generates a shock wave behind the nozzle with high pressure loading on the surrounding structures. If such a missile is launched from a flying military aircraft like a fighter jet or helicopter it has to be assured that the structural loading mainly of the rear part of the aircraft does not exceed design limits. In the past, common means of compliance were mainly based on experimental live fire tests with adequate data acquisition equipment on the aircraft structure being expensive in both time and cost.

Aircraft designers seek for more efficient approaches to analyse the structural loading under missile launch blast shock wave. Most attention was recently paid to computational fluid dynamics (CFD) methods to simulate the missile launch procedure and missile exhaust plume [1]. Such CFD analyses can be extended to a subsequent structural loading in a fluid-structure-interaction approach. For example, Kidd and Cloke [2] and Okamoto et al. [3] used CFD modelling to simulate the blast wave development due to the launch of a missile from a helicopter. Also the studies by Lee et al. [4, 5] and Fu et al. [6] present CFD models to simulate the missile launch from a canister in general or from a helicopter with the focus on the missile movement. Simulations of the exhaust plume from launcher rockets using CFD modelling were performed by Gusman et al. [7], Alexeenko et al. [8], Ebrahimi et al. [9], and Vitkin et al. [10]. The complex flow fields generated by missile propulsion systems are challenging for CFD simulations with Navier-Stokes flow solvers though.

Another approach to model the missile exhaust plume, based on molecular gas dynamics, is the direct simulation Monte Carlo (DSMC) method, which was used for the modelling of the plume flow field of a space orbiter by Rault [11] and of a satellite thruster by Park et al. [12], Lee and Choi [13], and Chung et al. [14]. The coupling of this method with a CFD analysis is also a common approach that is often used in space sciences, which can be found, for example, in the numerical studies of rocket plumes of spacecraft thrusters by Feng et al. [15], Gatsonis et al. [16], and Lumpkin et al. [17]. Markelov et al. [18] compute the thermal and mechanical loads from a rocket plume acting on a satellite using this combined approach. Also Giordano et al. [19] used this method to optimise a satellite design based on the effects of the thruster plumes. Morris et al. [20] used this approach for the analysis of a lunar lander on a dusty surface. However, these methods require numerous details and parameters of the specific missile, its solid/liquid propellant, and chemically reacting kinetics, which are typically beyond the knowledge of the structural engineer. The computational efforts for such detailed analyses are immense, too.

Therefore, it is desired to handle the whole structural analysis in the framework of a conventional explicit finite element (FE) simulation with a more efficient method for the complex blast shock wave loading of the aircraft. First of all, this paper presents an assessment of potential methods for this load case in commercial explicit FE software with an evaluation of the representation of the real physical phenomena. A new empirical approach is then derived from an experimental test campaign and finally applied to a helicopter structure for verification and validation. Merits and limits of this new method are discussed and guidelines for the structural engineer are given on how to use it for other missile launch load cases. It has to be pointed out that due to the sensitivity and confidentiality of this topic no quantitative values will be published in this paper.

#### 2. Experimental Work

In order to develop a numerical method for the representation of the shock wave during a missile launch it is necessary to have detailed data on the time- and position-dependent pressure distribution behind a missile during an actual launch. Therefore, an experimental campaign of two different test series was conducted.

The first test series was performed under laboratory conditions with the missile being fixed to a solid test rig and being launched under this condition. Four pressure sensors were mounted behind the missile’s nozzle with different lateral and axial distances to record the shock wave pressure loading (Figure 1(a)). The illustration in Figure 1(b) depicts the qualitative results of the different sensors in terms of overpressure versus time. The curves show a very sudden overpressure increase up to a peak value, followed by a decay even down to negative values and a final convergence to atmospheric pressure. This general behaviour matches exactly the typical overpressure versus time curves of explosive blast shock waves, also referred to as the Friedlander curve [21], suggesting that the theory of explosive blast can also be applied here. Taking a closer look at the results, the pressure peak of sensors 1 and 2 occurs earlier compared to sensors 3 and 4, because they are located closer to the missile nozzle. However, it will be emphasised that the peak value of sensors 1 and 2 is higher than the value of sensors 3 and 4. This violates the general blast theory, which expects the pressure peak to reduce with increasing distance from the detonation point [21]. In contrast to a singular explosive detonation, the launch of a missile is a continuous, more complex process, resulting from the combustion of the propellant and the resulting hypersonic jet behind the missile’s nozzle. The physical flow phenomena occurring during a missile launch are well described and verified by CFD analyses in [2]. After ignition of the missile the compressed air is pushed out of the nozzle and generates a pressure wave with spherical expansion with the speed of sound. In addition, a hypersonic jet is created by the combustion with a spatial orientation in backward direction. This jet interacts with the aforementioned pressure wave in a complex manner and leads to higher overpressure values in the vicinity of the jet flow axis, explaining the higher peak value for sensors 3 and 4 in Figure 1(b).

**(a)**

**(b)**

The second test series was performed by launching the missile directly from a helicopter, both on ground and in flight (Figure 2(a)). For this purpose, six pressure sensors were mounted onto the surface of the helicopter tail boom at various relevant positions. In addition, four strain gauges were used to measure the structural deformation and to compare it with design limits. An Acra Control D120 data acquisition system was used to link all the measurement equipment. A total of eight test repetitions were performed in order to obtain reliable test results. The overpressure versus time plots in Figure 2(b) show similar trends as in the previous test series. First of all, the general Friedlander-like pressure curve can be observed for each position. The order of peak time corresponds to the distance of the respective sensors to the missile nozzle. The maximum peak value, again, strongly depends on the distance to the jet flow axis (-axis) though and not only on the overall distance to the nozzle.

**(a)**

**(b)**

#### 3. Assessment of Numerical Methods

Keeping in mind the target of deriving an efficient numerical approach for the simulation of this shock wave loading in commercial explicit FE software, a short assessment of potential numerical strategies will be given.

One promising option, based on conventional blast theory and the Friedlander pressure curve, is the ConWep method, which is available in current releases of LS-Dyna [22] and Abaqus/Explicit [23]. ConWep, initially developed by the US Army to calculate the blast effects of conventional weapons, is an empirical approach to calculate blast overpressures and their temporal and spatial expansion, just based on the input of an equivalent mass of TNT and the coordinates of the detonation point. It was applied to air blast simulations of composite panels in [24, 25], of composite vessels in [26], and of honeycomb sandwich panels in [27]. Although this method is very efficient, it is only suitable for singular explosions and does not represent in a realistic way the observed pressure loading under missile launch with the higher overpressure close to the jet flow axis.

Another popular option to model blast shock wave loading with realistic fluid-structure interaction is the coupled Eulerian-Lagrangian approach (abbreviated as CEL or ALE in different FE software). With this technique, a fixed Eulerian mesh is generated that is initially filled with air, normally being modelled as an ideal gas. The shock wave loading results from the detonation modelling of an explosive charge, typically using a Jones-Wilkins-Lee equation of state [28]. The loaded structure is modelled as a Lagrangian body within this Eulerian mesh, coupled by a contact formulation [29]. This approach was used, for example, for blast simulations of composite plates in [30] and for simulations of explosive blast loads inside metallic aircraft fuselages with LS-Dyna [31, 32] or Dytran [33] and inside composite fuselages in [34]. Further applications to blast analyses of armoured ground vehicle are documented in [35] efforts, thus making it an inefficient technique. However, in order to achieve accurate results with this method, the Eulerian mesh needs to be as small as possible, typically leading to millions of elements and immense computational efforts, thus making it inefficient technique. Furthermore, it is again only suitable for singular or multiple explosions, but not for modelling the combustion jet flow during missile launch.

Besides this complex fluid-structure-interaction method, the simplest approach to achieve a time-dependent shock wave pressure loading is to apply the pressure based on tabular input as a function of time, either with a simple triangular pulse or by approximating the Friedlander curve. This method was applied, for example, in [36], for explosive blast simulations of stringer-stiffened aluminium plates and in a similar manner for composite plates in [37]. Comparable blast simulations of foam core sandwich [38], honeycomb sandwich [39], and corrugated aluminium sandwich structures [40] have also been reported. A direct application for helicopter structures is reported by McCarthy [41], who used this approach to model the gun blast pressure wave acting on the aircraft structure. However, it can be stated that a uniform loading of the whole aircraft structure would not be realistic and dividing the structure into several subregions with different pressure-time functions to approximate the test results makes the approach inefficient and still not accurate enough.

The method of choice should be not only time-dependent but also position-dependent, as described in [37, 42]. By analytical equations that were derived from physical considerations and experimental test data, this approach should calculate and apply the accurate pressure at each time step for each element of the structural model. The basis for this approach is available in Abaqus/Explicit as a user-defined subroutine VDLOAD, generating nonuniform distributed loads as a function of time and the position coordinates , , and . Within this subroutine the FE solver provides for each time step and each integration point the input values of the current time (value totalTime), the position coordinates (value curCoords), and the orientation of the element (value dirCos). Based on this information, the user is requested to calculate the current pressure, which is then returned to the solver as the output value of the subroutine.

#### 4. Development of Empirical Approach

The major challenge is the development of accurate and physically meaningful analytical equations to calculate this current overpressure value as a function of time and position coordinates. As a starting point, it can be observed from the experimental results in Figures 1(b) and 2(b) that the shock wave loading at each spatial position can be approximated by a Friedlander curve (see Figure 3), which is expressed asIn these equations is the atmospheric pressure, is the peak overpressure, is the arrival time of the shock front, is the duration of the positive pressure phase after the peak overpressure, and is a wave form parameter describing the pressure decay.

The three fundamental values that are needed for each spatial position are basically , , and . Taking into account the aforementioned experimental test results, these values can be derived for all pressure sensor positions. The empirical approach is then based on the interpolation between these known values in order to obtain the whole time- and position-dependent overpressure. Since the test results of the missile launch from the helicopter appear to be more realistic compared to the laboratory test launch with a fixed missile, those data from one representative test were used for the following study. Another important factor that needs to be taken into account is the angular orientation of the pressure sensors to the shock front during the test. Since the pressure sensors were directly mounted onto the helicopter tail boom surface, they were not oriented normal to the shock front. Under normal orientation they would have recorded a higher maximum value, which can be compared to the well-known example of putting a hand out of the window of a driving car, where the pressure on the hand (=sensor) depends on the angle to the driving direction. Therefore, the maximum value had to be recalculated using the angle between the sensor normal direction (taken from CAD model data of the helicopter structure) and the straight line between sensor and missile nozzle.

The arrival time of the shock wave was then extracted from the test results and plotted against the distance of the sensor position to the missile nozzle, which is defined asTo interpolate between those experimental values in Figure 4(a), a trend line or analytical function was determined. From a simplified point of view, the arrival time is in linear relation to the distance just based on the speed of sound for the wave propagation. Hence, a linear approach was adopted with the parameter :The experimentally recorded peak overpressure values were not plotted as a function of the spatial distance but of a weighted distance , which is defined asThis approach, based on a function similar to a Mach cone, allows capturing the experimentally observed effect of higher shock wave pressure loads in the vicinity of the jet flow axis (-axis). For interpolation a power function was used with the parameters and (Figure 4(b)):It has to be mentioned that this approach is a simplification and only valid for near field pressure prediction. The effect of decreasing peak pressures in large distances on -axis is not covered in this equation.

**(a)**

**(b)**

**(c)**

The duration of the positive pressure phase was finally extracted for all sensors from the experiment and plotted against the distance (Figure 4(c)). An exponential function was adopted for interpolation, using the two parameters and :With these empirical equations derived from the test data the final overpressure can be calculated for each time step and for every position. For this purpose, (3), (4), (5), (6), and (7) are integrated into (2), leading to the combined equationThe parameter values , , , , and were derived according to the trend lines shown in Figure 4 but are not presented here due to confidentiality. The wave form parameter was set to 1.2 as taken from literature [21].

As a final step, the orientation of the loaded elements with respect to the shock front needs to be taken into account again, since the current value represents the maximum pressure if the element is oriented perpendicularly to the shock front. For this purpose, the final pressure value is recalculated using the angle between the element’s normal (which is provided by the solver) and the connecting line of element and missile nozzle. Furthermore, the pressure is only applied to an element, if this angle is below 90°. In case of angles larger than 90°, the element would be in the shadow and, for reasons of simplification, would not be loaded.

#### 5. Application to Helicopter Tail Boom Analysis

Finally, this approach was adopted to simulate the missile launch shock wave pressure loading on a helicopter structure. The helicopter model is similar to the one of the live fire test campaign and limited to the tail boom, which was fully clamped at the nodes of the forward edge. The structure consists mainly of lightweight composite and sandwich structures, which were modelled using layered shell elements of type S4 with the Hashin composite damage model and an average size of 15 mm, leading to a total number of 137.000 elements for the whole model. Equation (8) was implemented into the user subroutine VDLOAD in Abaqus/Explicit to efficiently simulate the distributed load, which was applied to the whole exterior surface of the helicopter tail boom. Based on this equation, the subroutine interpolates to apply the correct pressure load at each time step and every position of the model. The average time step size was 0.0002 ms.

The only available values to validate the simulation results are the pressure recordings and strain gauge data. The numerically predicted pressure curves in Figure 5 are in good agreement with the test data and enable a realistic loading of the aircraft structure. Hence, the strain values, illustrated in Figures 6(b) and 6(c) for two exemplary strain gauges, are also in good agreement and allow for a valuable prediction of the maximum strains to be expected under this loading. Implementation of structural damping in the simulation is recommended in order to reduce the oscillations from the sudden shock loading of the structure. An illustration of the maximum deflection of the tail boom during the shock wave loading is presented in Figure 7, using a deformation scale factor of 2. Although the deformation appears to be significant, the structural loading is uncritical as the strains are far below the design limit.

**(a) Test**

**(b) Simulation**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

#### 6. Conclusions

The development of an empirical approach to represent the blast shock wave loading resulting from the launch of a missile from an aircraft based on an experimental test campaign was presented. Similarities to the blast shock wave of a conventional singular detonation of an explosive charge were highlighted in terms of the possibility to apply the Friedlander curve to approximate the pressure-time response. However, large differences to such a conventional explosion also occur resulting from the physical flow phenomena of the continuous jet flow from the missile nozzle, making approaches like ConWep or coupled Eulerian-Lagrangian simulations inappropriate. Hence, an empirical equation was developed to represent the overpressure at each time and spatial position in a very efficient way. Although a test series of a pressure sensor-instrumented missile launch is required to determine the necessary parameters for the equation, this loading model can afterwards be applied by the aircraft engineer to various aircraft types or configurations. An enhancement of the equation or different interpolation functions could be assessed if a larger test database compared to this study was available, using more pressure sensors at further positions.

#### Nomenclature

: | Constant in equation for peak overpressure, Pa |

: | Constant in equation for positive pressure duration, s |

: | Speed of sound, m·s^{−1} |

: | Distance of sensor to nozzle, m |

: | Weighted distance |

: | Constant in equation for positive pressure duration, s·m^{−1} |

: | Constant in equation for peak overpressure |

: | Pressure, Pa |

: | Atmospheric pressure, Pa |

: | Peak overpressure, Pa |

: | Time, s |

: | Duration of positive pressure phase, s |

: | Arrival time of shock front, s |

: | Position coordinates, m |

: | Wave form parameter. |

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.