Real-time simulation of a fragmenting explosion for cylindrical warheads

Abstract

Explosion models based on finite element analysis (FEA) can be used to simulate how a warhead fragments. However, their execution times are extensive. Active protection systems need to make very fast predictions, before a fast attacking weapon hits the target. Fast execution times are also needed in real-time simulations where the impact of many different computer models is being assessed. Hence, FEA explosion models are not appropriate for these real-time systems. As a trade-off between accuracy and execution time, this paper creates a simulation of fragments from a warhead’s explosion, using simple analytical equations. The results are verified against explosion experimental data and FEA results. The developed model then can be made available for real-time simulation and fast computation.

1. Introduction

A modern active protection system such as the Defence Aid System (DAS) relies on a computer system with a quick reaction time. An active protection system identifies a threat and can quickly assemble a response before the threat strikes, thereby reducing the potential damage. The computer system incorporated into such a scheme must be able to assess multiple simultaneous threats, within fractions of a second.

Fast computer systems are also required to gage the effect of various threats on a vehicle under design, and allow a designer to determine the optimum armor and/or vehicle architecture required to increase a vehicle’s survivability. ¹ Faster models allow more threats to be assessed early in the design cycle, which allows designers to consider a larger range of design options. ²

One of the main threats to a military asset is an explosion, a complex physical event involving many variables. Many of these variables were first successfully analyzed during World War II by people such as Gurney, ³ Taylor, ⁴ and Mott. ⁵

Since World War II, with the progress of computers, more complicated models have been developed using methods such as finite element analysis (FEA) and particle hydrodynamics. These methods provide accurate results, but they are computationally intensive and are therefore unsuitable for models that require a fast computer execution time.

As a trade-off between accuracy and execution time, this paper creates a simulation of fragments created by a warhead’s explosion using simple analytical equations. These equations include improvements to the equations developed by Gurney ³ and Taylor ⁴ and modifications to Mott’s distribution ⁵ of the number and weight of fragments.

Mott’s distribution ⁵ is an essential part of the simulation; it is very accurate^6,7 and the weight of fragments is vital in determining their ability to penetrate a target. Several improvements to this distribution have been provided by researchers such as Elek and Jaramaz, ⁸ Victor, ⁹ and Hokanson et al. ¹⁰ However, these results only deliver small advances, so the simulation includes changes to the original Mott distribution.

The accuracy of the new model is validated against experimental data and simulations that use FEA models. The model’s overall execution speed is compared with the execution speed of FEA models.

The simulation consists of linked steps with each step dependent on a number of parameters and variables. As the steps are linked, the code for any step may be replaced if and when it can be improved.

The steps taken to set-up the simulation are as follows.

1. Define the dimensions of the cylindrical explosion. The target vehicle is at a given distance from the warhead with one side of the vehicle armored.

2. Use voxels to define the target vehicle. The probability that a voxel is hit by fragments will be calculated. Large sized voxels will execute quickly but provide less detailed results, while small voxels will take longer to execute but provide more detailed information.

3. Identify the circumferential fragments. This step provides the width of fragments. The size of circumferential fragments is set out in Felix and Harris’s paper. ¹¹ This is a random process, depending on the Poisson distribution, and each time the simulation is run a new set of fragments is created.

4. Identify the axial fragments. Their thickness will be the thickness of the casing and their length is set out in Felix and Harris’s paper ¹¹ (their width is set out in item 3 above). This is another random process, depending on another Poisson distribution. Each time the simulation is run a new set of fragments is created.

5. Define the shape of the expanding explosion cylinder just before it breaks up (i.e., just before the fragments leave the casing). This is set out in Felix et al.’s paper. ¹² As the casing expands, so do the fragments.

6. Calculate the initial angle of projection of fragments. This is set out in Felix et al.’s paper. ¹²

7. Calculate the initial velocity of fragments. This is set out in Felix’s paper. ¹³

8. Calculate the trajectory of fragments.

9. Determine which fragments hit the target vehicle.

10. Calculate where fragments hit the target.

11. Determine how fragments penetrate the target armor.

12. Determine the probability a voxel within the vehicle is hit by fragments.

The above steps are discussed in more detail as follows.

2. Approach

2.1. Define the cylindrical explosion and the position of the target

Figure 1 shows the overall layout of the simulation. The warhead’s cylinder stands on the ground (plane X = 0), positioned vertically, with the detonation point at coordinate (0,0,0). The length of the warhead is LL and its diameter is WW. The distance between the detonation point and the closest side of the target vehicle is D.

Figure 1.

The overall layout of the simulation.

The target vehicle is also positioned on the ground with the closest side of the vehicle to the warhead situated on the Y = D plane. The middle of the vehicle’s length is at position Z = 0.

