Abstract

Reactive material (RM) is a new type of energetic material, which is widely used in the military technology fields such as fragmentation warheads and shaped charge warheads. Violent chemical reactions take place in the impact process of reactive materials, and how to realize the macro numerical simulation of shock-induced energy release behavior of reactive materials is one of the most urgent problems to be solved for its future military applications. In this study, a numerical simulation approach and procedure is proposed, which can simulate the shock-induced energy release behavior of reactive materials on a macro scale. Firstly, program implementation of the mechanical-thermal-chemical coupled effect model for RM is realized in the second-development interface of LS-DYNA software. Then, the adaptive simulated annealing algorithm is used to fit the chemical reaction kinetic parameters of RM using the direct ballistics test data. Finally, the simulation calculation of the fragment penetrating upon steel plate is carried out to expand the applicability of the numerical simulation approach proposed in this study. The results show that the numerical simulation approach proposed in this study can reproduce the results of the direct ballistics test more accurately, which assumes practical significance for the engineering application of reactive materials in the military field in the future.

1. Introduction

Reactive material (RM) is a special energetic material, which is inert under normal temperature and pressure, while upon strong impact it can react dramatically and release a lot of chemical energy. Reactive materials are usually prepared by hot pressing of metal/nonmetal powder (Al, Ni, W, PTFE, Fe2O3, etc.) [1]. Warhead fragments made of reactive material can release a large amount of chemical energy on the target in addition to the kinetic energy damage [2, 3]. The inner core of PELE bombs made of reactive materials can enhance the cracking of shells after penetration, resulting in both physical and chemical damage to the target [4]. The metal shape charge of armor penetration warhead made of reactive material can form a reactive jet containing chemical reactions and thus enhance the penetration effect of the jet [5, 6]. In addition, reactive materials have many other applications in antibiological war agents, underwater explosion, and so on [7].

Compared with the inert materials, the biggest advantage of reactive materials is their shock-induced energy release. To accurately describe the shock-induced energy release behavior of reactive materials, it is necessary to probe the chemical reaction mechanism and identify the chemical kinetic parameters. In theory, the chemical kinetic parameters can be obtained through differential scanning calorimetry (DSC). However, under impact loading, the temperature in reactive material rises very fast (in microseconds), accompanied by strong mechanical effects such as plastic deformation, shear, and friction between the particles [8, 9], bringing about different reaction conditions, which are completely different from the boundary conditions of DSC experiments. Taking Al/PTFE reactive material as an example, Miao [10] finds at the heating rate of 30 K/min that micron Al powder cannot react with PTFE due to the coating of oxidation film on the surface while according to [11], the oxidation film on the surface of Al particles can be destroyed by the mechanical effect of the shock wave, causing the temperature threshold of impact reaction between Al and PTFE to be only 411k. Therefore, it is necessary to find other ways to obtain reaction parameters. Xiong [12] obtained the shock-induced temperature rise and chemical reaction threshold of reactive materials through theoretical calculation employing results of direct ballistic tests; he also derived the chemical kinetic parameters by fitting test data, which proved the effectiveness of this method. Ren [13] used Xiong’s theoretical method to fit the chemical kinetic parameters of three Al/Ni reactive materials and used a self-developed constitutive program to calculate the shock-induced energy release behavior of reactive materials. The simulation results (calculated with Xiong’s method) showed a certain difference from the direct ballistic test results.

In this study, a numerical simulation scheme and procedure is proposed, which can simulate the shock-induced energy release behavior of reactive materials on a macro scale. In the first step of this study, the chemical kinetics equation of reactive materials is built into LS-DYNA in the form of a user-defined equation of state (EOS), which realizes the program implementation of the mechanical-thermal-chemical coupled effect model. Then, the adaptive simulated annealing algorithm is used to fit direct ballistics test data, which yields chemical reaction kinetic parameters of the reactive material. Finally, the simulation calculation of fragment penetrating upon steel plate is carried out, and a comparison is made between the damage effect of reactive material fragments and that of inert material fragments, which further expands the applicability of the numerical simulation method proposed in this study.

2. Mechanical-Thermal-Chemical Coupled Effect Model for Reactive Materials and Its Program Implementation

In the impact process of reactive materials, violent chemical reactions take place, which is called shock-induced energy release behavior. With complex chemical reaction chain and diversity of reaction products, the shock-induced energy release behavior is hard to be represented by accurate chemical reaction equations. For this reason, we consider that the reactive material during the partial reaction is a mixture of reactants (unreacted reactive material) and reaction products (reactive material with complete reaction). A parameter y, standing for the extent of chemical reaction, is introduced to represent the mixing of reactant substances and reaction products. The relative volume and internal energy meet the following mixing rules:where represents the relative volume and E represents the internal energy per element of the initial volume. The subscript m stands for the mixture, r for the reactants, and p for the reaction products. Meanwhile, it is assumed that the pressure P and temperature T of reactants and products are in equilibrium:

