Numerical Simulation of Conical and Linear-Shaped Charges Using an Eulerian Elasto-Plastic Multi-Material Multi-Phase Flow Model with Detonation.

This study developed a hydrocode to numerically simulate both conical and linear-shaped charges using an Eulerian multi-material and multi-phase flow model. Elasto-plastic solids and the detonation of a high explosive charge were modeled using a Johnson-Cook material model and the programmed burn model, respectively. Further, the plasticity of the solids was calculated using a radial return mapping algorithm. The model was solved using a high-resolution computational fluid dynamics (CFD) technique on Cartesian grids. Material interfaces were tracked using the level-set method, and the boundary conditions were imposed using the ghost fluid method. The developed hydrocode was validated using high-speed impact problems. Consequently, the developed hydrocode was used to successfully simulate the evolution and penetration of metal jets in shaped charges after a detonation.

1. Introduction

A shaped charge is a highly explosive charge for focusing explosive energy in one direction using the Monroe effect (or the Neumann effect) [1]. Generally, a shaped charge consists of a high explosive with a hollow cavity, a thin metal liner, a detonator, and a metal casing. When an explosive is ignited, a high-pressure and high-temperature explosive gas, whose pressure is concentrated in the center of the metal liner, is generated, which results in large deformation. Consequently, this results in the plastic deformation of the liner beyond its elastic limit and the formation of a hypervelocity metal jet with a typical speed and temperature of 8–12 km/s and approximately 1200 K, respectively. Shaped charges are used in numerous military and industrial applications, such as high-explosive anti-tank warheads [2], emergency escape devices for aircraft pilots, stage and fairing separation devices for space rockets [3], and for the removal and destruction of engineering facilities. Depending on their shape, shaped charges are mainly divided into conical- and linear-shaped charges, wherein the conical-shaped charge is used for penetration-related applications, and the linear-shaped charge is used for cutting objects.

Accordingly, numerous studies have investigated the types of shaped charges. Walters and Zukas [1] summarized the general theory of a shaped charge. In addition, Molinari [4] investigated the conical-shaped charge using the finite element method. Further, Liu et al. [5,6] compared two types of conical-shaped charges with conic and semi-elliptic cavities using the meshfree particle simulation. However, no liner and casting materials were added to the model. Ma et al. [7] employed a hydrocode to analyze conical-shaped charges consisting of steel and copper liners using the operator splitting method, and compared the penetration performance of the shaped charges at various liner angles. Kim and Yoh [8] presented the boundary conditions for the reactive ghost fluid method in energetic multi-material flows. Zhang et al. [9] compared a conical-shaped charge made of a copper liner to that made of a tungsten-copper liner using a commercial software (AUTODYN). Chen and Liu [10] simulated a conical-shaped charge using an elasto-plastic solid model in the Eulerian framework and solved it using the CE/SE scheme. However, because no casting material was added to the configuration, the accuracy of the computed results was limited. Sambasivan et al. [11,12] calculated the large deformation and penetration problems of high-speed elasto-plastic solids using the Convex Essentially Non-Oscillatory (CENO) scheme and the ghost fluid method. However, since their model did not consider a high explosive, it cannot simulate the shaped charge problems. Feng et al. [13] computed a linear-shaped charge using the smoothed particle hydrodynamics (SPH) technique. Recently, Yi et al. [14] analyzed a shaped charge made of a polymer liner rather than metal. Pyka et al. [15] simulated a shaped charge using a commercial software (ABAQUS), and modeled the metal jet using the smoothed particle hydrodynamics method. Du et al. [16] numerically simulated conical-shaped charges to investigate the penetration performance of a jet using a commercial software (LS-DYNA), and compared the numerical results to the experimental results. Peng et al. [17] experimentally analyzed the penetration–explosion effects of conical-shaped charges, and performed numerical simulation using a commercial software (AUTODYN).