In creating the simulation, the “standard” warhead used is a cylindrical warhead, whose length LL is 10 cm, diameter WW is 5 cm, and casing thickness t is 0.4 cm. The casing is made from hard homogeneous steel and the explosive charge is TNT. To test the validity of the model, a warhead whose length is 20 cm and diameter is 10 cm is also created. The target vehicle has a length of L, a height of H, and a width of W. In the presented simulation, L is equal to 2 m, H to 1 m, and W to 0.4 m.

Only fragments facing the target vehicle are analyzed, as only these fragments can hit the vehicle. Four coordinates define a fragment and the fragment’s thickness is equal to the thickness of the casing (Figure 2 shows the four coordinates).

1. When the case expands, fragments maintain their shape. This will result in small gaps between the long circumferential strings.

2. When the case expands, fragments expand. This will result in the fragments becoming isosceles trapezoids rather than cuboids.

Figure 2.

A fragment and its coordinates.

The simulation initially assembles the fragments on the warhead’s cylindrical casing. As the casing expands there are two different ways to provide these four coordinates.

The advantages of these approaches are as follows.

1. In the first option the mathematics is slightly more complex, and the gaps created by the expanding casing can be used to hold dust or small fragments.

2. In the second option the mathematics is less complex. As the fragments form isosceles trapezoids, the thickness of the fragment is not constant for the whole of the fragment. As the fragment’s thickness is small compared with the length or width of the warhead, the change from a cuboid to a trapezoid is small and the model can assume the fragment’s thickness is constant and equal to the thickness of the casing.

A combination of the two options is probably the most likely scenario; in the first option the casing is expanding so a fragment should change its shape and, in the second option, fragments will not have gaps between the axial lengths. This paper will develop equations that use option 2.

2.2. Use voxels to define the target vehicle

The target vehicle is divided into small cubes called voxels with the length of the sides of a voxel being variable. A simulation with large voxel lengths can quickly determine if the results are worthy of more detail and, if this is the case, smaller voxels can be used to provide more detail, but they would require a longer time to execute the simulation. Whichever size of voxel is chosen, it is possible to determine the number of times a voxel is touched by fragments.

2.2.1. The target vehicle

Figure 3 shows the orthographic projection of a target vehicle. The white space, which is external to the vehicle, is also included in the simulation. The part of the vehicle on the side elevation A, B, C, and D contains armor.

Figure 3.

The shape of the vehicle in the simulation.

Originally, the vehicle’s critical components were going to be included. However, in calculating the probability that voxels are hit by fragments, it is possible to identify positions in the vehicle that have a low and a high chance of being hit. It is then feasible to identify optimum positions for components. Consequently, the step of including components in the vehicle is not required.

2.3. Identify the circumferential fragments

This section determines the width of fragments. Fragments are assumed to be cuboids with a given width and length, generated using Poisson distributions. These assumptions are taken from Mott’s work on explosion fragments. ¹⁴ How circumferential and axial fragments are created is outlined in Felix and Harris’s paper. ¹¹

The simulation model assumes the first circumferential fragment touches the Z-axis. This limits the initial position of the first fragment but makes the simulation easier to execute, having a fixed starting point. However, this approach may have an impact on which fragments penetrate the target vehicle. As the fragments are generated randomly, the Monte Carlo method of repeating the experiment many times and then taking an average is a suggested future enhancement.

Each circumferential fragment is formed from part of the cylindrical casing (see Figure 4). It has a given angle to the Z-axis and a given length. Figure 4 shows two fragments: fragment AB, shown in red, whose length is L₁ and angle to the Z-axis is equal to $2{\theta }_1$ ; and fragment BC, also shown in red, whose length is L₂ and angle to the Z-axis is equal to $2{\theta }_1+2{\theta }_2$ .

Figure 4.

Cylindrical fragments and their dimensions. (Color online only.)

2.4. Identify the axial fragments

Each circumferential fragment creates a long string of metal and from this string several fragments are then created. ¹⁵ To define the fragment number, the circumferential fragment number forms the first part of the fragment number and the string number forms the second part of the of the fragment number (see Figure 5).

Figure 5.

The two-dimensional naming convention.

To calculate the coordinates of the fragments around the circumference of the cylinder, the value of ${\theta }_i$ (in radians) is as follows:

$${\theta }_i=L_i/2r$$ (1)

where $2{\theta }_i$ is the angle formed by the arc length of the ith width (circumferential) fragment, $L_i$ is its length, and r is the radius of the cylinder. The width of the fragment is bisected to create the midpoints of the fragments, which are the points used to determine the trajectory of the fragment.

Four coordinates that define the fragment and the midpoints A, B are shown in Figure 2. For the first circumferential fragment (1,1), as shown in Figure 5, the coordinates are denoted as (1,1,1), (2,1,1), (3,1,1), and (4,1,1).

