## Abstract

The ignition threshold of an energetic material (EM) quantifies the macroscopic conditions for the onset of self-sustaining chemical reactions. The threshold is an important theoretical and practical measure of material attributes that relate to safety and reliability. Historically, the thresholds are measured experimentally. Here, we present a new Lagrangian computational framework for establishing the probabilistic ignition thresholds of heterogeneous EM out of the evolutions of coupled mechanical-thermal-chemical processes using mesoscale simulations. The simulations explicitly account for microstructural heterogeneities, constituent properties, and interfacial processes and capture processes responsible for the development of material damage and the formation of hotspots in which chemical reactions initiate. The specific mechanisms tracked include viscoelasticity, viscoplasticity, fracture, post-fracture contact, frictional heating, heat conduction, reactive chemical heating, gaseous product generation, and convective heat transfer. To determine the ignition threshold, the minimum macroscopic loading required to achieve self-sustaining chemical reactions with a rate of reactive heat generation exceeding the rate of heat loss due to conduction and other dissipative mechanisms is determined. Probabilistic quantification of the processes and the thresholds are obtained via the use of statistically equivalent microstructure sample sets (SEMSS). The predictions are in agreement with available experimental data.

## 1 Introduction

Under rapid mechanical loading, energetic materials (EM) undergo thermal-mechanical deformation that can lead to the formation of localized heating due to bulk inelasticity, void collapse, fracture, friction, and volumetric compression. Material heterogeneities enhance the localization of heating into small regions that are called hotspots. Inside the hotspots, the elevated temperatures can lead to the initiation of a chemical reaction which releases heat that further exacerbates the coupled mechanical-thermal-chemical process. On the other hand, thermal conduction, convection, and other transfer processes tend to lower the temperatures inside the hotspots, thereby contributing to “quenching” the reactions and at the same time increasing the temperature on the periphery of the hotspots. The outcome of the competing effects depends on the relative rates of local heat generation and heat loss that are determined by the microstructure, constituent attributes, and overall applied loading. For a given material and microstructure, higher load intensity and more rapid load application correspond to a higher likelihood of sufficient chemical heating rate to overcome the competing heat loss and self-sustainment of reaction and, therefore, higher likelihood of ignition.

The ignition threshold of an energetic material (EM) quantifies the macroscopic conditions for the onset of self-sustaining chemical reactions. The threshold is a theoretically important and practically used measure of material attributes related to safety and reliability. Historically, ignition thresholds are measured experimentally. Weingart et al. [1] performed shock initiation experiments in the pressure range between 3.1 and 28 GPa using a thin flyer launched by an electric gun (E-gun). The flyer impacts a polymer-bonded explosive (PBX) pellet, and the ignition threshold is reported in the Walker–Wasley form [2] using loading pressure and time duration of loading. James [3,4] added specific kinetic energy to the Walker–Wasley threshold and extended the range of impact loading. The James threshold relation is expressed between the specific kinetic energy and the accumulative energy imparted into the material. The relation provides good descriptions of the results of experiments on PBXs with aluminum and plastic flyers [1,5,6]. Welle et al. [7] modified the James relation by replacing the specific kinetic energy with the rate of energy input into the material (which they called power flux). The results for different types of HMX (Octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine) obtained from their electric gun-driven plate impact experiments are well described by the modified James ignition threshold. Although the experiments provide direct quantification, they do not lend themselves to the systematic exploration of the effects of heterogeneities such as microstructural attributes and defects [8], owing to the fact that such attributes cannot be systematically varied [9] or controlled [10] in the experiments. Furthermore, experiments do not provide insight into the mechanisms underlying the manifestation of the ignition behavior nor can they allow the potential behavior of materials not yet available in the laboratory to be analyzed. Experiments are also expensive and require large numbers of tests to pinpoint the precise threshold condition. Consequently, only limited ignition threshold data are available, even for commonly used and otherwise well-studied EM. Computational simulations and prediction provide a means to address these issues and also enable the design and exploration of new materials not yet in existence.

In recent years, Barua et al. [11] developed a Lagrangian cohesive finite element method (CFEM) for computationally predicting the probabilistic ignition behavior of PBX under relatively low-intensity loading. Kim et al. [12] later extended the framework to the regime of shock pulse loading with load pressures up to 13.8 GPa. The predictions are in good agreement with experimental measurements. Subsequent analyses have concerned the effects of initial interfacial defects in microstructure [8] and aluminization [13] on the probabilistic ignition behavior of PBX. The framework allowed, for the first time, ignition thresholds to be predicted out of microstructure and fundamental material properties. The approach uses microstructure-explicit simulations that account for microstructural effects, constituent properties, and interfacial responses. To enable the determination of the ignition threshold, a phenomenological ignition criterion that relates the critical hotspot sizes and critical hotspot temperatures is used. This extrinsic criterion is determined independently from chemical kinetics calculations [14], outside the CFEM simulations. In essence, the criterion allows the effects of chemical heating to be accounted for in a purely thermal-mechanical computational framework. As such, it represents an indirect means of incorporating the effects of chemistry. Although good agreements are obtained with measured ignition thresholds [12,15], there is a desire and need to directly incorporate the effects of chemistry in the simulations in a coupled mechanical-thermal-chemical manner. Such a framework, which tracks the ignition as a natural outcome of the coupled mechanical-thermal-chemical processes, would also allow assessment of the accuracy of the purely thermal-mechanical approach with the independent ignition criterion.

