Simulation of shock initiation in explosives using a model combining high computational efficiency with a free choice of mixture rules

Models for shock initiation in explosives must consider how energy transfers from products to reactants. This is based on different energy-apportionment assumptions, which affect the results for shock initiation. This study proposes a robust model of shock initiation in explosives using a free choice of energy-apportionment assumptions. The reacting explosive is treated as a mixture of reactants and products under different energy-apportionment assumptions. The equations of state of the mixture are efficiently solved by refining the Cochran–Chan concept of the real volume fraction and introducing a real energy fraction term. The validity, efficiency, and universality of the proposed model are verified in one-dimensional numerical simulations of the shock initiation of homogeneous (nitromethane) and heterogeneous (PBX 9404) explosives. Compared to the conventional Cochran–Chan and Johnson–Tang–Forest models, the proposed model has a better balance of computational efficiency and universality.Models for shock initiation in explosives must consider how energy transfers from products to reactants. This is based on different energy-apportionment assumptions, which affect the results for shock initiation. This study proposes a robust model of shock initiation in explosives using a free choice of energy-apportionment assumptions. The reacting explosive is treated as a mixture of reactants and products under different energy-apportionment assumptions. The equations of state of the mixture are efficiently solved by refining the Cochran–Chan concept of the real volume fraction and introducing a real energy fraction term. The validity, efficiency, and universality of the proposed model are verified in one-dimensional numerical simulations of the shock initiation of homogeneous (nitromethane) and heterogeneous (PBX 9404) explosives. Compared to the conventional Cochran–Chan and Johnson–Tang–Forest models, the proposed model has a better balance of computational efficiency and universality.