The coordinates of fragment (i,j) where i >1 and j >1 are as follows:

$$\begin{array}{c}x(1,i,j)=\sum _{n=1}^{n=j-1}{L}_{i,n};y(1,i,j)=r\sin \left(2\sum _{n=1}^{n=i-1}\left({\theta }_n\right)\right);\hfill \\ z(1,i,j)=r\cos \left(2\sum _{n=1}^{n=i-1}\left({\theta }_n\right)\right)\hfill \\ x(2,i,j)=\sum _{n=1}^{n=j}{L}_{i,n};y(2,i,j)=r\sin \left(2\sum _{n=1}^{n=i-1}\left({\theta }_n\right)\right);\hfill \\ z(2,i,j)=r\cos \left(2\sum _{n=1}^{n=i-1}\left({\theta }_n\right)\right)\hfill \\ x(3,i,j)=\sum _{n=1}^{n=j}{L}_{i,n};y(3,i,j)=r\sin \left(2\sum _{n=1}^{n=i}\left({\theta }_n\right)\right);\hfill \\ z(1,i,j)=r\cos \left(2\sum _{n=1}^{n=i}\left({\theta }_n\right)\right)\hfill \\ x(4,i,j)=\sum _{n=1}^{n=j-1}{L}_{i,n};y(4,i,j)=r\sin \left(2\sum _{n=1}^{n=i}\left({\theta }_n\right)\right);\hfill \\ z(4,i,j)=r\cos \left(2\sum _{n=1}^{n=i}\left({\theta }_n\right)\right)\hfill \end{array}$$ (2)

$\mathrm{when}j=1;\sum _{n=1}^{n=j-1}{L}_{i,n}\mathrm{equals}0$ ,

$\mathrm{when}i=1;\sin \left(2\sum _{n=1}^{n=i-1}\left({\theta }_n\right)\right)\mathrm{equals}0\mathrm{and}\cos \left(2\sum _{n=1}^{n=i-1}\left({\theta }_n\right)\right)\mathrm{equals}1$ ,

$\mathrm{and}\mathrm{where}{L}_{i,n}$ is the length of the nth fragment along the cylinder’s axis.

2.5. Define the shape of the expanding explosion cylinder just before it breaks up

The warhead expands in the YZ plane (see Figure 6). The expansion of the warhead is symmetrical about the X-axis and the expansion at a given x value can be calculated for fragments with midpoints along an axis denoted by $Y{Z}_{\theta }$ , and 0 ≤θ≤ 360 degrees. The equation of the expanded warhead is a simplified version of the quartic equation in Felix and Harris’s paper ¹¹ and assumes the maximum expansion is at the midpoint of the warhead:

$$Y{Z}_{\theta }=-(16/3){(x-0.5)}^4+.33$$ (3)

Figure 6.

The casing expands and the fragment increases in length.

and the general formula of this equation where LL is the length of the warhead and WW its width is as follows:

$$Y{Z}_{\theta }=-{(16/3)(\mathrm{WWx}/\mathrm{LL}-0.5)}^4+.33$$ (4)

Set Equation (4) to the following:

$$Y{Z}_{\theta }=f(x,\mathrm{WW},\mathrm{LL})$$ (5)

The fragment on the cylindrical casing increases its length as the casing expands (see Figure 6).

The expanded length of a fragment is referred to as $\mathrm{new}\mathrm{length}\left(L_L(i,j)\right)$ , which satisfies the following equation:

