Research on the Crushing Process of PELE Casing Material Based on the Crack-Softening Algorithm and Stochastic Failure Algorithm

In order to more realistically reflect the penetrating and crushing process of a PELE (Penetration with Enhanced Lateral Efficiency) projectile, the stochastic failure algorithm and crack-softening algorithm were added to the corresponding material in this paper. According to the theoretical analysis of the two algorithms, the material failure parameters (stochastic constant γ, fracture energy Gf, and tensile strength σT) were determined. Then, four sets of simulation conditions ((a) no crack softening, (b) no stochastic failure, (c) no crack softening and no stochastic failure, and (d) crack softening and stochastic failure) were designed to qualitatively describe the influences of the failure algorithms, which were simulated by the finite element analysis software AUTODYN. The qualitative comparison results indicate that the simulation results after adding the two algorithms were closer to the actual situation. Finally, ten groups of simulation conditions were designed to quantitatively analyze the coincidence degree between the simulation results and the experimental results by means of two parameters: the residual velocity of the projectile and the maximum radial velocity of fragments. The results show that the simulation results coincide well with the experimental results and the errors were small. Therefore, the ideas proposed in this paper are scientific, and the conclusions obtained can provide guidance for engineering research.


Introduction
The PELE (Penetration with Enhanced Lateral Efficiency) projectile [1,2] is a new type of armor-piercing warhead that can transform part of the axial kinetic energy into radial kinetic energy by using the differences in material properties between the projectile casing and the inner core. The PELE projectile was originally developed jointly by the French-German Research Institute Saint Louis (ISL), the Diehl Munitionssysteme, and the GEKE Technology [1][2][3][4][5]. This new type of projectile is mainly composed of the outer casing and an inner core, and its structural diagram is shown in Figure 1. The outer casing is generally made of a dense metal, such as tungsten alloy and steel, and the inner core is generally made of a material with a relatively low density such as aluminum and plastic. Compared with the conventional armor-piercing projectiles, a PELE projectile can form a larger transverse reaming during penetration and produce fragment effects after perforating the target plate (as shown in Figure 2), which makes it significantly better than the conventional armor-piercing projectiles, especially in terms of damage efficiency. In order to investigate the damage efficiency and the influencing factors of PELE projectiles, a lot of research has been carried out from the following three aspects: experiment, theoretical analysis, and numerical simulation. Rheinmetall GmbH [3] conducted a series of tests of large-caliber PELE projectiles on reinforced concrete targets, brick walls, and sandbag walls, and the reaming of target plates under different conditions were measured. Paulus et al. [5] reported a large number of experimental results of the PELE projectile impacting the target plate, which were realized with the powder gun and the light gas gun. In the experiments, the axial residual velocity of the projectile and the radial velocity of fragments were measured by using flash X-ray photography. Paulus' test results are the most representative at present, and this paper will also take his test results as a reference for comparison. Zhu et al. [6,7] used the ø12.7 mm PELE to impact a 2 mm thick armored steel target plate, and the post-effect target was a 1 mm thick aluminum plate. The effective dispersion area of fragments on the post-effect targets and the number of fragments were measured. Jiang et al. [8] used the ø30 mm PELE to impact the spaced target plates, which were composed of three 3 mm thick A3 steel and a 15 mm thick armored steel. Tu et al. [9] carried out experiments on the impact of a ø30 mm PELE on different thickness target plates, where the dispersion of fragments and the residual velocity of projectiles were measured. To sum up, the theoretical analysis model of a PELE projectile usually includes two parts: the axial penetration of the projectile and the fragment radial dispersion. The axial penetration process of a projectile is basically described by the plug model, and there is little difference between the various theoretical models. Therefore, many scholars mostly focused on the fragment dispersion models. Paulus [5], Zhu [6], Du [10], and Fan [11] simplified the crushing process of a PELE projectile, and the calculation methods of the maximum radial velocity of fragments were put forward. Based on Mott fragmentation theory, Verreault et al. [12] proposed a model to describe the axial crushing of a PELE projectile and the fragment dispersion. Although the experiment can give the true crushing state of PELE casing material, the test cost and period are longer, and the intermediate crushing process of PELE casing material cannot be given. Therefore, many scholars have studied the fragment effects of PELE casing material by numerical simulation. Paulus et al. [2,5] used LS-DYNA-3D code to simulate the crushing process of PELE casing material under different working conditions and compared it with the theoretical model. Wang et al. [13] used the three-dimensional SPH (Smoothed Particle Hydrodynamics) dynamics algorithm to simulate the fragment effects of PELE casing material at a normal projectile velocity, and the effects of target thickness, target material, and projectile velocity on the projectile crushing were compared. Verreault In order to investigate the damage efficiency and the influencing factors of PELE projectiles, a lot of research has been carried out from the following three aspects: experiment, theoretical analysis, and numerical simulation. Rheinmetall GmbH [3] conducted a series of tests of large-caliber PELE projectiles on reinforced concrete targets, brick walls, and sandbag walls, and the reaming of target plates under different conditions were measured. Paulus et al. [5] reported a large number of experimental results of the PELE projectile impacting the target plate, which were realized with the powder gun and the light gas gun. In the experiments, the axial residual velocity of the projectile and the radial velocity of fragments were measured by using flash X-ray photography. Paulus' test results are the most representative at present, and this paper will also take his test results as a reference for comparison. Zhu et al. [6,7] used the ø12.7 mm PELE to impact a 2 mm thick armored steel target plate, and the post-effect target was a 1 mm thick aluminum plate. The effective dispersion area of fragments on the post-effect targets and the number of fragments were measured. Jiang et al. [8] used the ø30 mm PELE to impact the spaced target plates, which were composed of three 3 mm thick A3 steel and a 15 mm thick armored steel. Tu et al. [9] carried out experiments on the impact of a ø30 mm PELE on different thickness target plates, where the dispersion of fragments and the residual velocity of projectiles were measured. To sum up, the theoretical analysis model of a PELE projectile usually includes two parts: the axial penetration of the projectile and the fragment radial dispersion. The axial penetration process of a projectile is basically described by the plug model, and there is little difference between the various theoretical models. Therefore, many scholars mostly focused on the fragment dispersion models. Paulus [5], Zhu [6], Du [10], and Fan [11] simplified the crushing process of a PELE projectile, and the calculation methods of the maximum radial velocity of fragments were put forward. Based on Mott fragmentation theory, Verreault et al. [12] proposed a model to describe the axial crushing of a PELE projectile and the fragment dispersion. Although the experiment can give the true crushing state of PELE casing material, the test cost and period are longer, and the intermediate crushing process of PELE casing material cannot be given. Therefore, many scholars have studied the fragment effects of PELE casing material by numerical simulation. Paulus et al. [2,5] used LS-DYNA-3D code to simulate the crushing process of PELE casing material under different working conditions and compared it with the theoretical model. Wang et al. [13] used the three-dimensional SPH (Smoothed Particle Hydrodynamics) dynamics algorithm to simulate the fragment effects of PELE casing material at a normal projectile velocity, and the effects of target thickness, target material, and projectile velocity on the projectile crushing were compared. Verreault In order to investigate the damage efficiency and the influencing factors of PELE projectiles, a lot of research has been carried out from the following three aspects: experiment, theoretical analysis, and numerical simulation. Rheinmetall GmbH [3] conducted a series of tests of large-caliber PELE projectiles on reinforced concrete targets, brick walls, and sandbag walls, and the reaming of target plates under different conditions were measured. Paulus et al. [5] reported a large number of experimental results of the PELE projectile impacting the target plate, which were realized with the powder gun and the light gas gun. In the experiments, the axial residual velocity of the projectile and the radial velocity of fragments were measured by using flash X-ray photography. Paulus' test results are the most representative at present, and this paper will also take his test results as a reference for comparison. Zhu et al. [6,7] used the ø12.7 mm PELE to impact a 2 mm thick armored steel target plate, and the post-effect target was a 1 mm thick aluminum plate. The effective dispersion area of fragments on the post-effect targets and the number of fragments were measured. Jiang et al. [8] used the ø30 mm PELE to impact the spaced target plates, which were composed of three 3 mm thick A3 steel and a 15 mm thick armored steel. Tu et al. [9] carried out experiments on the impact of a ø30 mm PELE on different thickness target plates, where the dispersion of fragments and the residual velocity of projectiles were measured. To sum up, the theoretical analysis model of a PELE projectile usually includes two parts: the axial penetration of the projectile and the fragment radial dispersion. The axial penetration process of a projectile is basically described by the plug model, and there is little difference between the various theoretical models. Therefore, many scholars mostly focused on the fragment dispersion models. Paulus [5], Zhu [6], Du [10], and Fan [11] simplified the crushing process of a PELE projectile, and the calculation methods of the maximum radial velocity of fragments were put forward. Based on Mott fragmentation theory, Verreault et al. [12] proposed a model to describe the axial crushing of a PELE projectile and the fragment dispersion. Although the experiment can give the true crushing state of PELE casing material, the test cost and period are longer, and the intermediate crushing process of PELE casing material cannot be given. Therefore, many scholars have studied the fragment effects of PELE casing material by numerical simulation. Paulus et al. [2,5] used LS-DYNA-3D code to simulate the crushing process of PELE casing material under different working conditions and compared it with the theoretical model. Wang et al. [13] used the three-dimensional SPH (Smoothed Particle Hydrodynamics) dynamics algorithm to simulate the fragment effects of PELE casing material at a normal projectile velocity, and the effects of target thickness, target material, and projectile velocity on the projectile crushing were compared. Verreault [14] carried out numerical simulations by using the AUTODYN software to verify his proposed model for describing the fragment effects of PELE casing material.
This paper aims to study the crushing process of PELE casing material from the aspect of numerical simulations. As seen from the extensive literature research, most existing numerical simulations usually adopt the relatively simple failure model to characterize the penetration and crushing process of PELE casing material, which leads to many problems such as it cannot reflect the random failure of materials. Therefore, this paper hopes to make the numerical simulation results closer to the actual situation by adding the stochastic failure algorithm and the crack-softening algorithm to the material. The test results in Reference [5] were used as the reference standard to verify the scientific and rationality of the numerical simulation method proposed in this paper by means of qualitative analysis and quantitative analysis.