I. INTRODUCTION
Shock initiation in explosives remains a hot topic in the detonation field. To improve the safety of munitions, explosives should be sufficiently insensitive to unpredictable stimuli. However, if the sensitivity is lowered, the output pressure of the explosive train needs to be increased to initiate the explosion of the main charge more reliably. Therefore, shock initiation models are used to assess the hazards and reliability of explosives.
The process of shock initiation in homogeneous explosives has been well studied. 1 In the process, the increase in bulk temperature due to the induced shock controls the explosion. At high temperatures, the explosive decomposes and releases some energy, which further increases the temperature. When the critical temperature is reached, the reactants explode and a compression wave forms. The compression wave overtakes the leading shock and finally builds up to a steady detonation. Thus, a detonation is the result of a thermally activated reaction that obeys the Arrhenius rate law. The Arrhenius rate has been successfully used to calculate the shock initiation of some liquid 1,2 and single-crystal 1-3 explosives.
However, shock initiation in heterogeneous explosives is more complicated. 4 The shock-todetonation transition (SDT) of heterogeneous explosives results from localized regions of high temperature, 4 known as hot spots. 5 The formation of hot spots is determined by a competition between the generation and loss of heat through heat transfer. Heat is generated in a variety of ways, such as through cavity collapse, 6-8 friction, [9][10][11] and shear. 12,13 Various mesoscale models that account for shock interactions with the inhomogeneities in explosives have been developed to model the formation of hot spots. These models involve one or more a Electronic mail: zhangxr@bit.edu.cn heat generation processes, that is, hot-spot formation mechanisms, of which cavity collapse may be the most widely studied. Due to the interaction of the shock with the cavity, cavity collapse involves multiple heat generation processes, such as the adiabatic compression of gas in the cavity, shear between explosive material that has softened and still hard explosive material, and the viscous or plastic flow of the explosives. The ignition of some liquid and solid explosives by the adiabatic compression of gas has been observed 14,15 under the low-velocity impact of a striker. Andersen and Gillespie 16 modeled the process and concluded that adiabatic compression is a viable mechanism in the timescale of a low-velocity impact. However, for the shock initiation of explosives, the ignition only occurs under suitable hot spot size 17 and adiabatic compression is not a dominated mechanism. 18, 19 Frey 12 proposed a model where rapid shear resulted in a temperature high enough to ignite the explosive in shear bands. Based on the shear band concept, Kipp 20 simulated the SDT of PBX 9404 (94 wt% HMX, 3 wt% nitrocellulose, and 1 wt% tris-beta-chloroethylphosphate) and successfully reproduced the experimental data. However, during the process of cavity collapse, viscoplastic work caused by the viscous flow of explosives can produce a hot spot with a much higher temperature. 7,21 Besides cavity collapse, experiments have shown that friction and shear are also possible ignition mechanisms. 22,23 Models to assess friction 24,25 and shear 7 have been developed and they reveal the corresponding physical mechanisms.
The above mesoscale models with specific mechanisms for generating hot spots give insights into the ignition of explosives. However, the SDT may comprise a parallel realization of several hot-spot mechanisms with the dominant one depending on the physical and chemical properties of the particular explosives. 26 Moreover, full-scale simulations of shock initiation via mesoscale models are quite costly.
In this respect, macroscale or phenomenological models, which do not consider hot-spot mechanisms, may lose some localized information but are more efficient. The models attempt to reproduce the experimental data using appropriate forms of the equation of state (EOS) and reaction rate. The frequently referenced models are the Forest Fire model, 27 the Ignition and Growth model (I&G), 28,29 and the Johnson-Tang-Forest model. 30 The I&G model, as the most widely used phenomenological model, has been incorporated into several hydrodynamic codes, such as LS-DYNA 31 and AUTODYN. 32 It has been used successfully to reproduce numerous experimental data sets on shock initiation and detonation. 33 In this model, a reaction progress parameter, F, divides the SDT into three phases. In the first phase (F = 0), the reactant is solid. The third phase (F = 1) is when the explosive has completely reacted, so that there are only gaseous detonation products. The intermediate phase (0 < F < 1) corresponds to the partially reacted state, and the explosive under this phase is called the reacting explosive. The reaction progress parameter F is derived via the reaction rate equation. In each phase, the material is characterized by an EOS. For reactants and products, the EOS is obtained from shock compression data and detonation experiments, respectively. The reacting explosive is treated as a mixture of reactants and products. Therefore, the mixture EOS 34 is used to characterize the material properties of the reacting explosive. Other phenomenological models 35-37 adopt similar procedures but use different forms of reaction rates or EOSs.
The mixture EOS should be based on certain mixture rules, that is, several physical assumptions. The first assumption is that the particles of reactants and products move with the same velocity, 30 which results in perfect mixing. 35 The second assumption is that there is mechanical equilibrium between reactants and products. 38 The majority of available phenomenological models [27][28][29][30][35][36][37] adopt these assumptions.
Moreover, since the mixture EOS contains four unknown variables, namely, the specific volumes and internal specific energies of both components (reactant and product), another assumption is needed to derive them. This assumption allows one to solve the mixture EOS in a closed form and predicts the energy transfer from products to reactants. 39 We call this the energy-apportionment assumption.
Noteworthy is that Johnson et al., 30 Cowperthwaite, 40 and Wackerle 41 reported that the choice of energy-apportionment assumption has no apparent effect on the numerical results for the shock initiation of PBX 9404. Kubota et al. 42 explains why the choice of energy-apportionment assumption was negligible. Based on an analysis of the slopes of various constant-pressure curves of the two components in specific internal energy versus specific volume coordinates, these authors concluded that the energy-apportionment assumption has little impact on the numerical results of shock initiation if the curves have the same slopes. Any difference in the slopes (especially in the low-pressure region) will yield different results.
Therefore, a robust energy-balance model that can utilize different energy-apportionment assumptions is quite topical. The Johnson-Tang-Forest model 30 partially satisfied this requirement but was found to be too time-consuming, insofar as it required the inversion of a 4 × 4 matrix in each iteration. In 1979, Cochran and Chan 38 introduced a procedure based on determining the real volume fraction (RVF) occupied by the reactants, which reduced the order of the iteration equations from 4 to 1 and, thus, significantly improved the computational efficiency. However, the temperature equilibrium condition used as the energy-apportionment assumption in this model strongly restricted its adaptability and wide application.
The present study aims to devise an efficient model with a free choice of energy-apportionment assumptions that can simulate shock initiation in explosives. Based on a number of physical assumptions, a general model, which allows the application of various forms of EOS and reaction rates and a flexible choice of energy-apportionment assumptions, is developed. To improve the efficiency of the model, the Cochran-Chan RVF 38 is extended to energy, with the introduction of the respective real energy fraction (REF). Using the RVF and the REF reduces the order of the iteration equations and, thus, greatly improves the computational efficiency of the model. The feasibility and effectiveness of the proposed model were proved in shock initiation simulations for nitromethane and PBX 9404.

