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.


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.

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]: Materials 2022, 15, 1700 where ρ is the density; u and v are the xand y-velocity, respectively; E is the specific total energy; p is the pressure; s ij is the deviatoric stress tensor; and S 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, σ ij , can be divided into a dilatational term, −pδ ij , and a deviatoric term, s ij , as where p is obtained using the Mie-Grüneisen equation-of-state (EOS) and s ij 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 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.

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 fifthorder 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: 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, ∆τ = 0.5∆x, 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. 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 The velocity vector and the stress tensor should be rotated using the local coordinates system on the interface. The unit normal vector, n, the unit tangent vector, t, and the rotation tensor, R, on the local coordinate system are defined as n = n x , n y (10) 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, R1 and R2, were calculated by interpolating from the values of the surrounding internal nodes. The interpolation function can be expressed as where a is the coefficients. If there are four neighboring nodes around R1 and R2, 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 R1 or R2, 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 φ 1 and φ 2 are the level-set functions of two materials and ε = CFL · ∆x. As the deviatoric stress tensor in the "ghost" node obtained using the ghost fluid method does not satisfy the J 1 invariant condition, J 1 ≡ s nn + s tt = 0, the stress tensor can be corrected as The overall procedure of the present method is summarized in Algorithm 1.

Algorithm 1
The numerical procedure for the present hydrocode.

1:
Generate mesh and setup initial conditions. 2: Moving level-set fields, φ i , using the level-set method.

3:
Populate the ghost cells around the interface by solving Equation (9). 4: Impose boundary conditions according to the type of boundary by using the ghost fluid method. 5: Solve the governing equation (Equation (1)) to obtain the elastic predictor step for solids. 6: Apply the radial return mapping algorithm to perform the plastic corrector step. 7: Solve the auxiliary equations for the temperature of the solids. 8: Apply the programmed burn model to the explosive charge. 9: Solve the governing equation (Equation (4)) for explosive gas. 10: Update the time step. 11: Go back to step 2 and repeat until the last time step is reached.

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.

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 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 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 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.

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

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 5 shows the comparison of the numerical results of the pressure and effective plastic strain in the copper rod obtained in this study at = 80 μs 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 previ-  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 t = 80 µs 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  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 ε p 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 shows the comparison of the numerical results of the pressure and effective plastic strain in the copper rod obtained in this study at = 80 μs 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.
(a) (b) Figure 5. Interface shape and the contours of the pressure and effective plastic strain in the copper rod after collision at t = 80 µs. Results of (a) this study and (b) Sambasivan et al. [12]. Reprinted with permission from Elsevier.

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 (v = 540 m/s) 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. 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 ( = 540 m/s) 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 7 shows the computed interface shapes and density contours at five time instants up to 100 µs 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 100 µs after impact. The image revealed that the current hydrocode can successfully simulate the solid penetration process by a high-speed metal jet. those of previous studies [12,29]. Figure 8 shows the interface shapes and density contours plotted at eight time instants up to 100 μs after impact. The image revealed that the current hydrocode can successfully simulate the solid penetration process by a high-speed metal jet.

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.

Materials 2022, 15, x FOR PEER REVIEW 12 of 25
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 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  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 t = 32 µs, 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 t = 52 µs 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 t = 52 µs, the maximum density, pressure, and temperature of the jet were 15, 342 kg/m 3 , 610 MPa, and 1260 K, 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 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  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. 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 linearshaped charge was significantly shorter than that of the conical charge owing to the twodimensional 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 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 = 5.6 μs. 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 = 7 μs 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  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 twodimensional 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 t = 5.6 µs. 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 t = 7 µs 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.

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 suc-

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 suc-

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.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. Governing Equations
Appendix A. 1