Basic Theory of Mott Ring
A Mott ring [15] is a one-dimensional structure with a uniform expansion. The material will break under the action of the circumferential tensile stress, and the location of generated crack is random, as shown in Figure 3. Once the cracks are generated, the unloading waves immediately propagate to both sides, resulting in stress unloading. The unloading wave is commonly known as a Mott wave. When the fracture is completed, the unloading area formed by the adjacent cracks becomes fragmented. In the Mott ring, the material is regarded as having rigid-ideal plasticity; that is, the wave front material is plastic and the post wave material is rigid. The tensile stress is the constant yield stress Y, and the strain rate . ε = V/r is regarded as a constant.  [14] carried out numerical simulations by using the AUTODYN software to verify his proposed model for describing the fragment effects of PELE casing material. This paper aims to study the crushing process of PELE casing material from the aspect of numerical simulations. As seen from the extensive literature research, most existing numerical simulations usually adopt the relatively simple failure model to characterize the penetration and crushing process of PELE casing material, which leads to many problems such as it cannot reflect the random failure of materials. Therefore, this paper hopes to make the numerical simulation results closer to the actual situation by adding the stochastic failure algorithm and the crack-softening algorithm to the material. The test results in Reference [5] were used as the reference standard to verify the scientific and rationality of the numerical simulation method proposed in this paper by means of qualitative analysis and quantitative analysis.