A. Governing equations
The process of shock initiation is governed by conservation equations. In the Lagrangian framework, the conservation equations can be written as follows: where u i is velocity, ρ is density, e is specific internal energy, and σ ij is the Cauchy stress. Here, the tensors are expressed using the indicial notation, and the Einstein summation convention is used. The Cauchy stress can be split into the deviatoric stresses s ij and pressure p, as follows: where δ ij is the Kronecker delta. The pressure is calculated using the EOS of the explosive. In shock initiation, the deviatoric stresses in the explosive are neglected since they have a much smaller yield strength compared to the incident shock wave. 34 In Eq. (3),Q is the heat source per unit mass that is liberated when the explosive reacts. Thus, the heat source can be expressed asQ = ∆HḞ, where ∆H is the heat of reaction andḞ is the reaction rate.

B. Reaction rate
The process of shock initiation in homogeneous explosives results from a thermal explosion. The decomposition can be characterized by the Arrhenius rate. 1 Shock initiation in heterogeneous explosives is controlled by hot spots. 4,5 The hot-spot mechanism is still not very well understood. 26 Moreover, experimentally measuring hot spots is difficult. However, the explosives that have hot spots can be treated as homogeneous, because the hot spots are much smaller than the cell size in the numerical calculations. Therefore, a phenomenological model can describe shock initiation in heterogeneous explosives. The I&G model 28,29 may be the most widely used and has been incorporated into some commercial software such as LS-DYNA and Autodyn. In the I&G model, the reaction rate isḞ where η = ρ/ρ 0 − 1 and subscript 0 denotes the initial state. I, b, a, x, G 1 , c, d, y, G 2 , e, g, z, F igmax , F G1max , and F G2min are constants. The three terms on the right-hand side represent three stages of hot spots: (1) formation and ignition, (2) inward or outward growth, and (3) coalescence at high temperature and pressure. The I&G model has been verified using experimental data on shock initiation in heterogeneous and homogeneous explosives. 43

C. Equations of state of reactants and products
The reacting explosive consists of reactants and products and thus, the EOSs of the reactants and products are needed in simulations of shock initiation. The Jones-Wilkins-Lee (JWL) EOS 44 is widely used for detonation products. The isentrope of the JWL EOS through the Chapman-Jouguet point can be expressed as where the subscript s denotes the isentrope, V = v/v 0 is the relative volume, and A, B, R 1 , R 2 , and C are constants. With the isentrope as a reference curve, the regular JWL EOS can be derived from the thermodynamic relationship (∂E/∂V ) s = −p. By substituting Eq. (6) into the relationship and then integrating it, the energy in the isentrope is reduced to where E = ρ 0 e is the energy per unit initial volume. With the Gruneisen coefficient ω held constant, substituting the pressure and energy for the isentrope into the general Gruneisen EOS 45 gives the regular JWL EOS as follows: The JWL EOS can be used for the reactants by replacing the isentrope with a Hugoniot curve as the reference. Moreover, if a temperature variation assessment is required, one can adopt the temperature-dependent JWL equation (hereinafter referred to as JWLT): 46 where C is the average heat capacity. 47