$$\begin{array}{c}\mathrm{new}\mathrm{length}\left(L_L(i,j)\right)\hfill \\ =\sqrt{{L_{i,j}}^2+({(f(x_2,\mathrm{WW},\mathrm{LL})-f(x_1,\mathrm{WW},\mathrm{LL}))}^2}\hfill \end{array}$$ (6)

where $L_{i,j}$ is the length of the fragment before expansion and x₁ and x₂ are the midpoint x values of the fragment edges closest to and further away from the detonation, respectively. So, as the casing expands replace r in Equation (2) with $f(x,\mathrm{WW},\mathrm{LL})$ and length $L_{i,n}$ with $L_L(i,n)$ .

2.6. Calculate the initial angle of projection of fragments

The initial angle of projection of a fragment depends on the Taylor angle, developed in Felix et al.’s paper ¹² and shown in the following equation:

$$\begin{array}{c}\mathrm{Angle}=\left(\frac{1}{\lambda }\right)\left(\frac{180}{\pi }\right)\mathrm{ta}{\mathrm{n}}^{-1}\hfill \\ (-(16/3)(4\mathrm{WW}/\mathrm{LL}){(\mathrm{WWx}/\mathrm{LL}-0.5)}^3)+\mathrm{MTA}\hfill \end{array}$$ (7)

where x is the distance from the detonation and $\lambda$ is approximately set equal to 3, and the maximum velocity occurs at the warhead’s midpoint. The angle is assumed to be in radians; if this is not the case then the term $180/\pi \mathrm{can}\mathrm{be}\mathrm{omitted}.$ The term MTA is a simplified version of the equation given in Felix et al.’s paper and is shown as follows:

$$\mathrm{MTA}=-(\mathrm{TA}(1-0.65)/{0.5}^2){((\mathrm{WWx}/\mathrm{LL})-0.5)}^2+\mathrm{TA}$$ (8)

where TA is the Taylor angle and $\mathrm{TA}={\sin }^{-1}(G/2U)$ , where G is the Gurney velocity and U is the detonation velocity.

Figure 7 shows two possible Taylor angles, denoted by ${\alpha }_1$ and ${\alpha }_2$ . Unfortunately, no research papers could be found that identify how fragments break away from the casing. There are many ways the fragments can be projected from the casing. For example, one edge is projected first or both edges are projected simultaneously. The fragment is quickly forced into a direction aligning with the Taylor angle by the escaping gases. Consequently, in this paper, it is assumed the edge of the fragment closest to the detonation point is the first edge to be projected from the casing. This assumption does not change the velocity or the angle of projection of the fragment, but has a small effect on the impact position of the fragment on the target.

Figure 7.

An expanding case and fragment about to be projected.

In Figure 7 the fragment with the Taylor angle ${\alpha }_2$ has its edge furthest away from the detonation, positioned on the casing. This edge acts as an axis of rotation for the edge closest to the detonation.

Knowing the initial projection angle of a fragment ${\alpha }_2$ , the length of the fragment $L_L$ , and the two coordinates of the fragment on the casing (x₁,yz₁) and (x₂,yz₂), the coordinates of the fragment as it is about to be projected from the casing can be calculated. Figure 7 identifies these coordinates as (x₂,yz₂) for the edge on the casing and (x₂+L_Lsin ${\alpha }_2$ , yz₂+L_Lcos ${\alpha }_2$ ) for the rotated edge.

2.7. Calculate the initial velocity of fragments

The initial velocity of a fragment is calculated from the enhanced Gurney equation. The equation used in this simulation is taken from Felix’s paper. ¹³

The three-dimensional coordinates of the velocity are given by the following equation:

$$\begin{array}{c}{V}_x=V\sin \propto \hfill \\ V_y=V\cos \propto \sin ({\theta }_n+2\sum _{i=1}^{i=n-1}{\theta }_i)\hfill \\ V_z=V\cos \propto \cos ({\theta }_n+2\sum _{i=1}^{i=n-1}{\theta }_i)\hfill \end{array}$$ (9)

where V is the velocity of a fragment, which depends on the enhanced Gurney velocity, as shown in Felix’s paper, ¹³ ∝ is the initial angle of projection, and θ_i is the angle the midpoint of the ith fragment makes with the circumference.

2.8. Calculate the trajectory of fragments

This simulation assumes the exploding warhead and target vehicle are close and, if this is the case, the trajectory of fragments is a straight line, as indicated by Fitzpatrick. ¹⁶ Therefore, both gravity and air resistance can be ignored. In addition, the rotation of fragments also is ignored.

2.9. Determine which fragments hit the target vehicle

In the simulation, the closest face of the target is on the plane Y = D. The warhead is placed central to this face at a distance D from the target face. To hit the target vehicle the maximum angles the fragments can have relative to the Y-axis are shown in Figure 8.

Figure 8.

Which fragments hit the target.

The maximum angle ${\mathrm{\Phi }}_1$ relative to the Y-axis for fragments to hit the target vehicle is as follows:

$${\mathrm{\Phi }}_1=\mathrm{ta}{\mathrm{n}}^{-1}\left(\frac{\pm L}{2(D-R)}\right)$$ (10)

where R is the expanded radius of the casing, D is the distance to the target from the warhead’s axis, and L is the length of the target.

The maximum angle ${\mathrm{\Phi }}_2$ relative to the Y-axis for fragments to hit the target vehicle is as follows:

$${\mathrm{\Phi }}_2=\mathrm{ta}{\mathrm{n}}^{-1}\left(\frac{H}{(D-R)}\right)$$ (11)

where H is the height of the target and R and D are defined in Equation (10).

2.10. Calculate where fragments hit the target

Fragments are trapezoids with their axis of rotation determined by the midpoints of their two parallel sides, as shown in Figure 9.

Figure 9.

Axis of rotation of fragments.

In Figure 9(a) edge DC and Figure 9(b) edge AB are the rotation edges of the fragment.

Assuming that $(x_1,y_1,z_1)$ is the midpoint of the fragment’s edge on the casing and $(x_2,y_2,z_2)$ is the midpoint of the fragment’s edge rotated about $(x_1,y_1,z_1).$

The equation of the line joining two points $(x_1,y_1,z_1)$ and $(x_2,y_2,z_2)$ in three dimensions is as follows:

$$\begin{array}{c}x\left(t\right)=x_1+(x_2-x_1)*t\hfill \\ y\left(t\right)=y_1+(y_2-y_1)*t\hfill \\ z\left(t\right)=z_1+(z_2-z_1)*t\hfill \end{array}$$ (12)

where t is a constant. Equation (12) has different values of t to determine different points on the straight line.

Then, using Figure 7, the coordinates of the fragment rotated about the edge on the casing are ( $x_1+L_L\sin \alpha ,$ $y_1+L_L\cos \alpha \sin \theta ,$ $z_1+L_L\cos \alpha \cos \theta ).$

Figure 10 shows the Y and Z coordinates of the fragment about to be projected. As only one view of the warhead is shown, only the coordinate x₁ can be seen.

Figure 10.

Coordinates of a fragment about to be projected.

The equations of the straight line joining the midpoints of the parallel edges of the fragment are as follows:

$$\begin{array}{c}x\left(t\right)=x_1+\left(L_L\sin \mathrm{\alpha }\right)t\hfill \\ y\left(t\right)=y_1+L_L\left(\cos \mathrm{\alpha }\right)\left(\sin \mathrm{\theta }\right)t\hfill \\ z\left(t\right)=z_1+L_L\left(\cos \mathrm{\alpha }\right)\left(\cos \mathrm{\theta }\right)t\hfill \end{array}$$ (13)

where $L_L$ is the length of the fragment, $\alpha$ is the Taylor angle of projection of the fragment, t is the value of a point on the line varying from 0 to 1 for the line located on the fragment, and $\theta$ is the angle the fragment makes with the Z-axis.

It is now necessary to extend the straight line and calculate the value of t so that the fragment hits the target vehicle.

The fragment hits the target when y(t) is equal to D. So:

$$\begin{array}{c}D=y_1+L_L\left(\cos \alpha \right)\left(\sin \theta \right)t\mathrm{or}\hfill \\ \frac{D-y_1}{L_L\cos \alpha \sin \theta }=t\hfill \end{array}$$ (14)

So, solving Equation (13) using the value of t obtained in Equation (14), the coordinates of the fragment hitting the target are as follows:

$$\begin{array}{c}{x}_1+\frac{(D-y_1)L_L\sin \alpha }{L_L\cos \alpha \sin \theta };\mathrm{i}.\mathrm{e}.,x\prime =x_1+\frac{(D-y_1)\tan \alpha }{\sin \theta }\hfill \\ y_1+\frac{(D-y_1)L_L\cos \alpha \sin \theta }{L_L\cos \alpha \sin \theta };\mathrm{i}.\mathrm{e}.,y\prime =D\hfill \\ z_1+L_L\cos \alpha \cos \theta \frac{(D-y_1)}{L_L\cos \alpha \sin \theta };\mathrm{i}.\mathrm{e}.,z\prime =z_1+(D-y_1)\cot \theta \hfill \end{array}$$ (15)

The values of the coordinate ( $x\prime ,y\prime ,z\prime )$ represent where the fragment hits the target vehicle, with the initial position of the fragment located on the warhead’s casing at coordinate (x₁,y₁,z₁).

2.11. Determine how fragments penetrate the target armor

The THOR equation to determine residual velocity is as follows:

$$V_r=V_s-{10}^c{\left(\mathrm{tA}\right)}^{\alpha }{m}^{\beta }{\left(\sec \varnothing \right)}^{\gamma }{V_s}^{\lambda }$$ (16)

where the constants $c,\alpha ,\beta ,\gamma$ , and $\lambda$ for different metal targets and fragments made from steel with “no particular fragment shape” are shown in Table 1. The variable $V_s$ is the fragment striking velocity in feet per second, $m$ is the weight of the original fragment in grains (7000 grains are equal to 1 lb), $\varnothing$ is the angle between the trajectory of the fragment and the normal to the target material, t is the target thickness in inches, and A is the average impact area of the fragment in square inches.

Table 1.

THOR equation constants for residual velocity for no particular fragment shape.

Metal	c	α	β	γ	λ
Face-hardened steel	4.356	0.674	−0.791	0.989	0.434
Mild homogeneous steel	6.399	0.889	−0.945	1.262	0.019
Hard homogeneous steel	6.475	0.889	−0.945	1.262	0.019

For full penetration, the value of $V_r$ (the residual velocity) is set to zero (the velocity of the fragment exiting the target armor is zero, so the fragment just penetrates the target material). Initially, the thickness of the armor (t) is set to 4 mm and the striking area of the fragment (A) is calculated from the width of the fragment multiplied by its thickness. Also, it is assumed the angle of incidence of the fragment ( $\varnothing$ ) will not ricochet if the angle ≥80 degrees. Below this value, it is assumed the fragment will ricochet.

As fragments have a high velocity when they impact the target vehicle, research suggests that they will lose mass in penetrating the armor (see THOR ¹⁷ ). This is because the target material resists being penetrated by fragments. THOR ¹⁷ indicates this resistance affects the fragment through a loss in velocity and weight. However, this paper is only interested in which voxels are touched by fragments, which is not dependent on the fragment’s mass. The fragment residual mass is therefore ignored.

The angle of impact is given by the following equation:

$$\begin{array}{c}\mathrm{si}{\mathrm{n}}^{-1}\frac{\left(D\left(\cos \alpha \right)\left(\cos \theta \right)\right)}{\sqrt{D^2}\sqrt{\mathrm{si}{\mathrm{n}}^2\alpha +\left(\mathrm{co}{\mathrm{s}}^2\alpha \right)\left(\mathrm{si}{\mathrm{n}}^2\theta \right)+\left(\mathrm{co}{\mathrm{s}}^2\alpha \right)\mathrm{co}{\mathrm{s}}^2\theta )}}\hfill \\ =\mathrm{si}{\mathrm{n}}^{-1}\left(\left(\cos \alpha \right)\left(\cos \theta \right)\right)\hfill \end{array}$$ (17)