Basic Theory of Mott Ring
A Mott ring [15] is a one-dimensional structure with a uniform expansion. The material will break under the action of the circumferential tensile stress, and the location of generated crack is random, as shown in Figure 3. Once the cracks are generated, the unloading waves immediately propagate to both sides, resulting in stress unloading. The unloading wave is commonly known as a Mott wave. When the fracture is completed, the unloading area formed by the adjacent cracks becomes fragmented. In the Mott ring, the material is regarded as having rigid-ideal plasticity; that is, the wave front material is plastic and the post wave material is rigid. The tensile stress is the constant yield stress Y, and the strain rate ε =  V r is regarded as a constant. For the generation of a single crack in the Mott ring, Mott took a simplified approach. It is assumed that there is no unloading process in the material where the crack is generated, and the tensile stress is always zero. Therefore, the crack generation process and its energy loss can be ignored. Under the mechanical assumption of the Mott ring, the circumferential velocity distribution of the material is shown in Figure 4. For the generation of a single crack in the Mott ring, Mott took a simplified approach. It is assumed that there is no unloading process in the material where the crack is generated, and the tensile stress is always zero. Therefore, the crack generation process and its energy loss can be ignored. Under the mechanical assumption of the Mott ring, the circumferential velocity distribution of the material is shown in Figure 4.  In Figure 4, the material position is h, the Mott wave position is x, and the circumferential velocity is u. The crack is generated at t = 0 and h = 0, and then the Mott waves propagate in both positive and negative directions. Since the material behind the wave is rigid, the circumferential velocity of the material is equal everywhere after the Mott wave, and the circumferential velocity u(t) can be expressed as: where h0 is the investigation distance. Within this distance, the total momentum can be expressed as follows: ( ) ρε ρε ρε Within the investigation distance, the change rate of total momentum over time is determined by the load. Since the material at the crack generation (h = 0) is not affected by the load, and the load at h = h0 is the constant yield stress Y, there is a momentum conservation relationship: By integrating Equation (3), the expression of a Mott wave position changing with time x(t) is obtained as follows:

Grady One-Dimensional Fracture Theory and the Crack-Softening Algorithm
Grady and Kipp [15][16][17] thought that the energy consumed by the crack generation cannot be ignored in some cases. Therefore, the following assumptions were proposed: the stress at the crack is not immediately unloaded, the tensile stress is gradually reduced from the yield stress σ = Y to σ = 0, and the unloading curve of the material is linear, as shown in Figure 5. In Figure 4, the material position is h, the Mott wave position is x, and the circumferential velocity is u. The crack is generated at t = 0 and h = 0, and then the Mott waves propagate in both positive and negative directions. Since the material behind the wave is rigid, the circumferential velocity of the material is equal everywhere after the Mott wave, and the circumferential velocity u(t) can be expressed as: where h 0 is the investigation distance. Within this distance, the total momentum can be expressed as follows: Within the investigation distance, the change rate of total momentum over time is determined by the load. Since the material at the crack generation (h = 0) is not affected by the load, and the load at h = h 0 is the constant yield stress Y, there is a momentum conservation relationship: By integrating Equation (3), the expression of a Mott wave position changing with time x(t) is obtained as follows:

Grady One-Dimensional Fracture Theory and the Crack-Softening Algorithm
Grady and Kipp [15][16][17] thought that the energy consumed by the crack generation cannot be ignored in some cases. Therefore, the following assumptions were proposed: the stress at the crack is not immediately unloaded, the tensile stress is gradually reduced from the yield stress σ = Y to σ = 0, and the unloading curve of the material is linear, as shown in Figure 5. In Figure 5, y(t) is the crack-opening displacement, σ(0) is the tensile stress of the material where the crack occurs (h = 0), and Γ is the fracture energy. When y(t) reaches the critical value yc, the unloading process is considered complete. Since σ(0) is not zero, the momentum conservation relationship of a Mott ring can be rewritten as: In the Grady one-dimensional fracture model, the material still satisfies the rigid-ideal plasticity assumption, and there is a relationship: Thus, the expression of the unloading length x(t) can be obtained as follows: The crack-softening algorithm refers to the Grady one-dimensional fracture model, and its action process can be described as follows: when the grid reaches the tensile failure condition, it must undergo an unloading process to achieve complete failure; the grid stress-strain curve is linear, and the stress cannot exceed a certain critical value during the unloading process; when the grid is no longer subjected to the tensile stress, the unloading process will terminate. The action process of the crack-softening algorithm is shown in Figure 6.  In Figure 5, y(t) is the crack-opening displacement, σ(0) is the tensile stress of the material where the crack occurs (h = 0), and Γ is the fracture energy. When y(t) reaches the critical value y c , the unloading process is considered complete. Since σ(0) is not zero, the momentum conservation relationship of a Mott ring can be rewritten as: In the Grady one-dimensional fracture model, the material still satisfies the rigid-ideal plasticity assumption, and there is a relationship: dy/dt = .
εx. Thus, the expression of the unloading length x(t) can be obtained as follows: The crack-softening algorithm refers to the Grady one-dimensional fracture model, and its action process can be described as follows: when the grid reaches the tensile failure condition, it must undergo an unloading process to achieve complete failure; the grid stress-strain curve is linear, and the stress cannot exceed a certain critical value during the unloading process; when the grid is no longer subjected to the tensile stress, the unloading process will terminate. The action process of the crack-softening algorithm is shown in Figure 6. In Figure 5, y(t) is the crack-opening displacement, σ(0) is the tensile stress of the material where the crack occurs (h = 0), and Γ is the fracture energy. When y(t) reaches the critical value yc, the unloading process is considered complete. Since σ(0) is not zero, the momentum conservation relationship of a Mott ring can be rewritten as: In the Grady one-dimensional fracture model, the material still satisfies the rigid-ideal plasticity assumption, and there is a relationship: Thus, the expression of the unloading length x(t) can be obtained as follows: The crack-softening algorithm refers to the Grady one-dimensional fracture model, and its action process can be described as follows: when the grid reaches the tensile failure condition, it must undergo an unloading process to achieve complete failure; the grid stress-strain curve is linear, and the stress cannot exceed a certain critical value during the unloading process; when the grid is no longer subjected to the tensile stress, the unloading process will terminate. The action process of the crack-softening algorithm is shown in Figure 6.  In Figure 6, σ fail represents the grid-failure stress; ε cr represents the crack strain, which is the strain generated by the grid during the unloading process; L is the dimension along the grid stretching direction; and G f is the fracture energy, corresponding to the above Γ. The unloading characteristics of the grids are determined by G f . When ε cr reaches the critical value ε U , the unloading process is completed. Thus, Figure 6 can be expressed as an analytical expression: The critical crack strain ε U can be expressed as: During the unloading process, the maximum tensile stress σ max that the grid can withstand is expressed as: where Dam is the damage degree of the grid. Here, the following definitions are made: Dam = 0: it represents the grid does not reach the failure condition. Dam = 1: it represents the grid has gone through the unloading process and completely invalid.
In addition, when the grid is completely invalid, it will no longer be subjected to the tensile shear stress.
The main control parameter G f of the crack-softening algorithm can be calculated from the value of K f (dynamic fracture toughness) [18]. The relationship between the two parameters is: where E is the Young's modulus of material, and K f can be determined using the spallation test. For different materials, such as steels with different chemical compositions or heat treatment processes, the value of K f usually differs greatly [16].