A. The equation of state for reacting explosives
In shock initiation, a modest shock is incident on the explosive. This builds up to a detonation wave, which consists of a leading shock and a reaction zone. At the shock front, there is no or very little reaction. 39,48 Therefore, the EOS of the reactants can be used to calculate the jump in the state variables from the front of the wave to the back. At the end of the reaction zone, the reaction progress parameter, F = 1, and the solid reactants have been converted entirely into detonation products. Thus, the EOS for the products can be used for the region behind the reaction zone.
However, it is quite problematic to resolve a specific reaction zone, where the products may be composed of multiple species. Besides, the composition of the products varies with the reaction. Therefore, it is hard to simulate the shock initiation in explosives with a detailed account of the composition of the products. To tackle this issue, it is assumed that the reacting explosive consists of reactants and permanent products. The latter term implies the composition of the products remains unchanged. Therefore, the reaction is treated as A → B and, thus, a two-phase mixture EOS is used to characterize the reacting explosive in the reaction zone. To solve the mixture EOS, physical assumptions or mixture rules are applied. The first assumption, which implies that the particles of reactants and products move with the same velocity, has the following formulation: 30 where subscripts 1, 2, and c denote reactants, products, and reacting explosives, respectively, while e and are the specific internal energies and volumes, respectively. The second assumption of mechanical equilibrium is achieved after a few collisions between particles of reactants and those of products. 48 Since Eqs. (10)-(12) contain four unknown variables ( 1 , e 1 , v 2 , e 2 ), another assumption, is used to solve this set of equations in closed form. This assumption is about how energy transfers from products to reactants. We call it the energy-apportionment assumption.

B. Energy-apportionment assumptions
The total specific internal energy shared by the reactants and products is controlled by energy conservation. Here the energy-apportionment assumption specifies how the total specific energy is apportioned between two components. Temperature equilibrium (TEQ) 38 and isentropic reactants (ISE) 40 are frequently used as energy-apportionment assumptions. The TEQ assumption implies there is sufficient time for energy to be transferred between reactants and products. The energy transfer is driven by the temperature gradient between the two components. As soon as the temperatures of the reactants and products become equal, i.e., δ 4 = T 1 − T 2 = 0, the energy transfer ceases.
The ISE approach presumes that there is not enough time for energy transfer and thus, the reactants and products are thermally isolated. Therefore, the products take up all the energy released by the reaction of the explosives and the reactants are adiabatic, i.e., δ 4 = e 1 − e S 1 = 0, where e S 1 is the specific internal energy of the corresponding isentropic line.
TEQ and ISE are the two limiting cases. The apportionment rules for the total specific energy in a real detonation could be any compromise solutions falling within the above range. 39 The choice of energy-apportionment assumptions can affect the numerical results for the shock initiation. 42

C. The two-iteration-variable model
The mixture EOS reduced to Eqs. (10)-(13) requires the solution of the following fourdimensional root-finding problem: The approximate solution can be obtained via Newton's method, which is quite robust and efficient: where  30 used this procedure to calculate the shock initiation of PBX 9404 with the TEQ and the ISE energy-apportionment assumptions. However, the procedure is costly because the 4 × 4 matrix B needs to be inverted during each iteration. Cochran and Chan 38 reduced the order of iteration equations from 4 to 1 through three critical steps: (1) replacement of the energy equation by the temperature equation, (2) representing the EOS in the form p = p( , T ), and (3) defining the RVF occupied by the reactants. This procedure significantly improved the computational efficiency. However, the Cochran-Chan model uses only the TEQ assumption and the form of the EOS is restricted to p = p( , T ).
This deficiency is tackled in this study by extending the RVF to energy and introducing the REF in the following form: Using the proposed definition of the REF and the RVF, the variables v 1 , e 1 , 2 , and e 2 can be expressed as and thus, the mixture EOS can be expressed as Eq. (22) is still solved iteratively. Here, we adopt Newton's method, and then the iteration equation reduces to where The definitions of the RVF and REF take advantage of the linearity of Eq. (10) and Eq. (11) and reduces the order of the iteration equations, compared to the Johnson-Tang-Forest model. 30 Therefore, Eq. (22) is more efficient than Eq. (15). In addition, our model does not restrict the choice of energy-apportionment assumption or the form of the EOS. Note that the Johnson-Tang-Forest model, the Cochran-Chan model, and our model are numerical models and thus, they can be used with different EOSs and reaction rates.

A. Dynamic1D
To test our model and compare it with the other two models, the three models were implemented in the Dynamic1D hydrocode, which is used to solve the one-dimensional reactive flow in the Lagrange framework. In one dimension, the conservation equations (1), (2), and (3) take the following form: where q is the artificial viscosity, which is used to smear the discontinuities in the shock front and detonation front. This set of equations is solved by the finite difference method in the Dynamic1D code using the same difference scheme as that applied in the well-known SIN code, 49 which was developed for solving the reactive-flow problem. The artificial viscosity q in Dynamic1D adopts the Wilkins form: 50 where c l and c q are constants, and c s is the speed of sound. Once the initial and boundary conditions are set, the variables can be advanced with each time step based on the stability condition.

