1 Online

RESEARCH ARTICLE | AUGUST 27 2018
Multi-scale shock-to-detonation simulation of pressed

energetic material: A meso-informed ignition and growth
model 
O. Sen; N. K. Rai; A. S. Diggs; D. B. Hardin; H. S. Udaykumar 
J. Appl. Phys. 124, 085110 (2018)

https://doi.org/10.1063/1.5046185

View Export
Online Citation
18 September 2024 16:11:18

JOURNAL OF APPLIED PHYSICS 124, 085110 (2018)
Multi-scale shock-to-detonation simulation of pressed energetic material:

A meso-informed ignition and growth model
O. Sen,1 N. K. Rai,1 A. S. Diggs,2 D. B. Hardin,2 and H. S. Udaykumar1,a)
1
Mechanical and Industrial Engineering, The University of Iowa, Iowa City, Iowa 52242, USA
2
Air Force Research Laboratory, Munitions Directorate (AFRL/RW), Eglin AFB, Eglin, Florida 32542, USA
(Received 26 June 2018; accepted 4 August 2018; published online 27 August 2018)
This work presents a multiscale modeling framework for predictive simulations of shock-to-
detonation transition (SDT) in pressed energetic (HMX) materials. The macro-scale computations
of SDT are performed using an ignition and growth (IG) model. However, unlike in the traditional
semi-empirical ignition-and-growth model, which relies on empirical fits, in this work meso-scale
void collapse simulations are used to supply the ignition and growth rates. This results in a macro-
scale model which is sensitive to the meso-structure of the energetic material. Energy localization
at the meso-scale due to hotspot ignition and growth is reflected in the shock response of the ener-
getic material via surrogate models for ignition and growth rates. Ensembles of meso-scale reactive
void collapse simulations are used to train the surrogate model using a Bayesian Kriging approach.
This meso-informed Ignition and Growth (MES-IG) model is applied to perform SDT simulations
of pressed HMXs with different porosity and void diameters. The computations are successfully
validated against experimental pop-plots. Additionally, the critical energy for SDT is computed
and the experimentally observed P2s ss ¼ constant relations are recovered using the MES-IG model.
While the multiscale framework in this paper is applied in the context of an ignition-and-growth
model, the overall surrogate model-based multiscale approach can be adapted to any macro-scale
model for predicting SDT in heterogeneous energetic materials. Published by AIP Publishing.
https://doi.org/10.1063/1.5046185
18 September 2024 16:11:18

I. INTRODUCTION can be performed on samples that approach the macro-scale
size. This calls for developing multi-scale modeling strate-
Predictive modeling of shock-to-detonation transitions
gies for solving shock-initiated detonation problems in real
(SDT) in heterogeneous energetic (HE) materials (plastic-
energetic materials. In such models, the meso-scale physics
bonded explosives, neat pressed HMX etc.) is useful for sev-
is transmitted to the macro-scale through appropriate closure
eral applications, including controlled energy release and
models that appear as source terms in the macro-scale gov-
safe handling of propellant and explosive materials.
erning equations. This work presents a framework for a
Initiation of detonations in heterogeneous energetic materials
meso-informed multi-scale model for pressed energetic
occurs at the meso-scale, i.e., at the scale of the HE grains or
materials, where the meso-scale information is derived from
particles.1 When shock loads are imposed, energy localiza-
high-fidelity reactive void collapse simulations.
tion occurs at defects in the material, which leads to hot-
The present multiscaling framework is developed by
spots; reaction fronts propagate from such hot spots and may
treating the energetic material as consisting of two primary,
strengthen the imposed shock, leading to detonation.2 The
microstructure of heterogeneous energetic (HE) materials is widely separated scales: viz., the macro-scale and meso-
replete with defects at which hotspots can form, including scale. At the macro-scale, all defects, grains, porosities, etc.,
voids, cracks, and crystal-crystal and crystal-binder interfa- are treated as sub-grid or unresolved features. At this scale,
ces. To simulate the macro-scale response of the HEs, the the homogenized material is typically treated as a single
meso-scale dynamics and localization of energy must be cap- phase and macro-models such as ignition-and-growth (IG),8
tured accurately. However, the size of HE samples that are Johnson–Tang–Forest (JTF),9 Scaled Unified Reactive Front
of interest to typical engineering applications is of the order (SURF),10 etc., are used to model the heat release due to
of centimeters or greater.3 Capturing localized meso-scale reactions in the materials induced by shock passage.
events like hotspot formation in macroscopic samples of Alternately, the material can be treated as a multiphase mix-
HEs is computationally infeasible due to stringent spatial ture and modeled using mixture theory as proposed by Baer
and temporal resolution requirements.4,5 For a typical meso- and Nunziato.11,12 In either case, such macro-scale
structure in a millimeter-sized sample of HE, shown in Fig. approaches do not explicitly account for the crucial subgrid-
1,6,7 billions of grid points would be required to adequately scale heterogeneities, but the effects of such heterogeneities
resolve the meso-scale physics. It is therefore unlikely, in the are modeled via closure terms8–12 in the homogenized
foreseeable future, that fully resolved meso-scale simulations macro-scale governing equations. In the original IG model
(as well as in JTF and SURF), the models for heat release are
a)
Author to whom correspondence should be addressed: ush@engineering. phenomenological in nature and contain empirical constants
uiowa.edu that are fit using physical experiments. Recently, there have
0021-8979/2018/124(8)/085110/21/$30.00 124, 085110-1 Published by AIP Publishing.

085110-2 Sen et al. J. Appl. Phys. 124, 085110 (2018)
Deflagration to Detonation (DDT) scenarios, by analyzing the

compaction of granular energetic solids (HMX). In this case,
hotspots are assumed to be generated at contact points
between individual grains of HMX during compaction, i.e.,
void collapse is not involved. Gambino et al.25 use this model
in recent work on a multi-scale simulation of compaction-
initiated detonation. Zhou and co-workers26,27 use meso-scale
simulations to develop a probabilistic approach to estimate
the response of pressed HMX and HMX-based PBX at the
macro-scale. This wide variety of approaches reflects analyti-
cal (Gonthier24), heuristic (Akiki20), and semi-analytical
(Horie,19 Massoni,15 Jackson2,17,18) approaches to link sub-
grid (i.e., meso-scale) physics to the resolved (i.e., macro-
scale) computational predictive models.
In the multiscale models mentioned above and in empir-
ical/phenomenological models, such as ignition-and-growth
FIG. 1. SEM image7 of a pressed HMX sample of size 0.5 mm 0.375 mm. (IG8) JTF,9 SURF,10 statistical hotspot28 models, etc., the
The void diameters in the sample are of the order of a few microns. key quantity of interest (QoI) at the macro-scale is a macro-
scale reaction progress variable k. When k ¼ 0, the material
in a macro-scale control volume (cv) is unreacted explosive
been efforts to develop closure models by utilizing meso- and when k ¼ 1, the reaction is complete and gaseous reac-
scale analyses/simulations;13,14 in silico experiments allevi- tion products are formed. The rate k_ determines the rate of
ate the cost and time demands of physical experiments. energy deposition. In some models,2,15,19 the energy deposi-
The driving objective of multi-scale modeling of ener- _ is placed in the energy equation as a
tion rate e_ meso / kQ
getic materials is to, in some way, incorporate physics at the source term, where Q is a measure of the heat of reaction. In
grain (i.e., meso-) scale into the overall macro-scale other models,8–10,20,29 the progress variable k itself is used to
18 September 2024 16:11:18

response. In general, the primary contribution of meso-scale transition from a solid reactant equation of state to a product
dynamics is energy localization by hotspot formation, equation of state. This carries the material in a control vol-
growth, and interaction. Provided a suitable model for ume from a cold reactant Hugoniot to a reaction product
energy deposition rate due to hotspot formation and growth Hugoniot. The pressure rises with an associated increase in
is delivered to the macro-model, the imposed shock load can internal energy. If the reaction rate k_ is sufficiently high, the
transition to a detonation wave. There have been several energy deposited will strengthen the shock and lead to a
routes taken to model such energy localization. Massoni propagating detonation wave. Thus, the energy deposition
et al.15 developed a model for meso-scale hotspot dynamics rate within each macro-scale control volume, when supplied
based on the spherical void collapse (Carroll-Holt) model of at the appropriate rate, can lead to a macro-scale simulation
Kang et al.16 Such collapse occurs in the situation where of shock-to-detonation transition (SDT).
plastic dissipation plays a major role but the model was also In the present meso-informed ignition and growth
applied to situations where hydrodynamic void collapse is (MES-IG) framework, high resolution meso-scale simula-
expected. Nevertheless, good predictions, i.e., comparisons tions4,5,21–23,30–33 are used to close the macro-scale model.
with experimentally measures pop-plots were obtained. The goal is to infuse meso-structural and meso-scale dynam-
Jackson et al.17 have modeled the coupling between the ics information into the macro-scale via the energy
meso- and macro-models and obtained an energy localiza- deposition rate e_ meso : This link between the meso- and
tion model that in Ref. 18 used a temperature-based macro-scales is provided in the form of surrogate or meta-
Arrhenius type model. This was later17 replaced by an igni- models,13,34–36 in the procedure shown in Fig. 2. The surro-
tion-and-growth type density based model. These models are gates are trained by machine learning approaches using data
based on earlier meso-scale analysis by Jackson et al.2 on from high-resolution meso-scale simulations. No assump-
spherical void collapse. The hotspots were created in the tions on the hotspot temperature, size, or criticality are
solid energetic material by lowering the density at the loca- needed at the macro-scale. Thus, the limitations of the meso-
tion of hotspots resulting in high local temperatures. Hamate scale model are confined to the material models, reaction
and Horie19 include a meso-mechanics model that fits param- rate laws, and numerical resolution and accuracy of the
eters to experimental pop-plots so that hot-spot generated meso-scale simulations. As these models improve, and using
energy deposition of the appropriate magnitude can be used larger ensembles of simulations or (if necessary) higher reso-
at the macro-scale. Akiki and Menon20 present a mechanistic lution simulations, the fidelity of the closure model (e_ meso )
model where the energy deposition rate in a homogeneous can be further improved; in other words, the epistemic
solid material was provided by fitting the meso-scale simula- (reducible) uncertainties37–40 can be reduced. The remaining
tion data of Tran and Udaykumar21,22 and Kapahi and uncertainties then will be of the aleatory (irreducible) type,
Udaykumar23 to extract the time scale and magnitude of heat stemming from inherent stochasticity/variability of the
release due to void collapse. Gonthier24 developed a hotspot meso-structure. Thus, the present work represents a general
generated energy deposition rate formula applicable in approach to modeling heterogeneous energetic materials
18 September 2024 16:11:18