Grady Spallation Theory and Spall Strength of Materials
In the 1980s, Grady [19] conducted extensive research on the spall strength of condensed matter. He considered that when the condensed matter impacts the target at high speed, fracture will occur inside the material due to the tensile stress exceeding its tensile strength, and this phenomenon is called spallation. According to the microscopic mechanism of material fracture, the spallation is usually divided into brittle spallation and ductile spallation. The Grady spallation model considers that the material undergoes elastic volume deformation under linear tensile stress, and the strain rate is constant. The tensile stress of a material P can be expressed as: where ρ is the material density, c 0 is the material bulk velocity, . ε is the volumetric strain rate, and t is the material tensile time. The above situation is illustrated in Figure 7.  Under the action of tensile stress, the elastic volume and kinetic energy of a material gradually increases, and the material will crack when the total energy is greater than the energy required for crushing. The relationship between energy and time in the brittle spallation model is shown in Figure  8, and the energy condition for brittle spallation is as follows: where P is the tensile stress, ρ is the material density, c0 is the material bulk velocity, ε is the volumetric strain rate, and s is the fragment size. Kc is the static fracture toughness of material, which is close to the value of Kf (dynamic fracture toughness) mentioned above and can be interchanged [15,20]. The first term on the left side of Inequality (12) is the elastic volume energy U, the second term is the kinetic energy T, and the right side of Inequality (12) represents the energy required for brittle spallation Γ. Since there is the following relationship: T ≤ U/15, T can be ignored. In addition, the fragment size s must satisfy the following inequality: Obviously, when the two Inequalities (12) and (13) take the equal sign at the same time, the energy required for spallation is the smallest, as shown in Figure 8. The spall strength Ps can be solved using the above two equations:  Under the action of tensile stress, the elastic volume and kinetic energy of a material gradually increases, and the material will crack when the total energy is greater than the energy required for crushing. The relationship between energy and time in the brittle spallation model is shown in Figure 8, and the energy condition for brittle spallation is as follows: where P is the tensile stress, ρ is the material density, c 0 is the material bulk velocity, . ε is the volumetric strain rate, and s is the fragment size. K c is the static fracture toughness of material, which is close to the value of K f (dynamic fracture toughness) mentioned above and can be interchanged [15,20].
The first term on the left side of Inequality (12) is the elastic volume energy U, the second term is the kinetic energy T, and the right side of Inequality (12) represents the energy required for brittle spallation Γ. Since there is the following relationship: T ≤ U/15, T can be ignored. In addition, the fragment size s must satisfy the following inequality: Obviously, when the two Inequalities (12) and (13) take the equal sign at the same time, the energy required for spallation is the smallest, as shown in Figure 8. The spall strength P s can be solved using the above two equations: Under the action of tensile stress, the elastic volume and kinetic energy of a material gradually increases, and the material will crack when the total energy is greater than the energy required for crushing. The relationship between energy and time in the brittle spallation model is shown in Figure  8, and the energy condition for brittle spallation is as follows: where P is the tensile stress, ρ is the material density, c0 is the material bulk velocity, ε is the volumetric strain rate, and s is the fragment size. Kc is the static fracture toughness of material, which is close to the value of Kf (dynamic fracture toughness) mentioned above and can be interchanged [15,20]. The first term on the left side of Inequality (12) is the elastic volume energy U, the second term is the kinetic energy T, and the right side of Inequality (12) represents the energy required for brittle spallation Γ. Since there is the following relationship: T ≤ U/15, T can be ignored. In addition, the fragment size s must satisfy the following inequality: Obviously, when the two Inequalities (12) and (13) take the equal sign at the same time, the energy required for spallation is the smallest, as shown in Figure 8. The spall strength Ps can be solved using the above two equations:  In the numerical simulation, the material has various tensile failure models, such as the hydrostatic pressure and principal stress failure models, taking a given hydrostatic pressure and principal stress as the material tensile failure criterion, respectively. Reference [21] shows that when the principal stress failure model is used, the tensile strength of a material can be estimated by the Grady spallation theory. In other words, when the actual tensile strength of a material is unknown, the principal tensile failure stress σ T in the numerical simulation can be replaced by the theoretical spall strength P s . Therefore, in the following numerical simulation, the principal stress failure model will be added to the metal core material of a PELE projectile.