. Eulerian Elasto-Plastic Solid Model
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 ρ is the density; u and v are the xand y-velocity, respectively; E is the specific total energy; p is the pressure; and s ij is the deviatoric stress tensor. The source terms in the vector S are given as follows: To calculate the temperature change caused by plastic deformation, auxiliary equations for effective plastic strain and temperature are required. As shockwave propagation and the plastic deformation process of a shaped charge occur within a short time, the temperature change in a solid could be mainly attributed to plastic deformation, and the other effects, such as heat transfer, are negligible. The resulting auxiliary equations are given as The effective plastic strain rate, . p , in the aforementioned equation can be expressed as where the consistency parameters Λ and ξ can be expressed as where ∆t is the time step, s ij,tr is the trial deviatoric stress, σ 0 Y is the yield stress at a previous time step, and h is the hardening coefficient. The stress power caused by the plastic work, . W p , is computed as where s e is the von Mises effective stress. The Taylor-Quinney (T-Q) parameter, β, which indicates the fraction of mechanical power converted to thermal power, is taken as Experimental results of Soares and Hokka [35] show that the T-Q parameter of some materials, such as Ti, Fe, Cu, and Sn, have values of approximately 0.8-0.9 under sufficient plastic strain. The stress tensor, σ ij , can be divided into a dilatational term, −pδ ij , and a deviatoric term, s ij , as where p is obtained using the Mie-Grüneisen equation-of-state (EOS) and s ij is calculated using the plastic flow theory.

Appendix A.2. Material Model
The Johnson-Cook model [20] is expressed as

≤
. p 0 and T 0 ≤ T < T m . T 0 and T m are the reference room and melting temperatures of materials, respectively. Compared to the Prandtl-Ruess model [36], which considers only the rate-independent strain hardening effect (the first term), the Johnson-Cook model considers the rate-dependent strain hardening effect (the second term) and thermal softening effect (the third term).

Appendix A.3. Radial Return Mapping Algorithm
The following yield function f s ij , σ Y was introduced to determine whether a solid is in a pure elastic ( f ≤ 0) or plastic ( f > 0) domain: To consider the plastic deformation of elasto-plastic solids, a radial return mapping algorithm [21,22] was applied. This algorithm can be divided into two steps. In the first step, the effective stress is calculated by assuming that the solid is a perfectly elastic material. In the second step, if the calculated effective stress is less than the yield stress, the effective stress is used; otherwise, it is corrected back to the yield surface using the stress. If the stress calculated in the first step is denoted as the "trial" stress, s ij,tr , and the stress corrected in the second step is denoted as the "corrected" stress, s ij,cor , the resulting stress can be obtained using The dilatational stress (pressure) of solids was modeled using the Mie-Grüneisen [37,38] EOS, which can be expressed as To prevent the square of the sound speed from becoming negative, Equation (A26) can be rewritten as follows: where p 0 is the reference pressure, e 0 is the reference internal energy, ρ 0 is the density of unstressed material, c 0 is the bulk sound speed at zero pressure, Γ 0 is the Grüneisen parameter, and s is the coefficient that relates the shock speed to the particle velocity. The sound speed of compressible solids, c, is given as Appendix A. 5

. Programmed Burn Model
The programmed burn model based on the theory of detonation shock dynamics assumes that a detonation wave propagates at a constant velocity after the detonation. This constant velocity is known as the detonation velocity, and its value is determined experimentally according to the type of explosive. As the detonation wave propagates, the explosive power is burned and turned into an explosive gas, thus releasing its internal energy. In the programmed burn model, the internal energy of the explosive was calculated by introducing the Heaviside function located at the front of the detonation wave, which can be expressed as e = e 0 (1 − H) + e product p F , ρ H where H is the Heaviside function and e 0 is the initial internal energy. The internal energy of the explosive product (gas), e product , was modeled using the Jones-Wilkins-Lee (JWL) EOS. The programmed burn model introduces a combustion function, F, which has a value of 0 when no explosive is burned and a value of 1 when completely combusted, and this value varies smoothly for several computational cells around the front of the detonation wave. The combustion function used in this study is where n b = 2 or 3, and ρ CJ is the density at the Chapman-Jouguet (C-J) state. The burn interval, ∆t b , is given as where m b = 3-6 and V d is the detonation velocity. The JWL EOS was used to model the thermodynamic state of the explosive gas, and can be expressed as where θ = ρ/ρ 0 , ρ 0 is the reference density at room temperature, A and B. are the linear coefficients, R 1 and R 2 are the nonlinear coefficients, and ω is the unitless coefficient.