To account for the effect of chemistry in analyses, reactive burn models are typically used to characterize the heat contribution of the chemical decomposition into reaction products of the energetic materials [16]. Hubbard and Johnson [17] first used an Arrhenius-type reaction rate model to analyze the initiation of homogeneous EM [18,19]. Later, McGuire and Tarver’s three-step kinetic model [20] accounts for both an endothermic decomposition process and a subsequent exothermic gas-phase decomposition [21]. Henson et al. [22] developed a decomposition model that captures component processes including the initial beta to the delta phase transition, solid to gas decompositions, and gas-phase ignition. Lee and Tarver proposed the ignition and growth model [18], which includes an ignition phase with hotspot formation, reaction, and thermal runaway and a growth phase for the propagation of reactive waves triggered by burned hotspots [23]. The History Variable Reactive Burn (HVRB) model [24] is another example; it is usually used in conjunction with the Grüneisen EOS for unreactive solids and the Jones–Wilkins–Lee (JWL) EOS for gaseous products [25]. Similarly, the surface burn model [26,27] that considers the reaction rate’s temperature dependence is also often used with the Grüneisen EOS and JWL EOS for the reactant and product, respectively. Most reactive burn models focus on the continuum characteristics of EM reaction [19]. Relatively speaking, the reaction in the context and at the scale of heterogeneous microstructures has only been studied recently and explicit account of heterogeneous microstructures is still a developing endeavor.

Here, we present a new Lagrangian computational framework for establishing the probabilistic ignition thresholds of heterogeneous EM using mesoscale simulations. The framework allows the ignition to evolve out of the coupled mechanical-thermal-chemical processes, with chemistry intrinsically accounted for as part of the simulations, thereby eliminating the need for the use of an extrinsic hotspot-based ignition criterion previously used [14]. The simulations explicitly account for microstructural heterogeneities, constituent properties, interfacial processes, and capture processes responsible for the development of damage and formation of hotspots in which chemical reactions initiate. The specific mechanisms tracked include viscoelasticity, viscoplasticity, fracture, post-fracture contact, frictional heating, heat conduction, reactive chemical heating, gaseous product generation, heat conduction, and convective heat transfer from reaction products to the surface of the unreacted material. The reactive heat release is characterized by an Arrhenius chemical kinetics model and a convective heat transfer model. To determine ignition thresholds, the minimum loading required to achieve self-sustaining chemical reactions with a rate of reactive heat generation exceeding the rate of heat loss due to conduction and other transport mechanisms are determined. Probabilistic quantification of the processes and the thresholds is obtained via the use of statistically equivalent microstructure sample sets (SEMSS) [28]. The coupled mechanical-thermal-chemical framework is used for the analyses of responses up to the attainment of criticality that leads to the determination of the ignition thresholds of energetic materials. As a Lagrangian framework, it is not useful for post-ignition analyses of significant reaction propagation, deflagration, transition to detonation, or the propagation of detonation. Such post-ignition can and should be analyzed using Eulerian approaches. However, the benefit of this Lagrangian framework is that it allows factors and mechanisms that otherwise cannot be tracked by Eulerian approaches to be resolved and explicitly analyzed. Such factors and mechanisms include microstructure evolution, material interfaces, fracture, post-fracture contact and friction along interfaces, and frictional heating. The coupled mechanical-thermal-chemical framework is used to replace the previously used phenomenological hotspot-based ignition criterion that embodies an indirect incorporation of the effects of chemistry.

This paper consists of two parts. The first part describes the computational framework and the materials analyzed. PBX 9501 that consists of HMX grains and Estane (thermoplastic polyurethane) is the material of choice here because of the availability of data [29]. The second part discusses the prediction of initiation thresholds using the novel framework. Particular emphasis is placed on quantifying the probabilistic nature of the ignition thresholds resulting from intrinsic material heterogeneities. For this purpose, sets of multiple samples with statistically equivalent microstructures are generated and used, in direct mimicking of experimental quantification of statistical distributions of material behavior with multiple specimens.

## 2 Framework of Analysis

### 2.1 Materials and Microstructures.

Polymer-bonded explosives (PBXs) are heterogeneous composites consisting of energetic particles and a polymer binder. Some of the most common energetic particles in practical use include HMX, RDX, TATB, and PETN. Among these, HMX holds the highest energy density; therefore, it has drawn intense interest for several decades [30–33]. HMX-based PBXs include PBX 9501, which consists of HMX (∼95 wt%), Estane (∼2.5 wt%), and Bis-2,2-dinitropropyl-acetal/formal (BDNPA/F) as a nitroplasticizer ingredient (2.5 wt%).