Lagrangian and arbitrary Lagrangian and Eulerian (ALE) methods, except the SPH method, involve considerable complexity due to the need for mesh management to handle large deformations of the boundary. Therefore, periodic re-meshing is required to maintain good mesh quality, which increases the computational cost. Although the SPH method has the advantage of handling the large deformation efficiently, it is difficult to accurately represent the boundary surface with this method compared to the Eulerian method. As the metal jet of a shaped charge undergoes large deformation within a very short time and behaves similarly to a compressible fluid flow, it is advantageous to simulate it in an Eulerian framework using modern CFD gas dynamics techniques. Therefore, this study developed a hydrocode to simulate conical- and linear-shaped charges in an Eulerian framework on Cartesian meshes. Elasto-plastic solids were modeled using the CENO scheme [18] and the ghost fluid method [11,12,19] to simulate the large deformation and penetration of a metal liner and a casing. Further, the detonation of a high explosive was simulated using the programmed burn model [10]. For the multi-material interface treatment, the precise boundary conditions between an explosive gas, solids, and void were presented in detail. The developed hydrocode was validated using previously reported numerical and experimental data for 1-D solid impact and the Taylor bar problems.

The main difference between the developed hydrocode and the previous ones is that the present hydrocode newly combines the high-resolution elasto-plastic solid model [11,12] with the detonation model of explosive gas [10] on an Eulerian framework to simulate shaped charge problems. In addition, unlike Chen and Liu [10], since the present hydrocode can be applied regardless of the number of materials, it can solve more practical shaped charge problems.

2. Governing Equations

The governing equations for elasto-plastic solids under high-speed and high-strain rates in two-dimensional or axisymmetric geometry in an Eulerian framework can be represented as follows [11,12]: where where is the density; and are the x- and y-velocity, respectively; is the specific total energy; is the pressure; is the deviatoric stress tensor; and is the source vector. To calculate the temperature change caused by plastic deformation, auxiliary equations for effective plastic strain and temperature are required. The stress tensor, , can be divided into a dilatational term, , and a deviatoric term, , as where is obtained using the Mie–Grüneisen equation-of-state (EOS) and is calculated using the plastic flow theory. This study utilized the Johnson–Cook model [20] as the material model for the elasto-plastic solids. The plastic deformation was calculated using the von Mises J2 flow theory. To consider the plastic deformation of elasto-plastic solids, a radial return mapping algorithm [21,22] was applied. Further details on the governing equations, the auxiliary equations, the material model, and the radial return mapping algorithm appear in Appendix A.

A high-speed explosive gas can be described using the hyperbolic conservation laws of gas dynamics. The governing equation for an explosive gas in two-dimensional or axisymmetric coordinates is where

To simulate the detonation process of an explosive charge, a pressure-based programmed burn model [10] was used in the hydrocode developed in this study. The internal energy of the explosive product (gas) was modeled using the Jones–Wilkins–Lee (JWL) EOS. Table 1 shows the parameters of the JWL EOS of TNT and Composition B explosives. More details on the programmed burn model are given in Appendix A.

Table 1

Jones–Wilkins–Lee (JWL) equation-of-state (EOS) parameters of TNT and Composition B.

	ρ_0 (kg/m^3)	p_CJ (GPa)	v_d (m/s)	A (GPa)	B (GPa)	R_1	R_2	ω	I_0 (J/kg)
TNT	1630	21	6930	371.2	3.210	4.15	0.95	0.30	4.29×10^6
Composition B	1717	29	7980	524.2	7.678	4.20	1.10	0.34	4.95×10^6

3. Numerical Method

The governing equation for elasto-plastic solids was discretized using the third-order Runge–Kutta method in time and the third-order CENO scheme [18] using the finite difference method and the Lax–Friedrichs flux in space, which is similar to that reported in a previous study [11]. As the CENO scheme uses a component-wise method that does not consider the complicated characteristic decomposition of the governing equations via eigenvalues and eigenvectors, the scheme is somewhat less accurate than the well-known weighted ENO (WENO) scheme [23]. However, this method is easy to implement and requires lower computational costs. The third-order Runge–Kutta scheme and the fifth-order WENO scheme (WENO5) [23] are used to solve the auxiliary equations and the explosive gas equations. The level-set method [24] is used to identify the boundary of each material in a multi-material environment. The details on the numerical methods are reproduced in Appendix B for completeness.

Ghost Fluid Method