Due to the assumption of temperature equilibrium, it is very important to obtain the temperature of reactants and reaction products in the calculation process. So, we need to take temperature as an independent variable in the EOS to facilitate temperature calculation of the mixture. The JWL equation [14] in the form of pressure-relative volume-temperature (P--T) takes temperature and relative volume as independent variables. Therefore, in this study, JWL equations are adopted as EOS for reactants and products during the impact process:where Cv represents the specific heat capacity of the material and A, B, ω, R1, and R2 are related parameters of the material. The extent of chemical reaction y is controlled by chemical kinetics equation. Therefore, if the JWL equation for the reactant and the reaction product and the chemical kinetics equation of the reactive material are known, the EOS that describe the reaction behavior of the reactive material can be constructed, so as to calculate the physical quantities of the reactive material under different impact conditions.

In LS-DYNA, a second-development interface is reserved, allowing users to perform second development on material models such as constitutive equations and EOS, so as to help the users define the updating modes of physical quantities such as stress-strain, pressure, and internal energy by themselves. In this study, the JWL EOS of the reactants and products and the chemical kinetics equation is embedded into LS-DYNA in the form of user-defined EOS, and then the numerical simulation is carried out. In this way, program implementation of the mechanical-thermal-chemical coupled effect model for reactive materials is realized. In the program flow of the user-defined EOS, which considers mechanical-thermal-chemical coupled effect, updating of each physical quantity is shown in Figure 1.

The calculation process for the user-defined EOS is divided into the following seven steps:Step 1: Enter Tt, Pt, yt, and so forth into LS-DYNA main program.Step 2: Calculate the predicted temperature value, in the current time step (i.e., t +1 step).In the iterative scheme of the procedure, the temperature can be calculated through the predictor-corrector method. The prediction value in t +1 step can be calculated based on the values in the previous time step (i.e., t step) with the following equation:where ΔTs represents the temperature rise due to plastic deformation. q = 0.5(qt + qt+1) represents the artificial viscous force. J represents the partial derivative of internal energy with volume. represents the increment of the relative volume of the material in the current time step. In equation (4), some variables cannot be directly obtained in the LS-DYNA main program, so the form of equation (4) needs to be changed. In LS-DYNA, the element’s specific internal energy, E=e/_v0, represents the internal energy per element's initial relative volume, where e is the element’s internal energy and 0 is the element’s initial relative volume. In t +1 step, E is calculated as follows:σ stands for deviatoric stress. The first term of equation (5) represents the plastic work, the second term represents the volume compression work, and the last term represents the element’s specific internal energy of t step. Equation (5) is rearranged aswhere is transferred into the user-defined EOS subroutine by the dummy argument, specen. Combining equations (4), (5), and (6) yields the temperature prediction value:Step 3: Calculate the predicted pressure and other physical quantities.According to [15], an intermediate variable β=(1 y)/ can be defined to assist the calculation, where β means the reactant volume fraction in the total volume of the mixture. Based on the definition of β, the following relationship can be derived:Firstly, the pressure of reactants and products is calculated using the initial value of β. Then, the Newton iteration method is employed to iterate β:where and are calculated based on . When iteration stops and =( + )/2 is taken as the predicted pressure value for the current time step, then , , and are used to calculate , , and .Step 4: Calculate the corrected temperature, Tt+1.When , , , and are obtained, the corrected value of temperature, Tt+1can calculated as follows:where and .Step 5: Refer to the process in step 3 to update other physical quantities.Step 6: Calculate yt+1.According to the reaction kinetics model given by Arrhenius [16], the chemical reaction rate (dy/dt) follows the formwhere t is the reaction time, A0 is the preexponential factor, Ea is the apparent activation energy, T is the thermodynamic temperature, and Ru is the universal gas constant. f (y) is the reaction mechanism function, and its form is determined by the specific reaction type. The n-dimensional control reaction model proposed by Avrami–Erofeev can be employed to describe the solid reaction under rapid temperature rise [16]:where n is the reaction index. Zhang et al. [12, 17] derived the first-order differential relationship between the thermodynamic temperature, T, and the extent of reaction, y, by assuming that the reaction rate is positively correlated with time:Given the reaction critical temperature, Tcr, and the critical extent of reaction, ycr (assuming a minimum, 10−4 in this study), the differential equation of equation (13) can be used to calculate y through numerical integration. Here, the fourth-order Runge–Kutta method is adopted:whereStep 7: Feed the updated variables to LS-DYNA main program.