FIG. 2. Schematic representation of the multiscale modeling mechanism. The framework involves performing high-fidelity resolved mesoscale simulations of
reactive void-collapse on sub-samples drawn from a real microstructure of an HE (a) and (b). The numerical experiments are reactive meso-scale computations
_
of collapse of voids as shown in (c). The ensemble of mesoscale computations are used to train a surrogate model to quantify the reaction progress variable (k)
in the macro-scale model using the Modified Bayesian Kriging Method. The macro-scale computational model probes the surrogate on the fly to determine
e_ meso or k_ and transitions the material from a shock to a detonation Hugonoit for predictive modeling of SDTs in HEs.
using a variety of macro-scale models,8–10 while also pre- Macro-scale computations are performed in such a
senting a route to further improve the closure models that homogenized sample, as shown in Fig. 2(e). When reac-
inform the macro-scale. tions ensue, gaseous products are formed and in such
In this paper, the overall multiscale approach is described regions, the homogenized material in a macro-scale con-
in the context of an ignition-and-growth model. The IG trol volume is a mixture of solid and gas.
model, introduced by Tarver8 and coworkers, is a semi- (2) A meso-scale computational model for reactive void-
empirical approach to provide energy deposition at the macro- collapse simulations: At the meso-scale, the energy locali-
scale. The goal of the present work is to incorporate the meso- zation at hotspots is fully resolved by accounting for the
scale physics, particularly hotspot formation and growth, into voids and defects in the energetic material. Calculations
the framework of the ignition-and-growth model. This meso- are performed to follow the reactive void collapse and the
informed ignition and growth (MES-IG) model is described in resulting hot spot evolution. An ensemble of such high-
detail in Sec. II. The MES-IG model is then applied to simu- resolution reactive void collapse calculations,4,5,21–23,30–32
late the shock response of a pressed HMX material in Sec. III. as illustrated in Fig. 2(c), is employed to extract informa-
Comparison of the simulation results with experimentally tion on the ignition and growth rates of hotspots.
obtained pop-plots15,41 validates the overall procedure and is (3) A surrogate model to bridge meso- and macro-scales:
shown in Secs. III B and III C. Conclusions and directions for The energy localization due to hot-spot formation in the
future work are discussed in Sec. IV. subgrid (meso-) scale is transmitted to the macro-scale
through surrogate models13,14,34,35,42 that link the meso-
II. MULTISCALE FORMULATION AND METHODS
scale hotspot physics to the macro-scale response to
The present multiscale modeling framework, illustrated shock loading.
in Fig. 2, comprises three main components to simulate the
The computational models at the meso- and the macro-
shock-to-detonation transition in a pressed energetic material:
scales are presented in Secs. II A and II B. The link between
(1) A homogenized macro-scale model of the porous reac- the meso- and macro-scales is the key contribution of this
tive material: At the macro-scale, the meso-structure of paper, resulting in the MES-IG framework. The procedures
the pressed material is considered to be at the “subgrid” to link the scales are described in Secs. II C–II E.
or unresolved scale. The macro-scale representation of The three components of the overall MES-IG algorithm
the material is a homogenized porous solid material. are discussed in sequence below.
A. Governing equations
(a)
The conservation laws for mass, momentum, and energy
that apply at both meso- and the macro-scales are cast in the
Eulerian form, viz.,
@q @ ðqui Þ
þ ¼ 0; (1)
@t @xi
@ ðqui Þ @ ðqui uj rij Þ
þ ¼ 0; (2)
@t @xj
and (b)
@ ðqEÞ @ ðqEuj rij ui Þ _
þ ¼ E; (3)
@t @xj
where q and ui , respectively, denote the density and the

velocity components, E ¼ e þ 12 ui ui is the specific total
(c)
energy, and e is the specific internal energy. The source term
E_ in Eq. (3) is the rise in specific internal energy of the sys-
tem due to heat released in the decomposition of solid HMX
into gaseous reaction products. The Cauchy stress tensor rij
is of the form
rij ¼ Sij pdij ; (4)
where Sij is the deviatoric stress tensor and p is the pressure.
18 September 2024 16:11:18

FIG. 3. (a) A schematic of the computational set-up for performing meso-
While the above conservation laws are solved at both macro- scale simulations on a sample of pressed HMX material.7 (b) The set up for
and meso-scales, the constitutive relations for Sij and p are performing the macro-scale computations for SDT. At the macro-scale, a
homogenized material, devoid of voids and interfaces, is subjected to a
different at the two scales. These differences in treatment are shock pressure Ps is applied from the west side of the domain boundary. (c)
described in the following. The setup for performing high-fidelity meso-scale simulations of reactive
scv scv
void-collapse to construct the surrogate F_ ignition and F_ growth . A void with
B. Meso-scale modeling of material deformation diameter Dvoid is subjected to a shock pressure Ps , applied from the west
side of the domain boundary.
and reactive mechanics of void collapse
At the meso-scale, high-resolution reactive void collapse
Shock heating can lead to the melting of HMX; there-
calculations are performed, as shown in the setup shown in fore, thermal softening of HMX
Fig. 3(a). The HMX and void spaces are delineated using a is modeled
using the Kraut-
sharp-interface Eulerian framework presented in previous Kennedy relation, Tm ¼ Tm0 1 þ a DV
V0 , with the reference
work.4,5,23,30–32 The collapse of voids due to shock loading melt-temperature, Tm0 ¼ 552 K, and a ¼ 2(C-1/3), where C is
and the formation of hot spots are modeled with the solid the Gruneisen coefficient.47 Following the works of
HMX represented as an elasto-plastic material. Therefore, at Menikoff and Sewell,47 the Gruneisen coefficient is given by
the meso-scale, in addition to the above Eqs. (1)–(3), the C ¼ 1:1 0:2q0 =q, where the reference density, q0 ; is
deviatoric stress tensor, Sij , is evolved in the solid HMX 1905 kg/m3. Once the temperature exceeds the melting point
using the following equation:43–45 of HMX, the deviatoric strength terms are set to zero.
Furthermore, the specific heat of HMX is known to change
@ ðqSij Þ @ ðqSij uk Þ 2 significantly with temperature. The variation of specific heat
þ þ qGDkk dij 2qGDij ¼ 0; (5)
@t @xk 3 is modeled as a function of temperature as suggested by
Sewell and Menikoff.49
where Dij is the strain rate tensor and G is the shear modulus
Since the meso-scale simulations calculate the decom-
of the material. First, the deviatoric stresses are evolved by position of solid HMX into the gaseous products, the individ-
the elastic update above and then mapped back to the yield ual chemical species need to be evolved in time by solving
surface using a radial return algorithm.46 The yield surface is
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the species conservation equation
given by the function f ¼ Se ry , where Se ¼ 32 ðSij Sij Þ.
@qYi
The yield strength ry is taken to be a constant and set to ~Yi ¼ Y_ i ;
þ div qV (6)
@t
260 MPa for HMX,47 i.e., hardening, and visco-plastic
effects are neglected in the meso-scale computational model. where Yi is the mass fraction of the ith species and Y_ i is the
A detailed description of the solution techniques for the production rate source term for the ith species. Chemical
elasto-plastic governing equations and the radial return algo- decomposition of HMX is modeled in the Tarver 3-equation
rithm is provided in previous works.45,48 model50 via four different species. Appendix B contains the
details of the equation of state and reaction kinetics models While at the meso-scale, the heat released due to the
used at the meso-scale. As seen in Appendix B, the solid HMX chemical reaction of HMX is modeled via the source term E_
(species 1, with mass fraction Y1 ) at high temperature decom- in Eq. (3), at the macro-scale, E_ is set to zero. Instead, fol-
poses to fragments (species 2, with mass fraction Y2 ). The frag- lowing the conventional treatment in macro-scale mixture
ments further decompose to intermediate gases (species 3 with formulations,8,54 the chemical heat release due to decompo-
mass fraction Y3 ) and finally, the intermediate gases undergo sition of HMX into gaseous products is accounted for by
strong exothermic reactions to form gaseous products (species transitioning the mixture from a cold, unreacted solid
4, Y4 ), which leads to an increase in temperature. Hugoniot to a product Hugoniot. The last two terms on the
The numerical stiffness in solving the reactive set of right hand side of the Jones-Wilkins-Lee (JWL) equation of
equations [Eq. (6)] is circumvented by using a Strang state [Appendix A, Eq. (A2)] for the hot gaseous products
operator-splitting approach,51 where first the advection of accounts for the chemical energy release due to detonation,
species is performed using the flow time step to obtain pre- as explained in details in Refs. 55 and 56. This standard
dicted species values macro-scale treatment is described in detail in the works of
@qYi n Tarver8 and Kapila,54 among others. The products and reac-
~ Yi ¼ 0:
þ div qV (7) tants in the mixture in each macro-scale control volume are
@t
assumed to be in pressure and temperature equilibrium, i.e.,
In a second step, the evolution of the species mass fraction ps ¼ pg ¼ p; Ts ¼ Tg ¼ T; (10)
due to chemical reactions is calculated
where ps and Ts are the solid pressure and temperature and
dYinþ1
pg and Tg are the pressure and temperature for the gas phase.
¼ Y_ i : (8)
dt The overall specific internal energy and the density for the
The species evolution in Eq. (8) is advanced in time using a control volume are also functions of the specific internal
5th-order Runge-Kutta Fehlberg52 method, which uses an energy and the density of the solid and gas phases, viz.,
internal adaptive time-stepping scheme to deal with the stiff-
e ¼ ð1 kÞes þ keg (11)
ness of the chemical kinetics. Further details on the meso-
18 September 2024 16:11:18

scale reactive mechanics calculations are available in previ- and
ous work.4 1 1k k
The heat released due to the decomposition of solid ¼ þ ; (12)
q qs qg
HMX into the gaseous products is modeled via the 3-equation
model of Tarver et al.,50 as described in Appendix B. The where es and qs are the specific internal energy and the den-
change in temperature because of the chemical decomposition sity for the solid phase and eg and qg are the specific internal
of HMX is calculated by solving the evolution equation energy and density for the gas phase; the equations of state
qCp T_ ¼ Q_ R þ kr2 T ; (9) relating ps , qs , and es as well as pg , qg , and eg are specified
in Appendix A. In Eqs. (11) and (12), k is a macro-scale
where q is the density of HMX, Cp is the specific heat of reaction progress variable which is quantified by the mass-
HMX at constant pressure, T is the temperature, k is the ther- fraction of the reaction products. In the present MES-IG
mal conductivity of HMX, and Q_ R is the total heat release model, k is the key quantity that bridges the meso-scale and
rate because of the chemical reaction. The values of Cp , k, macro-scale energetic response of the material. The transi-
and Q_ R are obtained from the work of Tarver et al.50 The tion of the material in a macro-scale control volume from
source term in Eq. (3) is computed by setting E_ ¼ Cv T_ , solid unreacted HMX (k ¼ 0) to gaseous products (k ¼ 1) is
where Cv is the specific heat of HMX at constant volume. orchestrated through this reaction progress variable. The
Further details on the calculation of reactive void collapse macro-scale system of equations [Eqs. (1)–(4), (10)–(12)] is
using the Tarver 3-equation chemistry model have been closed with an evolution equation for the reaction progress
described in detail in previous works.4,31,53 variable, k. In this work, the reaction progress variable is
evolved using the ignition-and-growth model,8 as follows:
C. Macro-scale material models and reactive
dk _
mechanics: Mixture formulation ¼ k ignition H kignition;max k þ k_ growth H k kignition;max ;
dt
At the macro-scale, the material is considered to be a (13)
homogenized mixture of solid HMX and gaseous reaction where H is the Heaviside function
products. Simulations are then performed in a 1-D setting
corresponding to a flyer impact experiment on a coupon of 0; for Dk < 0
the pressed energetic material. The configuration is illus- H ðDkÞ ¼ (14)
1; for Dk 0:
trated in Fig. 3(b). The mixture is assumed to behave hydro-
dynamically,8,54 i.e., Sij is neglected in comparison to pdij in As in the standard ignition-and-growth model,8 the reaction
Eq. (4). The equations of state for the solid and the gas progress variable in the MES-IG model consists of an initial-
phases as well as the mixture laws for obtaining the homoge- phase ignition component, kignition , and a later-phase growth
nized pressure p are described in Appendix A. component, kgrowth . The switching constant, kignition;max

(selected to be 0.02 in this paper), dictates the transition from _k ignition ¼ Nvoid dMreacted; single void
the ignition to the growth phase. A variant of the ignition-and- (17a)
Mcv dt ignition
growth model57 also incorporates a completion term, given by
kcompletion , which is neglected as a first approximation in Eq. and
(13).
Once k is computed from Eq. (13), the macro-scale sys- _k growth ¼ Nvoid dMreacted; single void : (17b)
tem of equations (1)–(4) can be closed by combining the Mcv dt growth
product and reactant equations of state (described in

Here, Mreacted; single void is the mass of the reaction product
Appendix A) with the mixture rules for the gas and solid
gases formed due to the collapse of a single representative
phases, given in Eqs. (11) and (12).
void. Following the work of Bastea et al.,58 meso-scale
In the MES-IG framework, k_ ignition and k_ growth are the
counterparts to k_ ignition and k_ growth ; i.e., the meso-scale reac-
key quantities that bridge the meso- and macro-scale ener-
tion rate variables F_ ignition and F_ growth , can be defined as
getic response of the material. These reaction rates must be

functions of the meso-scale morphology (grain and defect 1 dMreacted; single void
sizes and shapes) to correctly predict mesostructured- F_ ignition ¼ (18a)
Mvoid dt ignition
dependent SDTs in macroscopic computations. In the present
MES-IG model, the values for k_ ignition and k_ growth are con- and
nected to the energy localization due to reactive void col-
lapse at the meso-scale, as explained next. 1 dMreacted; single void
F_ growth ¼ : (18b)
Mvoid dt growth
D. Obtaining k_ ignition and k_ growth - The MES-IG model Mvoid ¼ p4 q0 D2void is the mass of HMX contained in a void
to connect macro- and meso-scales of diameter Dvoid and q0 is the reference density of solid
HMX.
The macro-scale reaction progress rates in Eq. (13) track Note that by the above definitions F_ ignition and F_ growth
the evolution of ratios of the mass of reaction products to the are product mass fraction generation rates due to the collapse
18 September 2024 16:11:18