In this study, the microstructures that are computationally generated have an HMX grain volume fraction of 81% and an Estane binder volume fraction of 19%. The theoretical volume fraction of HMX in PBX 9501 is 92.6%. However, lower HMX volume fractions are typically observed in actual microstructures because some HMX particles are too small to be resolved and some are absorbed in the binder (resulting in what is often referred to as dirty binder). For example, Benson and Conley [34] observed a binder volume fraction of 26% from a micrograph of PBX 9501 whose theoretical binder volume fraction is only around 8%. Mas et al. [35] observed a binder volume fraction of 23% for PBX 9501 and reproduced the stress-strain behavior using an explicit finite element framework. Barua and Zhou [36] used the microstructures of PBX 9501 with an HMX volume fraction of 81–82% in their numerical study and obtained stress-strain curves that match experimental data. This paper uses the same HMX volume fraction and the same numerical framework as those used in Ref. [36].

The HMX particles typically have random polygonal shapes [34,37,38]. To design microstructures with similar attributes, a library of grains extracted from microstructures generated by the Voronoi tessellation is first established, with the sizes of the grains systematically tabulated. The grains are then used to compose the PBX microstructures with prescribed grain volume fractions and grain size distributions. Details of the method and the attributes of the microstructures are described in Refs. [28,39]. This approach allows large numbers of microstructure samples with prescribed grain size distributions and other attributes to be obtained efficiently. In this section and the following discussions, we use the term “samples” to indicate computationally generated microstructures for finite element simulations and “specimens” to indicate experimentally used ones. Sets of multiple samples (five for each set) with statistically equivalent attributes but different random distributions of the constituents are generated and used. Computationally analyzing the behavior of a set of multiple statistically similar samples is equivalent to carrying out experiments on multiple specimens of the same material, both allowing the statistical variations and probabilistic distributions of material behavior to be quantified. To illustrate the random variations in microstructure morphology and the statistical consistency among multiple samples, Fig. 1 shows five samples with the same HMX volume fraction of *η* = 0.81 and the corresponding variations in grain size distributions among the samples. The grains have an average size of 210 *μ*m and a monomodal size distribution with a standard deviation of 66 *μ*m.

### 2.2 Loading Configuration.

The computational model emulates the experiments with a thin flyer launched by an electric gun impacting on a PBX specimen as carried out by Weingart et al. [1] and Welle et al. [7,40]. The impact generates a loading pulse that propagates through the specimen, as seen in Fig. 2 of Ref. [41]. The boundary conditions and the loading conditions imitate the condition of such experiments. Specifically, impact loading is effected by applying a prescribed particle velocity at the impact face (left boundary of the sample) for a specified time duration, as shown in Fig. 2. The top and bottom boundaries are constrained such that lateral expansion does not occur. This is a 2D model and the conditions of plane-strain prevail. This configuration approximates the planar shock pulse loading of a sample under conditions of approximate macroscopic uniaxial strain. The imposed particle velocity and duration are chosen to correspond to the loading characteristics in the shock experiments [1,7,40]. In Weingart et al. [1], flyer velocities range from 1–5 km/s, and the flyer thickness varies from 1.27 mm to 25 *μ*m. The resulting shock pressure (or input stress) depends on the impedances $(\rho Us)$ of the flyer and the PBX sample. The range of loading conditions analyzed in the experiment corresponds to the imposed particle velocity range of *U _{p}* = 371–1960 m/s. In the calculations, the range of

*U*= 200–1200 m/s is considered. The pulse duration is determined by the flyer thickness and is equal to the time it takes the longitudinal wave to traverse a round trip in the flyer. The range of flyer thickness in the experiment of Welle et al. [7] corresponds to a pulse duration range of

_{p}*τ*= 10–90 ns. The pulse durations used in this study range from 20 ns to 4.5

*μ*s. To precisely determine the threshold (load) pulse duration for ignition at each load intensity (imposed particle velocity at impact face

*U*), multiple simulations are carried out for each sample with the same load intensity but successively longer durations. For example, for a given load intensity of

_{p}*U*= 200 m/s, the load duration of

_{p}*t*

_{pulse}= 1.50

*µ*s is used for all microstructures. Then, with the same load intensity, a slightly increased load duration is imposed on all microstructures. This process is repeated until the load duration reaches

*t*

_{pulse}= 4.50

*µ*s with the minimum increment of Δ

*t*

_{pulse}= 2 ns for all load velocities. This procedure is performed for all load intensities (

*U*= 200–1200 m/s). The loading conditions used in this study are listed in Table 1, including the imposed velocity and the range of pulse durations. The profile of the imposed load pulse at the impact face is shown in Fig. 2(b). The velocity increases rapidly from zero to

_{p}*U*over the ramp time of

_{p}*t*

_{ramp}= 10 ns to reduce the likelihood of numerical instability. This velocity is kept constant until the end of the pulse duration $\tau $. After the pulse ends $(t\u2265\tau )$, the impact face (left boundary) is released and no external loading is applied, while the boundaries on the top, bottom, and right remain constrained in their respective normal directions.

### 2.3 Lagrangian Computational Framework.