3. Fitting of Chemical Reaction Kinetic Parameters Using the Adaptive Simulated Annealing Algorithm

To fully model the mechanical-thermal-chemical coupled behavior of the reactive material, it is necessary to identify the chemical kinetic parameters of the material first. In the chemical kinetics equation (equation (13)), there are three parameters to be determined by experimental fitting. They are the activation energy Ea, reaction index n, and reaction critical temperature, Tcr. The fitting process of the above chemical reaction kinetic parameters is as follows in this study.

Firstly, the experimental conditions are simulated by using the solver obtained in the second section; then adaptive simulated annealing algorithm is used to adjust the parameters so that the calculated value and the experimental value of the extent of the reaction for the reactive material will match. In this way, the exact values of the above kinetic parameters can be obtained. Exploiting the existing data of the direct ballistic test (DBT), the application process of the above-mentioned parameter fitting method is described in detail below.

The direct ballistic test data and material parameters used in this section are taken from [13, 18]. References [13, 18] chose the Al/Ni reactive material. This material has good mechanical properties and is easy to prepare, so it is often used to make reactive warhead fragments. The direct ballistic test results for some Al/Ni reactive material are shown in Table 1. In the direct ballistics test setting, the front end of the sealed chamber (volume: 27 L) is made of steel sheet, and the back end is steel target plate. The projectile (cylindrical, Φ10  mm × 10  mm) made of Al/Ni reactive material is launched through a ballistic gun. The projectile first penetrates the steel sheet and the impacts on the steel target. In the process of projectile impact, a large amount of heat will be released and the air pressure in the chamber will be increased. By measuring the quasi-static pressure of the air in the chamber, the chemical energy released by the reactive material projectile can be calculated; and then y can be calculated by comparing the complete reaction heat of the reactive material. Applying polynomial fitting to the direct ballistic test data shown in Table 1, we find that the impact velocity threshold for the reaction of Al/Ni reactive material is 1054 m s−1.

3.1. Finite Element Model

The finite element model (Figure 2(a)) to simulate the direct ballistic test is constructed in LS-DYNA software (version 971 smp d R8.0.0, Livermore Software Technology Corp., Livermore, CA, U.S.). The model is a two-dimensional axisymmetric model, which simplifies calculation by using two-dimensional axisymmetric elements. The blue area represents the projectile, with a radius of 0.5 cm, a height of 1 cm, and a side length of 0.05 cm of each element; the red area stands for the steel target with a thickness of 2.5 cm, a radius of 2 cm, and a side length of 0.05 cm of each element. The velocity of the projectile is along the positive direction of the Y-axis, and the upper end of the steel target is a fixed boundary, which constrains the displacement of all nodes on this boundary in all degrees of freedom.

The steel target plate is modeled with Johnson–Cook (J–C) constitutive model and Gruneisen EOS (see Appendix A), whose parameters are taken from Tables 2 and 3. Meanwhile, the projectile made of Al/Ni mixture is modeled with J–C constitutive model and user-defined EOS. The parameters (Table 4) for the J–C constitutive model of the Al/Ni mixture are taken from [13].

3.2. Parameters for JWL Equation Obtainment

The parameters for the user-defined EOS of reactive materials fall into two parts. The first part is the parameters for the JWL equation of reactants and products, and the other part is the parameters for equation (13). The JWL equation parameters (A, B, R1, R2, ω, and Cv, altogether six parameters) of reactants and products can be obtained by fitting Hugoniot data of reactive materials. According to. [14], the parameters for JWL equations follow these two constraint relations:where represents material specific heat, which can be calculated with the mass ratios of the components. At the same time, it is generally assumed that R2 = R1/10. Thus, we can assume a set of values of R1 and ω to calculate the pressure in the object at each relative volume, and then compare it with the impact adiabat based on the Hugoniot data. Fit the parameters with genetic algorithms, and finally obtain the parameters for the JWL equation [14, 19].

The projectiles used in the direct ballistic test are made of Al/Ni mixtures, and the mass ratio is approximately 1 : 1. Because there are some voids in the mixture, the EOS of the porous mixture with voids is taken as the EOS of reactants. On the other hand, because there is no gas production and the overall y in the direct ballistic is not high, the EOS of the dense mixture is employed as the EOS of the reaction products.

Hugoniot data of dense and loose Al/Ni reactive materials [13] can be fitted with genetic algorithms to obtain the JWL equation parameters corresponding to reactants and products, as shown in Table 5.