mass of the solid HMX in the control volume (cv) and are of a single void in an otherwise homogeneous matrix of
given by HMX, a system that is commonly adopted in meso-scale
simulations of void collapse4,30,32,59,60 The quantities F_ ignition
_k ignition ¼ 1 dMreacted (15a) and F_ growth are obtained at distinct stages of the void collapse
Mcv dt ignition process; in the early post-collapse stage, F_ ignition is obtained,
and at a later stage of mature hotspot evolution, the F_ growth is
and
obtained. The calculation of both these meso-scale rates is
explained in detail in Sec. II E 2.
1 dMreacted
k_ growth ¼ ; (15b) Noting that the porosity / in a macro-scale control vol-
Mcv dt growth
ume is given by the quantity NvoidMMcv void , Eqs. (17) and (18) can
where Mcv is the mass of a macro-scale control volume, and be combined to obtain the following relation between the
Mreacted is the mass of the reaction products formed in the macro-scale rate k_ ignition=growth and the meso-scale rate
macro-scale control volume. The macro-scale rates in Eq. F_ ignition=growth
(15) are connected to the meso-scale reactive dynamics of
k_ ignition ¼ /F_ ignition (19a)
hotspot formation and growth as follows.
For a porous energetic material, such as pressed HMX, and
if the total number of pores or voids in a macroscopic control
volume is Nvoid, then the mass of the product formed from k_ growth ¼ /F_ growth : (19b)
the collapse of all the voids in the control volume is
Thus, k_ ignition and k_ growth in the macro-scale set of equations
X are connected to meso-scale product formation rates F_ ignition
Mreacted ¼ Mreacted; single void ¼ Nvoid Mreacted; single void :
N void and F_ growth , using Eq. (19).
(16) The modeling assumptions and the parameter space for
constructing the surrogate model functions for F_ ignition and
Here, it is assumed that all the voids in the control volume F_ growth are described next.
contribute equally to the local energy deposition rate, i.e.,
that the resulting hotspots are identical. This assumption can E. Calculating the meso-scale reaction progress rates
be relaxed if meso-structural information, such as the distri- F_ ignition and F_ growth
bution of void sizes and shapes, is known, for example from
imaged data. Such a general approach using imaged meso- 1. The model complexity and parameter space
for F_ ignition and F_ growth
structures6 is being pursued in an ongoing work.
Substituting Eq. (16) into Eq. (15), k_ ignition and k_ growth In general, the F_ ignition and F_ growth in Eq. (19) are func-
are obtained from tions of several input/design variables that are applicable at
the meso-scale, such as the loading and the voids/defects uncertainties resulting from this decoupling will need to be
present in the pressed HMX. In the present work, F_ ignition investigated in future work.
and F_ growth are assumed to be of the following form: Detailed descriptions of individual components of the
surrogate models in Eq. (21) are provided by Nassar et al.62
F_ ignition ¼ fignition ðPs ; ss ; Dvoid ; /; AR; hÞ (20a) and Roy et al.,61 where the requisite 584 void collapse simu-
lations have been performed to construct the multi-parameter
and surrogate models shown in Eq. (39). In this paper, however,
F_ growth ¼ fgrowth ðPs ; ss ; Dvoid ; /; AR; hÞ; (20b) a limited surrogate model is constructed with the goal of val-
idating the MES-IG model against the experiments of
i.e., the meso-scale product formation rates are functions of Vanpoperynghe et al.41 The pressed HMX composition in
the applied pressure Ps , the duration of the loading ss , the the experiments41 has a Theoretical Medical Density (TMD)
diameter of the voids Dvoid , the porosity /, the void aspect of 96%, and the global porosity of the sample is reported to
ratio AR, and the void orientation h. The last two quantities be 1.7%.15 Because of the small /, the effect of f ð/Þ in Eq.
apply to non-circular voids aligned at arbitrary orientations (21) can be neglected, an assumption that is supported by the
with respect to the incoming shock. The above functions are void field simulations of Roy et al.61 In addition, in the
each constructed in the form of a surrogate model, using absence of morphological information about the voids in the
methods developed in Sen et al.13,14 with input data obtained sample of pressed HMX in Massoni et al.,15 the present
from ensembles of numerical simulations of void collapse. work assumes circular voids in the material. Therefore, the
As shown in Eq. (20), F_ ignition and F_ growth depend on 6 effects of orientation and aspect ratio of voids in Eq. (21) are
design parameters. Because of the high-dimensionality of the neglected in the present multiscale computations. With these
scv
parameter space in Eq. (20), training a surrogate model using assumptions, Eq. (21) simplifies to F_ ignition ¼ F_ ignition ðPs ; ss ;
meso-scale simulations is computationally prohibitive. If the scv
Dvoid Þ and F_ growth ¼ F_ ðPs ; ss ; Dvoid Þ.
guideline in Sen et al.13 is followed, then at least 8 computa- growth
tions per parameter dimension are necessary to obtain a con- The surrogate construction process, in particular to the
verged surrogate model for F_ ignition and F_ growth . This implies present paper, is further simplified by noting that the dura-
that 86, i.e., more than 200,000 resolved meso-scale compu- tions of the shock loading, ss , are large compared to the
18 September 2024 16:11:18

tations of reactive void collapse are necessary. To overcome mesoscopic time-scales. This is because even a finite ss in
this curse of dimensionality, the complexity of the model in the macro-scale can be approximated as a sustained loading
Eq. (20) is reduced by assuming the following form: duration in the meso-scale. Typical macroscopic ss from
flyer-plate experiments is observed to be in the range of
scv
F_ ignition ¼ F_ ignition ðPs ; ss ; Dvoid Þ fignition
vv shape
ð/Þ fignition ð AR; hÞ 10–200 ns,27 which corresponds to shock pulses of length
25–550 lm (assuming a shock velocity of 2740 m/s for solid
(21a)
HMX). For Dvoid ¼ 0.1–1.0 lm, as reported in Massoni
and et al.,15 such macroscopic pulse lengths are at least an order
of magnitude higher than the meso-scale void-diameters pur-
scv shape
F_ growth ¼ F_ growth ðPs ; ss ; Dvoid Þ fgrowth
vv
ð/Þ fgrowth ð AR; hÞ: sued in the present work. Therefore, the effect of finite pulses
(21b) in the meso-scale can be neglected and it is sufficient to com-
pute F_ ignition and F_ growth by using sustained shocks in the
scv scv
In Eq. (21), F_ ignition and F_ growth are the meso-scale rate for meso-scale numerical experiments. Thus, the meso-scale
the product mass-fraction [Eq. (18)] due to the collapse of a reaction rates of interest to this paper are taken to be func-
single cylindrical void (hence the superscript scv) of diame- tions of Ps and Dvoid only, i.e.,
ter Dvoid , with the subscripts denoting the ignition and the scv
vv vv F_ ignition ¼ F_ ignition ðPs ; Dvoid Þ (22a)
growth phases. The functions fignition ð/Þ and fgrowth ð/Þ are
scv scv
modifiers to F _ and F _ , respectively, and account and
ignition growth
for void-void interactions in a material of porosity /. scv
shape
Similarly, fignition shape
ðAR; hÞ and fgrowth ðAR; hÞ account for the F_ growth ¼ F_ growth ðPs ; Dvoid Þ: (22b)
aspect ratio and the orientation of non-circular voids with
Hereinafter, we consider only single circular voids in calcu-
respect to the shock and are another set of modifiers to
scv scv lating the F_ ignition and F_ growth as expressed in Eq. (22). In par-
F_ ignition and F_ growth . The decomposition of F_ ignition and allel works (Nassar et al.,62 Roy et al.,61 and Rai and
F_ growth in Eq. (21) reduces the number of meso-scale compu- Udaykumar63) a wider parameter range of loading as well as
tations for training the surrogate from 200 000 to 83 þ 8 the void size, distribution, and shape is explored in detail.
þ 82 ¼ 584 simulations. This makes the training process for Those papers describe the physics and quantify the effects of
obtaining F_ ignition and F_ growth computationally tractable. The shock strength, pulse width, arbitrary void shapes, and orien-
price to be paid for the reduction in the number of simula- tations and void-void interactions; they also provide surro-
tions is that the effect of /, i.e., the void-void interactions, gate models for the meso-scale rates F_ ignition and F_ growth as
and that of AR; h, i.e., the void shape, are decoupled from the expressed in Eq. (22). In the following, we restrict attention
effects of the parameters Ps ; ss ; and Dvoid : The epistemic to the hotspot formation and energy deposition rates due to
collapse of single circular voids. The procedure for comput- Figure 5(a) shows that the initial value of Avoid is approxi-
ing F_ ignition and F_ growth from meso-scale simulations of col- mately 0.2 lm2 (Dvoid ¼ 0.5 lm). This value decreases when the
lapse of single circular voids is described next. shock reaches the void boundary and the void starts to collapse.
Once the void has completely collapsed, i.e., at t 0.18 ns, the
2. Computation of F_ ignition and F_ growth from Avoid value reaches zero. Due to the transformation of kinetic
void-collapse simulations energy to thermal energy, a rapid increase in temperature is seen
To obtain F_ ignition and F_ growth from reactive simulations in Fig. 5(b) and a hotspot of area Ahs is formed, as shown in Fig.
at the meso-scale, a sample of solid HMX containing a single 5(c). The reaction is initiated once the hotspot is formed, leading
circular void of diameter Dvoid is loaded with a shock of to accumulation of reaction products.
pressure Ps [shown in Fig. 3(c)]. Figures 4(a)–4(d) show the As shown in Eq. (19), the quantity of interest (QoI) from
temperature distribution in the bulk of a solid HMX, result- the meso-scale is the time derivative F_ of the reacted mass
ing from the collapse of a void of diameter Dvoid ¼ 0.5 lm fraction. By definition, F is the ratio of the mass of reaction
due to Ps ¼ 22.1 GPa. The time variations of the void area products (Mreacted ) to the mass of HMX that would nominally
Avoid , average hotspot temperature Ths , hotspot area Ahs , and fill the void (Mvoid )
ð
product mass fraction F in Fig. 5. The hotspot is defined as
qY4 dA
the region where the local temperature is greater than the Mreacted
F¼ ¼ A ; (24)
bulk temperature that results from shock heating. The aver- Mvoid qHMX Avoid
age hotspot temperature, Ths ; is the Favre average tempera-
ture in the hotspot region, defined as where q is the local density, Y4 is the mass fraction of final
ð gaseous species, qHMX is the density of solid HMX, and Avoid
qTdA is the volume of the void. This way of defining the mass frac-
Ahs tion of reaction products in a hotspot is similar to the one in
Ths ¼ ð ; (23)
Bastea et al.58 who also tracked the hotspot evolution after
qdA
Ahs
void collapse.
The time evolution of the product mass fraction F is
where Ahs is the hotspot area, q is the density, and T is the shown in Fig. 5(c). In the meso-scale computations, the QoI,
18 September 2024 16:11:18

local temperature in the hotspot region. F_ ignition is calculated in the initial phase of the hotspot
FIG. 4. Temperature contours (in K) from a representative meso-scale computation of the collapse of a 0.5 lm void. The simulation performed by loading the
solid HMX material with a pressure of 22.1 GPa. (a), (b), (c) and (d) are the contours at 0.06 ns, 0.18 ns, 2.05 ns and 2.15 ns, respectively.
18 September 2024 16:11:18