2.12. Determine the probability a voxel within the vehicle is hit by fragments

The steps taken to determine the probability a voxel within the vehicle is hit by fragments are as follows.

Firstly, calculate the impact position of the fragments on the vehicle. As the fragments are assumed to travel in straight lines, the initial coordinates, velocity, and angle of projection of a fragment enable the simulation to determine where the fragments hit the target vehicle and their angle of incidence.

Secondly, the length, width, and height of the vehicle space is divided into voxels. The fragments that hit the target vehicle and penetrate the armor then pass through the vehicle voxel space. Some of the voxels in the vehicle voxel space are touched by the fragments. Each voxel contains a count of the number of fragments touching it, and at the end of the simulation the probability a voxel is hit is equal to the count in the voxel divided by the sum of all the voxel counts.

Having empty voxels inside the vehicle is a simplistic view because a vehicle contains many components. A future improvement could include these components, contained in one or more voxels. These voxels then need to include information of how they are affected by a fast-moving fragment.

Finally, calculate which voxels are touched by fragments. Fragments travel in a straight line (i.e., like a ray of light) and the method to calculate their path in the vehicle voxel space is based on the ray tracing paper by Amanatides and Woo. ¹⁸ The Amanatides and Woo model, in which light rays touch voxels, is modified to calculate how fragments, which travel in straight lines, traverse the voxel space.