3.3. Chemical Reaction Extent Extraction

For comparison with the direct ballistic test data, it is necessary to extract the calculation results of y of the projectile. In this study, the data for y in 20 selected elements in the projectile model are extracted. Then, the data are processed through the weighted mean method according to the radius scale at each element, resulting in the overall extent of reaction in the projectile. Elements selection is shown in Figure 2(b).

3.4. Chemical Reaction Kinetic Parameters Fitting

When the calculated y of the projectile is obtained, it may be different from the direct ballistic test data, so it is necessary to continuously optimize the chemical kinetic parameters to make the calculated results consistent with the test data. In this study, the adaptive simulated annealing (ASA) algorithm is used to fit the parameters. It is a global random optimization algorithm based on the Monte Carlo iterative strategy. The ASA algorithm is built on the similarity between the annealing process of solid material and the general combinatorial optimization problem: the internal energy of the material is replaced by the objective function, and the annealing time is replaced by the control parameters. A simulated annealing algorithm can jump over the trap of the local optimal solution with a certain probability to find the global optimal solution. This is a mature method to solve nonlinear optimization problems [20, 21]. LS-OPT (version 6.0.0, Livermore Software Technology Corp., Livermore, CA, USA) is a simulation-based optimization design tool [22], which integrates the ASA algorithm and can solve complex problems such as parameter identification and system optimization. In addition, LS-OPT can directly extract numerical simulation data from LS-DYNA, so as to realize the fitting of material parameters efficiently and conveniently.

In equation (13), there are three parameters to be identified, namely, Ea, n, and Tcr. Tcr can be obtained through numerical simulation for the case that the impact velocity threshold,  = 1054 m s−1, and the corresponding simulation result is, Tcr = 406 K. Moreover, the other two parameters, Ea and n, can be fitted by the ASA algorithm using the data of direct ballistics tests.