Mott Fracture Probability Density Function and Stochastic Failure Algorithm
In order to solve the problem of multiple cracks, Mott proposed to use the fracture probability density function λ(ε) to control the occurrence probability of cracks, and these functions are called the risk function and the conditional failure function. The expression λ(ε)dεdl is the probability that a crack will occur in the element with a length dl and strain ε when the strain increases by dε. In 1943, Mott [15] proposed three kinds of fracture probability functions λ(ε), which are: where n and σ are the Weibull distribution constants, and C and γ are the Gumbel distribution constants.
In the above three forms of functions, the exponential expression is the most commonly used. It is assumed that the sample is composed of N 0 unit length segments, each segment is stretched independently at the same strain rate, and the number of surviving segments is N. As a result, there is the following relationship: The above formula can also be rewritten as an expression in the integral form: Hence, the cumulative probability distribution of the whole sample before or at the strain ε can be obtained: By differentiating the above formula, the complementary cumulative probability density can also be obtained: The fracture frequency function λ(ε) is different for different materials, resulting in different expressions for the cumulative probability distribution function F(ε) and the cumulative probability density function f (ε).
The stochastic failure algorithm refers to the Mott stochastic failure theory, which discretizes the Mott theory. Each grid is automatically assigned a stochastic coefficient, and the size of the coefficient is the ratio of the true failure strain of grid to the rupture strain of the material given by the user. The relationship between the grid number and the stochastic factor is determined by the given stochastic constant γ.
In the AUTODYN simulation software, the stochastic failure algorithm is embedded in order to more realistically describe the fracture and failure of a material [18]. The expression of cumulative failure probability is: where ε N is the normalized strain, i.e., ε N = ε/ε c ; ε c is the fracture strain; and γ and C are the material constants. For Equation (22), there is a very important premise condition: when the material strain reaches the specified fracture strain, the failure probability of the grid is 50%. That is, there exists the following relationship: Thus, the relationship between γ and C can be obtained: Substituting Equation (24) into Equation (22), we can get: Similarly, by differentiating the above formula, the cumulative failure probability density expression can be obtained: Since the γ-C relationship is certain, it is only necessary to give a specific γ value to determine the stochastic dispersion of grids in actual operation. However, for materials with unknown fracture properties, the γ value can only be selected empirically. Under the condition that the given fracture strain is ε c = 0.035, the cumulative failure probability density f (ε) and the cumulative fracture probability P(ε) under different γ values are shown in Figures 9 and 10, respectively. relationship between the grid number and the stochastic factor is determined by the given stochastic constant γ.
In the AUTODYN simulation software, the stochastic failure algorithm is embedded in order to more realistically describe the fracture and failure of a material [18]. The expression of cumulative failure probability is: (22) where ε N is the normalized strain, i.e., ε ε ε = N c ; ε c is the fracture strain; and γ and C are the material constants. For Equation (22), there is a very important premise condition: when the material strain reaches the specified fracture strain, the failure probability of the grid is 50%. That is, there exists the following relationship: Thus, the relationship between γ and C can be obtained: Substituting Equation (24) into Equation (22), we can get: Similarly, by differentiating the above formula, the cumulative failure probability density expression can be obtained: Since the γ-C relationship is certain, it is only necessary to give a specific γ value to determine the stochastic dispersion of grids in actual operation. However, for materials with unknown fracture properties, the γ value can only be selected empirically. Under the condition that the given fracture strain is ε =    As can be seen from Figure 9, the material fracture probability density curve f(ε) becomes steeper as the γ value increases. The failure strain error is defined as the error between the strain at the center of the fracture probability density distribution f(ε) and the given fracture strain ε c , and the failure strain error decreases as the γ value increases. From this point of view, it is desirable that the γ value should be as large as possible. However, it can also be seen from Figure 10 that if the γ value is too large, the failure strain of the grid will become more concentrated, that is to say, the material crushing performance will tend to be uniform. In this respect, it is not desirable that the γ value is too large. Therefore, for a specific material, there is an optimal range of the γ value. In the following numerical simulation, the optimal γ value will be determined by trial calculation.

Numerical Simulation
In this paper, the nonlinear dynamics software AUTODYN (Century Dynamics, Fort Worth, TX, USA) [22] was used to simulate the crushing process of PELE casing material. AUTODYN is an explicit finite element analysis program, which is used to solve the highly nonlinear dynamic problems of solids, fluids, gases, and their interactions. More importantly, the software has a unique stochastic failure model, which can reflect the randomness and non-uniformity exhibited by the material.

Finite Element Model
The entire model was divided into three parts: the outer casing of the PELE projectile, the inner core of the PELE projectile, and the target plate. The PELE projectile length was 50 mm, the outer and inner diameters of the projectile were 10 mm and 6 mm, respectively, and the thickness of the projectile rear was 5 mm. The inner core was a ø6 mm× 45 mm cylinder. The length and width of the target plate were both 120 mm, and the thickness was selected according to different working conditions. The unstructured hexahedral mesh was used in both the projectile and target plate, and the average size of the grid was 0.25 mm. The whole model was established using HyperMesh software (Altair Engineering, Detroit, MI, USA), then the finite element model was imported into AUTODYN for the solution and the Lagrange algorithm was adopted. In order to reduce the number of grids and improve the computational efficiency, the meshing method of the target plate adopted the variable-step size method. The grids were denser in the center area, which was about two times the diameter of projectile, and the minimum grid size was 0.25 mm. The transmissive boundary condition was applied to the edge of the target plate. In the case of considering the vertical impact, the model adopted a 1/4 simplification; when considering the oblique penetration, the model adopted a 1/2 simplification; when considering the projectile rotation, the model was not simplified. The schematic diagram of the finite element model is shown in Figure 11. As can be seen from Figure 9, the material fracture probability density curve f (ε) becomes steeper as the γ value increases. The failure strain error is defined as the error between the strain at the center of the fracture probability density distribution f (ε) and the given fracture strain ε c , and the failure strain error decreases as the γ value increases. From this point of view, it is desirable that the γ value should be as large as possible. However, it can also be seen from Figure 10 that if the γ value is too large, the failure strain of the grid will become more concentrated, that is to say, the material crushing performance will tend to be uniform. In this respect, it is not desirable that the γ value is too large. Therefore, for a specific material, there is an optimal range of the γ value. In the following numerical simulation, the optimal γ value will be determined by trial calculation.