B. Shock initiation in nitromethane and PBX 9404
To verify the validity and computational efficiency of our model, we simulated shock initiation in the well-known nitromethane and PBX 9404. In the simulation, these explosives were subjected to sustained shocks by applying a shock pressure as a boundary condition.
Some Lagrangian positions in the explosive disc were labeled to record the histories of the state variables. The histories calculated by different models underwent a comparative analysis. In shock initiation experiments, embedded gauges can be used to record the pressure or particle velocities over time at different Lagrangian positions. These histories adequately describe the evolution of the flow field during shock initiation. The pressure histories constructed with the Cochran-Chan model and the Johnson-Tang-Forest model were used as robust reference data, which we compared with our model. This approach is considered quite reasonable, given that both models have been successfully used for more than 30 years.
To compare the computational efficiencies, we carried out simulations using the three models with different mesh densities. A comparative analysis of our model and the Cochran-Chan model was performed to identify the most efficient model, although the latter model cannot adopt different energy-apportionment assumptions. For the same preset mesh density, the three models used the same constant time step determined by the stability conditions to avoid the effect of cycle number on the computation time. The calculations were carried out with one CPU (Intel(R) Core(TM) i7-6700), and the computation times required for the calculation of the mixture EOS were measured.
The input data for the simulations included the EOS formulations of the reactants and products, the reaction rate, and the energy-apportionment assumption. Since the I&G model has been widely used to simulate shock initiation in heterogeneous explosives, it was adopted in the simulations for PBX 9404. This model uses the JWLT EOSs for the reactants and products and a three-term reaction rate. The parameters of the I&G model for PBX 9404 can be found in Lee and Tarver 28 and Tarver, Hallquist, and Erickson. 29 Although the original I&G model used the TEQ energy-apportionment assumption, in this study, the ISE energy-apportionment assumption was also used to illustrate that our model does not restrict the choice of energy-apportionment assumptions. When all three models were compared, only the TEQ energy-apportionment assumption was used, since it was the only choice for the Cochran-Chan model.
Shock initiation in a homogeneous explosive results from a thermal explosion. The decomposition of the explosive obeys the temperature-dependent Arrhenius law. However, the Arrhenius law does not allow one to estimate shock initiation adequately due to the temperature uncertainty in the reaction zone. 43 However, Tarver and Urtiew 43 recently reported the successful simulation of experimental data for homogenous explosives via the I&G model. Therefore, the I&G model was also used to simulate shock initiation in nitromethane.