scv scv
FIG. 5. Illustration of the procedure for obtaining F_ ignition and F_ growth from the mesoscale numerical experiments. The illustration is shown for Dvoid ¼ 0.5 lm,
subjected to a shock pressure of Ps ¼ 22.1 GPa. The figures show the variation of (a) void area, (b) hotspot average temperature, (c) hotspot area, and (d) the
products mass fraction F with time. The numbers 1 through 4 in (a)–(d) indicate different instances of void and hot-spot formation/growth. In (d), the time
scv scv
instances selected for computing the slopes to obtain F_ ignition and F_ growth are shown.
formation, while F_ growth is computed in the later phase of hot- reaches twice its value at point 2. The differences in F and
spot growth. Figure 5 shows the procedure to delineate the ini- time values for these points (3 and 4) are used to calculate
tial ignition and later growth phases to obtain F_ ignition and the growth rate F_ growth as
F_ growth . First, the void area variation with time [Fig. 5(a)] is
used to locate two points in the evolution of the void. The first F4 F3
F_ growth ¼ : (26)
point (instant 1 in Fig. 5) is when Avoid ¼ 90% of its original t4 t3
(undeformed) area. The second point (instant 2 in the figure)
The calculation of the ignition and growth rates at the meso-
is when the void area equals zero—i.e., when the void has
scale, as described above, is performed for a range of shock
totally collapsed. Points 1 and 2 are used to demarcate the
loading conditions and void diameters relevant to the experi-
ignition phase. The difference between the corresponding val-
mental pop-plots in Massoni et al.15 The resulting functions
ues of the product mass fraction F and the time values for
F_ ignition ðPs ; Dvoid Þ and F_ growth ðPs ; Dvoid Þ are used to calculate
those two points are used to find the ignition rate F_ ignition as
the macro-scale ignition and growth rates k_ ignition and k_ growth
F2 F1 in Eq. (39). The overall solution algorithm for meso-
F_ ignition ¼ : (25) informed SDT simulations is described next.
t2 t1
F. Numerical procedure for performing meso-informed
The hotspot growth rate is calculated using points 3 and 4.
SDTs in pressed HMX
Point 3 is when the hotspot area equals 1.8 its value at point
2 (recall that point 2 is when the void has totally collapsed Meso-informed computations of SDT in pressed HMXs
and the hotspot is ignited). Point 4 is when the hotspot area are performed in two stages in the MES-IG framework, in a
hierarchical multi-scale setting. The first stage involves con- described in Sen et al.13,14 The surrogates are constructed
struction of surrogate models for F_ ignition ðPs ; Dvoid Þ and for both Dvoid ¼ 0.5 and 1 lm and are supplied to the
F_ growth ðPs ; Dvoid Þ. Once the construction of the surrogate is macro-scale computational code to perform SDT in
complete, in the second stage, the surrogate models are used pressed HMX. The surrogate functions are plotted in Fig.
to compute k_ ignition and k_ growth to perform macro-scale simu- 8 for the values of Dvoid ¼ 0.5 and 1 lm.
lations of SDT. The two steps are described below.
2. Numerical procedure for performing SDT simulation
1. Construction of surrogate models for F_ ignition and
at the macro-scale
F_ growth
At the macro-scale, 1-D simulations of shock-loaded
The procedure for constructing the surrogate models for
pressed HMX samples are performed to construct a pop-plot
F_ ignition and F_ growth is as follows.
to compare with experimental data.41 To perform SDT simu-
1. First, the parameter range for Dvoid is estimated. Because the lations using the MES-IG model, the macro-scale flow field
goal of the present paper is to validate the MES-IG model is initialized with the flow variables [q,ui,e,k] ¼ [q0,0,e0,0]
against the experiments of Vanpoperynghe et al.,41 Dvoid is and loaded with an input pressure of Ps . A set-up for a typi-
selected to be comparable to those of Vanpoperynghe cal 1-dimensional macroscopic simulation is shown in Fig.
et al.41 The samples of pressed HMXs in the experiments41 3(b). The governing equations are solved in the following
are reported to have voids with a maximum diameter of 5 sequence to evolve the flow variables with time.
lm. However, based on electron microscopy images, it was
1. First, the conservation equations for mass, momentum,
shown that the samples have a large number of smaller
and energy, Eqs. (1)–(3), are solved to compute the flow
voids, with sizes below 1 lm. Massoni et al.15 validated the
variables, [q,ui,e]. The conservation laws are spatially dis-
computational model against the experiments of
Vanpoperynghe et al.41 using void diameters lying in the cretized using the 3rd-order essentially non-oscillatory
range of 0.2–0.5 lm, with a few 5 lm voids. In the present (ENO) scheme64–66 and numerically integrated using the
paper, the void sizes are assumed to lie in the range identi- 3rd-order Runge-Kutta time stepping scheme.67
fied by Massoni et al.15 and the surrogates are constructed 2. Following this, the mixture pressure p, and the densities
18 September 2024 16:11:18

for Dvoid ¼ 0.5 and 1 lm. Therefore, for the purpose of vali- of the gas and solid phases, qg and qs , are obtained by sat-
dation against experiments, instead of the two-dimensional isfying the P-T equilibrium conditions and the mixture
laws for internal energy and phase densities, by solving
surrogates, F_ ignition ðPs ; Dvoid ) and F_ growth ðPs ; Dvoid ), four
Eqs. (12), (A13), and (A14) in Appendix A iteratively
one-dimensional surrogates, viz., F_ ignition ðPs ÞjDvoid ¼0:5; 1:0lm using the Newton-Raphson method.67
and F_ growth ðPs ÞjD ¼0:5; 1:0lm are used in the present work.
void
3. Finally, the updated mixture pressure p is used to compute
2. Second, the parameter range for Ps is selected. The mini- the values of k_ ignition and k_ growth . For this, each macro-
mum value of Ps is selected to be 5 GPa. This is because, scopic control volume is assumed to correspond to a reali-
for pressures below 5 GPa, F_ ignition and F_ growth are found zation of a meso-scale numerical experiment. Therefore,
to be zero for Dvoid ¼ 0.5 and 1 lm. For lower pressures, the macro-scale mixture pressure p is the pressure for
the hot-spot formed due to void-collapse is quenched by which the energy localization at the meso-scale is to be
diffusion and does not lead to sustained hot-spot growth calculated. The surrogate model [Eq. (22)] is probed at
and species formation in the meso-scale. Therefore, below ðPs ¼ p; Dvoid Þ, to obtain the values of F_ ignition and
5 GPa the meso-scale product formation rate is zero and F_ growth , which are then used to calculate k_ ignition and
the lower limit of Ps is set to be 5 GPa for the surrogate k_ growth using Eq. (19). The computed k_ ignition and k_ growth
model. The upper limit for Ps is selected to be the Von- are employed in the reaction evolution model, Eq. (13), to
Neumann pressure, which is estimated to be 60 GPa from update the reaction progress variable k. The equation is
the simulations of Massoni et al.15 solved using a forward Euler method,67 with a sub-time
3. Next, ensembles of reactive mesoscale computations of step of Dt ¼ 0:001Dtflow , where Dtflow is the time-step for
void-collapse are performed for Dvoid ¼ 0.5 and 1 lm the Runge-Kutta scheme for solving Eqs. (1)–(3).
using the computational model described in Secs. II A and 4. The updated values of k and p are used in the next time
II B. The simulations are performed for varying values of step to evolve the mixture variables with time, returning
Ps , ranging from 5 through 60 GPa, at 5 GPa intervals. to step 1.
This amounts to performing 12 reactive void collapse
The above procedure is carried out to evolve the macro-
simulations for each of Dvoid ¼ 0.5 and 1 lm. A represen-
scale fields for a desired time duration, tend, which is set to
tative meso-scale simulation of the various stages of col-
be the time when a steady Chapman-Jouguet (CJ) structure
lapse of a void of 0.5 lm is shown in Figs. 4(a)–4(d). The
is observed in the material.
quantities F_ ignition and F_ growth are computed from each of
the simulations using Eqs. (25) and (26), as described in
III. RESULTS AND DISCUSSION
Sec. II D and also shown in Fig. 5.
4. The F_ ignition and F_ growth computed from the meso-scale The methods discussed in Sec. II are used to perform
simulations are used as inputs to train a Modified SDT simulations in a macro-scale sample of pressed HMX;
Bayesian Kriging method (MBKG) using the methods 1D simulations are performed assuming a coupon of material
loaded with a sustained planar shock of strength Ps , as shown experiments of Vanpoperynghe et al.?15,41 Second, is it possi-
in Fig. 3(b). The goal is to validate the MES-IG calculation ble to recover the Walker-Wasley68 relations, given by (27),
of sensitivity of pressed HMXs to shocks at the macro-scale. in the MES-IG framework? To systematically address the
Experimentally, sensitivity of HEs is characterized in two above questions, first the macro-scale computational code is
ways: verified using a standard IG model against the shock-to-deto-
nation simulations of LX-17 of Kapila et al.54 Following this,
(1) Measuring run-to-detonation distances.3,8
the standard IG model is replaced by the MES-IG model for
(2) Measuring the minimum energy needed for SDT in a HE
pressed HMX and is used to validate the run-to-detonation
material.27,68
computations against the experiments on pressed HMX.15,41
In the first approach, i.e., experiments that measure run- Finally, instead of sustained loadings, macroscopic pulses of
to-detonation distances, the HE material is loaded with a sus- finite duration are used to compute k and recover the Walker-
tained shock of a given strength. As the shock travels through Wasley relations for pressed HMX. The effect of void sizes
the material, it transitions to a detonation wave. The run-to- and porosities on k is demonstrated in the calculations.
detonation distance h traveled by the shock before transition-
ing to a detonation wave is measured for different input pres-
A. Verification of the macro-scale computational code
sures, Ps. In general for a given Ps, h will be different for
different HE materials. The plot of Ps vs h (i.e., the pop-plot3) To verify the current macro-scale computational code
is used to characterize the sensitivity of the HE material. against numerical simulations of Kapila et al.,54 one-dimensional
In the second approach, i.e., experiments that measure simulations of shock-to-detonation transition in LX-17 are per-
the minimum energy needed for SDTs, the HE material is formed. The material is loaded with a sustained velocity Up, at
impacted by a flyer of a thickness w. This produces a shock its west end. Similar to Kapila et al.,54 the solid reactant and the
of a finite duration, ss in the HE material. The loading dura- gaseous products are modeled using their respective JWL equa-
tion, ss , depends on the thickness of the flyer w. Based on ss , tions of state, while the reaction-progress variable is evolved
and the shock pressure Ps, James69 introduced a critical using the standard ignition-and-growth model. The constants for
energy threshold, Ecr defined as follows: the JWL equations of state and the IG models are the same as
those reported in the literature.54 The temporal evolution of the
18 September 2024 16:11:18

l
P2s ss ¼ Ecr ¼ k; (27) pressure in the material is calculated for mesh sizes of 10, 25,
A
and 50 lm to establish grid-independence.
where l and A are the impedance and the area of the flyer, For Up ¼ 1.45 km/s, the long-term CJ structure in the
respectively. If the energy of the loading is above Ecr, then material is shown in Fig. 6(a). The figure shows that the von
the input shock will transition to a steady detonation wave. Neumann pressure obtained from the present meso-scale
In general, the critical energy Ecr and therefore k depends on computations is 34.8 GPa, which is the same as that reported
the meso-structure of the HE material27 and on the specific in the literature.54 Figure 6(a) compares the CJ structure of
species of HE material70 as sensitivities of HE can vary over the detonation wave with those in the works of Kapila
a wide range. et al.54 at 1.1 ls, while Fig. 6(b) shows the evolution of the
To predict the sensitivity of a HE subject to shock- pressure distribution along the sample for t ¼ 0.4 ls to 1.2 ls.
loading using the MES-IG model, two questions are posed. As shown in Fig. 6(a), the results converge for mesh sizes
First, can the MES-IG framework be used to validate the run- below 25 lm and are in good agreement with those by
to-detonation distances of pressed HMXs against the Kapila et al.54
FIG. 6. (a) Comparison of the pressure distribution over a 1-Dimensional LX-17 sample at t ¼ 1.1 ls with the results obtained from Ref. 54. A sustained veloc-
ity of 1.45 km/s is maintained at the left end of the sample (x ¼ 0). The computations are performed for three different mesh sizes, as indicated in the figure. (b)
Evolution of pressure with time in a 1 dimensional sample of LX-17. A sustained velocity of 1.45 km/s is maintained at the left end of the sample (i.e., at
x ¼ 0). The pressure distributions are plotted for t ¼ 0.4 ls to 1.2 ls at an interval of 0.2 ls. The mesh size is 10 lm.
FIG. 7. Evolution of pressure with time in a 1 dimensional sample of LX-17. A sustained velocity of 1.07 km/s is maintained at the left end of the sample (i.e.,
at x ¼ 0). For (a) the pressure distributions are for t ¼ 0.05 ls to 1.8 ls, at an interval of 0.05 ls until t ¼ 0.4 ls and 0.3 ls for t ¼ 0.6 ls to 1.8 ls. In (b), the
pressure distributions are shown for t ¼ 2.5 ls to 4.5 ls at an interval of 0.5 ls. The mesh size is 10 lm for all the cases shown.
To further validate the macro-scale computational code, from resolved meso-scale simulations of reactive void-
a case with lower input velocity of Up ¼ 1.07 km/s is simu- collapse in HMX. The variations of the F_ ignition and F_ growth
lated. As reported in the literature,54 two distinct stages of with Ps and Dvoid are discussed next.
evolution of the incident shock-front are noted in Fig. 7. In
the first stage, shown in Fig. 7(a), the lead shock travels
1. Surrogate models for ignition and growth rates
through the material without being strengthened by a reac- from meso-scale computational experiments
tion front. Instead, a secondary wave builds up behind the
lead shock and collides with the shock at about x ¼ 5 mm. As shown in Eq. (22), F_ ignition and F_ growth are functions of
18 September 2024 16:11:18

