Modelling on Shock-Induced Energy Release Behavior of Reactive Materials considering Mechanical-Thermal-Chemical Coupled Effect

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 seconddevelopment interface of LS-DYNA software.+en, 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.+e 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.


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, Fe 2 O 3 , 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]. e 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]. e 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. erefore, it is necessary to find other ways to obtain reaction parameters. Xiong [12] obtained the shockinduced 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. e 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 mechanicalthermal-chemical coupled effect model. en, 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.

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. e relative volume and internal energy meet the following mixing rules: where v represents the relative volume and E represents the internal energy per element of the initial volume. e 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. e JWL equation [14] in the form of pressure-relative volumetemperature (P-v-T) takes temperature and relative volume as independent variables. erefore, in this study, JWL equations are adopted as EOS for reactants and products during the impact process: where C v represents the specific heat capacity of the material and A, B, ω, R 1 , and R 2 are related parameters of the material. e extent of chemical reaction y is controlled by chemical kinetics equation. erefore, 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 mechanicalthermal-chemical coupled effect, updating of each physical quantity is shown in Figure 1. e calculation process for the user-defined EOS is divided into the following seven steps: Step 1: Enter T t , P t , y t , and so forth into LS-DYNA main program.
Step 2: Calculate the predicted temperature value, T * t+1 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.
e 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 ΔT s represents the temperature rise due to plastic deformation. q � 0.5(q t + q t + 1 ) represents the artificial viscous force. J represents the partial derivative of internal energy with volume. Δv 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/_v 0 , represents the internal energy per element's initial relative volume, where e is the element's internal energy and v( ) 0 is the element's initial relative volume. In t + 1 step, E is calculated as follows: σ stands for deviatoric stress. e first term of equation (5) represents the plastic work, the second term Calculate y t+1 , using Eq. (13) Output: Step

LS-DYNA main program
Continue the next circulation Δ < 10 -6 Use T * t+1 and β to calculate Δ Δ = |P r -P p | represents the volume compression work, and the last term represents the element's specific internal energy of t step. Equation (5) is rearranged as where E 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 P * t+1 and other physical quantities. According to [15], an intermediate variable β�(1 − y)v r / v m 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 β. en, the Newton iteration method is employed to iterate β: where P * t+1 r and P * t+1 p are calculated based on β * t+1 . When |P * t+1 r − P * t+1 p | ≤ 10 − 6 iteration stops and P * t+1 r �(P * t+1 + P * t+1 p )/2 is taken as the predicted pressure value for the current time step, then T * t+1 , P * t+1 , and β * t+1 are used to calculate J * t+1 , C * t+1 v , and E * t+1 .
Step 5: Refer to the process in step 3 to update other physical quantities.
Step 6: Calculate y t + 1 . According to the reaction kinetics model given by Arrhenius [16], the chemical reaction rate (dy/dt) follows the form where t is the reaction time, A 0 is the preexponential factor, E a is the apparent activation energy, T is the thermodynamic temperature, and R u is the universal gas constant. f (y) is the reaction mechanism function, and its form is determined by the specific reaction type. e 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, T cr , and the critical extent of reaction, y cr (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: where Step 7: Feed the updated variables to LS-DYNA main program.

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. ey are the activation energy E a , reaction index n, and reaction critical temperature, T cr . e 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. e 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. is material has good mechanical properties and is easy to prepare, so it is often used to make reactive warhead fragments. e 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. e projectile (cylindrical, Φ10 mm × 10 mm) made of Al/Ni reactive material is launched through a ballistic gun. e 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 .

Finite Element Model.
e 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.). e model is a two-dimensional axisymmetric model, which simplifies calculation by using two-dimensional axisymmetric elements. e 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. e 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. e 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. e parameters (Table 4) for the J-C constitutive model of the Al/Ni mixture are taken from [13].

Parameters for JWL Equation Obtainment.
e parameters for the user-defined EOS of reactive materials fall into two parts. e first part is the parameters for the JWL equation of reactants and products, and the other part is the parameters for equation (13). e JWL equation parameters (A, B, R 1 , R 2 , ω, and C v , 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 C v represents material specific heat, which can be calculated with the mass ratios of the components. At the same time, it is generally assumed that R 2 � R 1 /10. us, we can assume a set of values of R 1 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]. e 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.

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. en, 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).

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. e 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      Table 4: Al/Ni reactive material parameters of the J-C constitutive model [13].  global optimal solution. is 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, E a , n, and T cr . T cr can be obtained through numerical simulation for the case that the impact velocity threshold, v cr � 1054 m s − 1 , and the corresponding simulation result is, T cr � 406 K. Moreover, the other two parameters, E a and n, can be fitted by the ASA algorithm using the data of direct ballistics tests. e calculation procedure of LS-OPT is designed to incorporate ASA to fit parameters E a and n, as shown in Figure 3. Firstly, enter the values of T cr , and set E a 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.
irdly, 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. e optimization objective S is where Y i represents the data of numerical results and y i represents the direct ballistic test data. e chemical kinetic parameters of the reactive material are obtained after 10 iterations (when S � 1.546 × 10 − 4 ) (refer to Table 6).

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 4-6), 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. is 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. e 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. e baseline operation function in LS-OPT is employed to calculate y for all test impact velocities. e 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.
e chemical parameters fitted in [13] are T cr � 422 K, E a � 20880 J mol − 1 , and n � 1.6. y calculated with these parameters is obviously too large for high impact velocities. e 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 onedimensional 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. erefore, 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. e 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.

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 V7.5 mm spherical reactive material fragment into a steel plate, as shown in Figure 8. e impact velocity is 2000 m s − 1 , and the velocity direction is along the negative direction of the Y-axis.
Shock and Vibration e parameters of the reactive material are taken from Tables 4-6, 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. e 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. e 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. is 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. erefore, compared with the inert fragment, the damage effect of the reactive fragment is greatly reinforced, which is consistent with [7].  Fit. curve of sim.

Conclusion
(1) e second-development interface of LS-DYNA is exploited to incorporate the user-defined EOS, including the chemical kinetics equation of reactive materials. e user-defined EOS can characterize the mechanical-thermal-chemical coupled behavior of reactive materials during its shock-induced energy release process. us, 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. e 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. e 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.