The ghost fluid method is used for boundary conditions at the material interface [11,12,19]. First, the physical properties on the interface are extrapolated around the interface using the following equation: where

Applying the forward Euler method in time and the first-order upwind method in space to Equation (6), the following equation is derived where the time step, , is used. Since the purpose of Equation (9) is to populate a thin layer of ghost cells around the interface needed for the numerical method by solving it for a few time steps, it is not necessary to obtain an accurate solution in time and space; the first-order accuracy is sufficient.

According to the physical variables, there are four boundary conditions: Neumann, Dirichlet, continuity, and extrapolation conditions. The value of the “ghost” node can be calculated from the values of the internal “real” nodes, as shown in Figure 1.

Figure 1

Boundary conditions of the ghost fluid method.

The velocity vector and the stress tensor should be rotated using the local coordinates system on the interface. The unit normal vector, , the unit tangent vector, , and the rotation tensor, , on the local coordinate system are defined as

Thus, the rotated component in the local coordinate can be calculated using

The values of the “ghost” nodes are imposed according to the type of boundary. In this study, the following four cases were considered.

Solid–solid interface:

Solid–void (or gas–void) interface:

Solid–gas interface:

Gas–solid interface:

The “real” nodes values, and , were calculated by interpolating from the values of the surrounding internal nodes. The interpolation function can be expressed as where is the coefficients. If there are four neighboring nodes around and , including the boundary, the bilinear interpolation is used. For the Dirichlet condition, the coefficients are calculated as

For the Neumann condition, the coefficients can be expressed as

If there are less than four neighboring nodes around or , the least squared interpolation is used, which can be expressed as

To determine whether two materials collide or not, the following condition is used [12]: where and are the level-set functions of two materials and .

As the deviatoric stress tensor in the “ghost” node obtained using the ghost fluid method does not satisfy the invariant condition, , the stress tensor can be corrected as

The overall procedure of the present method is summarized in Algorithm 1.

Generate mesh and setup initial conditions.

Moving level-set fields, , using the level-set method.

Populate the ghost cells around the interface by solving Equation (9).

Impose boundary conditions according to the type of boundary by using the ghost fluid method.

Solve the governing equation (Equation (1)) to obtain the elastic predictor step for solids.

Apply the radial return mapping algorithm to perform the plastic corrector step.

Solve the auxiliary equations for the temperature of the solids.

Apply the programmed burn model to the explosive charge.

Solve the governing equation (Equation (4)) for explosive gas.

Update the time step.

Go back to step 2 and repeat until the last time step is reached.

4. Results and Discussion

4.1. One-Dimensional Shockwave Propagation in Elasto-Plastic Solid

This problem considered the propagation of shockwaves in a solid after a 20 m/s impact on one side in one dimension (Figure 2). Table 2 shows the parameters of the Mie–Grüneisen EOS of the solid used in this problem.

Figure 2

Shockwave propagation in solids by impact.

Table 2

Parameters of the Mie–Grüneisen EOS.

Γ	s	c_0 (m/s)	ρ_0 (kg/m^3)	e_0 (J/kg)	G (GPa)	Y (MPa)
2.0	1.49	3940	8930	110,920	45	90

Figure 3 shows the computed shockwave propagations of the elastic and elasto-plastic solids. As shown in Figure 3a, only one shockwave propagates inside the solid in a pure elastic solid; however, in an elasto-plastic solid, first, a weak elastic wave known as an elastic precursor propagates, after which a strong plastic wave propagates, as shown in Figure 3b. The comparison of these results to those by Udaykumar et al. [25] in Table 3 indicates that the two results are consistent.

Figure 3

Computed shockwave propagation in an (a) elastic and (b) elasto-plastic solid by a 20 m/s impact at the left end.

Table 3

Comparison of the results of this study to those of Udaykumar et al. [25].

	Elastic Precursor		Plastic Shock
	Udaykumar et al. [25]	Present	Udaykumar et al. [25]	Present
Density (kg/m^3)	8938.9	8939.8	8973.5	8973.3

4.2. Taylor Bar Impact Problem

The Taylor bar problem, which has been simulated by numerous researchers, considers that a conical copper rod collides with a rigid wall at high speed (227 m/s), as shown in Figure 4. The copper rod has an initial radius of 3.2 mm and a length of 32.4 mm, and is considered an elasto-plastic material.

