## 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 [3033]. 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.

Fig. 1
Fig. 1
Close modal

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 $(ρ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 Up = 371–1960 m/s. In the calculations, the range of Up = 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 τ = 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 Up), 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 Up = 200 m/s, the load duration of tpulse = 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 tpulse = 4.50 µs with the minimum increment of Δtpulse = 2 ns for all load velocities. This procedure is performed for all load intensities (Up = 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 Up over the ramp time of tramp = 10 ns to reduce the likelihood of numerical instability. This velocity is kept constant until the end of the pulse duration $τ$. After the pulse ends $(t≥τ)$, 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.

Fig. 2
Fig. 2
Close modal
Table 1

Up (m/s)20040060080010001200
Range of pulse duration tpulse (ns)1500–4500200–50060–14040–9020–8020–50
Up (m/s)20040060080010001200
Range of pulse duration tpulse (ns)1500–4500200–50060–14040–9020–8020–50

### 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.

The computational analysis involves solving a set of Lagrangian multi-physics equations for the solids that are given by equations for a balance of mass, momentum, and energy in the form of
${detF=J=ρ0ρ,∫Vτ:δDdV=∫St⋅δu˙dS+∫Vρ0b⋅δu˙dV−∫Vρ0∂2u∂t2⋅δu˙dV,and∫Vρ0cpT˙δTdV=∫Vχhτ:DpδTdV+∫Sfμfn⋅τ⋅n[vΓ]δTdS+∫Skc(F−1F−T∇T)⋅nδTdS−∫Vkc(F−1F−T∇T)⋅∇δTdV+∫Vh(F−1F−T∇T)δTdV+∫VH˙chemδTdV$
(1)
where $F$, J, $ρ0$, $ρ$, $τ$, $D$, $Dp$, $t$, $b$, $μf$, $u$, $[vΓ]$, cp, T, $χh$, kc, h, and $H˙chem$ 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.

The thermo-mechanical constitutive relations for the constituents are as follows. The stress tensor is decomposed into a hydrostatic part and a deviatoric part, i.e.,
$σij=−Pδij+σij′$
(2)
where $σij$ is the Cauchy stress tensor and P is the hydrostatic pressure that is obtained by
$P=−13(σ11+σ22+σ33)=−13σii$
(3)
In the discussions to follow, the hydrostatic stress is often referred to as “pressure,” with a negative sign referring to compressive stress state. The Cauchy stress relates to the Kirchhoff stress via
$τij=Jσij$
(4)
where $J=det(F)$ is the Jacobian, with $F$ being the deformation gradient. The deviatoric part of the stress tensor carried by the HMX grains follows an elasto-viscoplastic constitutive law and the hydrostatic part follows the Mie–Grüneisen (MG) equation of state (EOS), which is discussed in Sec. 2.3.3.
The deviatoric constitutive behavior of the grains is described by
$τ^′=L:(D′−Dp′)$
(5)
where $L$ is the tensor of elastic moduli and $τ^′$ is the deviatoric part of the Jaumann rate of the Kirchhoff stress. For isotropic elastic response,
$L=2μI~+λI⊗I$
(6)
Here, $I$ is the second-order identity tensor, $I~$ is the fourth-order identity tensor, λ and μ are Lamé’s first and second constants. $D′$ 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
$D′=De′+Dp′$
(7)
where $Dp′$ is the viscoplastic part of $D′$ in the form of
(8)
Here, $σ¯$ is the Misses equivalent stress, $τ′$ is the deviatoric part of the Kirchoff stress, and $ε¯˙$ is the equivalent plastic strain rate which has the form of
${ε¯˙=ε¯˙1ε¯˙2ε¯˙1+ε¯˙2,ε¯˙1=ε¯˙0[σ¯g(ε¯,T)]m,ε¯˙2=ε¯˙mexp[−ag(ε¯,T)],andg(ε¯,T)=σ0(1+ε¯ε0)N{1−β[(TT0)κ−1]}$
(9)
where $ε¯=∫0tε¯˙dt$ is the equivalent plastic strain, $ε¯˙0$ and $ε¯˙m$ are reference strain rates, 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, $ε¯˙0$ and $ε¯˙m$. The parameter $σ0$ is the quasi-static yield stress, $ε0$ is a reference strain, N is the strain hardening exponent, T0 is a reference temperature, and $β$ and $κ$ are thermal softening parameters. The function $g(ε¯,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].
Table 2

Parameters in viscoplastic constitutive model of HMX

$σ0(MPa)$$ε0$NT0 (K)$β$
2605.88 × 10−40.02930.0
$ε¯˙0(s−1)$m$ε¯˙m(s−1)$a (1/MPa)κ
1 × 10−4100.08.0 × 101222.50.0
$σ0(MPa)$$ε0$NT0 (K)$β$
2605.88 × 10−40.02930.0
$ε¯˙0(s−1)$m$ε¯˙m(s−1)$a (1/MPa)κ
1 × 10−4100.08.0 × 101222.50.0
The binder considered is Estane, which follows the Generalized Maxwell model (GMM) [46] in the form of
$σ(τ)=∫0τ2G(t−t′)∂εD∂t′dt′+∫0τK0(t−t′)∂εH∂t′dt′$
(10)
where $σ$ represents the Cauchy stress, $εD$ and $εH$ refer to the deviatoric and hydrostatic portions of the Eulerian strain tensor, and $τ$ and t refer to physical and reduced times, respectively. The bulk modulus K0 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 $τr$ in the form of
$G(t)=Ge+∑i=1nGiexp(−tτir)$
(11)
where Ge is the long-term modulus when the binder is fully relaxed, and $τir$ and Gi are the relaxation time and the modulus of ith 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
$log(aT)=−C1(T−T0)C2+(T−T0)$
(12)
where aT is the superposition parameter, C1 and C2 are empirical parameters, T is the temperature, and T0 is the reference temperature. The value of the parameters is given in Table 3. Mas et al. [48] showed that the storage modulus $G′(ω)$ can be represented in terms of the Prony series parameters by using the equation of
$G′(ω)=Ge+∑i=1nGiω2(τir)21+ω2(τir)2$
(13)
where $τir=1.5aT/ω$ as in Mas et al. [48]. They provided the storage modulus $G′(ω)$ and the Prony series parameters as listed in Table 3.
Table 3

Parameters of Prony series for Estane binder in PBX 9501 from Ref. [48]

Frequency (Hz)Gi (MPa)Frequency (Hz)Gi (MPa)
10−60.004171052.6182
10−50.0074110612.882
10−40.0158510752.481
10−30.03802108223.87
10−20.06761109436.52
10−10.089131010457.09
10.11561011346.74
1010.14221012251.19
1020.16221013177.83
1030.22181014117.49
1040.4753101575.858
Ge = 0
Frequency (Hz)Gi (MPa)Frequency (Hz)Gi (MPa)
10−60.004171052.6182
10−50.0074110612.882
10−40.0158510752.481
10−30.03802108223.87
10−20.06761109436.52
10−10.089131010457.09
10.11561011346.74
1010.14221012251.19
1020.16221013177.83
1030.22181014117.49
1040.4753101575.858
Ge = 0

Note: WLF parameters: C1 = 6.5, C2 = 120, T0 = 292 K.

#### 2.3.3 Equation of State.

The volumetric part of the stress tensor for the Estane binder is described by the Birch–Murnaghan (BM) EOS. The specific form of the equation is
$τh=32K0J(J−(7/3)−J−(5/3))[1+34(K0′−4)(J−(2/3)−1)]$
(14)
where $τh=τii=τ11+τ22+τ33$ is the hydrostatic part of the Kirchoff stress, which is the product of the Jacobian and the negative of the hydrostatic pressure. K0 is the bulk modulus, and $K0′=(∂K0/∂P)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
$∂J∂t=J⋅tr(D)$
(15)
the rate of change of the hydrostatic Kirchhoff stress is only a function of the Jacobian and rate of deformation, i.e.,
$∂τh∂t=f(dVdV0,tr(D))$
(16)

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 K0, decreases as the temperature increases. However, the first derivative $K0′$ shows consistent values with the average of $K0′=12.95$. We chose the entire set of K0 values in Ref. [37] with $K0′=12.95$. The parameters of the BM EOS for Estane are listed in Table 4.

Table 4

Parameters of Birch-Murnaghan (BM) equation of state (EOS) for Estane binder in PBX 9501

ParametersEstane
K0 (GPa)T dependent [37]
$K0′$12.95
ParametersEstane
K0 (GPa)T dependent [37]
$K0′$12.95
The volumetric part of the stress tensor for solid HMX granules is described by the MG EOS. The MG EOS relates the hydrostatic pressure to the volume of a solid at a given temperature by
$P−P0=Γv(E−E0)$
(17)
where P, v, and $Γ$ are pressure, specific volume, and Grüneisen parameter. Internal energy E is a function of temperature. P0 and E0 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.,
$PH−P0=Γv(EH−E0)$
(18)
where PH and EH are reference points on the Hugoniot curve. From the Hugoniot equations for the conservation of mass,
$ρ0Us=ρ(Us−Up)$
(19)
where $ρ$ is the density after shock compression, $ρ0$ is the density of uninfluenced material, Us is the shock velocity, and Up is the particle velocity. Rearrangement of Eq. (19) yields
$UpUs=1−ρ0ρ=1−vv0≜χ$
(20)
where v0 is the specific volume of uninfluenced material. For HMX, Us and Up has a linear relation given as
$Us=C0+sUp$
(21)
where C0 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
$Us=C01−χs$
(22)
From the Hugoniot equation for the conservation of momentum
$PH=ρ0UsUp$
(23)
and noting Eqs. (20)(22), we obtain
$PH=ρ0χC02(1−χs)2$
(24)
Similarly, from the Hugoniot equation for the conservation of energy
$EH=12PH(v0−v)$
(25)
Substituting Eqs. (24) and (25) into Eq. (18) and using the relation of
$Γ=Γ0vv0$
(26)
(in which $Γ0$ is Grüneisen parameter at reference state) and the thermodynamic identity [49]
$P0=−∂e0∂v$
(27)
we obtain the ordinary differential equation for E0 in the form of
$ρ0χC02(1−χs)2(1−Γ0χ2)+∂E0∂v+Γ0E0v0=0$
(28)
Gust [49] uses a power series expansion of $χ$ to estimate the expression for E0 with less than 10% error, i.e.,
$E0=E00+E01χ+E02χ2+E03χ3+⋯$
(29)
where E00, E01, E02, E03… are coefficients of the power series. Substituting Eq. (29) into Eq. (27) yields
$P0=1v0(E01+2E02χ+3E03χ2+⋯)$
(30)
Substituting Eqs. (26), (29), and (30) into Eq. (17) and neglecting the higher-order terms, we obtain
$P=1v0[−E00Γ0+E01(1−Γ0χ)+2E02χ(1−Γ02χ)]+Γ0eI0$
(31)
where eI0 = E0/v0 is the specific internal energy at reference state. Considering the boundary conditions E0 = 0 and P0 = 0 if v = v0 or $χ=0$, we have E00 = 0 and E01 = 0. Equation (31) is further simplified to
$P=2E02χv0(1−Γ02χ)+Γ0eI0$
(32)
where the coefficient E02 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
$τh=ρ0χC02(1−χs)2(1−Γ0χ2)+Γ0eI0$
(33)
The bulk modulus is determined using Eq. (33) through
$K0=−vdτhdχdχdv$
(34)

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/v0 < 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.

Fig. 3
Fig. 3
Close modal
Table 5

Parameters of the Mie-Grüneisen (MG) equation of state (EOS) for solid HMX in PBX 9501

$Γ0$$C0(km/s)$$ρ0(kg/m3)$s
1.12.74019102.6
$Γ0$$C0(km/s)$$ρ0(kg/m3)$s
1.12.74019102.6
As the reaction of solid HMX takes place inside hotspot regions, hot gaseous reaction products are formed and confined temporarily in the unreactive granular bed, leading to local high-pressure regions. An equation of state (EOS) that describes the state of matter of the reaction products is required. It has been found that an EOS derived from a thermodynamic approach for the chemical reaction processes can over predict the energy release in the system by 20% [52]. Instead, Wilkins et al. [53] studied the one-dimension flow of PBX-9404 and LX04-01 driving the motion of a spherical metal shell to obtain the Chapman–Jouguet (CJ) adiabat (i.e., adiabatic curve) both experimentally and numerically. Using this more mechanical approach [52], an EOS that better describes the behavior of the chemical reaction products in a one-dimensional sphere test is given as [53]
$P=av^−Q+B(1−ω′Rv^)exp(−Rv^)+ω′Ed0v^$
(35)
where a, Q, B, $ω′$, and R are material-dependent EOS parameters. The detonation energy Ed0 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^=ρ0ρ=vv0$
(36)
with $ρ$ and $ρ0$ being the current and initial densities and v and v0 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
$Pgas=A(1−ω′R1v^)exp(−R1v^)+B(1−ω′R2v^)exp(−R2v^)+ω′Ed0v^$
(37)
where A and B are linear pressure coefficients, $ω′$ is the fractional part of the normal Tait equation adiabatic part, R1 and R2 are principal and secondary eigenvalues [54], Ed0 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.
Fig. 4
Fig. 4
Close modal
Table 6

Jones–Wilkins–Lee (JWL) equation of state (EOS) parameters for HMX reaction products [58]

Chapman–Jouguet (CJ) parameterEquation of State (EOS) parameters
$Ed0(100GPa)$$A(100GPa)$B (100 GPa)R1R2$ω′$
0.10507.7830.070714.201.000.30
Chapman–Jouguet (CJ) parameterEquation of State (EOS) parameters
$Ed0(100GPa)$$A(100GPa)$B (100 GPa)R1R2$ω′$
0.10507.7830.070714.201.000.30

For simplicity in the calculation of the pressure in the HMX gaseous reaction products Pgas, 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.,

$P=(1−r)[ρ0′χ′C02(1−χ′s)2(1−Γ0χ′2)+Γ0eI0]+r[A(1−ω′R1v^)exp(−R1v^)+B(1−ω′R2v^)exp(−R2v^)+ω′Ed0v^]$
(38)
where 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.

The rate of chemical decomposition of HMX can be determined by an Arrhenius-type chemical kinetics model. Henson et al.’s model [22] accounts for three decomposition processes of HMX, including a nucleation and growth stage describing the phase transition from β-HMX to δ-HMX, the decomposition of δ-HMX into NO2 and CH2O, 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
$k=A0exp(−Ea/RT)$
(39)
where A0 is the pre-exponential factor, Ea 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 T0 = 680 K. Among other models [5966] plotted in Fig. 5, the activation energy Ea is considered temperature-dependent in the model proposed by German et al. [67]. A0 = 4.71 × 1013 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.

Fig. 5
Fig. 5
Close modal
For numerical implementation, the reaction rate can be written as [59]
$−dmdt=kmx$
(40)
where 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]
$mm0=exp(−kt)$
(41)
where m0 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.

The temperature in the material under dynamic loading rises locally due to plastic dissipation, viscoelastic dissipation, frictional dissipation along interfaces, heat conduction, heat due to bulk compression, chemical energy release, and conductive/convective heat transfer. We assume there exists a density-temperature equilibrium [68] in the HMX system such that the mass densities and temperatures are the same for the unreacted solids and the reaction products. Therefore, the local internal energy and the hydrostatic pressure of the solid and the gas mixture satisfy the linear rule of mixture [69], i.e.,
${E=(1−r)Es+rEgP=(1−r)Ps+rPg$
(42)
where Es, Eg, Ps, and Pg are internal energies and hydrostatic pressures for the solid and the gas, respectively. r is the reaction ratio. Under this condition,
$ρdEdt=ρ[(1−r)dEsdt+rdEgdt]+ρdrdt(Eg−Es)$
(43)
Note that
$ρdrdt(Eg−Es)=−(1−r)H˙chem$
(44)
where $H˙chem$ is the chemical heat release rate per unit volume of solid due to the decomposition of the reactant. Equation (43) can be rewritten as
$ρdEdt=ρ[(1−r)cp+rcv]∂T∂t−(1−r)H˙chem,$
(45)
where cp and cv are specific heats at constant pressure and constant volume, respectively. Thus, the specific form of the heat equation can be obtained as
$ρ[(1−r)cp+rcv]∂T∂t=(1−r)(kc∇2T+χhW˙p+H˙ve+H˙fric+H˙bulk+H˙chem+H˙inconv)+r(−Pg∇⋅u−H˙outconv)$
(46)
where kc is thermal conductivity, $χh$ is the fraction of plastic work that is converted into heat, $W˙p$ is the rate of plastic work per unit volume, $H˙ve$ is the rate of viscoelastic dissipation per unit volume, $H˙fric$ is the rate of frictional dissipation per unit volume, $H˙bulk$ is the rate of heat release per unit volume due to volumetric compression, $H˙inconv$ and $H˙outconv$ are input and output convective heat transfer rates per unit volume, Pg is the gas pressure, and $u$ is the velocity vector.
When r = 0, only purely unreacted solid HMX exists and Eq. (46) simplifies down to the specific form of the solid energy balance equation
$ρcp∂T∂t=kc∇2T+ηW˙p+H˙ve+H˙fric+H˙bulk+H˙chem+H˙inconv$
(47)
When 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
$ρcv∂T∂t=−Pg∇⋅u−H˙outconv$
(48)
For a volume element ΔV that contains a total area of ΔS of internal frictional surface pairs,
$H˙fric=1ΔV∫ΔSμσnvreldS$
(49)
where $μ$ is the coefficient of friction, $σn$ is the normal force on the contact surface pair, and vrel is the relative sliding velocity at the contact surface pair. Table 7 shows the parameters of the ratio of plastic work converted into heat $(η)$, density $(ρ)$, specific heat (cp), and thermal conductivity (kc) for HMX and the Estane binder of PBX 9501.
Table 7

Material properties of HMX and Estane binder

PropertiesHMXEstane
Density ρ (kg m−3)1910 [70]1190 [71]
Specific Heat cp (kJ kg−1 K−1)1.254 [72]1.500 [36]
Thermal conductivity kc (W m−1 K−1)0.52 [73,74]0.2
Ratio of heat dissipation from plastic deformation $η$0.9 [75]N/A
PropertiesHMXEstane
Density ρ (kg m−3)1910 [70]1190 [71]
Specific Heat cp (kJ kg−1 K−1)1.254 [72]1.500 [36]
Thermal conductivity kc (W m−1 K−1)0.52 [73,74]0.2
Ratio of heat dissipation from plastic deformation $η$0.9 [75]N/A

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.

In Fig. 6, the Hugoniot is illustrated as an assembly of the initial state $(v0,τh0)$ in front of the shock wave and all possible states $(vi,τhi)$ behind the shock front. The Hugoniot curve should be distinguished from the loading path. The Hugoniot is the locus of all states that a shock transition can reach [79], whereas the loading path is represented by a Rayleigh line—the secant of the Hugoniot curve linking the initial and final states in an instantaneous loading process. Unloading process follows the isentropic curve, which is very close to the Hugoniot curve. Practically, the Hugoniot curve is considered the unloading path. The area between a Rayleigh line and the corresponding piece of the Hugoniot curve is considered the energy release in the form of heat due to the bulk loading. The specific form of heat release rate between two adjacent states $(vi,τhi)$ and $(vi+1,τh(i+1))$ over time Δt can be evaluated using the area between the two curves as
$H˙bulk=(vi−vi+1)τh(vi)+τh(vi+1)2−∫vi+1viτh(v)dvΔt$
(50)
where $τh(v)$ is given by Eq. (33) and the integral can be evaluated via integration of $τh(v)$ which is given as
$∫vi+1viτh(v)dv=−C02ρ0s3[v0(2s−Γ0)2(1−χs)+v0(s−Γ0)ln(1−χs)+Γ0vs2]|vi+1vi$
(51)
Fig. 6
Fig. 6
Close modal
Overall, the decomposition process of HMX is an exothermic process. The associated heat release rate $H˙chem$ in Eq. (46) is
$H˙chem=dn^HMXdtΔHc∘,HMX$
(52)
where $n^HMX$ is the amount of unreacted HMX per unit volume, and the enthalpy of combustion $ΔHc∘,HMX=−2820kJ/mol$ [80]. The temperature history can be calculated and is shown in Fig. 7 for an initial temperature of T0 = 680 K.
Fig. 7
Fig. 7
Close modal

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.

The heat flowrate across the unit area due to advection is
$H˙conv=h∇T=ρgascp,gasugas∇T$
(53)
where h is the heat transfer coefficient, $ρgas$ and cp,gas are the density and specific heat of the gaseous reaction products, and $∇T$ is the temperature difference between the gas and unreacted solid per unit length. ugas 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]
$ugas=−κϕμ′∇Pgas$
(54)
where $∇Pgas$ is the pressure gradient in gaseous reaction products, $ϕ$ is the porosity of the unreacted solid HMX, $μ′$ is the viscosity of the gas, and $κ$ is the permeability. For implementation simplicity, the material-dependent parameters in Eqs. (53) and (54) are considered constant and are listed in Table 8. These parameters are in the appropriate ranges given in Refs. [8589] and are further calibrated based on independent simulations using CTH [90], a hydrocode from the Sandia National Laboratories.
Table 8

Convective heat transfer parameters

$μ′(105Pa⋅s)$$ϕ$$κ(10−14m2)$cp,gas (kJ/kg −K)$ρgas(kg/m3)$
1.00.011.8014401910
$μ′(105Pa⋅s)$$ϕ$$κ(10−14m2)$cp,gas (kJ/kg −K)$ρgas(kg/m3)$
1.00.011.8014401910

As illustrated later in Fig. 8 in Sec. 3, convective heating rapidly increases the temperature of the adjacent HMX particles toward the activation of chemical reaction, and consequently can lead to a self-sustained combustion process across the material.

Fig. 8
Fig. 8
Close modal

## 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 Up = 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 $τ=112ns$ (corresponding energy fluence E = 24.24 J/cm2). 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 $τ=114ns$ (E = 24.74 J/cm2) 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/cm2. 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 $(Π=0.02−1.49GW/cm2)$. The difference in the input energy fluence between the two (the minimum increment in $τ$ 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 $τ$ (and therefore E) increases, there is a distinct threshold E level (24.49 kJ/cm2) 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.

Fig. 9
Fig. 9
Close modal

### 3.1 Analytical Quantification of Ignition Threshold.

The ignition threshold is quantified using the modified James threshold relation [3,7,15]
$ΠcΠ+EcE=1$
(55)
where $Π$ is the power flux $(Π=PUp)$ that measures the rate at which work is imparted to the material per unit area of the impact face. The energy fluence (or input energy, E) is the time integration of the power flux $(Π)$ over the pulse duration. Ec and $Π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−Π$ 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].

Fig. 10
Fig. 10
Close modal

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 $Π$ 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 Ec and $Πc$ are listed in Table 9.

Table 9

Parameters in James initiation threshold relation

PBX 9501Πc (GW/cm2)Ec (kJ/cm2)
Hotspot-based ignition criterion0.01190.0260
Coupled mechanical-thermal-chemical framework0.01380.0195
PBX 9501Πc (GW/cm2)Ec (kJ/cm2)
Hotspot-based ignition criterion0.01190.0260
Coupled mechanical-thermal-chemical framework0.01380.0195

### 3.2 Ignition Probability Map.

Each data point in Fig. 10 represents one load instance of one sample. Multiple samples with statistically similar microstructure attributes are used. This process mimics the use of multiple specimens of the same material batch in experiments to obtain a statistical and probabilistic characterization of a material’s behavior. Here, the scatter among the data points in Fig. 10 is attributed to statistical variations in heterogeneities in the microstructure [94] and can be used to quantify the uncertainties in the ignition behavior of the material arising from the variations among the microstructures. For this purpose, the James number J is introduced to the threshold equation [Eq. (55)] [94] via,
$ΠcΠ+EcE=1J$
(56)
where 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].
The data in Fig. 10 conforms to the Gaussian distribution function $P(J)$, as shown in probability density function (PDE) form in Fig. 11. This distribution is symmetric about the threshold line for J = 1. The cumulative probability can be expressed in terms of J or E and $Π$ [12] as
$P(J)=1σ2π∫−∞Jexp[−(x−1)22σ2]dx=12[1+erf(J−12σ)]$
(57)
and
$P(E,Π)=12+12erf[12σ(EΠEcΠ+EΠc−1)]$
(58)
Fig. 11
Fig. 11
Close modal