The numerical approach explicitly captures hotspot evolution due to thermo-mechanical energy dissipation inside the microstructures. The simulations are performed using a recently developed Lagrangian cohesive finite element (CFEM) framework [11,36,42]. This framework allows quantification of the effects of microstructure and thermal-mechanical-chemical processes, including bulk deformation, interfacial debonding, fracture of grains, subsequent frictional heating, temperature-dependent material decomposition, and subsequent convective/conductive heat transfer to be considered. The constitutive relations for the grains are those of a hydrostatic stress-dependent, elasto-viscoplastic material. The binder follows an elastic-viscoelastic constitutive law. Arbitrary fracture patterns along grain-grain and grain-binder boundaries and inside individual constituents are explicitly captured by the use of cohesive elements embedded throughout the microstructure, along all element boundaries. The cohesive elements follow a bilinear traction-separation law which consists of an initial reversible separation process within a certain separation limit, followed by irreversible damage and separation beyond the limit. A contact detection algorithm and a subsequent contact force model are used for surfaces after a fracture. The Coulomb friction damping model is used for surfaces that are in contact. Fourier’s heat conduction model is coupled with mechanical deformation and failure models to account for thermal conduction in the material. The reaction rate of the solid is based on an Arrhenius-type reaction kinetics. Convective heat transfer is used to track the energy. Details of the cohesive element framework and algorithm and the finite element numerical scheme can be found in Ref. [36].

#### 2.3.1 Governing Equations.

*J*, $\rho 0$, $\rho $, $\tau $, $D$, $Dp$, $t$, $b$, $\mu f$, $u$, $[v\Gamma ]$,

*c*

_{p},

*T*, $\chi h$,

*k*

_{c},

*h*, and $H\u02d9chem$ denote, respectively, deformation gradient, Jacobian, mass density in the reference configuration, mass density in the current configuration, Kirchhoff stress, rate of deformation, plastic part of rate of deformation, traction on a surface with normal $n$, body force, friction coefficient, displacement, relative velocity at discontinuity interface, specific heat, temperature, fraction of plastic work converted to heat, heat conductivity, convective heat transfer coefficient, and heat rate per unit volume due to chemical reaction.

#### 2.3.2 Constitutive Laws.

*P*is the hydrostatic pressure that is obtained by

*λ*and

*μ*are Lamé’s first and second constants. $D\u2032$ in Eq. (5) is the deviatoric part of the rate of deformation, which can be decomposed into an elastic part and a viscoplastic part as

*m*and

*a*are rate sensitivity parameters for a low regime of strain rate and a high regime of strain rate, respectively. The specific range of the regime varies depending on the reference strain rates, $\epsilon \xaf\u02d90$ and $\epsilon \xaf\u02d9m$. The parameter $\sigma 0$ is the quasi-static yield stress, $\epsilon 0$ is a reference strain,

*N*is the strain hardening exponent,

*T*

_{0}is a reference temperature, and $\beta $ and $\kappa $ are thermal softening parameters. The function $g(\epsilon \xaf,T)$ represents the quasi-static stress-strain response at ambient temperature. The above relations consider strain hardening and strain-rate dependence of plasticity. The details of the above constitutive relations and descriptions of the parameters can be found in Ref. [43]. The parameters of the plasticity model for HMX used in this study are listed in Table 2. The parameters are calibrated to match the experimental wave profile obtained by Ref. [44]. The verification of the calibrated parameters is described in Ref. [45].

*t*refer to physical and reduced times, respectively. The bulk modulus

*K*

_{0}of the polymer is assumed to be a constant as in Refs. [46,47]. The deviatoric part of the constitutive behavior of the polymer binder is described by a Prony series. The shear modulus

*G*is assumed to vary with the relaxation time $\tau r$ in the form of

*G*

_{e}is the long-term modulus when the binder is fully relaxed, and $\tau ir$ and

*G*

_{i}are the relaxation time and the modulus of

*i*th mode, respectively. The modulus of the binder is highly dependent on the temperature and the loading rate. The frequency-modulus relations for different temperatures are shifted and superposed by using the Williams–Landell–Ferry (WLF) shift function with the form of

*a*

_{T}is the superposition parameter,

*C*

_{1}and

*C*

_{2}are empirical parameters,

*T*is the temperature, and

*T*

_{0}is the reference temperature. The value of the parameters is given in Table 3. Mas et al. [48] showed that the storage modulus $G\u2032(\omega )$ can be represented in terms of the Prony series parameters by using the equation of

#### 2.3.3 Equation of State.

*K*

_{0}is the bulk modulus, and $K0\u2032=(\u2202K0/\u2202P)P=0$. For the implementation of the BM EOS, a time incremental form of the above is used. Since the time rate of change of the Jacobian is

This relation, details of which are obtained by differentiating Eq. (14), is needed in evaluating the responses of the constituents.

Dattelbaum and Stevens [37] obtained the parameters of the BM EOS for the Estane binder of PBX 9501 at various temperatures based on their experiments. The initial bulk modulus *K*_{0}, decreases as the temperature increases. However, the first derivative $K0\u2032$ shows consistent values with the average of $K0\u2032=12.95$. We chose the entire set of *K*_{0} values in Ref. [37] with $K0\u2032=12.95$. The parameters of the BM EOS for Estane are listed in Table 4.

*P*,