Figure 4

Schematic of the Taylor bar impact problem.

Figure 5 shows the comparison of the numerical results of the pressure and effective plastic strain in the copper rod obtained in this study at to those obtained by Sambasivan et al. [12], where the deformed shape of the rod and the inner contour plots are similar to each other. Table 4 summarizes the comparison of the computed results of the final geometry of the copper rod obtained in previous studies to those obtained in this study. Although the final shape may exhibit slightly different computational results depending on the physical model and numerical technique used, the computed results of this study were consistent with those of previous studies within a reasonable error range. The computed results for the maximum show slight differences from those of previous researchers. As an analysis of the causes of these differences is beyond the scope of this study, further comparison and analysis of the different models and numerical methods may improve the present method in future studies.

Figure 5

Interface shape and the contours of the pressure and effective plastic strain in the copper rod after collision at . Results of (a) this study and (b) Sambasivan et al. [12]. Reprinted with permission from Elsevier.

Table 4

Comparison of the results of the final geometry of the copper rod in this study to those of previous studies.

	Final Height (mm)	Final Base Radius (mm)	Maximum ε¯_p
Present	21.44	6.65	2.672
Sambasivan et al. [12]	21.53	7.05	3.169
Mehmandoust and Pishevar [26]	21.80	6.53	2.178
Tran and Udaykumar [27]	21.15	7.15	2.86
Udaykumar et al. [25]	21.4	6.97–7.24	-
Camacho and Ortiz [28]	21.42–21.44	7.21–7.24	2.97–3.25

4.3. Generation and Penetration of a High-Speed Metal Jet

This problem was designed to verify whether the present hydrocode can successfully simulate the generation of a metal jet and the penetration of the target material. Figure 6 shows a schematic of the problem, where a strong impact () was applied to the bottom of the hemispherical groove of the copper liner with a thickness of 29 mm and a radius of 15 mm, and a 50 mm thick steel plate was placed 71 mm apart from the copper plate. The high-speed metal jet was created by a shockwave after a strong impact to the opposite side of the copper liner with a hemispherical groove, which is a similar mechanism to that of the jet generation of shaped charges.

Figure 6

Schematic of the high-speed copper jet penetrating steel problem.

Figure 7 shows the computed interface shapes and density contours at five time instants up to after impact, where the generation and growth of the copper jet were well simulated. The results of this study on the evolution of the interface were similar to those of previous studies [12,29]. Figure 8 shows the interface shapes and density contours plotted at eight time instants up to after impact. The image revealed that the current hydrocode can successfully simulate the solid penetration process by a high-speed metal jet.

Figure 7

Evolution of the metal jet by the impact of a hemispherical groove in copper. The interface shapes and density contours at five time instants up to after impact are shown.

Figure 8

Computed results of the penetration of steel by the copper jet generated from the hemispherical groove. The interface shapes and density contours at eight time instants up to after impact are shown.

4.4. Conical-Shaped Charge

A conical-shaped charge problem was numerically simulated using the developed hydrocode to examine the detonation of a high explosive and the resulting jet of metal liner. Figure 9 shows the schematic of this problem. Composition B, copper, and aluminum were used for the explosive charge, liner, and casing, respectively. The explosive was assumed to detonate from the bottom plane.

Figure 9

Schematic of a conical shaped charge problem.

Figure 10 shows the expansion of the explosive gas and the evolution of a metal jet by the detonation of an explosive charge at seven different time instants. As the detonation wave propagated, the pressure was concentrated in the center of the copper liner, and a jet began to form. As the wave propagation proceeded, a high-speed thin jet was formed in the front of the liner, the lateral wing of the liner contracted towards the center, and a very slow slug that occupied most of the mass emerged in the rear part of the liner. In addition, the angle of the lateral liner wing gradually increased, and became horizontal (90°) at approximately , after which it is formed in the reverse direction. In addition, the aluminum case expanded owing to the pressure of the detonation wave, and the thickness of the case reduced at a specific part. Figure 11 shows the density, pressure, and temperature contours at after detonation, where the density of the liner was high at the front of the jet and the center of the slug, most of the pressure was concentrated in the central part of the liner, and the temperature was high along the central axis of the liner. At , the maximum density, pressure, and temperature of the jet were , , and , respectively. The evolution of the jet demonstrated in this study was consistent to those of previous experiments and simulations [4,5,6,7,9,10,30,31,32], indicating that the hydrocode developed in this study can successfully simulate the conical-shaped charge.