The calculation procedure of LS-OPT is designed to incorporate ASA to fit parameters Ea and n, as shown in Figure 3. Firstly, enter the values of Tcr, and set Ea and n as the variables to be fitted. Secondly, the space-filling sampling method in Radial Basis Function Method is employed to conduct sampling in the variable space, which can select a large number of experimental points for simulation calculation, and avoid errors caused by excessive assumptions. Thirdly, three sets of tests are simulated (2#, 3#, 7# are selected in this study) to extract calculation results of y. Fourthly, the optimization objective S is constructed, which needs to be minimized. Fifthly, an adaptive simulated annealing algorithm (ASA) is used to iterate and finally obtain the reaction kinetic parameters. The optimization objective S iswhere Yi represents the data of numerical results and yi represents the direct ballistic test data.

The chemical kinetic parameters of the reactive material are obtained after 10 iterations (when S = 1.546 × 10−4) (refer to Table 6).

4. Results Analysis

In order to verify the validity of the fitting approach, the mechanical parameters and chemical kinetic parameters of the material (listed in Tables 46), which are calculated based on the test data, are used to perform numerical simulation calculation for #6 and #9 tests. y of the projectile calculated by simulation is shown in Figure 4.

It can be seen from Figure 4 that the chemical reaction is obviously enhanced with increased impact velocity. And the area with high y gradually shrinks along the radial direction of the projectile. This is due to the attenuation of shock wave propagation and the influence of lateral rarefaction wave, which lowers the temperature of elements on the edge, leading to a reduced extent of reaction.

For the case that the impact velocity is 1532 m s−1, the pressure-time curves (Figure 5) of elements 42#, 36#, 101# (red circle in Figure 2(b)) are extracted, and the peak value of the shock pressure is obtained as 24 GPa. The shock pressure is calculated with the one-dimensional shock wave theory to be 21 GPa [13]. Considering the error of JWL equation fitting results and the error of the numerical method, we hold that the pressure calculation of numerical simulation is effective.

The baseline operation function in LS-OPT is employed to calculate y for all test impact velocities. The calculation results are shown in Table 7 and Figure 6. In Figure 6, the numerical simulation results in this study (triangles in Figure 6) are compared with the direct ballistic test data (rectangles in Figure 6) and Ren’s numerical simulation results [13] (circles in Figure 6). It can be found that the results obtained with the numerical approach in this study are closer to the test data, which gains an obvious advantage over Ren’s numerical approach. Table 7 shows the relative error (RE) between the slope and interception of the three fitting curves.

The chemical parameters fitted in [13] are Tcr = 422 K, Ea = 20880 J mol−1, and n = 1.6. y calculated with these parameters is obviously too large for high impact velocities. The reason for the difference is that the temperature used in the fitting is the overall temperature (473 K at the current impact velocity), which is calculated by employing one-dimensional shock wave theory and shock Hugoniot, without considering the influence of shock wave attenuation and rarefaction wave. By contrast, in finite element simulation, the temperature is expressed for individual elements, so the temperature inside the projectile is not uniform (in Figure 7, red represents temperatures greater than 473 K, pure blue represents 298 K, and the selected time is the moment when the initial shock wave is unloaded). When the temperature reaches a certain threshold, the extent of reaction increases rapidly with the increase of temperature. Therefore, the overall extent of reaction is mainly contributed by the elements at high temperature. Simultaneously, when the bullet hits the target plate, strong plastic deformation occurs in front of the bullet and it expands to the radial direction, so the problem cannot be simply reduced to a one-dimensional problem. The parameters fitted with a numerical method based on LS-DYNA and LS-OPT in this study can be readily applied to numerical simulation of reactive materials.

5. Further Discussion

In order to further verify the ability of the numerical approach constructed in this study to simulate the shock-induced energy release behavior of reactive materials, a numerical simulation is carried out on the penetration of a Ф7.5 mm spherical reactive material fragment into a steel plate, as shown in Figure 8. The impact velocity is 2000 m s−1, and the velocity direction is along the negative direction of the Y-axis. The parameters of the reactive material are taken from Tables 46, and the parameters of the steel plate are taken from Tables 2 and 3. Considering the relatively large deformation of the fragment, we adopt the Smoothed Particle Hydrodynamics (SPH) method [23] for the simulation. The SPH numerical model is shown in Figure 8, where the blue particles stand for the steel plate and the red particles for the spherical fragment.

On the basis of the SPH model above, inert material fragment and reactive material fragment with the same mechanical properties are, respectively, used for numerical simulation. The simulation results are shown in Figure 9. Figure 9(a) is the changing temperature contours of both two fragments during penetration. Figure 9(b) shows the average bulk temperature-time curves of both two fragments during penetration (red for reactive material, blue for inert material). According to the simulation results, it can be seen that the reactive material can release a large amount of chemical energy upon impact, which makes the average bulk temperature of the reactive fragment as high as about 1000K, compared with a temperature of 400K for the inert fragment. This means the surface temperature (red area in Figure 9(a)) of the reactive fragment exceeds the melting point of steel (∼1750 K), meaning a stronger ignition and detonation effect on military targets such as fuel tanks and explosives. Therefore, compared with the inert fragment, the damage effect of the reactive fragment is greatly reinforced, which is consistent with [7].

6. Conclusion

(1)The second-development interface of LS-DYNA is exploited to incorporate the user-defined EOS, including the chemical kinetics equation of reactive materials. The user-defined EOS can characterize the mechanical-thermal-chemical coupled behavior of reactive materials during its shock-induced energy release process. Thus, the constructed solver with the user-defined EOS incorporated can be used to perform finite element simulation on shock-induced energy release behavior of reactive materials. The simulation results can provide guidance for relevant experimental design.(2)A method based on the finite element method and ASA algorithm is proposed to fit the chemical kinetic parameters of reactive materials. The parameters fitted with this method can effectively reproduce the results of direct ballistic tests. Compared with the parameters obtained by theoretic formula fitting, the parameters obtained with this method can reproduce the direct ballistic tests more accurately, which assumes practical significance for the engineering application of reactive materials in the military field in the future.

Appendix

In general, for metal materials upon high speed impact, the J–C constitutive model and Gruneisen EOS can be used to describe their mechanical behavior.

The J–C constitutive model [24] describes the relationship between the flow stress σy and the effective plastic strain :where is the normalized effective plastic strain rate; A is the initial yield stress of the material under reference strain rate, , and reference temperature, Troom; and B and N are the strain hardening modulus and the hardening exponent of the material, respectively. C is the strain rate hardening parameter and M is the thermal softening parameter. If the room temperature is Troom and the melting point is Tmelt, the relative temperature TH is defined as

The pressure change of common materials under pressure can be described by Gruneisen EOS [24] in the form of pressure-specific volume-internal energy (P--E):where C0 is sound velocity; S1,S2, S3 are the material parameters; γ0 is the Gruneisen coefficient; a is the first-order volume-correction coefficient of γ0; and u = / − 1. The J–C constitutive model parameters and Gruneisen EOS parameters of steel 45 are listed in Tables 2 and 3, respectively.

Data Availability

The datasets generated and analyzed in the current study may be obtained from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was financed and fully supported by the National Science Foundation of China (Grant no.11672328), which is greatly appreciated by the authors.