In Eq. (58), the standard deviation $σ$, cutoff energy fluence Ec and cutoff power flux $Π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 $Π$ 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.

Fig. 12
Fig. 12
Close modal

## 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

1.
Weingart
,
R. C.
,
Jackson
,
R. K.
,
Honodel
,
C. A.
, and
Lee
,
R. S.
,
1980
, “
Shock Initiation of Pbx-9404 by Electrically Driven Flyer Plates
,”
Propellants Explos.
,
5
(
6
), pp.
158
162
.
2.
Walker
,
F. E.
, and
Wasley
,
R. J.
,
1969
, “
Critical Energy for Shock Initiation of Heterogeneous Explosives
,”
Explosivstoffe
,
17
(
1
), pp.
9
13
.
3.
James
,
H. R.
,
1996
, “
An Extension to the Critical Energy Criterion Used to Predict Shock Initiation Thresholds
,”
Propellants, Explos., Pyrotech.
,
21
(
1
), pp.
8
13
.
4.
James
,
H.
, “
Links Between Macroscopic Behaviour and Explosive Morphology in Shock to Detonation Transitions
,”
Proceedings of the 13th International Detonation Symposium
,
Norfolk, VA
,
July 23
, pp.
952
961
.
5.
Gittings
,
E. F.
,
1965
, “
Initiation of a Solid Explosive by a Short-Duration Shock
,”
Proceedings of the Fourth Symposium (International) on Detonation
,
White Oak, MD
,
Oct. 12
, pp.
373
380
.
6.
Trott
,
B. D.
, and
Jung
,
R. G.
,
1970
, “
Effect of Pulse Duration on the Impact Sensitivity of Solid Explosives
,”
Proceedings of the Fifth Symposium (International) on Detonation
,
,
Aug. 18
, pp.
191
205
.
7.
Welle
,
E. J.
,
Molek
,
C. D.
,
Wixom
,
R. R.
, and
Samuels
,
P.
,
2014
, “
Microstructural Effects on the Ignition Behavior of HMX
,”
J. Phys.: Conf. Ser.
,
500
(
5
), p.
052049
.
8.
Wei
,
Y.
,
Kim
,
S.
,
Horie
,
Y.
, and
Zhou
,
M.
,
2018
, “
Quantification of Probabilistic Ignition Thresholds of Polymer-Bonded Explosives With Microstructure Defects
,”
J. Appl. Phys.
,
124
(
16
), p.
165110
.
9.
Doherty
,
R. M.
, and
Watt
,
D. S.
,
2008
, “
Relationship Between RDX Properties and Sensitivity
,”
Propellants, Explos., Pyrotech.
,
33
(
1
), pp.
4
13
.
10.
Teipel
,
U.
,
2004
,
Energetic Materials: Particle Processing and Characterization
,
Wiley-VCH
,
Weinheim, Germany
.
11.
Barua
,
A.
,
Kim
,
S.
,
Horie
,
Y.
, and
Zhou
,
M.
,
2013
, “
Prediction of Probabilistic Ignition Behavior of Polymer-Bonded Explosives From Microstructural Stochasticity
,”
J. Appl. Phys.
,
113
(
18
), p.
184907
.
12.
Kim
,
S.
,
Wei
,
Y.
,
Horie
,
Y.
, and
Zhou
,
M.
,
2018
, “
Prediction of Shock Initiation Thresholds and Ignition Probability of Polymer-Bonded Explosives Using Mesoscale Simulations
,”
J. Mech. Phys. Solids
,
114
, pp.
97
116
.
13.
Miller
,
C.
,
Kim
,
S.
,
Horie
,
Y.
, and
Zhou
,
M.
,
2019
, “
Ignition Thresholds of Aluminized HMX-Based Polymer-Bonded Explosives
,”
,
9
(
4
), p.
045103
.
14.
Tarver
,
C. M.
,
Chidester
,
S. K.
, and
Nichols
,
A. L.
,
1996
, “
Critical Conditions for Impact- and Shock-Induced hot Spots in Solid Explosives
,”
J. Phys. Chem.
,
100
(
14
), pp.
5794
5799
.
15.
Kim
,
S.
,
Miller
,
C.
,
Horie
,
Y.
,
Molek
,
C.
,
Welle
,
E.
, and
Zhou
,
M.
,
2016
, “
Computational Prediction of Probabilistic Ignition Threshold of Pressed Granular Octahydro-1,3,5,7-Tetranitro-1,2,3,5-Tetrazocine (HMX) Under Shock Loading
,”
J. Appl. Phys.
,
120
(
11
), p.
115902
.
16.
Reynaud
,
M.
,
Sorin
,
R.
,
Dubois
,
V.
, and
Desbiens
,
N.
,
2020
, “
WGT: A Mesoscale-Informed Reactive Burn Model
,”
J. Appl. Phys.
,
127
(
6
), p.
065901
.
17.
Hubbard
,
H. W.
, and
Johnson
,
M. H.
,
1959
, “
Initiation of Detonations
,”
J. Appl. Phys.
,
30
(
5
), pp.
765
769
.
18.
Lee
,
E. L.
, and
Tarver
,
C. M.
,
1980
, “
Phenomenological Model of Shock Initiation in Heterogeneous Explosives
,”
Phys. Fluids
,
23
(
12
), pp.
2362
2372
.
19.
Handley
,
C. A.
,
Lambourn
,
B. D.
,
Whitworth
,
N. J.
,
James
,
H. R.
, and
Belfield
,
W. J.
,
2018
, “
Understanding the Shock and Detonation Response of High Explosives at the Continuum and Meso Scales
,”
Appl. Phys. Rev.
,
5
(
1
), p.
011303
.
20.
McGuire
,
R.
, and
Tarver
,
C.
,
1981
,
Chemical-decomposition Models for the Thermal Explosion of Confined HMX, TATB, RDX, and TNT Explosives
,
Lawrence Livermore National Lab.
,
CA
.
21.
Tarver
,
C. M.
, and
Tran
,
T. D.
,
2004
, “
Thermal Decomposition Models for HMX-Based Plastic Bonded Explosives
,”
Combust. Flame
,
137
(
1–2
), pp.
50
62
.
22.
Henson
,
B. F.
,
Asay
,
B. W.
,
Smilowitz
,
L. B.
, and
Dickson
,
P. M.
,
2002
, “
Ignition Chemistry in HMX From Thermal Explosion to Detonation
,”
Proceedings of the Shock Compression of Condensed Matter-2001, Pts 1 and 2
,
Atlanta, GA
,
June 24
, Vol.
620
, pp.
1069
1072
.
23.
Menikoff
,
R.
, and
Shaw
,
M. S.
,
2010
, “
Reactive Burn Models and Ignition & Growth Concept
,”
EPJ Web Conf.
,
10
, p.
00003
.
24.
Starkenberg
,
J.
, and
Dorsey
,
T. M.
,
1998
,
An Assessment of the Performance of the History Variable Reactive bum Explosive Initiation Model in the CTH Code
,
Army Research Lab Aberdeen Proving Ground MD Weapons Technology Directorate
.
25.
Todd
,
S. N.
,
Anderson
,
M. U.
,
Caipen
,
T. L.
, and
,
D. E.
,
2009
, “
Non-Shock Initiation Model for Explosive Families: Numerical Results
,”
Proceedings of the Shock Compression of Condensed Matter—2009, Pts 1 and 2
,
Nashville, TN
,
June 28
, Vol. 1195, p.
361
.
26.
Partom
,
Y.
,
1998
, “
Predicting PBX-9404 Initiation and Detonation Data with a Calibrated Reaction Model
,”
Proceedings of the Eleventh International Detonation Symposium
,
Snowmass Village, CO
,
Aug. 31
, pp.
33300
33305
.
27.
Partom
,
Y.
,
2001
, “
Hydro-reactive Computations with a Temperature Dependent Reaction Rate
,”
Proceedings of the Shock Compression of Condensed Matter-2001, Pts 1 and 2
,
Atlanta, GA
,
June 24
, Vol.
620
, pp.
460
463
.
28.
Wei
,
Y.
,
Olsen
,
D. H.
,
Miller
,
C. M.
,
Wagner
,
K. B.
,
Keyhani
,
A.
,
,
N.
, and
Zhou
,
M.
,
2020
, “
Computational Design of Three-Dimensional Multi-Constituent Material Microstructure Sets with Prescribed Statistical Constituent and Geometric Attributes
,”
Multiscale Sci. Eng.
,
2
(
1
), pp.
1
13
.
29.
Mulford
,
R. N.
, and
Swift
,
D. C.
,
2001
, “
Mesoscale Modelling of Shock Initiation in HMX-Based Explosives
,”
Proceedings of the Shock Compression of Condensed Matter-2001, Pts 1 and 2
,
Atlanta, GA
,
June 24–29
, Vol.
620
, pp.
415
418
.
30.
Christiansen
,
D. E.
, and
Taylor
,
J. W.
,
1973
, “
HE Sensitivity Study
,”
Los Alamos Scientific Lab, Los Alamos, NM, Report No. LA-5440-MS
.
31.
Hayes
,
D. B.
,
1976
, “
A Pnt Detonation Criterion From Thermal Explosion Theory
,”
Proceedings of the 6th International Detonation Symposium
,
,
Aug. 24
, pp.
95
100
.
32.
Bennett
,
J. G.
,
Haberman
,
K. S.
,
Johnson
,
J. N.
, and
Asay
,
B. W.
,
1998
, “
A Constitutive Model for the non-Shock Ignition and Mechanical Response of High Explosives
,”
J. Mech. Phys. Solids
,
46
(
12
), pp.
2303
2322
.
33.
Austin
,
R. A.
,
Barton
,
N. R.
,
Reaugh
,
J. E.
, and
Fried
,
L. E.
,
2015
, “
Direct Numerical Simulation of Shear Localization and Decomposition Reactions in Shock-Loaded HMX Crystal
,”
J. Appl. Phys.
,
117
(
18
), p.
185902
.
34.
Benson
,
D. J.
, and
Conley
,
P.
,
1999
, “
Eulerian Finite-Element Simulations of Experimentally Acquired HMX Microstructures
,”
Modell. Simul. Mater. Sci. Eng.
,
7
(
3
), pp.
333
354
.
35.
Mas
,
E. M.
,
Clements
,
B. E.
,
Ionita
,
A.
, and
Peterson
,
P.
,
2006
, “
Finite Element Method Calculations on Statistically Consistent Microstructures of PBX 9501
,”
AIP Conf. Proc.
,
845
(
1
), pp.
487
490
.
36.
Barua
,
A.
, and
Zhou
,
M.
,
2011
, “
A Lagrangian Framework for Analyzing Microstructural Level Response of Polymer-Bonded Explosives
,”
Modell. Simul. Mater. Sci. Eng.
,
19
(
5
), p.
055001
.
37.
Dattelbaum
,
D. M.
, and
Stevens
,
L. L.
,
2008
, “Equations of State of Binders and Related Polymers,”
Static Compression of Energetic Materials
,
S. M.
Peiris
, and
G. J.
Piermarini
, ed.,
Springer
,
Virginia
, pp.
127
202
.
38.
Liu
,
C.
,
2003
, “
Specific Surface: A Missing Parameter in High-Explosive Modeling
,”
Los Alamos National Laboratory (LANL), Los Alamos, NM, USA, Report No. LA-UR-14-20512
.
39.
Kim
,
S.
,
Barua
,
A.
,
Horie
,
Y.
, and
Zhou
,
M.
,
2014
, “
Ignition Probability of Polymer-Bonded Explosives Accounting for Multiple Sources of Material Stochasticity
,”
J. Appl. Phys.
,
115
(
17
), p.
174902
.
40.
Welle
,
E. J.
,
Molek
,
C. D.
,
Wixom
,
R. R.
,
Samuels
,
P.
, and
Langhals
,
J.
, “
Microstructure Effects on the Initiation Threshold Behavior of HMX and PBXN-5
,”
Proceedings of the 15th International Detonation Symposium
,
San Francisco, CA
,
July 13
.
41.
Springer
,
H. K.
,
Tarver
,
C. M.
,
Reaugh
,
J. E.
, and
May
,
C. M.
,
2014
, “
Investigating Short-Pulse Shock Initiation in HMX-Based Explosives with Reactive Meso-Scale Simulations
,”
J. Phys.: Conf. Ser.
,
500
(
5
), p.
052041
.
42.
Barua
,
A.
,
Kim
,
S.
,
Horie
,
Y.
, and
Zhou
,
M.
,
2013
, “
Ignition Criterion for Heterogeneous Energetic Materials Based on Hotspot Size-Temperature Threshold
,”
J. Appl. Phys.
,
113
(
6
), p.
064906
.
43.
Zhou
,
M.
,
Needleman
,
A.
, and
Clifton
,
R. J.
,
1994
, “
Finite-element Simulations of Shear Localization in Plate Impact
,”
J. Mech. Phys. Solids
,
42
(
3
), pp.
423
458
.
44.
Dick
,
J. J.
,
Hooks
,
D. E.
,
Menikoff
,
R.
, and
Martinez
,
A. R.
,
2004
, “
Elastic-plastic Wave Profiles in Cyclotetramethylene Tetranitramine Crystals
,”
J. Appl. Phys.
,
96
(
1
), pp.
374
379
.
45.
Hardin
,
D. B.
,
2015
, “
The Role of Viscoplasticity in the Deformation and Ignition Response of Polymer Bonded Explosives
,”
Doctor of Philosophy, Georgia Institute of Technology
.
46.
Mas
,
E. M.
, and
Clements
,
B. E.
,
1996
, “
A Viscoelastic Model for PBX Binders
,” http://lib-www.lanl.gov/la-pubs/00818442.pdf,
Report No. LA-UR-01-3492
.
47.
Wu
,
Y.-Q.
, and
Huang
,
F.-L.
,
2009
, “
A Micromechanical Model for Predicting Combined Damage of Particles and Interface Debonding in PBX Explosives
,”
Mech. Mater.
,
41
(
1
), pp.
27
47
.
48.
Mas
,
E. M.
,
Clements
,
B. E.
,
Blumenthal
,
B.
,
,
C. M.
,
Gray
,
G. T.
, and
Liu
,
C.
,
2002
, “
A Viscoelastic Model for PBX Binders
,”
AIP Conf. Proc.
,
620
(
1
), pp.
661
664
.
49.
Gust
,
W. H.
,
1982
, “
High-Impact Deformation of Metal Cylinders at Elevated-Temperatures
,”
J. Appl. Phys.
,
53
(
5
), pp.
3566
3575
.
50.
Menikoff
,
R.
, and
Sewell
,
T. D.
,
2001
, “
Constituent Properties of HMX Needed for Meso-Scale Simulations
,”
Los Alamos National Lab., Report No. LA-UR-00-3804-rev
.
51.
Gibbs
,
T. R.
, and
Popolato
,
A.
,
1980
,
LASL Explosive Property Data
,
University of California
,
Berkeley
.
52.
Wilkins
,
M. L.
,
1999
,
Computer Simulation of Dynamic Phenomena
,
Springer
,
Berlin, NY
.
53.
Wilkins
,
M. L.
,
Squier
,
B.
, and
Halperin
,
B.
,
1964
, “
Equation of State for Detonation Products of PBX 9404 and LX04-01
,”
Proceedings of the Symposium (International) on Combustion
,
Cambridge, UK
,
Aug. 17
,
Elsevier
, pp.
769
778
.
54.
Urtiew
,
P. A.
, and
Hayes
,
B.
,
1991
, “
Parametric Study of the Dynamic Jwl-Eos for Detonation Products
,”
Combust., Explos. Shock Waves
,
27
(
4
), pp.
505
514
.
55.
Jones
,
H.
, and
Miller
,
A. R.
,
1948
, “
The Detonation of Solid Explosives—the Equilibrium Conditions in the Detonation Wave-Front and the Adiabatic Expansion of the Products of Detonation
,”
Proc. R. Soc. London, Ser. A
,
194
(
1039
), pp.
480
507
.
56.
Kury
,
J.
,
Hornig
,
H.
,
Lee
,
E.
,
McDonnel
,
J.
,
Ornellas
,
D.
,
Finger
,
M.
,
Strange
,
F.
, and
Wilkins
,
M.
,
1965
, “
Metal Acceleration by Chemical Explosives
,”
Proceedings of the Fourth Symposium (International) on Detonation
,
White Oak, MD
,
Oct. 12
, pp.
3
13
.
57.
Lee
,
E.
,
Hornig
,
H.
, and
Kury
,
J.
,
1968
,
Adiabatic Expansion of High Explosive Detonation Products
,
University of California Radiation Laboratory at Livermore
,
Livermore, CA
.
58.
Lee
,
E.
,
Finger
,
M.
, and
Collins
,
W.
,
1973
,
JWL Equation of State Coefficients for High Explosives
,
Lawrence Livermore National Lab. (LLNL)
,
Livermore, CA
.
59.
Skinner
,
D.
,
Olson
,
D.
, and
Block-Bolten
,
A.
,
1998
, “
Electrostatic Discharge Ignition of Energetic Materials
,”
Propellants, Explos., Pyrotech.
,
23
(
1
), pp.
34
42
.
60.
Brill
,
T. B.
,
Gongwer
,
P. E.
, and
Williams
,
G. K.
,
1994
, “
Thermal-Decomposition of Energetic Materials. 66. Kinetic Compensation Effects in Hmx, Rdx, and Nto
,”
J. Phys. Chem.
,
98
(
47
), pp.
12242
12247
.
61.
Suryanarayana
,
B.
, and
Graybush
,
R.
,
2010
, “
Thermal Decomposition of 1, 3, 5, 7-Tetranitro-l, 3, 5, 7-Tetrazacyclooctane (HMX): A Mass Spectrometric Study of the Products From 3-HMX
,”
Proceedings of 39th Congress on Industrial Chemistry, Gr. XS
,
,
July 25
, pp.
24
591
.
62.
Kimura
,
J.
, and
Kubota
,
N.
,
1980
, “
Thermal-Decomposition Process of Hmx
,”
Propellants Explos.
,
5
(
1
), pp.
1
8
.
63.
Kraeutle
,
K.
,
1981
,
18th JANNAF Combustion Meeting
,
, pp.
19
23
.
64.
Medvedev
,
A. I.
,
Sakovich
,
G. V.
, and
Kostantinov
,
V. V.
,
1977
,
Meeting on Kinetics and Mechanism of Chemical Reactions in Solids (Abstracts of Reports)
,
Novosibirsk: Institute of Physical Chemistry of Foundation of Processing of Mineral, Siberian Branch of Academy of Sciences of the USSR
.
65.
Klimenko
,
G. K.
,
1975
, “
Combustion and Explosion
,”
Proceedings of the Fourth All-Union Symposium on Combustion and Explosion [in Russian]
,
Nauka, Moscow
.
66.
Belyayeva
,
M. S.
,
Klimenko
,
G. K.
,
Babaytseva
,
L. T.
, and
Stolyarov
,
P. N.
,
1977
,”
Fifth All Union Symposium on Combustion and Detonation
.
67.
German
,
V.
,
Grebennikova
,
S.
,
Kornilova
,
L. Y.
,
Lobanova
,
S.
, and
Fomicheva
,
L.
,
2002
, “
Thermal Decomposition of PETN and HMX Over a Wide Temperature Range
,”
Proceedings of the 12th Symposium (International) on Detonation, Office of Naval Research
,
San Diego, CA
,
Aug. 11
.
68.
,
C. L.
,
2008
,
Numerical Modeling of Explosives and Propellants
,
CRC Press
,
Boca Raton
.
69.
Kittell
,
D. E.
,
Schmitt
,
R. G.
,
Tuttle
,
L. W.
, and
,
E. N.
,
2018
,
Implementation of a CREST Multistate Reactive Burn Model in CTH for two Solid High Explosives
,
Sandia National Lab. (SNL-NM)
,
Albuquerque, NM, USA
.
70.
Saw
,
C. K.
,
2002
, “
Kinetics of HMX and Phase Transitions: Effects of Grain Size at Elevated Temperature
,”
Proceedings of the 12th International Detonation Symposium
,
San Diego, CA
,
Aug. 11
.
71.
Bourne
,
N. K.
, and
Gray
,
G. T.
,
2005
, “
Dynamic Response of Binders; Teflon, Estane and Kel-F-800
,”
J. Appl. Phys.
,
98
(
12
), p.
123503
.
72.
Shoemaker
,
R. L.
,
Stark
,
J. A.
,
Koshigoe
,
L. G.
, and
Taylor
,
R. E.
,
1985
, “Thermophysical Properties of Propellants,”
Thermal Conductivity 18
,
T.
Ashworth
and
D. R.
Smith
, ed.,
Springer US
,
Boston, MA
, pp.
199
211
.
73.
Gonthier
,
K. A.
,
2003
, “
Modeling and Analysis of Reactive Compaction for Granular Energetic Solids
,”
Combust. Sci. Technol.
,
175
(
9
), pp.
1679
1709
.
74.
Wemhoff
,
A. P.
,
Burnham
,
A. K.
,
Nichols
,
A. L.
, and
Knap
,
J.
,
2007
, “
Calibration Methods for the Extended Prout-Tompkins Chemical Kinetics Model and Derived Cookoff Parameters for RDX, HMX, LX-10 and PBXN-109
,”
Proceedings of the ASME-JSME Thermal Engineering Summer Heat Transfer Conference, American Society of Mechanical Engineers
,
,
July 8
,
American Society of Mechanical Engineers
, pp.
625
632
.
75.
Hodowany
,
J.
,
Ravichandran
,
G.
,
Rosakis
,
A. J.
, and
Rosakis
,
P.
,
2000
, “
Partition of Plastic Work Into Heat and Stored Energy in Metals
,”
Exp. Mech.
,
40
(
2
), pp.
113
123
.
76.
Zhai
,
J.
,
Tomar
,
V.
, and
Zhou
,
M.
,
2004
, “
Micromechanical Simulation of Dynamic Fracture Using the Cohesive Finite Element Method
,”
ASME J. Eng. Mater. Technol.
,
126
(
2
), pp.
179
191
.
77.
Green
,
L.
,
Weston
,
A.
, and
Van Velkinburg
,
J.
,
1971
,
Mechanical and Frictional Behavior of Skid Test Hemispherical Billets
,
California Univ., Livermore, Lawrence Livermore Lab
.
78.
Chidester
,
S.
,
Green
,
L.
, and
Lee
,
C.
,
1993
,
A Frictional Work Predictive Method for the Initiation of Solid High Explosives From low-Pressure Impacts
,
Lawrence Livermore National Lab.
,
CA
.
79.
Zukas
,
J. A.
,
2004
,
Introduction to Hydrocodes
,
Elsevier
,
Amsterdam; Boston
.
80.
Krien
,
G.
,
Licht
,
H.
, and
Zierath
,
J.
,
1973
, “
Thermochemical Investigation of Nitramines
,”
Thermochim. Acta
,
6
(
5
), pp.
465
472
.
81.
,
M. W.
,
Peterson
,
N. L.
,
Pilcher
,
D. T.
,
Hopkins
,
B. D.
, and
Krier
,
H.
,
1977
, “
Convective Combustion Modeling Applied to Deflagration-to-Detonation Transition of Hmx
,”
Combust. Flame
,
30
(
3
), pp.
231
241
.
82.
Kim
,
B.
,
Choi
,
S.
, and
Yoh
,
J. J.
,
2019
, “
Modeling the Shock-Induced Multiple Reactions in a Random bed of Metallic Granules in an Energetic Material
,”
Combust. Flame
,
210
, pp.
54
70
.
83.
Liu
,
W.-G.
,
2014
,
First-principle Studies of the Initiation Mechanism of Energetic Materials
,
California Institute of Technology
.
84.
Weismiller
,
M. R.
,
Malchi
,
J. Y.
,
Yetter
,
R. A.
, and
Foley
,
T. J.
,
2009
, “
Dependence of Flame Propagation on Pressure and Pressurizing gas for an Al/CuO Nanoscale Thermite
,”
Proc. Combust. Inst.
,
32
(
2
), pp.
1895
1903
.
85.
Bastea
,
S.
, “
Transport Properties of Fluid Mixtures at High Pressures and Temperatures. Application to the Detonation Products of HMX
,”
Proceedings of the 12th International Detonation Symposium
,
Wyndham San Diego at Emerald Plaza
.
86.
Parker
,
G. R.
,
Dickson
,
P. M.
,
Asay
,
B. W.
,
Smilowitz
,
L. B.
,
Henson
,
B. F.
, and
Perry
,
W. L.
,
2006
, “
Understanding the Mechanisms Leading to gas Permeation in Thermally Damaged PBX 9501
,”
Shock Compression of Condensed Matter—2005, Pts 1 and 2
,
845
, pp.
1101
1104
.
87.
Asay
,
B.
,
Parker
,
G.
,
Dickson
,
P.
,
Henson
,
B.
, and
Smilowitz
,
L.
,
2004
, “
Dynamic Measurement of the Permeability of an Explosive Undergoing Thermal Damage
,”
Energ. Mater.
,
22
(
1
), pp.
25
39
.
88.
Koshigoe
,
L.
,
Shoemaker
,
R.
, and
Taylor
,
R.
,
1983
,
Specific Heat of Octahydro-1, 3, 5, 7-Tetranitro-1, 3, 5, 7-Tetrazocine (HMX)
,
Purdue Univ Lafayette in Thermophysical Properties Research Lab
.
89.
Koshigoe
,
L. G.
,
Shoemaker
,
R. L.
, and
Taylor
,
R. E.
,
1984
, “
Specific-Heat of Hmx
,”
AIAA J.
,
22
(
11
), pp.
1600
1601
.
90.
Mcglaun
,
J. M.
,
Thompson
,
S. L.
, and
Elrick
,
M. G.
,
1990
, “
CTH—a Three-Dimensional Shock-Wave Physics Code
,”
Int. J. Impact Eng.
,
10
(
1–4
), pp.
351
360
.
91.
Chidester
,
S. K.
,
Tarver
,
C. M.
, and
Garza
,
R.
,
1999
, “
Low Amplitude Impact Testing and Analysis of Pristine and Aged Solid High Explosives
,”
Proceedings of the Eleventh (International) Symposium on Detonation
,
Snowbird, UT
,
June 27
.
92.
Idar
,
D. J.
,
Straight
,
J. W.
,
Osborn
,
M. A.
,
Coulter
,
W. L.
, and
Buntain
,
G. A.
,
2000
, “
Low Amplitude Impact of Damaged PBX 9501
,”
Shock Compression of Condensed Matter-1999, Pts 1 and 2
,
505
, pp.
655
658
.
93.
Gustavsen
,
R. L.
,
Dattelbaum
,
D. M.
,
Johnson
,
C. E.
, and
Bartram
,
B. D.
,
2013
, “
Experimental Studies of Rod Impact on Bare/Uncovered PBX 9501 Explosive
,”
Procedia Eng.
,
58
, pp.
147
156
.
94.
Gresshoff
,
M.
, and
Hrousis
,
C. A.
,
2010
, “
Probabilistic Shock Threshold Criterion
,”
Proceedings of the 14th International Detonation Symposium
,
Coeur d'Alene, ID
,
Apr. 11
.
95.
Keyhani
,
A.
,
Kim
,
S.
,
Horie
,
Y.
, and
Zhou
,
M.
,
2019
, “
Energy Dissipation in Polymer-Bonded Explosives with Various Levels of Constituent Plasticity and Internal Friction
,”
Comput. Mater. Sci.
,
159
, pp.
136
149
.