A. Shock initiation in nitromethane
A sustained shock with a pressure of 7.5 GPa was induced in nitromethane at t = 0 µs. The propagation of the shock in nitromethane is shown using a time sequence of pressure profiles in Fig. 1. The pressure profiles depict a typical shock initiation in a homogeneous explosive. The shock induced in the nitromethane results in a sharp rise of the bulk temperature and pressure near the shock front. When the temperature is high enough, the explosive decomposes and releases energy, which further increases the temperature. When the critical temperature of the explosive is reached, it explodes to form a compression wave. The compression wave overtakes the leading shock at x = 0.39 cm as shown in Fig. 2. This results in an overdriven detonation, which quickly decays to a steady-state detonation, with a pressure reduction at the leading shock front.
The calculation results depicted in Figs  were chosen because they represent the three typical phases of detonation: (i) the compression wave overtaking the shock, (ii) overdriven detonation, and (iii) the steady-state detonation. As expected, the pressure histories obtained using our model were identical to those of the other two models, as shown in Fig. 3. As expected from our theoretical predictions, we found that our model and the Cochran-Chan model required less computation time than the Johnson-Tang-Forest model for the same simulations. To compare the computational efficiency quantitatively, the three models were used to simulate the same case with different mesh densities. The relative efficiency ratio of each model was defined via its normalization with the output of the Johnson-Tang-Forest model, i.e., RE = t J /t T or t J /t C , where t is the computation time. Subscripts T, C, and J denote our model, the Cochran-Chan model, and the Johnson-Tang-Forest model, respectively.
As shown in Fig. 4 energy-apportionment assumptions, though its efficiency is 2.5 times higher. This improvement can be attributed to the definition of the RVF and REF, which reduced the order of the iteration equations, according to the analysis in Section III C. Although the relative efficiency ratio (RE = 2.5) of our model was lower by 19% than that of the Cochran-Chan model (RE = 3.1), this deficiency is compensated for by the free choice of energy-apportionment assumption, which the Cochran-Chan model cannot do. Figure 5 depicts the shock initiation in PBX 9404 under a sustained shock of 5.0 GPa. After the shock has been induced into the explosive, numerous hot spots form due to cavity collapse, friction, and other factors. The formation of hot spots results in an increase of the pressure at the leading shock, and then the shock rapidly transforms into a steady-state detonation. We validated our model and compared the computational efficiency of the three models as we did for nitromethane in Section V A. Figure 6 shows the pressure histories of PBX 9404 at the three  Figure 7 shows that the respective RE values are consistent with the results for shock initiation in nitromethane. This strongly indicates that our model has a high computational efficiency and that the definition of RE is reasonable. Compared to the Johnson-Tang-Forest model, our model is not only more efficient but also more robust. Figures 8 and 9 show the histories of the specific internal energy of reactants and products at the above three Lagrangian positions. The histories indicate that the specific internal energies, e 1 and e 2 , are discontinuous when the leading shock reaches the Lagrangian positions. The discontinuity of the iteration variables in the Johnson-Tang-Forest model increases the possibility of non-convergence and reduces the rate of convergence.

B. Shock initiation of PBX 9404
In our model, the REF introduced in Section III C plays the role of specific energies e 1 and e 2 in the Johnson-Tang-Forest model. Figure 10  To compare the three models, the TEQ energy-apportionment assumption, which is the only choice in the Cochran-Chan model, was used in the simulations. We also tested the performance of our model and the Johnson-Tang-Forest model using the ISE energy-apportionment assumption for the same shock initiation of PBX 9404. Figure 11 illustrates the excellent agreement between the pressure histories calculated by the two models at the three Lagrangian positions. This close agreement confirms that our model can adopt different energy-apportionment assumptions and accurately estimate the shock initiation. With the ISE energy-apportionment assumption, the computational efficiency of our model is about 2.5 times higher than that of the Johnson-Tang-Forest model. Figure 12 shows the simulation results calculated by our model using the TEQ and ISE energyapportionment assumptions. Although the small difference between the pressure histories obtained with these two different assumptions increases as the distance increases, it is negligible compared to the experimental errors. In other words, the shock initiation in PBX 9404 is insensitive to the choice of energy-apportionment assumption, which supports previous findings in the literature. 40,41,51 However, this insensitivity of PBX 9404 to the energy-apportionment assumption is attributed to the slopes of the constant-pressure lines (in the e versus plane) of the reactants and the products being the same. 42 If the slopes of the two components were different, the choice of energyapportionment assumption may have a reasonable effect on the shock initiation results. Therefore, a model that can adopt different energy-apportionment assumptions is quite important in optimizing this process.

VI. CONCLUSION
A model for shock initiation in explosives was developed with a free choice of energy apportionment assumptions. The model was used to simulate shock initiation in nitromethane and PBX 9404 and two conventional models were used to produce reference data. Compared to the conventional models, the proposed model has a better balance of computational efficiency and universality or a wider choice of energy-apportionment assumptions. Although the computational efficiency of the proposed model was lower by 19% than that of the Cochran-Chan model, 38 it surpasses the latter because it allows a free choice of energy-apportionment assumption. The computation results obtained by the proposed model were identical to those of the Johnson-Tang-Forest model 30 but the computational efficiency was higher by 2.5 times. This improvement was achieved by using both the REF and RVF of the reactants. The proposed procedure can be extended to other mixtures, including alloys. Based on the successful application of our model to one-dimensional shock initiation, follow-up studies are envisaged: (i) to implement the model in two-and three-dimensional codes to further verify its validity and (ii) to examine multiple explosives to determine where the choice of energy-apportionment assumption has a significant effect on shock initiation.