This forms a steady CJ-wave in the second stage, which the shock pressures and the void diameters. In the present paper,
thereafter propagates downstream [as shown in Fig. 7(b)]. the void sizes are assumed to lie in similar ranges as Massoni
Thus, the two stages leading to the formation of a steady CJ et al.15 and the following simulations for validating the MES-IG
wave are clearly identified in the present macro-scale com- model are limited to voids of diameters 0.5 and 1 lm.
putations. Therefore, the current computations not only pre- Surrogates for F_ ignition and F_ growth covering a wider range of
dict the correct von Neumann pressure and the CJ structure, loading conditions and void diameters are shown in the works
but also capture the dynamics leading to the formation of a of Nassar et al.,62 Roy et al.,61 and Rai and Udaykumar.63
steady CJ wave in SDT scenarios. The surrogates for F_ ignition and F_ growth with pressure are
shown in Fig. 8 for two void diameters, Dvoid ¼ 0.5 and
1.0 lm. As expected, F_ ignition and F_ growth increase monotoni-
B. Validation of the MES-IG model – run-to-detonation cally with pressure. Figure 8 shows that below 25 GPa,
distances in pressed HMX
F_ ignition and F_ growth are higher for Dvoid ¼ 1.0 lm than for
In the MES-IG model, the ignition and growth rates for Dvoid ¼ 0.5 lm. However, above 25 GPa, the F_ ignition
pressed HMXs are obtained using surrogates constructed and F_ growth for Dvoid ¼ 0.5 lm voids exceed those for
FIG. 8. Ignition and growth rates computed from mesoscale numerical experiments of the collapse of single voids of diameters 1.0 and 0.5 lm over a range of
input pressures. The pressure values at which the mesoscale computational experiments were performed are denoted by solid (for Dvoid ¼ 0.5 lm) and hollow
(for Dvoid ¼ 1 lm) circles.
Dvoid ¼ 1.0 lm. Thus, the larger (1 lm) void leads to a faster detonation transitions in pressed HMXs. The first computa-
growing hotspot at lower pressures, while the smaller tion is performed for a material with / ¼ 3% and
(0.5 lm) void hotspot grows faster at higher pressures. This Dvoid ¼ 0.5 lm. The F_ ignition and F_ growth for Dvoid ¼ 0.5 lm
behavior is aligned with the experimental and theoretical are used to obtain the macro-scale reaction rates k_ ignition and
observations of Price,71 explained in further detail by Nassar k_ growth , using Eq. (19). A sustained loading of 12 GPa is
et al.,62 and supported with scaling arguments by Rai and maintained at the west end of the domain. The flow variables
Udaykumar.63 Briefly, the observed dependency of F_ ignition are allowed to evolve to t ¼ 1.0ls when a steady CJ structure
and F_ growth on the void size is because, for collapse at lower is observed.
pressures, the hotspot temperature is lower and thermal dif- The temporal evolution of the pressure distribution in
fusion is more effective in containing the growth of the the material is shown in Fig. 9. As shown in the figure, there
smaller (Dvoid ¼ 0.5 lm) hotspot. At higher pressures, the are two stages of the shock-to-detonation transition in the
temperatures of the hotspots are high for both the void diam- material. In the initial stage, i.e., until the shock travels a dis-
eters; therefore, the growth of the hotspot due to chemical tance of 0.5 mm, the incident shock strengthens slowly. This
reactions overwhelms the effects of thermal diffusion, even is because the F_ growth term is low up to a pressure of 25 GPa
for the smaller void. There is a cross-over in the sensitivity (as discussed in Sec. III B 1). At later stages, after the pres-
of voids as the imposed shock pressure increases, as pointed sure in the material exceeds 25 GPa, the shock strengthens
out in previous work.71 The implications of this crossover faster and a von Neumann spike is developed at x ¼ 2.5 mm.
for the run-to-detonation and critical energy behavior are The transition from shock to detonation is also observed in
examined below. the x-t plots in Fig. 9(b). Until a steady CJ structure is
formed, the slope of the shock locus in the x-t plot progres-
2. Meso-informed macro-scale computations for sively changes with time. After x ¼ 2.5 mm (or equivalently
pressed HMXs with Dvoid 5 0.5 lm and / 5 3%
t ¼ 0.4 ls), the von Neumann pressure is reached and a
The F_ ignition and F_ growth surrogates of Sec. III B 1 are steady CJ structure is formed. The slope of the shock-locus
used to perform macro-scale simulations of shock-to- in the x-t plot [Fig. 9(b)] remains constant thereafter,
18 September 2024 16:11:18
FIG. 9. (a) Evolution of the pressure with time for a 1-dimensional sample of pressed HMX comprising voids of diameter 0.5 lm with porosities of / ¼ 3%
and 1.7% (shown in solid and dashed-dotted lines, respectively) for an input pressure of 12 GPa. The pressure distributions are for t ¼ 0.02 ls to 0.68 ls
recorded at an interval of 0.08 ls until t ¼ 0.29 ls and an interval of 0.06 ls thereafter. Figures (b) and (c) show the x-t plots for the pressure in GPa for the
samples of pressed HMX comprising voids of diameters of 0.5 lm and with porosities of 3% and 1.7%, respectively. The mesh size is 10 lm in the simulations.
In (b) and (c), h denotes the run-to-detonation distances of the respective samples.
indicating that a steady detonation front velocity is reached. the pressed HMX material was also observed in Ref. 3,
Since the von Neumann pressure is first attained when the where the sensitivity of the material was shown to increase
shock-front reaches x ¼ 2.5 mm, the run-to-detonation dis- with /.
tance h is taken to be 2.5 mm in this case. In subsequent
cases, the x-t plots are used to locate the value of h in the b. Effect of Dvoid. To study the effect of void diameter
same way. Dvoid on the run-to-detonation distances, the porosity of the
material is fixed at / ¼ 3%, and two values of Dvoid are stud-
3. Effect of porosity /, and void-diameter Dvoid , ied, viz., 1 lm and 0.5 lm. The material is loaded with a sus-
on run-to-detonation distances in pressed HMX tained pressure of 12 GPa and the run-to-detonation
a. Effect of /. To study the effect of porosity on the run- distances are compared for the two void sizes.
to-detonation distance, two values of / ¼ 3% and 1.7% are Figure 10 compares the SDT in HMX with Dvoid
studied, while maintaining Dvoid ¼ 0.5 lm. Similar to Sec. ¼ 0.5 lm and 1.0 lm. Figure 10(a) shows that there are two
III B 2, a sustained pressure of 12 GPa is applied at x ¼ 0. distinct stages in SDTs for the two void-diameters. The ear-
The flow variables are allowed to evolve to t ¼ 1.0ls, i.e., lier stage is until the shock-front reaches x ¼ 2 mm. In this
until a steady CJ wave is observed. stage, the lead shock strengthens faster for Dvoid ¼ 1.0 lm
The process of transition of the shock to a detonation than for Dvoid ¼ 0.5 lm. In the later stages, i.e., once the pres-
wave is compared for / ¼ 3% and 1.7% in Fig. 9. The figure sure exceeds 25 GPa, the growth rate is higher for
shows that run-to-detonation distance h decreases as porosity Dvoid ¼ 0.5 lm than for Dvoid ¼ 1.0 lm. The explanation for
increases; h ¼ 2.5 mm for / ¼ 3%, while it is 4.1 mm for increased shock strengthening for Dvoid ¼ 0.5 lm after
/ ¼ 1.7% This is because, for a constant Dvoid , k_ ignition and 25 GPa follows from the discussion in Sec. III B 1. As was
k_ growth are directly proportional to /, as shown in Eq. (19). discussed, at lower pressures, F_ ignition and F_ growth are higher
Therefore, as observed in Fig. 9(a), the reaction front grows for Dvoid ¼ 1.0 lm than for Dvoid ¼ 0.5 lm. After 25 GPa,
faster and the von Neumann pressure is reached earlier for F_ ignition and F_ growth for Dvoid ¼ 0.5 lm exceed those of
/ ¼ 3% than for / ¼ 1.7%. This can also be observed in the Dvoid ¼ 1.0 lm. Therefore, the run-to-detonation distance is
x-t plots in Figs. 9(b) and 9(c). The decrease in h with / in smaller for Dvoid ¼ 0.5 lm than that for Dvoid ¼ 1.0 lm.
18 September 2024 16:11:18
FIG. 10. (a) Evolution of the pressure with time for a 1-dimensional sample of pressed HMX with porosity, / ¼ 3% and void diameters of 1 and 0.5 lm (shown
in dashed and solid lines, respectively) for an input pressure of 12 GPa. The pressure distributions are for t ¼ 0.02 ls to 0.68 ls recorded at an interval of 0.08
ls until t ¼ 0.29 ls and an interval of 0.06 ls thereafter. Figures (b) and (c) show the x-t plots for the pressure in GPa for the samples of pressed HMX compris-
ing voids of diameters of 1 and 0.5 lm, respectively, and a porosity of 3%. The mesh size is 10 lm for the simulations. In (b) and (c), h denotes the run-to-deto-
nation distances of the respective samples.
4. Comparison of the computed run-to-detonation experimental sample,15 the current MES-IG predictions of
distances with experimental observations sensitivity provide good agreement with experiments.
To compare the run-to-detonation distances with the
C. Characterizing the sensitivity of pressed HMXs
experiments,41 a sustained pressure is applied to the material
for pulses of finite duration
at x ¼ 0. The pressure is varied from 8 GPa to 20 GPa. Three
combinations of void diameters, Dvoid , and porosity / are In this section of the paper, the current MES-IG frame-
chosen: (/, Dvoid ) ¼ (3%, 0.5 lm), (3%, 1.0 lm) and (1.7%, work is used to predict the constant k in Eq. (27) to obtain
0.5 lm). The run-to-detonation distance h is taken to be the the critical energy criterion for pressed HMX. For this pur-
distance at which the von Neumann pressure is first recorded pose, pulses of finite duration (ss ) are applied to the material.
in the material, as depicted in the x-t plots. The pop-plot First, the effect of a subcritical ss (i.e., a loading duration
obtained from the computations is compared with experi- which does not lead to detonation) is compared with that of a
mental data.15,41 The pop-plot in Fig. 8 also shows the exper- super-critical ss for a fixed Ps. Following this, Ps and ss are
imentally obtained fitted line.15,41 both varied to construct the criticality envelopes and recover
Figure 11 shows the run-to-detonation distances com- the Walker-Wasley relations for pressed HMX.
puted from the MES-IG model for the three parameter sets.
As observed in the figure, the present computations are in 1. Critical loading conditions for pulses of finite
good agreement with experiments.41 Furthermore, the slope thickness for a given pressed HMX sample
of the h vs Ps plot matches the slope obtained from experi-
To study SDTs under loads of finite duration, a shock of
ments for (/, Dvoid ) ¼ (3%,0.5 lm). The slope is also similar
pressure 12 GPa is applied at the west end of the pressed
for the lower porosity case for the 0.5 lm void, i.e., for (/,
HMX sample with (/, Dvoid ) ¼ (3%, 0.5 lm). The time evo-
Dvoid ) ¼ (1.7%,0.5 lm). However, because h increases with
lution of the shock in the material with two values of pulse
decreasing porosity, the (/, Dvoid ) ¼ (1.7%,0.5 lm) cases
duration, ss ¼ 70 and 80 ns, is compared with that for a sus-
have higher run-to-detonation distances than the experi-
ments. In contrast to Dvoid ¼ 0.5 lm, the slope of the h vs Ps tained loading (i.e., ss ! 1). As shown in Figs. 12(a)–12(c),
plot is different from the experimental slope for the short pulse shock (ss ¼ 70 ns) results in a no-go condi-
tion, i.e., the shock does not transition into detonation.
18 September 2024 16:11:18