Figure 11 provides a flowchart showing how fragments hitting the voxel space of the vehicle are counted.

Figure 11.

Flowchart of fragments hitting a vehicle.

3. Validation

The steps outlined in Section 2 were programmed in MATLAB. This means the execution time probably can be improved by using another programming language.

The simulation is validated by creating witness plates of fragments hitting the target vehicle. These can then be compared with experimental witness plates and witness plates created by alternative simulations.

The timing of the complete simulation can be compared with other simulations’ timings. The objective of this paper is to have a model that is faster than models using FEA models, so the comparison should be against FEA models.

4. Results and discussion

4.1. Results of simulation

Figure 12 shows examples of where fragments in the simulation first impact the vehicle. In Figure 12, “*” represents where the center of fragments hit the vehicle.

Figure 12.

Examples of patterns made by fragments.

Each picture in Figure 12 represents a different execution of the simulation. The size of the vehicle and the warhead are constant, but the creation of fragments depends on a random distribution, so the witness plates vary for each simulation. The size of fragment impacts is not shown in the simulation. Many fragments form a crater when they hit the vehicle’s external wall, and this can be included in the simulation by plotting circular impact craters with the width of the fragment equal to the diameter of the circle. This was not included in the simulation.

4.2. Comparison with other researchers’ experiments

The only simulations that provide witness plates are for FEA simulations, which is because FEA simulations are assumed to be more accurate than this paper’s witness plates.

Chen et al. ¹⁹ carried out experiments and simulations using the FEA program LS-Dyna, for cylindrical warhead explosions. An example of one of their simulated witness plates is shown in Figure 13. In this simulation the warhead contained JHL-3 explosive, and the casing had an outer diameter of 180 mm, height of 200 mm, and unknown thickness. The steel target was 1.5 m long and 1.26 m wide.

Figure 13.

Witness plate provided by Chen et al.’s simulation. ¹⁹

The Shapiro equation was used to determine the initial angle of projection of fragments and the Gurney equations to determine the initial velocity.