*v*, and $\Gamma $ are pressure, specific volume, and Grüneisen parameter. Internal energy

*E*is a function of temperature.

*P*

_{0}and

*E*

_{0}are 0 K states of pressure and specific internal energy, respectively. One can reference the Hugoniot curve to evaluate the 0 K curve, i.e.,

*P*

_{H}and

*E*

_{H}are reference points on the Hugoniot curve. From the Hugoniot equations for the conservation of mass,

*U*

_{s}is the shock velocity, and

*U*

_{p}is the particle velocity. Rearrangement of Eq. (19) yields

*v*

_{0}is the specific volume of uninfluenced material. For HMX,

*U*

_{s}and

*U*

_{p}has a linear relation given as

*C*

_{0}is referred to as the bulk speed of sound and

*s*is linear Hugoniot slope coefficient. These are material-dependent parameters. Combining Eqs. (20) and (21) yields

*E*

_{0}in the form of

*e*

_{I0}=

*E*

_{0}/

*v*

_{0}is the specific internal energy at reference state. Considering the boundary conditions

*E*

_{0}= 0 and

*P*

_{0}= 0 if

*v*=

*v*

_{0}or $\chi =0$, we have

*E*

_{00}= 0 and

*E*

_{01}= 0. Equation (31) is further simplified to

*E*

_{02}can be evaluated by substituting the second-order term of Eq. (29) into Eq. (28). Thus, the first-order MG EOS for evaluating the hydrostatic pressure can be given as