Dvoid ¼ 1 lm. As seen in Fig. 11, at lower pressures (Ps less
than 12 GPa), h is lower than the experimental values. In Although at the initial stages, the lead shock is slightly
contrast, at higher pressures, h is higher than the experi- strengthened [Fig. 12(a)] because of increasing k [Fig.
ments. At 12 GPa, the computational and the experimental 12(c)], the rarefaction front catches up and quenches the
values agree with each other. These observations follow lead-shock before the von Neumann pressure is reached.
from the discussions in Sec. III B 3, where it was shown that This is also observed in the x-t plots shown in Fig. 12(b).
the reaction rates are higher at lower pressure for Therefore, ss ¼ 70 ns results in a no-go condition, i.e., no det-
Dvoid ¼ 1 lm, but do not increase significantly at higher pres- onation is observed in the material.
sures. Therefore, h is lower than the experimental values for In contrast to ss ¼ 70 ns, the somewhat longer pulse of
d ¼ 1 lm. Finally, it is noted that in Fig. 11, for all three ss ¼ 80 ns results in a steady detonation wave, as observed in
cases of (/, Dvoid ) ¼ (3%, 0.5 lm), (3%, 1.0 lm), and (1.7%, Figs. 12(d)–12(f). A comparison of Figs. 12(c) and 12(f)
0.5 lm), the computed run-to-detonation distances are in shows that k increases faster for ss ¼ 80 ns than for
close agreement with the experimental values.41 Given the ss ¼ 70 ns. The build-up of the reaction front overcomes the
uncertainty in the porosity and void size values in the rarefaction effect in Fig. 12(d), and the von Neumann pres-
sure is reached at x ¼ 4.3 mm. As expected, the run-to-deto-
nation distances are lower for sustained pulses than for finite
pulse thicknesses. In the limit of a sustained loading (i.e., ss
! 1), h ¼ 2.5 mm, as seen in Fig. 12(g). For a sustained
loading, there is no rarefaction to compete with the reactions,
so that the reaction front grows faster than for ss ¼ 80 ns, as
seen by comparing Figs. 12(f) and 12(i).
A summary of the above observations is shown in Fig.
12(j), which compares the locus of the shock in the material
for ss ¼ 70 and 80 ns and for sustained loadings. As observed
in Fig. 12(j), the macro-scale go and no-go conditions are
successfully delineated using the present MES-IG model for
pressed HMX.
2. Critical energy thresholds for pressed HMX

FIG. 11. Comparison of the run-to-detonation distances of meso-scale
informed macroscopic simulations of pressed HMX samples, with those
for pulses of finite thickness
obtained from experiments.15,41 The experimental fit is shown as a solid To obtain the critical energy envelope, pulses of finite
line. The diamonds denote pressed HMXs with / ¼ 3% and Dvoid ¼ 1 lm.
The circles denote samples with Dvoid ¼ 0.5 lm and / ¼ 3% (solid circles) duration are applied at the west end of the 1D sample of
and 1.7% (hollow circles). pressed HMX. The length of the samples is taken to be
18 September 2024 16:11:18
FIG. 12. The effect of the pulse duration (ss ) on shock-to-detonation transition in samples of pressed HMX with Dvoid ¼ 0.5 lm and / ¼ 3%, for an input pres-
sure of 12 GPa. Figures (a)–(c) are for ss ¼ 70 ns, (d)–(f) are for ss ¼ 80 ns, while Figures (g)–(i) are for a sustained shock. Figures (a), (d), and (g) show the
pressure distribution, and Figures (c), (f), and (i) show the reaction progress variable, k, along the samples until t ¼ 1 ls, recorded at an interval of 0.05–0.06
ls. Figures (b), (e), and (h) are the x-t diagrams for the pressure distribution (in GPa) in the respective samples. Figure (h) is the time-history of the locus of
the shock front for the three values of ss as indicated in the legend. The mesh size for all the simulations is 10 lm.
10 mm. Similar to Sec. III B 4, three combinations of void critical condition is plotted against the input pressure to
diameters, Dvoid , and porosity / are chosen: (/, obtain the criticality envelope for pressed HMX.
Dvoid ) ¼ (3%, 0.5 lm), (3%, 1.0 lm) and (1.7%, 0.5 lm). The Figure 13 shows the criticality envelope obtained using
shock pressure is varied from 10 to 20 GPa, in 1 GPa incre- the MES-IG model. For all three combinations of (/, Dvoid ),
ments, while ss is varied from 10 to 200 ns in 10 ns incre- the Walker-Wasley relationship given by (27) is obtained
ments. For a fixed pressure, the minimum ss , at which the using the MES-IG model. The constants k are 0.48, 1.22, and
von Neumann pressure is attained in the material, is taken as 2.03 GW/cm2 for (/, Dvoid ) ¼ (3%, 1.0 lm), (3%, 0.5 lm),
the critical loading condition. The ss corresponding to the and (1.7%, 0.5 lm), respectively. As shown in Fig. 13, for
Dvoid ¼ 0.5 lm, k increases with decreasing porosity. In other F_ ignition and F_ growth values are lower than for the 1 lm
words, the sensitivity of the material increases with increas- void (see Fig. 8). Therefore, in the earlier stages of
ing /, which was also seen in the pop-plots discussed in Sec. SDT, the rarefaction effects suppress the reaction
III B 4. growth rates and the reaction front grows only gradually
As discussed earlier, both the pop-plot and the critical
energy criteria show that the sensitivity of the material
increases with /. However, the pop-plots and the critical
energy criteria show opposite trends in sensitivity with
respect to Dvoid . While Fig. 13 shows that k (or equivalently
Ecr ) decreases with Dvoid , the pop-plot (Fig. 11) shows that h
increases with Dvoid . This contradiction arises because k is
computed using pulses of finite duration, whereas h is com-
puted for sustained shock-loadings. Once the durations of the
pulses are taken into account, both pop-plots and the critical
energy criteria predict similar trends of sensitivity with
respect to Dvoid —as will be shown next.
To study the effect of ss on h for voids of different diam-
eters, the porosity of the material is fixed at 3% and the cal-
culations are performed for the two values of Dvoid ¼ 1 and
0.5 lm. The samples are loaded with 12 GPa shocks and
ss ¼ 80, 90, and 100 ns. The shocks are allowed to evolve
and the run-to-detonation distances (h) are computed for
each case. The results from the simulations are shown in Fig.
14, where the following three scenarios are observed:
(a) For ss ¼ 80 ns, (i.e., for the short pulse case), the von
18 September 2024 16:11:18

Neumann pressure is attained earlier for Dvoid ¼ 1 lm
case than for Dvoid ¼ 0.5 lm case, as observed in Fig.
14(a). The figure shows that the reaction front grows
uniformly for Dvoid ¼ 1 lm; in contrast, the growth of
the reaction front is slower for Dvoid ¼ 0.5 lm and is
non-uniform. In the earlier stages, the shock grows
slowly, while in the later stages, the growth is faster.
This is because, below 25 GPa, the F_ ignition and F_ growth
are higher for Dvoid ¼ 1 lm, as discussed in Sec. III B 1
(and seen from Fig. 8). The chemical energy release
rate overcomes the rarefaction effects, even when the
shock pressure is below 25 GPa. Therefore, the von
Neumann pressure is achieved earlier for Dvoid ¼ 1
lm. In contrast, for Dvoid ¼ 0.5 lm, below 25 GPa the
FIG. 13. Criticality envelope for pressed HMX samples for shock-to-detona-
tion, showing the transition from no SDT conditions to SDT conditions for a
values of pressure and ss . The diamonds denote pressed HMXs with / ¼ 3% FIG. 14. Evolution of the pressure with time along a specimen of pressed
and Dvoid ¼ 1 lm. The circles denote samples with Dvoid ¼ 0.5 lm and HMX with porosity ¼ 3%, comprising voids of diameter 1 lm (dashed lines)
/ ¼ 3% (solid circles) and 1.7% (hollow circles). The constants k1, k2, and and 0.5 lm (solid lines) for an input pressure of 12 GPa, and sS of (a) 80 ns,
k3 are 2.03, 1.22, and 0.48 GW/cm2, respectively. (b) 100 ns, and (c) 90 ns. The mesh size is 10 lm for the computations.
until 25 GPa. It is only beyond 25 GPa that the reaction MES-IG model is used to perform SDT simulations in pressed
front overcomes the rarefaction effects and the von HMX with different void diameters and porosity. The simula-
Neumann pressure is achieved. Therefore, h is higher tions are run to calculate run-to-detonation distances (h) and
for Dvoid ¼ 0.5 lm for the short-duration pulse, com- are validated against pop-plots obtained from experiments.41
pared to Dvoid ¼ 1.0 lm. Additionally, the critical energy for SDTs are also computed
(b) For ss ¼ 100 ns, i.e., the long pulse, h is higher for the to recover Walker-Wasley relations for pressed HMX.
larger diameter void, i.e., the von Neumann pressure is To validate the pop-plots, three combinations of the poros-
attained later in the Dvoid ¼ 1 lm case than in the Dvoid ity (/) and the void-diameter (Dvoid ) are used in the macro-
¼ 0.5 lm case. This is shown in Fig. 14(b). Similar to scale model. The values of (/, Dvoid ) are selected to be (/,
Fig. 14(a), for ss ¼ 100 ns, the reaction rates are higher Dvoid ) ¼ (3%, 0.5 lm), (3%, 1.0 lm) and (1.7%, 0.5 lm). For
initially for Dvoid ¼ 1 lm, whereas, after 25 GPa pres- all three combinations of (/, Dvoid ), the computed values of h
sure is reached, the reaction rates are higher for are in good agreement with the experimental observations.
Dvoid ¼ 0.5 lm. However, because the pulse is suffi- To simulate flyer-plate experiments, the MES-IG model
ciently long, rarefaction is not significant in either is used to capture SDTs using pulses of finite time duration.
case, i.e., for Dvoid ¼ 1 lm or 0.5 lm [as seen in Fig. It is found that the multiscale model naturally recovers the
14(b)]. Rarefaction effects do not compete with the Walker-Wasley relations for the aforementioned combina-
growing reaction front. The SDT is similar to the sus- tions of / and Dvoid . The meso-informed computations suc-
tained loading case described in Sec. III B 3, i.e., for cessfully demarcate the go-no-go boundaries, identifying the
the longer pulse of ss ¼ 100 ns, h is smaller for critical energy for SDT in pressed HMX.
Dvoid ¼ 0.5 lm. The present framework does not rely on a priori
(c) Between ss ¼ 80 and 100 ns, there exists a value of ss , assumptions about the hotspot distributions or void collapse
where h is the same for both Dvoid ¼ 1 lm and mechanisms at the meso-scale. Furthermore, although the
Dvoid ¼ 0.5 lm cases. This is shown in Fig. 14(c), for framework is demonstrated in the context of an IG model at
ss ¼ 90 ns. At this intermediate ss , the effects of rare- the macro-scale, the surrogate model based approach can be
faction and the respective chemical heat release are used in any general macro-scale model for predicting SDT in
such that for both Dvoid ¼ 1 and 0.5 lm, the von- HE materials. Therefore, the framework presented in this
18 September 2024 16:11:18

Neumann pressures are achieved at the same location paper can provide a useful and reliable tool for sensitivity
in the material. In other words, the run-to-detonation characterization of HE materials.
distances are the same for both the diameters.
However, as seen in Fig. 11(c), the route through which
the von Neumann pressure is reached is different for ACKNOWLEDGMENTS
the two diameters: the 0.5 lm void case starts off lag- We gratefully acknowledge the financial support from the
ging behind the 1 lm void but in the later stages the Air Force Office of Scientific Research under Grant No.
0.5 lm case catches up because of a steeper rise in the FA9550-15-1-0332 (Dynamic Materials Program, program
pressure to the von Neumann value. manager: Dr. Martin Schmidt) and under Grant No. SA0000506
In summary, the results in this section show that the present (Computational Mathematics Program, program manager: Dr.
MES-IG model is useful for sensitivity characterization and pre- Jean-Luc Cambier). The authors are also thankful to Dr. K. K.
dictive modeling of SDTs in pressed HMXs. The pop-plots Choi at the University of Iowa and Dr. Nicholas J. Gaul at
obtained for conditions pertaining to experimental meso- RAMDO LLC, Iowa City, for providing the computational code
structure and loading conditions15,41 agree with the experimental for the Modified Bayesian Kriging Method.
observations and the Walker-Wasley relations are recovered
naturally using the present MES-IG model. The model also cap- APPENDIX A: EQUATIONS OF STATE FOR THE
tures finer details of the effect of void sizes on the sensitivity, in GASEOUS REACTION PRODUCTS AND SOLID HMX
agreement with such effects noted by Price.71 Further discus- AT THE MACROSCALE
sions on the surrogate-model in extended parameter spaces and
For the gas phase, the Jones-Wilkins-Lee (JWL) equa-
detailed examination of the physics of void collapse connected
tion of state is used where the specific internal energy of the
to the behavior of F_ ignition and F_ growth are pursued in Nassar
et al.,62 Roy et al.61 and Rai and Udaykumar.63 gas is given by15,72,73