Ding et al. ²⁰ performed similar experiments and simulations, also using the FEA program LS-Dyna program, on cylindrical warhead explosions. An example of one of his witness plates is shown in Figure 14. In the simulation, the warhead contained TNT explosive, and the casing had an outer diameter of 100 mm, height of 80 mm, and unknown thickness. The steel target was 3.5 m from the warhead.

Figure 14.

Witness plate provided by Ding et al.’s simulation. ²⁰

The Shapiro equation was used to determine the initial angle of projection of fragments and the Gurney equations to determine the initial velocity.

There are similarities and differences between the results shown in Figure 12 and the simulations from Chen et al. ¹⁹ and Ding et al. ²⁰ The three simulations approximately form straight lines of fragment hits. The Chen et al. simulation shows more fragment hits than Ding et al.’s and this paper’s simulation.

In Figure 12, the “*” form straight lines because the simulation excludes yaw angle rotation in the fragment’s flight. Including yaw angle rotation causes the fragment to be subjected to a small amount of movement in the Z direction (this is the vehicle length axis). This will make the fragments represented on the witness plate more realistic. It is unknown if Chen et al. ¹⁹ and Ding et al. ²⁰ included a yaw angle, but their witness plates suggests they do not.

In addition to the above omission, Figure 12 does not include small fragments that are created as the casing expands and the long strings of circumferential fragments are ripped apart from one another. Again, it is unknown if Chen et al. ¹⁹ and Ding et al. ²⁰ include small fragments, but their witness plates suggest they do.

4.3. Comparison with experimental data

This section examines experimental results provided by three researchers and compares their results with this paper’s and two FEA simulation results.

Arnold and Rottenkolber’s paper ²¹ shows physical explosion witness plates created by fragment indentations. An example of a witness plate is shown in Figure 15.

Figure 15.

Example of an Arnold and Rottenkolber ²¹ witness plate.

Figure 16 also shows an image of fragment impacts on a witness plate. This figure is obtained using an infrared thermographic lock-in method that enhances the view of the fragment impacts, and the damage to the surface is visible around the fragment indentations. ²²

Figure 16.

Swiderski et al. ²² image of fragments hitting the plate.

4.3.1. Number of impacts for experimental simulations compared with FEA simulations

The number of impacts shown in the Chen et al. ¹⁹ and Ding et al. ²⁰ FEA simulation witness plates are greater than the number of impacts in the experimental witness plates. This is probably because the FEA simulation witness plates contain all the impacts, including small particle impacts, and it is probably difficult to see small impacts in the experimental witness plates.

4.3.2. Number of impacts for experimental simulation compared with this paper’s simulation

The number of impacts in this paper’s witness plates is similar to the number of impacts in the experimental witness plates. Although the number of impacts depend on factors such as the size of the warhead and the type of explosive, the number of fragments in this paper’s simulation correlates well with Mott’s accurate distribution of the number and weight of fragments. This suggests the number of impacts in this paper’s simulation is probably acceptable.

4.3.3. Are fragments grouped in the simulation?

Witness plates created by this paper’s simulation show many impacts occurring in a straight line with wide gaps between the straight lines. The simulation’s width fragments create long strings, which without fragment rotations form a straight line of fragment impacts. The gaps between these lines increases as the distance between the explosion and the target increases and, as the fragments do not rotate, these gaps are more apparent.

These observations suggest that small rotations should be included in the simulation. The work by De Vuyst et al. ²³ looks promising and a suggested partitioning of the casing’s surface, as shown in Figure 17, is worth investigating.

Figure 17.

An example of a split casing.

The experimental witness plates show a “random” group of fragment impacts. Figure 17 shows witness plates of fragment impacts from Arnold and Rottenkolber, ²¹ An et al., ²⁴ and Swiderski et al. ²² Although these witness plates indicate “random” fragment hits, the inserted lines, which “connect” three or more impacts, show that fragments hitting a target may be linked.

This indicates this paper’s simulation produces witness plates that are not totally inaccurate. Including fragment small rotations and dividing the length of the casing into smaller lengths, as shown by De Vuyst et al., ²³ would make the simulation more realistic. The simulation also can be improved by changing the “*” impacts into impact craters that depend on the known width of fragments.

4.4. Execution time

No similar simulations could be found to provide alternative execution times. The execution time of the simulation is 0.22 seconds, for 10 executions. This gives an average execution time of about 0.02 seconds. The code has not been optimized so this average time can probably be reduced.

The timings for FEA simulation are provided in papers by Babu et al. ²⁵ and The Lawrence Livermore National Laboratory and show that the execution times of FEA simulations are many times slower than the timing of this paper’s simulation.

Table 2 shows the breakdown of the average execution time for the steps in this simulation.

Table 2.

Approximate execution times of simulation modules.