Numerical Simulation
In this paper, the nonlinear dynamics software AUTODYN (Century Dynamics, Fort Worth, TX, USA) [22] was used to simulate the crushing process of PELE casing material. AUTODYN is an explicit finite element analysis program, which is used to solve the highly nonlinear dynamic problems of solids, fluids, gases, and their interactions. More importantly, the software has a unique stochastic failure model, which can reflect the randomness and non-uniformity exhibited by the material.

Finite Element Model
The entire model was divided into three parts: the outer casing of the PELE projectile, the inner core of the PELE projectile, and the target plate. The PELE projectile length was 50 mm, the outer and inner diameters of the projectile were 10 mm and 6 mm, respectively, and the thickness of the projectile rear was 5 mm. The inner core was a ø6 mm× 45 mm cylinder. The length and width of the target plate were both 120 mm, and the thickness was selected according to different working conditions. The unstructured hexahedral mesh was used in both the projectile and target plate, and the average size of the grid was 0.25 mm. The whole model was established using HyperMesh software (Altair Engineering, Detroit, MI, USA), then the finite element model was imported into AUTODYN for the solution and the Lagrange algorithm was adopted. In order to reduce the number of grids and improve the computational efficiency, the meshing method of the target plate adopted the variable-step size method. The grids were denser in the center area, which was about two times the diameter of projectile, and the minimum grid size was 0.25 mm. The transmissive boundary condition was applied to the edge of the target plate. In the case of considering the vertical impact, the model adopted a 1/4 simplification; when considering the oblique penetration, the model adopted a 1/2 simplification; when considering the projectile rotation, the model was not simplified. The schematic diagram of the finite element model is shown in Figure 11.

Determination of Key Parameters in the Material Failure Model
In this paper, the experiments carried out by Paulus [5] were used as the reference for comparison. In the reference experiments, the projectile casing was made of tungsten alloy (D-180 K), the inner core was made of aluminum (A-G3) and polyethylene (PE), and the target plate material was made of steel (XC48) and aluminum (A-U4G), but several experimental materials were not identical in the material library of AUTODYN. In order to solve this key problem, the materials were of the same type as the experimental materials, but different parameters were found in the material library of AUTODYN, which have the same equation of state and constitutive equation as the experimental materials. Then, all material parameters in the substitute material selected in AUTODYN were replaced with the parameters of the experimental material, which could ensure the consistency of the materials in the experiment and numerical simulation.
(1) Determination of the key parameter Gf of the crack-softening algorithm of the outer casing tungsten: The fracture energy Gf was determined based on the dynamic fracture toughness Kf and the Young's modulus E. Reference [23] provided the range of the dynamic fracture toughness of the fragile tungsten alloy Kf = 3-5 MPa⋅m 1/2 , and this paper takes the intermediate value Kf = 4 MPa⋅m 1/2 . Then, according to Gf = Kf/E (see Equation (10)) and E = 360 GPa, Gf = 45 J/m 2 was obtained.
(2) Determination of the key parameter γ of the stochastic failure algorithm of the outer casing tungsten: The stochastic constant γ was determined by the failure strain error analysis and numerical trial calculations. The test condition (the inner core material was polyethylene, the target plate was 3 mm thick steel, and the projectile velocity was 936 m/s) in Reference [5] was used as the reference standards. A series of γ values were substituted into the numerical simulation, and it was found that the length of the residual projectile varied greatly when the γ value was different, and the trial result is shown in Figure 12.
As shown in Figure 12, with the increase of the γ value, the changing law of the residual projectile length was rather complicated, and the simulation results at γ = 36.5 were in good agreement with the experimental results. Therefore, the stochastic constant in this paper was identified as γ = 36.5.

Determination of Key Parameters in the Material Failure Model
In this paper, the experiments carried out by Paulus [5] were used as the reference for comparison. In the reference experiments, the projectile casing was made of tungsten alloy (D-180 K), the inner core was made of aluminum (A-G3) and polyethylene (PE), and the target plate material was made of steel (XC48) and aluminum (A-U4G), but several experimental materials were not identical in the material library of AUTODYN. In order to solve this key problem, the materials were of the same type as the experimental materials, but different parameters were found in the material library of AUTODYN, which have the same equation of state and constitutive equation as the experimental materials. Then, all material parameters in the substitute material selected in AUTODYN were replaced with the parameters of the experimental material, which could ensure the consistency of the materials in the experiment and numerical simulation.
(1) Determination of the key parameter G f of the crack-softening algorithm of the outer casing tungsten: The fracture energy G f was determined based on the dynamic fracture toughness K f and the Young's modulus E. Reference [23] provided the range of the dynamic fracture toughness of the fragile tungsten alloy K f = 3-5 MPa·m 1/2 , and this paper takes the intermediate value K f = 4 MPa·m 1/2 . Then, according to G f = K f /E (see Equation (10)) and E = 360 GPa, G f = 45 J/m 2 was obtained.
(2) Determination of the key parameter γ of the stochastic failure algorithm of the outer casing tungsten: The stochastic constant γ was determined by the failure strain error analysis and numerical trial calculations. The test condition (the inner core material was polyethylene, the target plate was 3 mm thick steel, and the projectile velocity was 936 m/s) in Reference [5] was used as the reference standards. A series of γ values were substituted into the numerical simulation, and it was found that the length of the residual projectile varied greatly when the γ value was different, and the trial result is shown in Figure 12.