The parameters in the MG EOS are tabulated in Table 5 for solid HMX [50,51]. Figure 3 shows a comparison between the MG EOS and BM EOS for HMX under bulk compression. The parameters used for the BM EOS for HMX are chosen from Ref. [12]. As shown in Fig. 3(b), when *v*/*v*_{0} < 0.98, the MG EOS provides a higher bulk modulus than the BM EOS, making the material following the MG EOS more resistant to bulk compression than that following the BM EOS. The MG-EOS is used to describe the volumetric part of the stress tensor for unreacted solid HMX granules.

*a*,

*Q*,

*B*, $\omega \u2032$, and

*R*are material-dependent EOS parameters. The detonation energy

*E*

_{d0}contains chemical bond energy and kinetic energy associated with the momentum aspect of the flow [54]. $v^$ is the relative specific volume in the forms of

*v*and

*v*

_{0}being the current and initial specific volumes, respectively. In Eq. (35), the exponential term that describes an increased effect of the repulsive forces on the internal energy at small volumes [52] shares a similar form as that in Jones and Miller [55], where a theoretical approach was used to obtain the EOS for 2,4,6-trinitrotoluene (TNT) reaction products. To extend the application of Eq. (35) to lower pressures with cylinder tests, in which the expansion rate of reacting explosives is measured in a metal cylinder, Kury et al. [56] and Lee et al. [57] added a second exponential term to Eq. (35) such that

*A*and

*B*are linear pressure coefficients, $\omega \u2032$ is the fractional part of the normal Tait equation adiabatic part,

*R*

_{1}and

*R*

_{2}are principal and secondary eigenvalues [54],

*E*

_{d0}is the detonation energy. The relative specific volume $v^$ has been defined in Eq. (36). Equations (36) and (37) together are named the Jones–Wilkins–Lee (JWL) EOS [57]. The parameters for HMX are taken from Lee et al. [58] and are listed in Table 6. The JWL EOS for HMX reaction products is plotted in Fig. 4.

For simplicity in the calculation of the pressure in the HMX gaseous reaction products *P*_{gas}, we assume that the gaseous products for a certain amount of reactant only fill the volume occupied by the reactant if no reaction occurred.

In a mixture of unreacted HMX solid and its gaseous reaction products, the hydrostatic part of the stress tensor or pressure *P* in Eq. (2) is calculated as a weighted average of the solid hydrostatic pressure obtained from the MG EOS in Eq. (33) and the gas pressure obtained from the JWL EOS in Eq. (37), i.e.,

*r*is the extent of reaction. The elastic modulus, flow strength, and fracture energy of the mixture are multiplied by a factor (1 −

*r*) to account for the weakening of the mechanical properties of the mixture relative to the original unreacted solid.

#### 2.3.4 Chemical Reaction.

_{2}and CH

_{2}O, and the further bimolecular reaction and the formation of the HCO product. A simple Arrhenius model that implicitly accounts for the whole processes of the HMX chemical decomposition is convenient to use if only the concentrations of the initial reactant and the final gaseous reaction products are to be tracked. The reaction rate coefficient

*k*in the Arrhenius chemical kinetics model has the form of

*A*

_{0}is the pre-exponential factor,

*E*

_{a}is the activation energy,

*R*is the gas constant, and

*T*is the temperature.

There are multiple Arrhenius models that track the chemical decomposition of HMX. Figure 5 shows the time history for HMX to fully react at an initial temperature of *T*_{0} = 680 K. Among other models [59–66] plotted in Fig. 5, the activation energy *E*_{a} is considered temperature-dependent in the model proposed by German et al. [67]. *A*_{0} = 4.71 × 10^{13} s^{−1}. *R* = 8.314 J · K^{−1} · mol^{−1}. German et al.’s model is also computationally advantageous due to the slower evolution of the initial reaction products, which is beneficial to maintaining numerical stability. In this paper, German et al.’s model is used.

*m*denotes the current mass of the reactant and the exponent

*x*is the reaction order. The negative sign on the left-hand side of Eq. (40) ensures the reaction rate is positive. For simplification, only the first-order process is considered, i.e.,

*x*= 1. Integrating Eq. (40) between two time-steps yields [59]

*m*

_{0}represents the mass of unreacted HMX at the onset of a time step,

*m*denotes the mass at the end of the time-step, and the reaction rate coefficient

*k*is determined by Eq. (39). Note that the final state of

*m*in each step is only a function of the initial concentration in the time-step and temperature.

#### 2.3.5 Heat Equations.

*E*

_{s},

*E*

_{g},

*P*

_{s}, and

*P*

_{g}are internal energies and hydrostatic pressures for the solid and the gas, respectively.

*r*is the reaction ratio. Under this condition,

*c*

_{p}and

*c*

_{v}are specific heats at constant pressure and constant volume, respectively. Thus, the specific form of the heat equation can be obtained as

*k*

_{c}is thermal conductivity, $\chi h$ is the fraction of plastic work that is converted into heat, $W\u02d9p$ is the rate of plastic work per unit volume, $H\u02d9ve$ is the rate of viscoelastic dissipation per unit volume, $H\u02d9fric$ is the rate of frictional dissipation per unit volume, $H\u02d9bulk$ is the rate of heat release per unit volume due to volumetric compression, $H\u02d9inconv$ and $H\u02d9outconv$ are input and output convective heat transfer rates per unit volume,

*P*

_{g}is the gas pressure, and $u$ is the velocity vector.

*r*= 0, only purely unreacted solid HMX exists and Eq. (46) simplifies down to the specific form of the solid energy balance equation

*r*= 1, only gaseous products from the completed reaction permeates the local space and Eq. (46) degenerates to the specific form of the gas energy balance equation

*V*that contains a total area of Δ

*S*of internal frictional surface pairs,

*v*

_{rel}is the relative sliding velocity at the contact surface pair. Table 7 shows the parameters of the ratio of plastic work converted into heat $(\eta )$, density $(\rho )$, specific heat (

*c*

_{p}), and thermal conductivity (

*k*

_{c}) for HMX and the Estane binder of PBX 9501.

Arbitrary fracture patterns inside the individual constituents and interfacial debonding between particles and the binder are explicitly captured by the use of cohesive elements embedded throughout the finite element model. The cohesive elements follow a bilinear traction separation law described by Zhai et al. [76]. The cohesive relation embodies initial reversible separation processes with a certain separation limit, followed by irreversible damage and separation beyond the limit. If the separation reaches a critical distance, the cohesive surface pair is considered failed, and no stress is needed for further separation. Details of the cohesive element framework are provided in Barua and Zhou [36].

The formation of a crack (inside a constituent or along a grain-binder boundary) results in the creation of two surfaces. At each computational time-step, the entire domain is scanned and such surfaces are identified. The corresponding nodal coordinates of all possible pairs of surfaces are compared to detect surface contact and overlap. Penalty forces are applied to strongly discourage interpenetration and maintain proper contact of the surfaces. Detailed descriptions of the multi-step contact and the penalty force algorithm are given in Ref. [45]. Frictional heating, due to sliding along surfaces in contact, is assessed using the Coulomb friction law. The stick-slip state is determined by the normal force between the contact surface pair. Green et al. [77] obtained the frictional coefficient of approximately 0.3–0.7 for PBX 9404, and Chidester et al. [78] used the value of 0.5 based on the work of Green et al. [77]. In this analysis, we take the frictional coefficient to be 0.5.

*t*can be evaluated using the area between the two curves as

*T*

_{0}= 680 K.

Convective heat transfer is a process that hot fluids pass heat to a cooler solid surface. Specifically in HMX, it is associated with hot gaseous reaction products penetrating into pores in the cooler solid unreacted granules [81] to enlarge the local hotspot. Following the conductive heat transfer, the convective burning is the key mechanism that leads to the propagation of the combustive wave. Generally, diffusion and advection are the two mechanisms associated with convective heat transfer. In the rapid and violent reaction of HMX, the strong compression waves govern the energy transfer, to which diffusion has limited contribution [82]. Therefore, the rate at which the combustion wave propagates is not determined by the diffusion of reactive HMX [83]. In the present framework, only forced convection (i.e., advection) driven by the pressure gradient in the gaseous products of HMX is considered.

*h*is the heat transfer coefficient, $\rho gas$ and

*c*

_{p,gas}are the density and specific heat of the gaseous reaction products, and $\u2207T$ is the temperature difference between the gas and unreacted solid per unit length.

*u*

_{gas}is the velocity at which the gaseous reaction products move and is a function of pressure gradient in the gaseous reaction products, i.e., [84]

## 3 Results and Discussion

In order to establish the “go” or “no-go” threshold of ignition which is defined as the sustained propagation of reaction in the microstructure as indicated by the sustained spatial spread of hotspots, calculations and analyses are performed in the following steps. First, each microstructure sample in a given SEMSS, as described in Sec. 2.1, is subject to the loading conditions discussed in Sec. 2.2. Second, at each load intensity, the evolutions of the temperature field at successively longer load durations are analyzed. Note that, although the temperature field depends on the interplay among a number of heating and transfer processes as stated earlier, the chemical reactive heating and the convective heat transfer are the ultimate manifestations of reaction and ignition. Tracking and analysis of the temperature field evolution allow the macroscopic condition [combination of energy flux (energy input rate due to loading) and energy fluence (total energy imparted into the sample by loading)] corresponding to the criticality of onset of ignition to be determined. This is how the ignition thresholds are determined through the precise determination of the onset of ignition as the energy fluence *E* is increased by increasing the load duration.

To illustrate this analysis, Fig. 8 shows the response of the same sample under two levels of energy fluence corresponding to two pulse durations at the same load intensity of *U*_{p} = 600 m/s or input pressure of *P* = 3.69 GPa. Figure 8(a) (top row of images) shows the temperature field evolution after pulse loading over a duration of $\tau =112ns$ (corresponding energy fluence *E* = 24.24 J/cm^{2}). No rapid spread of hotspots is observed, indicating no ignition. On the other hand, Fig. 8(b) shows the result for a slightly longer duration of $\tau =114ns$ (*E* = 24.74 J/cm^{2}) but otherwise the same conditions. Hotspots with very high temperatures (>700 K) form and spread rapidly in a self-sustained manner, indicating ignition. The energy fluence for ignition at this loading intensity of *P* = 3.69 GPa is thus taken as *E* = (24.24 + 24.74)/2 = 24.49 J/cm^{2}. This power flux-energy fluence pair represents one data point in the macroscopic condition space for characterizing the ignition threshold of the same sample over a range of power fluxes analyzed $(\Pi =0.02\u22121.49GW/cm2)$. The difference in the input energy fluence between the two (the minimum increment in $\tau $ is 2 ns) trans-threshold loading conditions is ∼2%. Figure 9 compares the fraction of HMX that has fully reacted at *t* = 0.95 *μ*s under different levels of energy fluence (a measure of input energy). Each data point represents a single loading instance in the power flux-energy fluence space that represents the macroscopic ignition threshold for the sample. It is found that as $\tau $ (and therefore *E*) increases, there is a distinct threshold *E* level (24.49 kJ/cm^{2}) above which the reaction progresses very rapidly (i.e., ignition will occur). Although significant scatter exists due to the intrinsically heterogeneous nature of the material, it is clear that the ignition threshold can be established this way.

### 3.1 Analytical Quantification of Ignition Threshold.

*E*) is the time integration of the power flux $(\Pi )$ over the pulse duration.

*E*

_{c}and $\Pi c$ are fitting parameters that represent asymptotic thresholds for the critical energy fluence and the critical power flux, respectively.

Figure 10 shows a comparison in the $E\u2212\Pi $ space of the thresholds for PBX 9501 predicted using the coupled mechanical-thermal-chemical model developed here (circle data points) and using the approach involving a hotspot-based initiation criterion developed by Barua et al. [42] (square data points). The ignition-criterion-based approach entails simulations that are mechanical-thermal in nature, with the effects of chemistry accounted for through the ignition criterion. The criterion is outside the simulation and is used to determine the onset of ignition based on the size and temperature states of hotspots. As is the case in Ref. [12], the specific threshold information regarding the hotspot criticality is taken from the results of the chemical kinetics calculations by Tarver et al. [14].

Each data point obtained from simulation in Fig. 10 represents a load instance for one of the five microstructures. The solid lines are fits of the data points to the James threshold relation [Eq. (55)] with 50% probability of ignition. The initiation threshold of PBX 9501 obtained by explicitly accounting for chemical reaction (red data points) shows the same trend as that obtained using the hotspot-based ignition criterion of Barua et al. [42] (blue data points). The average energy fluence *E* needed for ignition decreases as the power flux $\Pi $ increases. The data points show less scatter at higher load intensities. The threshold relation obtained by the coupled thermal-mechanical-chemical simulations (red) is generally lower than the threshold line obtained using the criterion (blue). The primary factor leading to this result may be the fact the coupled thermal-mechanical-chemical simulations account for the effects of chemical reaction on the thermal-mechanical processes while the hotspot-based criterion approach does not. The experimental data points shown are estimates based on reports of “go” and “no-go” observations by Chidester et al. [91] and Idar et al. [92] for conditions of long duration impact events. Overall, the initiation thresholds predicted using the coupled thermal-mechanical-chemical simulations and the hotspot-based criterion approach are both in good agreement with the experiments of Idar et al. [92] and Chidester et al. [91]. The data point from Gustuvsen et al. [93] in Fig. 10 appears to be an outlier. This can be at least partly attributed to the fact Gustuvsen et al. [93] used the critical energy of another material (PBX 9404) in interpreting their experiment due to a lack of data for the actual material (PBX 9501). Factors contributing to the deviation may also include the use of the one-step Arrhenius chemical kinetics here. A more detailed reaction model is preferable when concentrations of intermediate decomposition products control the reaction rate. The values of *E*_{c} and $\Pi c$ are listed in Table 9.

### 3.2 Ignition Probability Map.

*J*is introduced to the threshold equation [Eq. (55)] [94] via,

*J*= 1 yields the relation in Eq. (55) which represents loading conditions leading to the “average” or 50% probability of ignition,

*J*> 1 corresponds to loading conditions resulting in ignition probabilities greater than 50%, and

*J*< 1 corresponds to conditions resulting in ignition probabilities less than 50% [12].

*J*= 1. The cumulative probability can be expressed in terms of

*J*or

*E*and $\Pi $ [12] as

In Eq. (58), the standard deviation $\sigma $, cutoff energy fluence *E*_{c} and cutoff power flux $\Pi c$ are material constants whose values are determined by fits to the data sets. Once these parameters are determined, the probability of ignition $P$ under any loading condition as measured by *E* and $\Pi $ can be calculated directly from Eq. (58). The ignition probability map obtained using the hotspot-based criterion approach is shown in Fig. 12(a), and the map obtained using the coupled mechanical-thermal-chemical framework is shown in Fig. 12(b). These maps are in general agreement, except that the probability predicted by the coupled mechanical-thermal-chemical framework is slightly lower than that predicted by the hotspot-based criterion approach.

## 4 Summary and Conclusion

The ignition thresholds of energetic materials that quantify the macroscopic conditions for the onset of self-sustained chemical reactions are theoretically important and practically used measures of material attributes. Historically, the thresholds are determined experimentally. A computational approach using coupled thermal-mechanical simulations and an extrinsic criterion accounting for the effects of chemical heat release has been developed and used in recent years. While the approach considers the mechanism of hotspot growth which underlies the criticality of ignition, it does not account for the effects of chemical heat on the thermal-mechanical processes. Here, we have presented a new Lagrangian computational framework that allows the probabilistic ignition thresholds of heterogeneous energetic materials (EM) to be determined as a direct outcome of the evolution of the coupled mechanical-thermal-chemical processes underlying the response of the materials. The new approach obviates the need for any extrinsic ignition criterion. The ignition threshold predicted using this new framework and the threshold predicted by the existing approach based on an extrinsic ignition criterion are in general agreement, with the threshold predicted by the new framework on average approximately 15.8% lower than that predicted by the criterion-based approach. The reason for this difference can be logically explained by the fact that the new framework accounts for the coupling between the chemical process and the thermal-mechanical process while the existing criterion-based approach does not. However, the difference is relatively small and both predictions are consistent with available experimental results reported in the literature, with the new approach providing a better match with the experimental data.

It is useful to point out that neither method involves curve-fitting with respect to experimental observations on initiation threshold, or requires prior information about the predicted behavior. Instead, the prediction is based on material microstructural attributes and fundamental constituent properties.

It is also important to point out that the coupled mechanical-thermal-chemical Lagrangian framework here is useful for analyses of responses up to the attainment of criticality that leads to the determination of the ignition thresholds of energetic materials. As a Lagrangian framework, it is not useful for post-ignition analyses of significant reaction propagation, deflagration, transition to detonation, or the propagation of detonation. However, the benefit of this Lagrangian framework is that it allows factors and mechanisms that otherwise cannot be tracked by Eulerian approaches to be resolved and explicitly analyzed. Such factors and mechanisms include microstructure evolution, material interfaces, fracture, post-fracture contact and friction along interfaces, and frictional heating. Previous studies [8,12,95] have shown that the frictional heating plays an important role in hotspot formation and thus significantly contributes to initiation under the shock conditions of impact up to 1000 m/s. Of course, like Eulerian and other approaches, the current framework allows elastic, viscoelastic, elastic-viscoplastic constituent response, pressure-dependence of the material response, bulk inelastic heating, heat conduction, chemical reaction, gaseous product generation, and convective heat transfer to be accounted for.

Finally, although only one material configuration is analyzed here to illustrate the use of this new framework, it can be used to study a wide range of material configurations and issues of high explosives, reactive materials, and propellants that contain not only energetic granules and polymers but also metallic fuels and oxidizers. In addition, microstructure defects such as microcracks, constituent clustering, and gradients can be explicitly tracked as well.

## Acknowledgment

The authors gratefully acknowledge the support from the Defense Threat Reduction Agency (DTRA) through project No. HDTRA1-18-1-0004 (Dr. Jeffery Davis) and the Air Force Office of Scientific Research through MURI project No. FA9550-19-1-0008 (Dr. Mitat Birkan). CM's contribution to this work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344. Lawrence Livermore National Security, LLC. LLNL-JRNL-820786. Some calculations are carried out on supercomputers at the AFRL DSRC of the U.S. DoD High Performance Computing Modernization Program. The authors also would like to thank Dr. David Kittell at Sandia National Laboratories for enlightening discussions that enhanced our understanding of some issues in the model development.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The authors attest that all data for this study are included in the paper. Data provided by a third party are listed in Acknowledgment.

## References

^{™}and Kel-F-800

^{™}