Step	Simulation module(s)	Time (milli seconds)
1	Create fragments on cylinder casing	12
2	Expand the casing and calculate the angle of projection and velocity for each fragment	0.5
3	Move fragments ready for projection	0.7
4	Calculate flight trajectories of fragments	0.8
5	Calculate impact of fragments on target and armor penetration	0.2
6	Move fragments through vehicle voxel space	8
7	Calculate probabilities of voxels being hit	0.1
	Total	22

4.4.1. The main variables that affect execution time

The main variables that affect execution timings in Table 2 are described in more detail as follows.

Step 1. For the described warhead the simulation creates about 2200 fragments. The main variable to determine the execution time is the number of created fragments. After all the fragments have been created about one third are available to hit the target vehicle (see Figure 8).

Steps 2–5. The execution times depend on the percentage of fragments that hit the target (see Figure 8).

Step 6. The execution time depends on how the fragments move through the voxel space and the size of a voxel. Based on the Amanatides and Woo model, ¹⁸ Figure 18 is used to provide an approximate value for the maximum number of voxels hit by a fragment.

Figure 18.

Estimated number of voxels hit by a fragment. (Color online only.)

The blue line represents the maximum fragment path and the gray shading represents fragment paths for fragments hitting the vehicle at point A. The vertical lines AB and AC together with the blue lines and gray shading can move horizontally, depending where the fragment initially hits the vehicle. The fragment is assumed to ricochet if the angle, in degrees, between a blue line and a vertical line is less than 90 – 80, where 80 represents the angle below which the fragment will ricochet.^26,27

For the simulation’s length/height combination, the maximum number of voxels hit by a fragment is equal to $10+\mathrm{roundup}(10/\left(\tan 80\right))=12$ , where 10 represents the number of voxels in the height of the vehicle.

For the simulation’s length/width combination, the maximum number of voxels hit by a fragment is equal to $4+\mathrm{roundup}(4/\left(\tan 80\right))=5$ , where 4 represents the number of voxels in the width of the vehicle.

The maximum number of voxels hit by a fragment is estimated by calculating which of the three velocity components will cause the fragment to exit the vehicle. Generally, the X velocity component is the smallest component and so most fragments exit the vehicle along the Y- or Z-axes. In our simulation, as the length of the vehicle is significantly greater than its width most fragments will exit the vehicle along the Y-axis. So, the maximum number of voxels hit by a fragment is estimated to be 5.

Step 7. The execution time depends on the number of voxels hit by fragments, which is estimated in step 6.

4.4.2. Reduce the simulation’s voxel size

If the voxel’s size is reduced, for example by a factor of 10, the execution time of step 6 will increase. Divide each side of a voxel into 10 parts then the number of small voxels (1/1000 the size of the original voxel) hit by a fragment is $10+\mathrm{roundup}(10/\left(\tan 80\right))=12,$ where 10 represents the reduction factor. So, reducing a voxel by a factor of 10 leads to multiplying the execution time of step 6 by a maximum of 12.

4.4.3. Estimate execution times for a different cylinder warhead

The estimated execution times for a different warhead can be calculated as follows.

Step 1. Use Mott’s equation to estimate the number of fragments. Calculate the ratio of the number of fragments divided by 2200. Multiply this ratio by the step 1 execution time in Table 2 to estimate the execution time for the different warhead’s step 1.

Steps 2–5. Calculate the fraction of fragments that can hit the target. Multiply the steps 2–5 execution times in Table 2 by this fraction to estimate the steps 2–5 execution times for the different warhead.

Step 6. For the different cylinder warhead and this paper’s warhead, calculate the average number of voxels hit by a fragment (as shown in Section 4.4) multiplied by the fraction of fragments that can hit the target (see Figure 19). Calculate the ratio of the different cylinder warhead’s result divided by this paper’s warhead’s result. Take this ratio and multiply it by the step 6 execution time in Table 2 to provide the estimated step 6 execution time for the different warhead.

Figure 19.

Witness plates from Arnold and Rottenkolber, ²¹ An et al., ²⁴ and Swiderski et al. ²²

5. Conclusions

This study presents the simulation of a fragmenting explosion.

The witness plates of the initial impact of fragments on the vehicle are compared with available witness plates from two simulated explosions by Chen et al. ¹⁹ and Ding et al. ²⁰ The witness plates are similar, but the comparison shows the simulation has fewer fragment hits. It is suspected this is because small fragments are not included in the simulation.

The simulation is also compared with experimental witness plates created by Arnold and Rottenkolber, ²¹ An et al., ²⁴ and Swiderski et al. ²² The witness plates are similar and there is an indication that fragment hits and the same angle of projection form a straight line. Delivering fragment impacts as craters rather than points and including small yaw angle rotations would enable the simulation to be more accurate.

In the computer execution time tests, the computation time of the simulation was about 22 milliseconds, many times faster than FEA simulations.