Material Model and Parameters
In this paper, the selection method of the equation of state and constitutive equations for all materials are referred to References [24,25]. The equation of state for all materials adopted Mie-Grüneisen, which is denoted as the "shock" equation of state [22] in AUTODYN, and its basic theoretical relationship is as follows: where Us is the shock wave velocity, u is the post-wave particle velocity, and c0 and s are the material Hugoniot parameters. The constitutive equation of the material is elastic-ideal plastic. The constitutive law is based on the von Mises yield criterion and the steady yield stress assumption, and the shear modulus of the material is regarded as a constant. This constitutive equation is described as the "von Mises" strength model in AUTODYN, which requires a given material shear modulus G and the flow stress Y. In addition, the dynamic yield strength is typically approximately 2-3 times the static value.
The projectile casing adopted the principal stress/strain failure model. It was considered that the material would fail when the tensile principal strain reached εT or the tensile principal stress reached σT. According to Reference [5], εT = 0.035 and σT = 2.8 GPa. In addition, the stochastic failure algorithm and crack-softening algorithm were added to the projectile casing material. According to the analysis above, the stochastic coefficients were γ = 36.5 and Gf = 45 J/m 2 . The aluminum core material adopted the principal stress failure model, and the tensile strength was σT = 0.5 GPa, the polyethylene core material was not added with failure. The material of the target plate adopted the plastic strain failure model, and it was considered that the material would fail when the plastic strain reached a certain value. All materials were added with an artificial erosion algorithm (Erosion) to ensure a normal calculation. The target plate material adopted the failure erosion algorithm (Failure), which was used to delete the invalid grid. The rest of the materials use the geometric strain erosion algorithm, which As shown in Figure 12, with the increase of the γ value, the changing law of the residual projectile length was rather complicated, and the simulation results at γ = 36.5 were in good agreement with the experimental results. Therefore, the stochastic constant in this paper was identified as γ = 36.5.
(3) Determination of the key parameter of principal stress failure model of the inner core material Al-6061: The tensile strength σ T , corresponding to P s above, could be calculated by the Grady spallation theory, where its specific expression is: σ T = 3ρc 0 K 2 f . ε 1/3 . According to Reference [15], the dynamic fracture toughness was taken as K f = 27.5 MPa·m 1/2 , and the strain rate was taken as . ε ≈ 10 3 −10 4 s −1 . Then, when the material density ρ 0 = 2.65 g/cm 3 and the volumetric sound velocity c 0 = 5176 m/s were substituted into the expression, we obtained σ T ≈ 0.31-0.68 GPa. Here, we take the intermediate value σ T = 0.5 GPa.

Material Model and Parameters
In this paper, the selection method of the equation of state and constitutive equations for all materials are referred to References [24,25]. The equation of state for all materials adopted Mie-Grüneisen, which is denoted as the "shock" equation of state [22] in AUTODYN, and its basic theoretical relationship is as follows: where U s is the shock wave velocity, u is the post-wave particle velocity, and c 0 and s are the material Hugoniot parameters. The constitutive equation of the material is elastic-ideal plastic. The constitutive law is based on the von Mises yield criterion and the steady yield stress assumption, and the shear modulus of the material is regarded as a constant. This constitutive equation is described as the "von Mises" strength model in AUTODYN, which requires a given material shear modulus G and the flow stress Y. In addition, the dynamic yield strength is typically approximately 2-3 times the static value.
The projectile casing adopted the principal stress/strain failure model. It was considered that the material would fail when the tensile principal strain reached ε T or the tensile principal stress reached σ T . According to Reference [5], ε T = 0.035 and σ T = 2.8 GPa. In addition, the stochastic failure algorithm and crack-softening algorithm were added to the projectile casing material. According to the analysis above, the stochastic coefficients were γ = 36.5 and G f = 45 J/m 2 . The aluminum core material adopted the principal stress failure model, and the tensile strength was σ T = 0.5 GPa, the polyethylene core material was not added with failure. The material of the target plate adopted the plastic strain failure model, and it was considered that the material would fail when the plastic strain reached a certain value. All materials were added with an artificial erosion algorithm (Erosion) to ensure a normal calculation. The target plate material adopted the failure erosion algorithm (Failure), which was used to delete the invalid grid. The rest of the materials use the geometric strain erosion algorithm, which removed the grid whose instantaneous geometric strain (Geometric Strain) was greater than the given value. In summary, the material parameters used in this paper are summarized in Table 1.

Influence of the Crack-Softening Algorithm and Stochastic Failure Algorithm on the Simulation Results
In order to more intuitively give the influence of the crack-softening algorithm and stochastic-failure algorithm on the simulation results, four groups of simulation conditions ((a) no crack softening, (b) no stochastic failure, (c) no crack softening and no stochastic failure, and (d) crack softening and stochastic failure) were designed in this paper. Moreover, the test conditions (the inner core material was polyethylene, the target plate was 8 mm thick aluminum, and the projectile velocity was 939 m/s) in Reference [5] were taken as the reference. The simulation results corresponding to each group of conditions were compared with the experimental results. The effects of adding the crack-softening algorithm and stochastic-failure algorithm before and after are shown in Figure 13.
As shown in Figure 13, when the crack-softening algorithm (a, c) was not added, the material failure was not in accordance with the actual situation. Moreover, the failure grids were usually distorted so that more grids were eroded, and the remaining failure area of material was usually smaller. When the stochastic failure algorithm was not added (b, c), the fragment length of projectile was larger and the shape of fragments was not ideal. After adding the stochastic-failure algorithm and crack-softening algorithm (d), the simulation results were more consistent with the actual conditions (e). The length and shape of fragments were ideal, and the projectile failure region was relatively uniform.
In summary, adding the crack-softening algorithm was beneficial to control the distortion of the failure grids, and adding the stochastic failure algorithm was beneficial to control the length and shape of fragments. Therefore, the crack-softening algorithm and stochastic-failure algorithm were conducive to simulating the real damage situation, and they were added in the subsequent simulation.

Quantitative Comparison of the Simulation Results and Experimental Results
Based on the analysis in Section 4.1, it can be qualitatively seen that the simulation results after adding the two failure algorithms were in good agreement with the experimental results. However, this chapter hopes to quantitatively analyze the coincidence degree between the simulation results and the experimental results by designing several sets of simulation conditions. The entire qualitative comparison process was based on two parameters: the residual velocity of the projectile and the maximum radial velocity of fragments. Referring to the experimental results in Reference [5], the specific setting of the simulation conditions in this paper were designed as shown in Table 2. The numerical simulation results were compared with the flash X-ray pictures [5] as shown in Figure 14. The maximum radial velocity of fragments and the residual velocity of the projectile obtained by the simulation are shown in Table 3.