IV. CONCLUSIONS eg qg ; Tg ¼ ekg ðqg Þ þ Cvg Tg : (A1)
This work validates a multiscale modeling framework for Cvg is the specific heat at constant volume, Tg is the tempera-
simulating shock-to-detonation transitions (SDT) in pressed ture of the product gases, and
HMX. To bridge scales, a surrogate model is constructed from
high fidelity meso-scale simulations of reactive void collapse. Ag R1g qq0g Bg R2g qq0g
ekg ðqg Þ ¼ e þ e þ cek
At the macro-scale, the meso-scale informed Ignition and q0 R1g q0 R2g
Growth (MES-IG) model is used. The reaction rates in the Cg0
k qg
MES-IG model are learned by training a Bayesian Kriging þ ; (A2)
code using a small ensemble of numerical experiments. The q0 Cg0 q0
where qg is the density of the gas and q0 is a reference den- TABLE II. HMX parameters for the JWL Equation of State for the product
phase.15
sity. The gas-phase pressure, pg , is a function of the internal
energy and density and is given by Ag (Pa) Bg (Pa) R1g R2g Cvg ðJ kg1 K1 Þ Cg0
pg ðqg ; eg Þ ¼ qg Cg0 ðeg ekg ðqg ÞÞ þ pkg ðqg Þ; (A3) 12:18 1011 12:3 109 4:766 1:102 815 0.408
where Cg0 is the Gruneisen coefficient and

dekg TABLE III. Variables at the CJ state and reference quantities for HMX.15
pkg ðqg Þ ¼ : (A4)
d q1
g PCJ (Pa) DCJ (m/s) TCJ ðKÞ q0 ðkg=m3 Þ e0 ðJ kg1 Þ T0 ðKÞ
9
The constants in Eq. (A1) are obtained from 38:8 10 9300 4113:7 1905 0 298
Cg0 þ1
q0
k ¼ PCJ pkg ðqCJ Þ qCJ Cg0 Cvg TCJ (A5) where Cs0 is the Gruneisen coefficient and
qCJ
deks
Ag R1g qq0 Bg R2g qq0 pks ðqs Þ ¼ : (A12)
cek ¼ e CJ e CJ 1
q0 R1g q0 R2g d
qs
1
PCJ pkg ðqCJ Þ þ eCJ (A6)
qCJ Cg0 The constants in the equations of state for the reactants and
2 products are given in Tables I–III.15
1 1 1 The equations of state for the solid and the gas
¼ PCJ (A7)
qCJ q0 q0 DCJ phases can be combined with the mixture laws [Eqs.
(10)–(12)] to compute the pressure in the macro-scale
and control volume. The mixture law for the specific internal

1 1 1 energy given by Eq. (11), and the pressure equilibrium
eCJ ¼ PCJ þ e0 : (A8)
18 September 2024 16:11:18

2 q0 qCJ condition, in Eq. (10) can be combined with Eqs. (A3)
and (A11) to give
In Eqs. (A5)–(A8), PCJ , TCJ , qCJ , eCJ; and DCJ are the pres-
p pks p pkg
sure, temperature, density, specific internal energy, and the ð1 k Þ þ eks ðqs Þ þ k þ ekg ðqg Þ ¼ e:
qs Cs0 qg Cg0
detonation speed at the CJ state respectively, while e0 is a
reference specific internal energy. (A13)
For the solid phase, the Cochran-Chan equation of state
Similarly, eliminating eg from Eqs. (A1) and (A3) and es
is used to model the thermomechanical behavior of the solid
from Eqs. (A9) and (A11) in favor of Ts and Tg , the tempera-
under shocks. In this equation of state (EOS), the specific
ture equilibrium condition [Eq. (10)] gives
internal energy of the solid reactant is given by15,72,74
p pks ðqs Þ p pkg ðqg Þ
es ðqs ; Ts Þ ¼ eks ðqs Þ þ Cvs Ts ; (A9) ¼ : (A14)
qs Cs0 Cvs qg Cg0 Cvg
where Cvs is the specific heat at constant volume, Ts is the
temperature of the solid, and Equation (12), (A13), and (A14) provide three equations in
" # terms of the mixture pressure p, and the densities of the gas
As q0 1R1s and solid phases, qg and qs and can be solved simultaneously
eks ðqs Þ ¼ 1
q0 ð1 R1s Þ qs for a specified value of k.
" #
Bs q0 1R2s APPENDIX B: EQUATIONS OF STATE AND REACTION
þ 1 Cvs T0 þ e0 ; KINETIC LAWS AT THE MESO-SCALE
q0 ð1 R2s Þ qs
The pressure at the meso-scale is obtained from a Birch-
(A10)
Murnaghan equation of state,47,49 which can be written in
where qs is the density of the solid phase and T0 is the refer- the general Mie-Gruneisen form as
ence temperature. The solid pressure, ps , is given by pðq; eÞ ¼ pk ðqÞ þ qCs ½e ek ðqÞ; (B1)
ps ðqs ; es Þ ¼ qs Cs0 ðes eks ðqs ÞÞ þ pks ðqs Þ; (A11) where