Figure 10

Snapshots of the jet evolution of a conical-shaped charge. The interface shapes and density contours are plotted at seven time instants up to after detonation.

Figure 11

Metal jet of a conical-shaped charge at : (a) density contour, (b) pressure contour, and (c) temperature contour.

4.5. Linear-Shaped Charge

Figure 12 shows a linear-shaped charge problem, where an aluminum liner, which was also used as the casing of the explosive charge with a thickness of 0.58 mm, surrounds a TNT explosive. As the linear-shaped charge has a two-dimensional geometry, compared to the conical-shaped charge, the cross-section of the explosive of the linear-shaped charge detonated at the same time. In addition, as there was no additional casing metal, which is usually thicker than the liner, in the linear-shaped charge, the case was destroyed earlier than in the conical-shaped charge, thus resulting in the rapid release of the pressure of the explosive gas to the external environment. Furthermore, the length of the jet of the linear-shaped charge was significantly shorter than that of the conical charge owing to the two-dimensional geometry and the pressure loss of the linear-shaped charge. Because of these characteristics, the linear-shaped charge is generally more difficult to numerically simulate than the conical-shaped charge.

Figure 12

Schematic of a linear-shaped charge problem.

Figure 13 shows the expansion of the explosive gas by the detonation, the deformation of the liner, and the formation of the jet. Most of the pressure generated by the two-dimensional detonation was concentrated in the center of the liner to create a jet, and some was used to expand the other part of the liner. Over time, a high-speed jet was formed in the front of the liner and a low-speed slug in the rear part of the liner. The pressure of the detonation wave reduced the thickness of a specific part of the case, and some parts were fragmented at approximately . In addition, the length of the jet was considerably shorter than that of the conical-shaped charge, and the surrounding casing moved forward together as the jet developed, and eventually, only the jet part was separated from the casing. Figure 14 shows the pressure, density, and temperature contours at after detonation. The casing was fragmented in several parts, through which the explosive gas escaped. In addition, most of the pressure was concentrated on the central wing portion of the jet, indicating a high temperature along the central axis of the jet. At this time, the jet was almost separated from the surrounding material and moved independently. These results are consistent with the experimental and numerical results of previous studies [13,33,34], indicating that the developed hydrocode can successfully simulate the linear-shaped charge.

Figure 13

Snapshots of the jet evolution of a linear-shaped charge. The interface shapes and pressure contours are plotted at six time instants up to after detonation.

Figure 14

Jet formation of a linear-shaped charge at : (a) pressure contour, (b) density contour, and (c) temperature contour.

5. Conclusions

In this study, a conical- and a linear-shaped charge were numerically simulated using an Eulerian multi-material and multi-phase flow model for elasto-plastic solids and a high-explosive gas with a detonation. The elasto-plastic solids and the detonation of a high explosive charge were modeled using the Johnson–Cook material model and the programmed burn model, respectively. In addition, the stress was corrected back to the yield surface using the radial return mapping algorithm. The model was solved on the Cartesian grids using a third-order CENO and the third-order Runge–Kutta method for the elasto-plastic solids, and using a fifth-order WENO and the third-order Runge–Kutta method for the constitutive equations and the explosive gas equations. The boundary conditions on the material interfaces were imposed depending on the two materials crossing the boundary using the level-set and ghost fluid methods. The developed hydrocode was validated using two high-speed impact problems, and the hydrocode was used to successfully simulate the evolution and penetration of metal jets in the shaped charges after the detonation of an explosive charge. The developed hydrocode is expected to help in the simulation of large deformation and penetration problems of elasto-plastic solids other than shaped charges. In addition, future studies will investigate shaped charges with a three-dimensional geometry, which may require adaptive mesh refinements and parallel processing techniques.