Quantitative Comparison of the Simulation Results and Experimental Results
Based on the analysis in Section 4.1, it can be qualitatively seen that the simulation results after adding the two failure algorithms were in good agreement with the experimental results. However, this chapter hopes to quantitatively analyze the coincidence degree between the simulation results and the experimental results by designing several sets of simulation conditions. The entire qualitative comparison process was based on two parameters: the residual velocity of the projectile and the maximum radial velocity of fragments. Referring to the experimental results in Reference [5], the specific setting of the simulation conditions in this paper were designed as shown in Table 2. The numerical simulation results were compared with the flash X-ray pictures [5] as shown in Figure 14. The maximum radial velocity of fragments and the residual velocity of the projectile obtained by the simulation are shown in Table 3.  It can be seen intuitively from Figure 14 and Table 3 that the simulation results were basically consistent with the experimental results and the errors were smaller. In addition, it was also found that the PELE projectile crushing performance was also greatly different due to the inner core material, and target plate material and thickness, which was embodied in the following aspects: (1) Comparing the following four sets of working conditions (#3-#9; #4-#10; #5-#7; #6-#8), the following conclusions can be obtained intuitively from Figure 14 and Table 3. When the core material was aluminum, the residual length of the projectile was larger, the size of fragments was smaller, and the radial velocity of fragments was higher than when the core material was polyethylene. (2) Comparing the following two sets of working conditions (#1-#5; #2-#6) from Table 3, for the same core material and the target plate thickness, when the target material changed from aluminum to steel, the residual velocity of projectile decreases and the radial velocity of fragments increased. Comparing the following two sets of working conditions (#1-#3; #2-#4) from Table 3, for the same core material and the target plate material, when the target plate thickness increased, the residual velocity of the projectile decreased and the radial velocity of the fragments increased.
Taking the #6 simulation condition as an example, the fragmentation process of the PELE projectile impacting the target plate is shown in Figure 15. In order to facilitate the observation of the crushing process of the outer casing, the core material is not shown. The core material was aluminum, the target plate was 3 mm steel, and the projectile impact velocity was 1262 m/s. As shown in Figure 15, at the moment of impacting the target plate, two shock waves propagating respectively toward the rear of the projectile and the back of the target plate were generated at the interface, and the plastic deformation of the projectile head caused the failure of the grids (a, e). While the plastic deformation region propagated toward the rear of projectile, the failure region of casing gradually expanded backward (b, f) under the combined action of the resistance of the target plate and the radial force of the inner core, and these failure regions were considered to be cracks. After the projectile perforated the target plate, the plug and the inner core continued to interact, and the casing cracks expanded further along the axial direction (c, g). When the projectile moved further forward, the fragments generated by the crack propagation gradually scattered under the radial velocity (d, h). In the numerical simulation, a series of Gauss observation points were set up on the outer casing and the target plug to record the velocity changes. The radial velocity-time history curve of the material at the Gauss point of projectile casing is shown in Figure 16a. The axial velocity-time history curves of the material at the Gauss point of the projectile casing and target plug are shown in Figure 16b. It can be seen from Figure 16a that the radial velocity of the fragments increased rapidly at ≈0-0.01 ms, because there existed a radial compression during this period. Then, some fluctuations occurred and it reached stability at about 0.05 ms, which indicates that the projectile had perforated the target plate. According to Figure 16b, the axial velocity of the projectile was gradually reduced under the resistance of the target plate. However, the axial velocity of the target plug was gradually increased under the interaction of the projectile and inner core. Finally, the two curves tend to be consistent around 0.115 ms, at which point the projectile and the target plug were no longer interacting. In addition, it was also found that the aluminum core projectile and the target plug have a shorter action time than the polyethylene core projectile. In the numerical simulation, a series of Gauss observation points were set up on the outer casing and the target plug to record the velocity changes. The radial velocity-time history curve of the material at the Gauss point of projectile casing is shown in Figure 16a. The axial velocity-time history curves of the material at the Gauss point of the projectile casing and target plug are shown in Figure 16b. In the numerical simulation, a series of Gauss observation points were set up on the outer casing and the target plug to record the velocity changes. The radial velocity-time history curve of the material at the Gauss point of projectile casing is shown in Figure 16a. The axial velocity-time history curves of the material at the Gauss point of the projectile casing and target plug are shown in Figure 16b. It can be seen from Figure 16a that the radial velocity of the fragments increased rapidly at ≈0-0.01 ms, because there existed a radial compression during this period. Then, some fluctuations occurred and it reached stability at about 0.05 ms, which indicates that the projectile had perforated the target plate. According to Figure 16b, the axial velocity of the projectile was gradually reduced under the resistance of the target plate. However, the axial velocity of the target plug was gradually increased under the interaction of the projectile and inner core. Finally, the two curves tend to be consistent around 0.115 ms, at which point the projectile and the target plug were no longer interacting. In addition, it was also found that the aluminum core projectile and the target plug have a shorter action time than the polyethylene core projectile. It can be seen from Figure 16a that the radial velocity of the fragments increased rapidly at ≈0-0.01 ms, because there existed a radial compression during this period. Then, some fluctuations occurred and it reached stability at about 0.05 ms, which indicates that the projectile had perforated the target plate. According to Figure 16b, the axial velocity of the projectile was gradually reduced under the resistance of the target plate. However, the axial velocity of the target plug was gradually increased under the interaction of the projectile and inner core. Finally, the two curves tend to be consistent around 0.115 ms, at which point the projectile and the target plug were no longer interacting. In addition, it was also found that the aluminum core projectile and the target plug have a shorter action time than the polyethylene core projectile.