" 5=3 #
TABLE I. HMX parameters for the Cochran-Chan Equation of State for the 3 q 7=3 q
solid phase.15 pk ðqÞ ¼ KT0
2 q0 q0
Cvs ðJ kg1 K1 Þ " " # #
q 2=3
As (Pa) Bs (Pa) R1s R2s Cs0
3 0
1:287 1010 1:342 1010 4:1 3:1 1087 0.93 1 þ KT0 4 1 ; (B2)
4 q0
TABLE IV. HMX chemical reaction parameters for Tarver 3-equation50 and secondary collapse on initiation thresholds,” Phys. Rev. Fluids 2(4),
model. 043202 (2017).
5
N. K. Rai, M. J. Schmidt, and H. S. Udaykumar, “Collapse of elongated
ln Z1 ðs1 Þ ln Z2 ðs1 Þ ln Z3 ðs1 Þ E1 ðkcal=molÞ E2 ðkcal=molÞ E3 ðkcal=molÞ voids in porous energetic materials: Effects of void orientation and aspect
ratio on initiation,” Phys. Rev. Fluids 2(4), 043201 (2017).
6
48:7 37:3 28:1 52:7 44.1 34.1 E. J. Welle, C. D. Molek, R. R. Wixom, and P. Samuels, Microstructural
Effects on the Ignition Behavior of HMX (IOP Publishing, 2014).
7
C. Molek and E. Welle, Image Obtained by Ryan Wixom (LANL), private
ð 1=q communication (2016).
1 8
E. L. Lee and C. M. Tarver, “Phenomenological model of shock initiation
e k ðq Þ ¼ e 0 pk ðqÞd : (B3)
1=q0 q in heterogeneous explosives,” Phys. Fluids 23(12), 2362–2372 (1980).
9
J. N. Johnson, P. K. Tang, and C. A. Forest, “Shock-wave initiation of het-
Chemical decomposition of HMX is modeled in the Tarver erogeneous reactive solids,” J. Appl. Phys. 57(9), 4323–4334 (1985).
10
R. Menikoff and M. S. Shaw, “The SURF model and the curvature effect
3-equation model50 via four different species. The three reac- for PBX 9502,” Combust. Theory Modell. 16(6), 1140–1169 (2012).
tions are 11
M. Baer and J. Nunziato, “A two-phase mixture theory for the deflagra-
tion-to-detonation transition (DDT) in reactive granular materials,” Int. J.
Reaction 1 : HMX ðC4 H8 N8 O8 Þ ! fragments ðCH2 NNO2 Þ Multiphase Flow 12(6), 861–889 (1986).
12
M. W. Crochet and K. A. Gonthier, “A multi-phase theory for the detona-
(B4) tion of granular explosives containing an arbitrary number of solid
components,” Int. J. Multiphase Flow 77, 76–89 (2015).
Reaction 2 : fragments ðCH2 NNO2 Þ 13
O. Sen, N. J. Gaul, K. K. Choi, G. Jacobs, and H. S. Udaykumar,
! intermediate gases ðCH2 O; N2 O; HCN; HNO2 Þ (B5) “Evaluation of kriging based surrogate models constructed from mesoscale
computations of shock interaction with particles,” J. Comput. Phys. 336,
Reaction 3 : 14
235–260 (2017).
O. Sen, S. Davis, G. Jacobs, and H. S. Udaykumar, “Evaluation of conver-
2 intermediate gases ðCH2 O; N2 O; HCN; HNO2 Þ gence behavior of metamodeling techniques for bridging scales in multi-
! final gases ðN2 ; H2 O; CO2 ; COÞ: (B6) scale multimaterial simulation,” J. Comput. Phys. 294, 585–604 (2015).
15
J. Massoni, R. Saurel, G. Baudin, and G. Demol, “A mechanistic model
for shock initiation of solid explosives,” Phys. Fluids 11(3), 710–736
The solid HMX (species 1, with mass fraction Y1 ) at high (1999).
temperature decomposes to fragments (species 2, with mass 16
J. Kang, P. Butler, and M. Baer, “A thermomechanical analysis of hot spot
18 September 2024 16:11:18

fraction Y2 ). The fragments further decompose to intermedi- formation in condensed-phase, energetic materials,” Combust. Flame
89(2), 117–139 (1992).
ate gases (species 3 with mass fraction Y3 ) and finally the 17
T. L. Jackson and Ju Zhang, “Density-based kinetics for mesoscale simula-
intermediate gases undergo strong exothermic reactions to tions of detonation initiation in energetic materials,” Combustion Theory
form gaseous products (species 4, Y4 ), which leads to an Modelling 21, 749–769 (2017).
18
increase in temperature. J. Zhang, T. L. Jackson, J. D. Buckmaster, and J. B. Freund, “Numerical
modeling of shock-to-detonation transition in energetic materials,”
The rate equations for all the species are given as Combust. Flame 159(4), 1769–1778 (2012).
19
Y. Hamate and Y. Horie, “Ignition and detonation of solid explosives: A
E1 micromechanical burn model,” Shock Waves 16(2), 125–147 (2006).
_
Y 1 ¼ Y1 Z1 exp (B7) 20
M. Akiki and S. Menon, “A model for hot spot formation in shocked ener-
RT getic materials,” Combust. Flame 162(5), 1759–1771 (2015).
21
L. Tran and H. S. Udaykumar, “Simulation of void collapse in an energetic
_ E1 E2 material, Part I: Inert case,” J. Propul. Power 22(5), 947–958 (2006).
Y 2 ¼ Y1 Z1 exp Y2 Z2 exp (B8) 22
L. Tran and H. S. Udaykumar, “Simulation of void collapse in an energetic
RT RT
material, Part 2: Reactive case,” J. Propul. Power 22(5), 959–974 (2006).
23
A. Kapahi and H. S. Udaykumar, “Dynamics of void collapse in shocked
_ E2 2 E3 energetic materials: Physics of void–void interactions,” Shock Waves
Y 3 ¼ Y2 Z2 exp Y3 Z3 exp (B9)
RT RT 23(6), 537–558 (2013).
24
P. T. Rao, K. A. Gonthier, and S. Chakravarthy, “Compaction shock dissi-
pation in low density granular explosive,” J. Appl. Phys. 119(22), 224904
E3
Y_ 4 ¼ Y32 Z3 exp ; (B10) (2016).
RT 25
J. R. Gambino, D. W. Schwendeman, and A. K. Kapila, “Numerical study
of multiscale compaction-initiated detonation,” Shock Waves 1–27
where Yi is the mass fraction of the ith species, Zi is the fre- 26
(2018).
quency factor for each reaction, Ej is the activation energy A. Barua and M. Zhou, “Computational analysis of temperature rises in micro-
structures of HMX-Estane PBXs,” Comput. Mech. 52(1), 151–159 (2013).
for each reaction, R is the universal gas constant, and T is 27
S. Kim, C. Miller, Y. Horie, C. Molek, E. Welle, and M. Zhou,
the temperature. The values for each of these constants are “Computational prediction of probabilistic ignition threshold of pressed
provided in Table IV.25 granular octahydro-1, 3, 5, 7-tetranitro-1, 2, 3, 5-tetrazocine (HMX) under
shock loading,” J. Appl. Phys. 120(11), 115902 (2016).
28
A. Nichols and C. M. Tarver, A Statistical Hot Spot Reactive Flow Model
1
J. E. Field, “Hot spot ignition mechanisms for explosives,” Acc. Chem. for Shock Initiation and Detonation of Solid High Explosives (Lawrence
Res. 25(11), 489–496 (1992). Livermore National Laboratory, CA, USA, 2002).
2 29
T. L. Jackson, J. D. Buckmaster, J. Zhang, and M. J. Anderson, “Pore col- A. K. Kapila, D. W. Schwendeman, J. R. Gambino, and W. D. Henshaw,
lapse in an energetic material from the micro-scale to the macro-scale,” “A numerical study of the dynamics of detonation initiated by cavity
Combust. Theory Modell. 19(3), 347–381 (2015). collapse,” Shock Waves 25(6), 545–572 (2015).
3 30
F. Garcia, K. S. Vandersall, and C. M. Tarver, “Shock initiation experi- N. K. Rai and H. Udaykumar, “Three-dimensional simulations of void col-
ments with ignition and growth modeling on low density HMX,” J. Phys.: lapse in energetic materials,” Phys. Rev. Fluids 3(3), 033201 (2018).
31
Conf. Ser. 500(5), 052048 (2014). N. K. Rai and H. S. Udaykumar, “Mesoscale simulation of reactive
4
N. K. Rai, M. J. Schmidt, and H. S. Udaykumar, “High-resolution simula- pressed energetic materials under shock loading,” J. Appl. Phys. 118(24),
tions of cylindrical void collapse in energetic materials: Effect of primary 245905 (2015).
32 53
A. Kapahi and H. Udaykumar, “Three-dimensional simulations of dynam- N. K. Rai, “Numerical framework for mesoscale simulation of heteroge-
ics of void collapse in energetic materials,” Shock Waves 25(2), 177–187 neous energetic materials,” Ph.D. dissertation (University of Iowa, 2015).
54
(2015). A. Kapila, D. Schwendeman, J. Bdzil, and W. Henshaw, “A study of deto-
33
R. Zhang, Y. Lan, G. B. Huang, and Z. B. Xu, “Universal approximation nation diffraction in the ignition-and-growth model,” Combust. Theory
of extreme learning machine with adaptive growth of hidden nodes,” IEEE Modell. 11(5), 781–822 (2007).
55
Trans. Neural Networks Learn. Syst. 23(2), 365–371 (2012). R. Menikoff, JWL Equation State (Los Alamos National Laboratory
34
N. J. Gaul, M. K. Cowles, H. Cho, K. Choi, and D. Lamb, “Modified (LANL), 2015).
56
Bayesian kriging for noisy response problems for reliability analysis,” in P. Souers, B. Wu, and L. Haselman, Jr., Detonation Equation of State at
ASME 2015 International Design Engineering Technical Conferences and LLNL, 1995, Revision 3 (Lawrence Livermore National Laboratory, CA,
Computers and Information in Engineering Conference (American Society USA, 1996).
57
of Mechanical Engineers, 2015). G. DeOliveira, A. Kapila, D. Schwendeman, J. Bdzil, W. Henshaw, and C.
35
L. Zhao, K. Choi, and I. Lee, “Metamodeling method using dynamic krig- Tarver, “Detonation diffraction, dead zones, and the ignition-and-growth
ing for design optimization,” AIAA J. 49(9), 2034–2046 (2011). model,” in The Thirteenth Symposium (International) on Detonation
36
S. Taverniers, T. S. Haut, K. Barros, F. J. Alexander, and T. Lookman, (2006).
58
“Physics-based statistical learning approach to mesoscopic model H. K. Springer, C. M. Tarver, and S. Bastea, “Effects of high shock pres-
selection,” Phys. Rev. E 92(5), 053301 (2015). sures and pore morphology on hot spot mechanisms in HMX,” AIP Conf.
37
D. Xiu and G. E. Karniadakis, “Modeling uncertainty in flow simulations Proc. 1793(1), 080002 (2017).
59
via generalized polynomial chaos,” J. Comput. Phys. 187(1), 137–167 A. Swantek and J. Austin, “Collapse of void arrays under stress wave
(2003). loading,” J. Fluid Mech. 649, 399–427 (2010).
38 60
R. Zhang and S. Mahadevan, “Model uncertainty and Bayesian updating L. Tran and H. S. Udaykumar, “A particle-level set-based sharp interface
in reliability-based inspection,” Struct. Saf. 22(2), 145–160 (2000). cartesian grid method for impact, penetration, and void collapse,”
39
B. D. Youn, K. K. Choi, L. Du, and D. Gorsich, “Integration of J. Comput. Phys. 193(2), 469–510 (2004).
61
possibility-based optimization and robust design for epistemic S. Roy, N. K. Rai, O. Sen, and H. S. Udaykumar, “Modeling meso-scale
uncertainty,” J. Mech. Des. 129(8), 876–882 (2007). energy localization in shocked HMX, Part II: Machine learned surrogate
40
I. Park and R. V. Grandhi, “A Bayesian statistical method for quantifying models for void shape and void-void interaction effects,” Shock Waves
model form uncertainty and two model combination methods,” Reliab. (submitted).
62
Eng. Syst. Saf. 129, 46–56 (2014). A. Nassar, N. K. Rai, O. Sen, and H. S. Udaykumar, “Modeling meso-
41
J. Vanpoperynghe, J. Sorel, J. Aveille, and J. Adenis, “Shock initiation of scale energy localization in shocked HMX, Part I: Machine-learned surro-
TATB and HMX explosive compositions,” in 8th Symposium gate model for effect of loading and void size,” Shock Waves (submitted).
63
(International) on Detonation (1985). N. K. Rai and H. S. Udaykumar, “Void collapse generated meso-scale
42
O. Sen, N. J. Gaul, K. Choi, G. Jacobs, and H. Udaykumar, “Evaluation of energy localization in shocked energetic materials: Non-dimensional
multifidelity surrogate modeling techniques to construct closure laws for parameters, regimes and criticality of hotspot,” Phys. Fluids
18 September 2024 16:11:18

drag in shock-particle interactions,” J. Comput. Phys. 371, 434–451 (unpublished).
64
(2018). C.-W. Shu and S. Osher, “Efficient implementation of essentially non-
43
A. Kapahi, J. Mousel, S. Sambasivan, and H. Udaykumar, “Parallel, sharp oscillatory shock-capturing schemes, II,” J. Comput. Phys. 83(1), 32–78
interface Eulerian approach to high-speed multi-material flows,” Comput. (1989).
65
Fluids 83, 144–156 (2013). C.-W. Shu and S. Osher, “Efficient implementation of essentially non-
44
N. K. Rai, A. Kapahi, and H. S. Udaykumar, “Treatment of contact separa- oscillatory shock-capturing schemes,” J. Comput. Phys. 77(2), 439–471
tion in Eulerian high-speed multimaterial dynamic simulations,” Int. J. (1988).
66
Numer. Methods Eng. 100(11), 793–813 (2014). H. Udaykumar, L. Tran, D. Belk, and K. Vanden, “An Eulerian method
45
S. Sambasivan, A. Kapahi, and H. S. Udaykumar, “Simulation of high for computation of multimaterial impact with ENO shock-capturing and
speed impact, penetration and fragmentation problems on locally refined sharp interfaces,” J. Comput. Phys. 186(1), 136–177 (2003).
67
Cartesian grids,” J. Comput. Phys. 235, 334–370 (2013). S. C. Chapra and R. P. Canale, Numerical Methods for Engineers
46
J.-P. Ponthot, “Unified stress update algorithms for the numerical simula- (McGraw-Hill, New York, 1998), Vol. 2.
68
tion of large deformation elasto-plastic and elasto-viscoplastic processes,” F. Walker and R. Wasley, “A general model for the shock initiation of
Int. J. Plast. 18(1), 91–126 (2002). explosives,” Propellants, Explos., Pyrotech. 1(4), 73–80 (1976).
47 69
R. Menikoff and T. D. Sewell, “Constituent properties of HMX needed for H. R. James, “An extension to the critical energy criterion used to predict
mesoscale simulations,” Combust. Theory Modell. 6(1), 103–125 (2002). shock initiation thresholds,” Propellants, Explos., Pyrotech. 21(1), 8–13
48
A. Kapahi, S. Sambasivan, and H. Udaykumar, “A three-dimensional (1996).
70
sharp interface Cartesian grid method for solving high speed multi- P. Howe, R. Frey, B. Taylor, and V. Boyle, in Sixth Symposium
material impact, penetration and fragmentation problems,” J. Comput. (International) on Detonation, California (1976).
71
Phys. 241, 308–332 (2013). D. Price, Effect of Particle Size on the Shock Sensitivity of Pure Porous
49
T. D. Sewell and R. Menikoff, Complete Equation of State for/3-HMX and HE (High Explosive) (Naval Surface Weapons Center, Silver Spring, MD,
Implications for Initiation (American Institute of Physics, 2003). 1986).
50 72
C. M. Tarver, S. K. Chidester, and A. L. Nichols, “Critical conditions for G. Baudin, “A Reaction Model for Aluminized PBX Applied to
impact-and shock-induced hot spots in solid explosives,” J. Phys. Chem. Underwater Explosive Calculations,” in Tenth Symposium (International)
100(14), 5794–5799 (1996). on Detonation (1993).
51 73
G. Strang, “On the construction and comparison of difference schemes,” E. Lee, H. Hornig, and J. Kury, Adiabatic Expansion of High Explosive
SIAM J. Numer. Anal. 5(3), 506–517 (1968). Detonation Products (University of California Radiation Laboratory,
52
E. Fehlberg, Classical Fifth-, Sixth-, Seventh-, and Eighth-Order Runge- Livermore, CA, USA, 1968).
74
Kutta Formulas with Stepsize Control (National Aeronautics and Space S. G. Cochran and J. Chan, “Shock initiation and detonation models in one
Administration, 1968). and two dimensions,” Technical Report (1979).

1 Online

Uploaded by

Copyright:

Available Formats

1 Online

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

1 Online

Uploaded by

Copyright:

Available Formats

RESEARCH ARTICLE | AUGUST 27 2018

Multi-scale shock-to-detonation simulation of pressed

J. Appl. Phys. 124, 085110 (2018)

18 September 2024 16:11:18

Multi-scale shock-to-detonation simulation of pressed energetic material:

18 September 2024 16:11:18

0021-8979/2018/124(8)/085110/21/$30.00 124, 085110-1 Published by AIP Publishing.

Deflagration to Detonation (DDT) scenarios, by analyzing the

18 September 2024 16:11:18

18 September 2024 16:11:18

where q and ui , respectively, denote the density and the

rij ¼ Sij pdij ; (4)

where Sij is the deviatoric stress tensor and p is the pressure.

18 September 2024 16:11:18

18 September 2024 16:11:18

product and reactant equations of state (described in

18 September 2024 16:11:18

18 September 2024 16:11:18

18 September 2024 16:11:18

18 September 2024 16:11:18

18 September 2024 16:11:18

18 September 2024 16:11:18

18 September 2024 16:11:18

18 September 2024 16:11:18

18 September 2024 16:11:18

18 September 2024 16:11:18

2. Critical energy thresholds for pressed HMX

18 September 2024 16:11:18

18 September 2024 16:11:18

18 September 2024 16:11:18

where Cg0 is the Gruneisen coefficient and

18 September 2024 16:11:18

ps ðqs ; es Þ ¼ qs Cs0 ðes eks ðqs ÞÞ þ pks ðqs Þ; (A11) where

18 September 2024 16:11:18

18 September 2024 16:11:18

You might also like