Abstract
This study uses the continuum-based discrete element method (CDEM) to build a new model for warhead safety distance, simulating the fragment fields over space and time for tungsten alloy and steel fragment focused warheads to derive respective safety distances for personnel. Regarding the effects of different warhead attitudes and fragment effectiveness definitions, it is found that, for single-shot focused warheads, the fragmentation risk is significantly higher for horizontal placement than for vertical placement; the fragment safety distance inneed increases nonlinearly as the threshold energy decreases. Meanwhile, the safety distance is less affected by the threshold energy for high-density-casing charges than for low density.
1. Introduction
The warhead is the key component determining the terminal effect of missiles. Most missile warheads currently use fragmentation as their main method [1]. Any accidental explosions during production or storage pose a great threat to lives and property, so safety zones must be set up around its production, storage, and transport areas. At the same time, as more houses and infrastructures are built in the vicinity of military enterprises with economic development, land shortage ensures that the safety zones must be moderate in size. Thus, research on the safety distance from fragments in the case of accidental explosion has great significance for the safety of military production.
The fragmentation safety distance for safety applications is based on hazard-free, as opposed to the fragmentation radius in military applications, based on destruction. It is generally believed that the safety distance can be demarcated into two parts: where the fragments are “effective” (definitely lethal), i.e., danger to people, and where the fragment field has a certain density, or a certain probability of harm to personnel [2, 3].
Fragments’ effectiveness is generally described by their kinetic energy or specific kinetic energy [1,4–6], with different threshold values based on the data, as shown in Table 1.
For density, the US Department of Defense’s (DOD’s) “1968 Ammunition and Explosives Safety Standards” were the first to propose that the density of fragments at the safety distance must be greater than 1000 ft2/fragment (56 m2/fragment). This indicator has been widely recognized for considering factors such as the average population density near the military facility and the public’s risk tolerance [7].
Other studies on fragment safety distance have been carried out on this basis. Empirical formulas were combined with experimental data at an early stage. For example, Oliver [8] regarded all fragments (with kinetic energy greater than 0 J) as effective to obtain the safety distance when detonating single-shell explosives as follows:where D is the safety distance (m) and is the total weight of the destroyed ammunition (kg), i.e., shell + explosive.
The 1992 edition of the DOD’s Ammunition and Explosives Safety Standards used the total mass of explosives in the warhead, taking 78 J as the effectiveness cutoff for kinetic energy to produce the D-Q criteria of safety distance:where D is the safety distance (m), Q is the total mass of the explosives (kg), and K is a risk factor determined empirically. Safety distances based on different requirements can be expressed directly using K values, i.e., K6, K7, and K8.
Similar criteria include Jarret [9], Ward [10] and China’s standard of 1990 [11] (effective kinetic energy cutoff of 98 J), etc. The accuracy of such methods is limited, only serving to qualitatively show how different damage criteria influence the distance.
As technology developed, warhead simulations began to appear. The commonly used finite element analysis generates a grid of the calculation area [12, 13], with high precision in the small-scale stage of warhead casing breakage and fragment acceleration [14]. For research on safety distance, however, the fragment dimensions may be millimeter-scale, and the displacement may reach hundreds of meters, so grid generation and finite element analysis become extremely difficult [15, 16].
Therefore, current fragment safety simulations mainly use iteration. Taking the ray tracing method developed by the Dutch TNO-PML [17] as an example, by inputting the velocity vector and resistance of each fragment, resistance and gravity are calculated iteratively to obtain the trajectory and safety distance of the fragment group. Moore et al. [18] added static explosion test data to a ray trace program to obtain the safe dynamic fragment distances from a horizontal warhead for both personnel and airborne aircraft. Wang et al. [19] estimated the fragment field of vertically placed 155 mm high explosives, using both the 78 J kinetic energy and 160 J/cm2 specific kinetic energy criteria to calculate their maximum kill radius, finding that the former was about 13% larger than the latter. Jiang et al. [20–22] used a similar method, taking 78 J–98 J as the kinetic energy criteria to calculate the fragment fields of various warheads, finding that the angle of the projectile at the time of detonation also affects the fragmentation radius.
Test conditions and fragment effectiveness criteria are the bases for research on the warhead fragmentation safety. Even testing the same warhead, changes in the two can significantly affect the definition of the safe fragmentation distance, thereby affecting data comparability and industry standards. Nevertheless, existing research in this area is not systematic.
The traditional finite element methodology has mature application in close-range warhead fragmentation fields, but it cannot be directly applied to calculate long-range fragmentation fields. Gridless iteration can describe long-range fragmentation fields, but its initial motion conditions (fragment velocity, direction, rotation, etc.) must be inputted through theoretical equations or static explosion testing. The accuracy of the former is poor, while the latter is expensive and has a long testing cycle, which restricts the accuracy and applicability of the iterative method.
This study uses the continuum-based discrete element (CDEM) method, which includes finite element and discreet element modules (calculated iteratively). The calculation results of the former are inherited by the latter as initial conditions, and a numerical method is used to completely describe the whole process from warhead casing breakup to fragments landing, improving the accuracy and efficiency of safety distance calculations.
Focused warheads concentrate the detonation energy and have higher fragment density and initial velocity than other warheads of the same weight [23], causing higher fragmentation risk. Thus, a certain type of focused warhead is selected to study the impact of projectile positioning and definition of fragment effectiveness on safety distance.
2. CDEM Method
Finite element and discrete particle element methods are combined to simulate the prefabricated fragment dispersal field with high efficiency. In the calculations, the explosive and casing are described using Lagrangian finite elements, structured as shown in Figure 1.
At the start of the calculation, the explosive, casing, and equivalent prefabricated fragment units are activated to perform the finite element calculation; the explosion is ignited at a given position in the explosive, driving the casing and equivalent prefabricated fragment layer to expand. Meanwhile, an equivalent detonation product escape algorithm is introduced to obtain accurate gas pressure. After the velocity of the prefabricated fragment equivalent layer reaches the constant stage, finite element calculation is stopped, all finite elements are inactivated, discrete element calculation is activated, and the velocity of the prefabricated elements is inherited from corresponding finite element calculations. After that, the flying process with drag force of the prefabricated fragment is simulated, and the dispersal field and target shooting is analyzed. The calculation process is shown in Figure 2.
2.1. Finite Element Portion
2.1.1. Fundamental Principles
An incremental method is used to calculate the stress and node deformation force of the finite elements [24]:where Bi, Δεi, Δσi, , and are, respectively, the strain matrix, incremental strain vector, incremental stress vector, integral coefficient, and Jacobian determinant at Gauss point i; and are the current and previous period stress vectors at Gauss point i; D, Δue, and Fe are, respectively, the element’s elasticity matrix, node displacement vector, and node deformation force vector; and N is the number of Gauss points.
Under external load, the finite element unit undergoes large translation and rotation, simulated here by updating the strain matrix (B) in real time.
After calculating the deformation force of the node, we calculate its resultant force:where F is the node’s resultant force, FE is its external force, Fe is its force contributed by deformation of the finite element, Fc is its force contributed through the contact interface, and Fd is the node’s damping force.
The node’s motion is calculated using Euler’s forward interpolation method:where a is the node’s acceleration, is its velocity, is its incremental displacement, is its total displacement, m is its mass, and is the time step. An explicit solution process can be realized based on alternate computation of equations (4) and (5).
2.1.2. Explosive Model
This study uses the Landau ignition model to describe the adiabatic constant detonation and gas detonation processes. The model uses the Landau-StanNewkovic formula (rate equation) (equation (4)) to calculate the gas expansion pressure, where and indicate the adiabatic indices in the first and second stages, respectively, (for general condensed explosives, and ); and , respectively, indicate the instantaneous pressure and volume of the high-pressure explosion frontier; and , respectively, indicate the pressure of the frontier in the initial period and the volume of the explosives; and and , respectively, indicate the pressure and volume at the high-pressure frontier at the boundary of the two-stage adiabatic process. and are given by equations (5) and (6), respectively, where is explosion heat (), is charge density (), and is detonation velocity ():
2.2. Discrete Element Computation
When the finite element calculation concludes, the prefabricated fragments inherit their velocity from the corresponding finite blocks. The fragments decelerate in the air under the combined action of gravity and drag. The dynamics calculations in equation (9), based on Newton’s laws, is to calculate the velocity and displacement of the fragments over time.
The resultant force on the fragments during flight iswhere F is the resultant force, G is gravity, and Fc is the drag on the fragments during flight.
The drag Fc is calculated as follows:where is air density, A is the fragments’ equivalent windward cross-section area, and is the coefficient of drag. In order to more accurately reflect the impact of drag on fragment flying velocity and trajectory, a velocity-dependent drag coefficient is used to solve for the drag Fc. Based on previous research results [25–27], for specific shape fragment, a correspondence relation exists between coefficient of drag and Mach number, as shown in Figure 3.
These curves are embedded in the calculation program, and the coefficient of drag is calculated using the resultant fragment speed at a given moment, finally solving for drag.
2.3. Experimental Verification
Wang and Guo [28] conducted static explosion tests on a certain focused fragment warhead with length 300 mm and outer diameter 230 mm, with 4 g heavy spherical tungsten alloy fragments placed inside. The charge was a high-energy explosive mainly composed of HMX. The detonation method was end cover center detonation. The warhead was placed upright at a height of 1 m, and a velocity measurement target and witness plate were placed 10 m from the ammunition to collect fragment velocity and density distribution. The test site lay out and witness plate after detonation are shown in Figure 4.
(a)
(b)
The CDEM described above was used to simulate the same warhead. The number density of fragments along a 360° annular vertical target plate at 10 m was obtained using 350 mm-wide frame in accord with the experiment, as shown in Figure 5. Comparison of the simulation results with the test results is shown in Table 2.
The initial velocity and number density of fragments 10 m from the explosion center obtained from simulation fit well with the experimental results, with deviation less than 10%, verifying the validity and accuracy of CDEM to calculate the fragment field.
3. Numerical Model
3.1. Modeling
For convenience, the above warhead is used here, building a finite element model as shown in Figure 6.
(a)
(b)
In Figure 6, the explosive charge is modeled with density 1780 kg/m3, detonation velocity 8070 m/s, and detonation pressure 32.6 × 109 Pa.
The aluminum inner casing, cap, and equivalent fragment layer are described using an elastic-plastic model with strain softening.
For comparative study, two types of fragments are loaded into the above frame: tungsten alloy balls, as above, and steel balls of equal volume. Each fragment uses a linear elastic model with the material parameters shown in Table 3:
The explosive, inner casing, end covers, and equivalent fragment layer are all modeled using a Lagrangian grid with interval 5 mm.
In the discrete element calculations, the coefficient of drag varies dynamically with speed, and the air density is 1.069 g/L (1000 m above sea level).
3.2. Analysis of Errors Accumulation in Calculation
Accumulation of errors takes place at each time step during simulation, which is of great significance to the accuracy and confidences of calculation result. Thus, a lot of studies have been made on this question, and the current research holds that the error itself mainly depends on accuracy of algorithm and grid, and on the number of time steps. A quantitatively estimate method considered the factors above has been proposed by Smirnov et al. [29, 30]:where is the relative error of integration in one-dimensional case. and are the sizes of one cell and domain, respectively. k is the order of accuracy of numerical scheme (k = 1 for CDEM). In the presence of several directions in integration, the errors are calculated as follows:
Considering the influence of time step, the accumulated error is calculated as follows:where n is the number of time steps in Navier–Stokes equations integration. Smax is defined as the allowable value of total error in simulation; then, the maximal allowable number of time steps for solving a problem for each code could be determined by the following formula:
The reliability index of results could be defined as follows:
Table 4 summarizes accumulation of error for CDEM for different grid resolutions and time steps. The allowable error was assumed 5%.
Accumulation of errors is negatively correlated with grid resolution and time step. Thus, Lagrangian grid with interval 5 mm (resolution: 60 × 46 × 46) and the time step of 1e − 8 s is selected to conduct calculation, considering both accuracy and calculation cost.
4. Analysis of Simulation Results
4.1. Effect of Warhead Attitude on Safety Distance
The current standard contains no requirements for the attitude of the warhead when calculating safety distance. The previous tests used both horizontal and vertical positioning, with different results.
Taking the tungsten fragments warhead as an example, a 2 m tall vertical annular plate is used, defining fragment effectiveness using the 78 J kinetic energy cutoff to calculate the safety distance for both vertically and horizontally placed warheads.
Figure 7 shows that the vertical tungsten ball warhead forms a distinct focal zone in the center, with significantly higher fragment speed, consistent with the design of the warhead. 150 μs after the detonation, the fragments have essentially finished accelerating. At this point, their maximum velocity is 2286 m/s. Projectile acceleration for a horizontally placed warhead is similar.
Then, the finite element module is inactivated and discrete element computation is activated. The results are shown in Figures 8 and 9.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
The fragments’ deceleration flight is calculated during discrete element computation. Due to the focus effect of warhead, most of the fragments are driven and concentrated into an extremely small angle during acceleration. The fragments group showed a narrow-annular shaped distribution in ballistic motion, as is shown in Figures 8 and 9.
There are significant differences in near-ground fragment field between vertical and horizontal placed warheads. The former one covers all extreme directions, and the latter one only affects a narrow angle. In order to illustrate the fragment risks for both attitudes, assess method from MIL-STD-2105D [4] is applied as follows.
Defining the y axis as vertical and the x and z axes as the horizontal planes, the warhead center is projected on the origin point of the horizontal plane, positioned along the x axis. Since the projection of the fragment trajectory on the financial plane is a ray, the polar angle θ of the flight can be represented as follows:where is the horizontal distance from the fragment to the origin at a given time, and is the z axis coordinate value of the fragment at that time, i.e.,
The fragmentation at a 10 m target plate is shown in Figure 10.
(a)
(b)
Within the angle of danger, the number of effective fragments is calculated using virtual target plates of different distance to finally calculate their number density. The result is shown in Figure 11, intersection between the effective fragments curve and the 56 m2/fragment threshold is the fragment safety distance of the warhead.
(a)
(b)
As shown in Figure 11, for this type of warhead, the number of fragments in both warhead attitudes decreases with distance. The required safety distance for the vertical position is only 15.8% of that in the horizontal position, for two reasons.
First, by definition, the scope of density used to calculate safety distance is that within the angle of danger (the angle with 90% of the target fragments). Due to the focusing effect shown in the figure, most of the fragments from the horizontal warhead are concentrated in the narrow direction of danger, whereas the angle of danger for the vertical warhead is 360°, implying a much larger area than with the horizontal position, correspondingly reducing fragment density and safety distance.
Second is kinetic energy. Fragments’ ballistic trace against a specific target is affected by their initial velocity and vertical projection angle, which are not evenly distributed along the projectile (especially for focused warheads). Thus, in vertical placement, the combination of fragments’ vertical projection angle and initial velocity is not “optimal”; in horizontal placement, since the projectile angle becomes radially distributed along the projectile body, it covers the whole area evenly, optimizing fragment trajectory to the greatest extent, increasing the kinetic energy at long distances.
Thus, calculations of safety distances for single warheads show that horizontal placement is more dangerous than vertical; research on the former is more meaningful.
4.2. Impact of Effectiveness Definitions on Safety Distances
The current standard uses different definitions of fragment effectiveness. The above two cases of focused warheads with tungsten and steel fragments are taken as continued examples. The warhead is placed horizontally 1 m from the ground, and the fragments are counted at 2 m high annular vertical target plates. The critical density is 56 m2 per effective fragment. The finite element calculation results for each warhead are shown in Table 5.
Fragment kinetic energy levels of 21 J (light injury), 39 J (French standard for lethality), 78 J (US standard for lethality), 98 J (Chinese standard for lethality) and specific kinetic energy levels of 78.6 J/cm2 (light injury) and 160 J/cm2 (lethality) are used to define fragment effectiveness and calculate safety distance.
Figure 12 shows several safety distances derived using different definitions of fragment effectiveness.
(a)
(b)
At a relatively close distance, the number density of effective fragments from the steel ball warhead is larger than that from the tungsten ball warhead of the same structure, because the material density of steel is much less than that of the tungsten alloy. The former is driven to a higher initial velocity by the detonation over a narrower angle of danger. Subsequently, the numeric density of the tungsten alloy fragments overtakes that of steel, due to their higher mass and material density and better speed maintenance during flight. The final safety distances appear after the intersection of the two curves. Qualitatively, speed retention contributes more to focused warheads’ safety distance than initial fragment speed.
As distance increases, at a certain point, the density of effective steel fragments drops off rapidly. Due to the focusing effect of the warhead, most steel balls have similar initial velocity, projection angle, and deceleration. A large number of target fragments fall below the threshold velocity at almost the same time, dramatically changing the density of effective fragments at a certain distance. For the tungsten balls, however, the density of effective fragments drops more steadily. They retain speed better, and a large number of target fragments remain above the threshold. The main reasons for the reduction in density are increasing statistical area and fragments landing. This result shows that the factors controlling safety distance are different for the two warheads. The former is kinetic energy, while the latter is fragment quantity.
As distance increases, the number of effective target fragments in the two warheads essentially decreases monotonously, demonstrating that, in this case, direct-shot fragments (where the highest ballistic point does not exceed the target height) are the main danger. As shown in Figure 4, some long-distance targets get more effective fragments than close ones: some fragments bypass the upper edge of the front target to hit the rear target, indicating a risk of overtopping. It is worth noting that as the energy threshold decreases, the portion of effective fragments increases.
Figure 13 shows the changes in safety distance with definition of kinetic energy and specific kinetic energy.
(a)
(b)
As the kinetic energy definition decreases, the safety distance increases continuously.
Longitudinally comparing the effects of different kinetic energy criteria, defining fragment effectiveness at 21 J (slight injury), 39 J (French lethality standard), 78 J (US lethality standard), and 98 J (Chinese lethality standard), the safety distance varies relatively little for tungsten balls: as the kinetic energy definition declines 78.6%, the distance only increases 2.15%. The influence on steel fragments is larger: the safety distance varies 30%.
Transversely comparing the differences for kinetic energy and specific kinetic energy, both safety distances defined using the lethality standard for specific kinetic energy (160 J/cm2) are close to the results obtained using the American standard (78 J), with less than 1% deviation. When judging by slight injury, the safety distances of tungsten fragments between kinetic and specific kinetic energy are close, with deviation of 2.2%. The safety distance of steel fragments for slight injury is however greatly affected by the definition of effectiveness, with the kinetic energy distance 15.5% than that for specific kinetic energy, consistent with the conclusions of Wang et al. [19].
Regardless of the criteria used, the required safety distance increases nonlinearly as the lethality standard is reduced, at a decreasing rate. This is because it is affected not only by the lethality criteria, but also just as much by the density of fragments. As the distance increases, the area of calculation likewise increases, and density becomes more salient as a controlling factor.
5. Conclusions
(1)This paper has applied the CDEM method to simulating warheads, building a cross-scale model to calculate safety distance in the whole space-time domain. Comparison with test results verifies the validity and accuracy of this model.(2)For single-shot focused fragment warheads, fragmentation risk from detonation is much higher in horizontal than in vertical placement.(3)Regardless of the criteria adopted, as the kinetic energy definition of fragment effectiveness decreases, the safety distance increases nonlinearly at a decreasing rate.(4)The safety distance for tungsten alloy fragments is less affected by the definition of effectiveness than for steel fragments.(5)For this type of warhead, the safety distance obtained using the lethal specific kinetic energy definition is close to that using kinetic energy; it is significantly farther than using the corresponding kinetic energy criteria. Further research should carefully select fragment effectiveness criteria.
Data Availability
The data used to support the findings of this study are available upon request to the corresponding author.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The research was funded by Major Special Projects of the National Defense Science and Technology Bureau. The associated grant number is not applicable due to confidentiality reasons.
Supplementary Materials
Several calculation results of safety distance with definition of kinetic energy and specific kinetic energy are attached to provide all the data in Figure 12. (Supplementary Materials)