Direct numerical simulation of shear localization and decomposition reactions in shock-loaded HMX crystal

Under relatively weak shock/impact loading conditions, the ignition of chemical reactions in solid (heterogeneous) explosives is controlled by localized regions of elevated temperature, i.e., hot spots.^1–3 These hot spots are formed by various dissipative processes (e.g., micro-jetting in collapsing pores, inelastic deformation, and frictional sliding), which convert mechanical work into thermal energy.⁴ A critical hot spot produces a small burning region that propagates in a self-sustained manner. It is the interaction of burn/pressure fronts from many hot spots that is responsible for the build-up to detonation.

Despite our knowledge of the importance of hot spots, a precise understanding of the ignition and growth process is lacking for a wide variety of shock-loaded solid high explosive (HE) formulations. As mentioned above, there is not a single universal mechanism of hot spot formation. The operative mechanisms depend, rather, on the properties of the material, the underlying micro/mesostructure, and the loading path/intensity. To improve our understanding of shock wave ignition in specific HE formulations, the mechanisms that give rise to hot spots should be identified and the hot spot sizes and temperatures should be quantified. The effects of chemical reactions in the hot spots should also be considered, as reaction products will interact with the surrounding material. Such information would undoubtedly find use in the development and calibration of reactive flow models that describe the ignition and growth of reactions from a macroscopic point of view.⁵ Given the time and length scales that are involved, the shock-induced hot spot dynamics are not amenable to experimental measurement. We have therefore adopted a computational approach.

This work seeks to develop a predictive model for the chemical reaction of porous octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX) crystal under shock wave loading. HMX has been selected for study because it serves as an energetic component in various HE formulations and there exists a reasonable base of experimental data from which a physically based model can be constructed. The continuum model developed in this work is unique in that it accounts for the crystal mechanics and temperature-driven decomposition reactions. This allows for the direct simulation of crystal-sensitive localization behavior and the chemical reaction of hot spots that are formed under shock/impact loading.

Pore collapse (both intra- and intergranular) is thought to play an important role in the generation of critical hot spots. We therefore considered, as a model problem, the shock-induced closure of a single intragranular pore in β-phase HMX crystal under a shock stress of ∼10 GPa. A systematic study was then undertaken to quantify the effects of loading intensity, pore size, conductive heat transfer, crystal strength, and liquid viscosity.

The most stable phase of HMX at room temperature and 1 atm is the β phase. The crystal structure of the β phase is monoclinic and can be represented in either of two equivalent space groups, P2₁/n or P2₁/c, with two molecules per unit cell.^6,7 A review of the known properties of HMX is available in the literature;⁸ however, the full characterization of thermal/mechanical properties under shock loading conditions remains a difficult task given the pressures and time scales involved and the occurrence of rapid exothermic reactions at high temperatures.

Chemical reaction pathways and kinetics for the decomposition of HMX-based formulations have been studied in experimental and modeling efforts (cf. Refs. 9–14 and work cited therein). These efforts have focused on thermal explosion experiments, where specimens are slowly heated-up (under various levels of confinement) until thermal runaway is achieved. This has led to the parameterization of a single-step reaction scheme¹² and various multi-step reaction networks.^11,12,15 Kinetic constants have also been measured for the decomposition of liquid-HMX over a limited range of temperature.¹⁶ At 1 atm, the β phase starts to transform to the δ phase (hexagonal crystal structure) at about 438 K.⁸ The decomposition reaction then proceeds from the δ phase provided the pressure is kept low enough.¹⁷ Accordingly, the β → δ transformation is usually taken as the first step in multi-step reaction paths. However, experiments performed at higher pressures (>1 GPa) showed that this phase transformation is inhibited and that decomposition reactions proceed directly from the β phase.¹⁷ Although the reaction path that is followed under shock wave loading has not been definitively identified, it seems reasonable to assume that some amount of β-HMX is melted by the shock deformation and that accelerated reactions proceed directly from the melt.⁸ 

Shock wave ignition (or initiation) is regarded as the first step of the shock-to-detonation transition. In this context, ignition refers to the production of critical hot spots, which, if left alone, burn in a self-sustaining manner. This work considers hot spots that are generated by the collapse of intragranular pores. These pores are formed during crystal growth and are distinct from intergranular pores that are also present in (polymer-bonded) crystal aggregates. The initial stages of reaction are the primary focus of this work. As such, the literature should be consulted for discussions of critical conditions for idealized hot spots,¹¹ the effectiveness of viscoplastic deformation as a means for igniting material about a collapsing pore,¹⁸ the sensitivity of dynamic pore collapse to loading rate and various material modeling assumptions,¹⁹ and the relative importance of intra- and intergranular pores with respect to ignition.²⁰ 

Previous efforts have considered the direct numerical simulation of pore collapse in various HE materials. For example, continuum-based approaches have been used to study intragranular pore collapse in HMX^21–23 and TATB²⁴ and the dynamic loading of HMX grain aggregates.^25–29 Discrete methods have also been employed. For example, reactive molecular dynamics (MD) simulations have been used to study pore collapse in PETN³⁰ and fuel oil inclusions in ammonium nitrate,³¹ whereas coarse-grained descriptions based on dissipative particle dynamics show promise for such calculations in RDX.^32,33 Previous continuum-based approaches have, for the most part, viewed the mechanical responses of crystalline HE materials as isotropic and independent of the rate of deformation. While this has been useful for studying certain effects and mapping out various model sensitivities, we have sought to include an improved description of the crystal mechanics. This is motivated by the fact that many HE crystals are of low symmetry and display strongly anisotropic thermal/mechanical responses. Furthermore, many materials exhibit a significant increase in flow stress at the high strain rates encountered under shock wave loading. This is particularly true when plastic deformation is mediated by slip (dislocation motion). Such considerations were, indeed, shown to be important when computing wave profiles in shock-loaded PETN crystal.³⁴ To progress towards improved descriptions of HE grain-scale responses, it appears necessary to account for the anisotropic thermo-elasto-viscoplastic behavior of the crystalline phases. This is important to our study of HMX crystal because inelastic deformation is an important source of localized heating.

In recent years, a number of frameworks have been developed to treat high-pressure crystal thermo-elasto-viscoplasticity with varying levels of thermodynamic consistency.^23,35–39 The crystal mechanics framework that was developed for β-HMX²³ is basis of this modeling effort. The previous model has been calibrated to elastic precursor wave decay measurements⁴⁰ and was used to simulate pore collapse in inert β-HMX crystal.²³ In this work, the crystal model is extended to account for the effects of exothermic decomposition reactions.

The material model provides a description of two phase transformations: crystal melting (β → liquid) and decomposition reactions that yield gaseous products. The decomposition chemistry is handled using the thermochemical code, Cheetah.⁴¹ Given the uncertainty in reaction path under shock wave loading, we have adopted a single-step (global) reaction instead of the more complex multi-step networks. The selected reaction addresses the decomposition of both β- and liquid-phase HMX, i.e.,

In this scheme, the chemical species to be tracked are HMX and the product gas mixture. The kinetics of this reaction are prescribed. In addition to the reaction given above, numerous fast reactions are allowed to occur among the product gas components. These reactions, which occur instantaneously, serve to maintain chemical equilibrium among the gas components and the HMX species. As a result, the product gas mixture is adjusted to include other gases (for example, C, H₂, CO₂, HCN, NO₂, N₂O, NH₃, and so on). The current product gas composition is therefore not given exactly by (1), but rather that which is obtained by minimizing the Gibbs free energy of system, subject to the kinetic constraint on HMX concentration.

Assuming a decomposition reaction that is first-order with respect to HMX, the species evolution rates are written as

where bracketed quantities refer to concentrations (in mol per volume) and k is the rate coefficient. An Arrhenius rate coefficient is utilized, i.e.,

where k₀ is the rate constant, E_a is the activation energy, R is the gas constant, and T is the temperature. Despite its simplicity, the single-step reaction provides a reasonable approximation of experimentally measured ignition times when coupled with the kinetic parameters of Henson.¹² The relative success of the single-step description points to the highly exothermic gas production step as limiting the overall rate of reaction. As noted above, this formalism assumes the decomposition kinetics for the β and liquid phases are the same. While this seems to be a reasonable starting point, the model might be improved by distinguishing the kinetics of reactions in the solid and liquid phases.¹⁶ 

The thermo-elasto-viscoplastic behavior of the β phase is described using a crystal model that was developed in previous work.²³ The crystal model incorporates an equation of state (EOS) and a strength law that is sensitive to the rate of deformation, the microstructure, and the thermodynamic state. Although β-HMX is generally regarded as brittle at ambient pressure,⁴² our working assumption is that the high-pressure state of the shock wave allows for larger extents of plastic deformation prior to cleaving (brittle fracture). A brief account of the model is provided in this section. The literature should be referenced for full details of the model and its parameterization.²³ In the following, all crystal quantities are referenced to the P2₁/c space group unless noted otherwise.

The crystal kinematics are written using a multiplicative decomposition of the deformation gradient, i.e.,

where $Fp$ describes the plastic shearing of the lattice, R is the lattice rotation, and V is the thermoelastic lattice stretch. Plastic deformation is considered to occur by crystallographic slip. As such, the velocity gradient in an intermediate configuration is written in terms of the crystallographic shearing rates, $γ̇α$⁠, and the slip system geometries, i.e.,

where $ŝα$ and $m̂α$ are unit vectors that describe the slip direction and slip plane normal of the α^th slip system, respectively. The slip kinetics transition from being dominated by thermal activation at lower shear stresses to phonon drag at higher shear stresses.

The slip systems utilized in this work are given in Table II. This listing includes the two slip systems that were identified in experiments, i.e., (001)[100] and $(1¯02)[201]$ (equivalently, (001)[100] and $(101)[101¯]$ in P2₁/n)^40,43 and additional slip systems that were identified in MD calculations.²³ These additional slip systems serve to close the flow surface. The slip system resistances are written as

where g₀ and s are constants, $h=ρ⊥/ρref⊥$ is the normalized dislocation density, and r^α are the slip system strength ratios, given in Table II. The dislocation density is an evolving internal state variable (ISV) in this description.

The thermoelastic relation is written as $τ¯=τ¯(ln(V¯),e)$⁠, where $τ¯$ and $V¯$ are the Kirchhoff stress and lattice stretch, respectively, which rotate with the crystal frame (e.g., $V¯=RTVR$⁠), and e is the specific internal energy. The functional form incorporates the second-order elastic constants, which have been computed by atomistics,⁴⁴ and a Murnaghan EOS, which expresses the pressure as follows:⁴⁵ 

In the above equation, K₀ is the reference bulk modulus, n = (∂K/∂p)₀, and α₀ is the reference volumetric coefficient of thermal expansion (CTE). The EOS parameters, which are given in Table I, were selected to reproduce Hugoniot data for solvent-pressed HMX grain aggregates that are close to fully density.⁴⁶ In Fig. 1, the model Hugoniot curve for the β-phase is plotted in pressure-volume space. This illustrates the agreement between the model and the aforementioned data as well as higher-pressure data for single crystals of unknown orientation.⁴⁷ 

To allow for crystal melting, a compression-dependent melt energy curve, e_m(ρ/ρ₀), was fit to a Lindemann law for HMX. The material is considered as fully melted when the internal energy exceeds the melt energy by an amount equal to the latent heat of melting (⁠$Δhm=h0fliq−h0fβ$⁠) and partially melted at intermediate energies. This represents a simplified treatment of the solid-liquid phase transformation, as melting kinetics are not considered. The main purpose of including the melting behavior is to account for the loss of strength associated with the production of liquid-HMX.

The distortional responses of the liquid phase are isotropic and described using a viscous Newtonian fluid law, i.e., $σ′=2ηD′$⁠, where $σ$ is the Cauchy stress, D is the symmetric part of the velocity gradient, η is the liquid viscosity, and the prime symbol refers to the deviatoric part of a tensor. The liquid viscosity is taken as constant, with a value of 5.5 cP. This corresponds to the value computed from atomistic simulations performed at 800 K and 1 atm.⁴⁸ The reference density and EOS of the liquid phase are taken to be identical to that of the β phase (cf. Table I). This is assumed purely for expedience and it is recognized that refined treatments of the liquid phase would improve the model.

The volumetric response of each product gas component (CO, H₂O, N₂, etc.) is described using a Buckingham exponential-6 potential.⁴¹ The distortional response of the product gas mixture is described using the viscous fluid law given above, where the gas mixture viscosity is assumed equal to that of the liquid phase for simplicity.

The inelastic work done on all phases (β, liquid, gas) is fully dissipated and converted to heat. The heat capacity of the β phase is described using an Einstein relation, i.e.,

where the constants have been selected as A₁ = 84 and θ₁ = 1000 K. This parameterization respects the experimental data at lower temperatures⁴⁹ and the classical limit at higher temperatures. The heat capacity of the liquid phase is assumed to be identical to that of the solid phase for simplicity. The calculations currently neglect thermal expansion effects in the solid and liquid phases and the heat of melting. The volumetric expansion of an intensely heated region is therefore not computed properly. By itself, this leads to temperatures that are too high, as the internal energy of an expanding volume under pressure should be reduced. However, the temperature correction is proportional to the product of the CTE (or Gruneisen parameter) and the volumetric strain increment,⁵⁰ which should be small compared to the temperature rise generated by plastic dissipation. To improve the fidelity of predicted thermal/mechanical responses, the measured anisotropic thermal expansion coefficients of the β phase should be incorporated,⁵¹ as well as more appropriate treatments of the liquid phase heat capacity and the heat of melting.

In non-adiabatic simulations, heat conduction is modeled using Fourier's law. The thermal conductivity of the β phase is assumed to be constant, isotropic, and independent of pressure and temperature. Analytical predictions suggest the thermal conductivity tensor of the β phase is close to isotropic.⁵² However, refined treatments might incorporate the effects of temperature and pressure on solid phase thermal conductivity.^11,52 The ambient state thermal conductivity of the β phase is given in Table I. The thermal conductivities of the liquid and gas phases are assumed to be equal to that of the solid phase for simplicity.

The material model considers a material point as a multi-phase entity (i.e., a container of solid, liquid, and gas), which resides at uniform stress and temperature. A scheme must therefore be defined to homogenize the stresses and temperatures of the constituent phases. The homogenization is carried out in a hierarchical fashion, considering first the species components (e.g., β- and liquid-HMX) and then the species (HMX and product gas) to compute the effective state.

The deviatoric stress of the HMX species is computed as a weighted average of stresses in the solid and liquid phases. The weighting factors are based on the fraction of melted material, which is taken as f_m = (e − e_m)/Δh_m (limited to the range of 0–1). The deviatoric stress of a material point is then computed as a weighted average of stresses in the HMX and product gas species. In this computation, the weighting factors are the mass fractions of reactants and products. The species are maintained in pressure/temperature equilibrium by the Cheetah code. The equilibrated pressure and temperature are computed by minimizing the Gibbs free energy while holding (ρ, e) and the species concentrations, [X_i], fixed.

The problem of interest is intragranular pore collapse under shock wave loading. To simulate wave propagation in porous crystal samples, a numerical model was constructed using the multi-physics arbitrary Lagrangian-Eulerian finite element (FE) code, ALE3D.⁵³ In this code, the material and mesh are permitted to undergo independent motions with algorithms accounting for advection of material among the computational zones (elements). In this way the mesh is incrementally relaxed, which allows for large deformations while avoiding excessive mesh distortion. The coupling of the mechanics, chemistry, and heat transport is handled using a sequential operator splitting scheme. The numerical methods employed here are similar to those outlined in a previous review.⁵⁴ 

The simulations were performed in 2D (plane strain) to allow for tractable computations. The computational domain is rendered by inserting a single pore near the center of a rectangular slab of β-HMX crystal. The pore is cylindrical in shape and filled with ambient air. The diameter of the pore is denoted as d and the planar dimensions of the crystal sample are about 25 d × 25 d. Two-dimensional axisymmetric simulations (for spherical pores) were not considered because the crystal lattice does not display this symmetry. Although full 3D calculations are certainly desirable, we believe such calculations are premature at this stage of model development.

The crystal sample, which is initially at rest, resides at an initial pressure and temperature of 1 atm and 298 K, respectively. A planar shock wave is generated by prescribing the axial component of velocity on the left-hand surface of the crystal. The prescribed velocity rises instantaneously, similar to conditions imposed at the impact face in a plate experiment. The top and bottom surfaces of the crystal sample are periodic and the right-hand surface is restrained by a rigid frictionless wall. A single shock wave transit of the slab is simulated.

The initial dislocation density of the full crystal sample is taken as 0.0307 μm^–2, which corresponds to the value determined in the model calibration.²³ For this nominal value of dislocation density, (g^α/r^α)_t=0 = 0.103 GPa and the individual slip system initial strength levels (g^α) may be found using the strength ratios given in Table II. However, the initial dislocation density field is not uniform but rather randomly distributed in space to provide a more realistic description of microstructure heterogeneity. This is done because the mesh size utilized in these calculations is small enough to expect fluctuations in dislocation density among the computational zones. Following arguments presented earlier,²³ the zonal quantity h is taken to be log-normally distributed with a mean of 0.0679 and a standard deviation of 0.0672. This prescription yields a coefficient of variation of r ∼ 1, which indicates strong fluctuations in dislocation density among individual zones.

The real viscosity of the material is augmented by an artificial viscosity for the purpose of shock capturing. This is done to spread out the shock front when the material model predicts close to zero dissipation. Such is the case for extremely sharp velocity gradients, where the viscoplastic strength law predicts negligible amounts of plastic shearing and, therefore, large elastic deformations. When a numerical viscosity is used as part of the shock capturing scheme, transients in the wave structure give rise to entropy errors.⁵⁵ Transients occur, for example, when a shock wave interacts with a material boundary or with other shock waves. The entropy errors take the form of an excess energy and a deficit in density.⁵⁶ These errors require consideration because we are interested in computing temperature fields about a collapsing pore. An estimate of this excess heating effect will be given in Sec. V.

There is a large parameter space that can be studied with regard to intragranular pore collapse and it is rather difficult to suppose which effects are the most important ahead of time. In this work, we constructed a reference case and then systematically modified that reference case to investigate the effects of stress wave amplitude, pore size, heat conduction, crystal strength, and liquid viscosity. This list, which is by no means exhaustive, was selected to map out basic responses over a range of loading conditions and to assess the effects of various material modeling assumptions. The model predictions are now examined.

The reference case involves a pore diameter of 1 μm, an imposed boundary velocity u = 1 μm/ns (which generates a peak axial stress σ₁₁ = 9.4 GPa and a shock wave velocity D = 5.057 μm/ns), and a crystal orientation that is obtained as follows: the normal of the $(1¯1¯1)$ plane is aligned with the shock direction (x-axis); the crystal is then rotated about the shock direction until the $[11¯0]$ direction forms an angle of θ with the y-axis. The reference angle was selected as θ = –38.8° to be consistent with previous work.²³ It is noted that model predictions are relatively insensitive to the value of θ when impacting normal to the $(1¯1¯1)$ plane.⁵⁷ Finally, the reference case is locally adiabatic.

For the reference case, a mesh study indicated that a mesh size of 8 nm was needed to converge the behaviors of interest. The details of this study are given below. This amounts to a fine meshing of the domain, as the pore diameter is spanned by 125 elements. To mesh the entire crystal sample, about 10⁷ elements were required. A single pore collapse simulation, run in parallel on 512 cores, required about 72 h of wall clock time.

To illustrate the nominal stresses in the crystalline phase under these loading conditions, the stress components of the wave are plotted in Fig. 2. These profiles illustrate the structure of the stress wave just prior to its arrival at the pore location (x = 14 μm). As shown in Fig. 2, the shock wave is mostly unsplit at this run distance, i.e., a distinct precursor wave has not yet formed. The kink in the wave structure (near x = 8 μm) is, however, due to the onset of precursor wave formation. If the shock is allowed to propagate a few hundred microns further into a monolithic (non-porous) sample, a distinct precursor wave would separate from the main shock wave.

The pressure and temperature fields that are generated when the porous crystal sample is shocked to 9.4 GPa are depicted in Fig. 3. In this figure, the observation windows are fixed in space and contain only the central portion of the sample. The times given for each snapshot are relative to the time at which the shock front reaches the left-hand side of the pore. Let us first consider the pressure fields in Fig. 3. Behind the front, the nominal pressure (p) and deviatoric stress (s₁₁) are approximately 6.6 GPa and 2.8 GPa, respectively. The high level of deviatoric stress is due to the rapid compression of the crystal, which produces large (elastic) strains that await relaxation by plastic flow. When the shock front crosses the pore, a release wave is emitted from the crystal-air interface (0.2 ns). The release wave is followed by a secondary stress wave that is generated when the pore is fully closed (0.5 ns). This disturbance propagates away from the pore collapse region and beyond the observation window (1.2 ns). The simulations are run until the secondary wave begins to interact with the boundaries, which allows for a post-collapse simulation time of approximately 2 ns.

The temperature fields that are given in Fig. 3 illustrate the nature of energy localization under these conditions. The narrow bands of elevated temperature that grow away from the pore are shear bands. An enlarged view of the shear band patterning is given in Fig. 4. In this figure, the temperature field is plotted later in time (2.6 ns) and for the full crystal sample. The shock wave generates large, highly non-equilibrium shear stresses in the crystal sample. These shear stresses drive dislocation motion on the slip planes, which serve to relax the crystal. As a result, the shear bands tend to grow along the crystallographic planes that are well-oriented for slip (e.g., the (011) and $(1¯02)$ planes, as shown in Fig. 4). The shear band tips are quickly heated to melting, at which point shear loads are transferred to the surrounding material. This allows for the propagation of shear bands that are filled with liquid-HMX. The shear bands may therefore be viewed as melt cracks.

The observed shear banding behavior indicates that the effects of material heterogeneities may not be as localized as one might assume. For example, the hot spots generated by the growing shear bands are considerably larger and more dispersed than the footprint of the original pore. The shear bands in the lower right-hand quadrant of Fig. 4 are approximately 50 nm in width, with temperatures ranging from 500 to 900 K. The bands propagate at a rate of ∼4.6 μm/ns and eventually reach the boundaries of the crystal sample. The simulation must be stopped at this point because severe element distortions develop when a shear band intersects a boundary. This occurs because the FE code is not equipped to advect material across the periodic boundaries. Such type of advection is needed to relax the mesh on a boundary.

The shear banding of low-symmetry HE crystals under shock wave loading has been documented in experimental work.^58,59 For example, RDX crystals that were shock-loaded up to ∼13 GPa exhibited shear bands and beaded-up volumes of material on the surfaces of recovered samples. The beads are thought to be molten material that was squeezed out of the shear bands during the shock deformation and which solidified on the surface upon release. Shear banding is also predicted in atomistic simulations of shock wave propagation in α-HMX⁶⁰ and α-RDX.⁶¹ In the atomistic simulations, quantitative analysis revealed that the shear banding regions are composed of liquid-like phases. These experimental and numerical results are encouraging to us, as this is precisely what is predicted by our HMX model: shear bands that are filled with liquid-HMX.

Previous analytical treatments have considered shear banding as a mechanism for igniting explosive materials. For example, stability analyses can be carried out to estimate the critical thickness of an explosive shear band.⁶² The role of crack growth in the ignition of mechanically loaded HE crystals has also been considered, where heating is derived from the frictional sliding of closed crack surfaces.⁶³ When coupled with a model of crack burning, the internal pressure of the crack can lead to unstable crack growth and large scale reaction. The shear banding behavior that is predicted in our simulations is somewhat connected to these ideas. Although the shear bands remain inert for the reference case, it will be shown that slightly higher values of liquid viscosity can lead to reactive shear bands.

The pore collapse process is illustrated in Fig. 5, where (a) stress concentrations drive small amounts of localized plasticity as the wave front begins to interact with the pore; (b) shear localization (melt-cracking) contributes to the formation of two fluid jets; (c) the jets impinge on the downstream side of the pore and are re-directed towards each other; (d) the jets interact and create a backward-moving jet, which collides with a small tertiary jet that formed on the upstream side; (e) the collision of these jets yields a pressure of ∼30 GPa, which drives up the temperature enough to react a small amount of material (the red-colored region in the figure); and (f) the pore is fully closed and the reaction zone is beginning to grow in size. The pore collapse time for the reference case is ∼0.5 ns.

In the classical view, the symmetry of a collapsing pore is determined by the dimensionless Reynolds number (cf. Ref. 18 and work cited therein). The Reynolds number for pore collapse may be written as

where p is the shock pressure and ρ and η are the density and viscosity of the medium, respectively. When Re ≪ 1, the process is dominated by viscous forces and the pore is expected to collapse in self-similar manner, maintaining spherical (cylindrical) symmetry. When Re ≫ 1, the collapse is dominated by inertial forces and jetting is expected.

Let us consider the Reynolds number for the reference case. The viscosity of the solid phase is derived from the rate-sensitivities of crystallographic shearing, $∂γ̇α/∂τα$⁠, where τ^α denotes the projection of the stress tensor onto the α^th slip system. When subjected to strong, sudden loading, the initial viscosity of the solid phase is estimated to be approximately 100 kP. The solid viscosity is, however, inversely proportional to the dislocation density and stands to be reduced by a large amount as the dislocation density multiplies up. Allowing for a 1000-fold increase in dislocation density, Re ≤ 0.4 for pore collapse in the crystal phase. This indicates that the pore should maintain cylindrical symmetry during its collapse. However, the closure process depicted in Fig. 5 is highly asymmetric. This asymmetry is caused by localized melting of the β phase around the periphery of the HMX/air interface. Jets of liquid-HMX are injected into the pore because the liquid phase is strengthless and of much lower viscosity (5.5 cP). Therefore, in the reference case, pore collapse is initially dominated by the high viscosity of the solid phase, which tends to maintain cylindrical symmetry. This symmetry is broken when the solid phase is partially melted.

Entropy errors are generated when shock interactions are captured using an artificial viscosity. This is relevant to the current problem because the liquid jets impinge on the downstream side of the pore at speeds of 2–3 μm/ns. To quantify this effect, additional simulations were performed to compute the excess heating at a driven boundary (impact face). Drive velocities of 1–1.5 μm/ns were simulated, as this corresponds to the impact of one crystal slab, traveling at 2–3 μm/ns, onto another crystal slab that is stationary. The boundary temperatures were 100–200 K higher than the plateau temperatures found a few mesh lengths away from boundary. Therefore, the temperatures generated by jet impingement during pore collapse are likely too high, by about 100–200 K. Although the entropy errors are unfortunate, they are somewhat unavoidable given the standard shock capturing methods employed in most hydrocodes.

The evolving morphology of the reaction zone is illustrated in Fig. 6. In this figure, zonal mass fractions of product gas are plotted at five time instances after full closure of the pore. The reaction zone corresponds to the red-colored region in Fig. 5(f). Gaseous products evolve out of the liquid phase, wherein the transition between unreacted and fully reacted regions is resolved over a distance of ∼35 nm, or ∼10 distorted elements. The growth of the reaction zone is dictated by the temperatures achieved by mechanical deformation and exothermic heat release, as these calculations are adiabatic. It is noted that the shear bands simulated in the reference case are not hot enough to achieve significant degrees of reaction on this time scale. Simulations run out to longer times might allow for reaction in the shear bands, although such calculations should consider the conduction of heat away from the bands. The relatively inert behavior of the shear bands is, of course, tied to model assumptions pertaining to the liquid decomposition kinetics and viscosity.

The sensitivity of model predictions to the mesh size is now considered. To quantify the overall reactivity of a sample, the relative mass of products is defined as ξ ≡ m_p/m_pore, where m_p is the mass of the product phase and m_pore is the pore mass (i.e., the reference mass of crystal that would fit inside the initial pore). In Fig. 7(a), the relative mass of products is plotted as a function of time for successively finer mesh sizes, Δx. Here, and in all subsequent product curves, the origin of the time coordinate corresponds to the arrival of the shock at the upstream side of the pore. The product curve is reasonably well converged at a mesh size of 8 nm. In this case, about 10% of the pore mass is reacted during a post-collapse compression time of 2 ns.

In the absence of sufficient regularization, the simulated localization behavior will always depend on the mesh size. Mesh-dependent behavior can be avoided by including, for example, dissipation (viscous flow, heat conduction) and/or non-local material behavior. Although the plastic response of each phase is viscous, the reference case simulations are adiabatic. To assess the mesh-dependence of the shear banding behavior, the average width of one of the larger shear bands (the forward-most band in the lower right-hand quadrant of Fig. 4) was computed over a range of mesh size. The extent of the shear band was identified by thresholding the temperature field. The average width was then computed by sampling multiple locations along the band length at the final simulation time. As shown in Fig. 7(b), the average width of the leading band is reasonably well converged at Δx = 8 nm. A similar procedure was used to compute the average width of the finer shear bands (crack branches) in the lower left-hand quadrant of Fig. 4. The average width of the finer bands (plotted in Fig. 7(b)) is again reasonably well converged at Δx = 8 nm. Therefore, a mesh size of 8 nm appears to be sufficient to resolve the behaviors of interest. This is the mesh size that has been used in the reference case and all simulations presented in Secs. V B–V F.

We now consider a challenging numerical issue: artificial heat transfer. The solution tends to be smeared out in space due to material advection among computational zones. For example, if part of a hot (reacted) zone is advected into a neighboring zone that contains cold (unreacted) material, the temperature of the initially cold zone will be uniformly raised due to the mixing of hot and cold volumes. There is not a physical basis for this heat transfer; it is simply an error introduced by the numerical treatment.

An attempt was made to circumvent this problem by shutting off advection among zones undergoing reaction. The idea was that advection among fully reacted zones and among unreacted zones is benign, whereas advection among partially reacted zones is problematic. This effort was not successful because the partially reacted zones are precisely those that require mesh relaxation, as they undergo large volume expansions during the decomposition reaction. Performing the calculations with the aforementioned advection constraint led to sliver elements and mesh tangling when zones started to react. The solution was then compromised and had to be stopped. The fine mesh size that is utilized in this work helps to mitigate against artificial heat transfer, as advection errors scale with mesh size. However, there is still some amount of artificial heat transfer that stems from material advection and this remains an open issue.

It may be useful to devise a treatment that distinguishes the temperatures of the species (HMX and product gas) instead of assuming an equilibrated mixture temperature. Mixing among computational zones could then be handled on a species-by-species basis. In this way, when product gas is advected it would only be allowed to mix with other product gas and not serve to heat up the cold unreacted phase. A scheme like this has not yet been implemented.

Model predictions are now examined over a range of axial shock stress (σ₁₁). All other aspects of the calculation are identical to the reference case. In Fig. 8, the shear banding (melt cracking) regions that are induced at 6.5, 9.4, and 10.7 GPa are illustrated by plotting the fraction of melted material, f_m. In each of these calculations, the stress wave arrives at the pore as a mostly unsplit wave with essentially no precursor separation (cf. Fig. 2). At each stress level, a pool of molten HMX is formed in the pore collapse region. At 6.5 GPa, a few shear bands grow out of the pore collapse region and branch off to form finer localization bands. When the stress is increased to 9.4 GPa, the shear bands grow in four general directions and exhibit greater amounts of branching. At 10.7 GPa, the spacing of molten sheets is yet finer. At the highest shock stress, numerous bands are nucleated away from the pore (e.g., 5–10 μm below and above it). The nucleated bands are not topologically connected to the fluidized pore collapse region or to shear bands that grew out of it. Attempts to simulate stress wave amplitudes greater than 11 GPa were not successful, as shear bands were immediately nucleated at the driven boundary (impact face). These cracks go on to intersect the lateral boundaries and induce element distortions that are too large to handle with the FE code.

The crystal sample reactivities, as measured by ξ(t), are plotted in Fig. 9. At 6.5 GPa, there is little reaction on this time scale. At 9.4 GPa, approximately 8% of the pore mass is reacted during a compression time of approximately 2 ns. Increasing the stress to 10.7 GPa does not appear to have a large effect on reactivity. It may be noted, however, that the product curve at 9.4 GPa is concave-down (decelerating) over most all of its trajectory, whereas the product curve at 10.7 GPa exhibits a slight up-turn at about 1.8 ns. It would be significant if this increased reaction rate is carried out later in time. Barring such speculation, the simulations indicate that the shock-loaded crystals exhibit modest extents of reaction for a shock stress of ∼10 GPa and a compression time of a few nanoseconds.

Micron-sized pores are of primary interest when considering the shock initiation of HMX-based formulations.¹¹ To begin scoping out effects related to pore size, the collapse of a smaller pore (d = 0.5 μm) was simulated. All other aspects of the calculation are identical to the reference case.

The temperature fields for the smaller pore are qualitatively similar to those shown earlier for the larger pore (Fig. 4) and are not repeated here. Again, for the smaller pore, products are concentrated in the pore collapse region and not in the shear bands. The product curves, ξ(t), for the smaller and larger pores are given in Fig. 10. The smaller and larger pores yield similar levels of reactivity on this time scale when the product masses are normalized. The smaller pore begins to form products earlier because the collapse time is shorter. To compare these measures of reactivity on an absolute basis, one must account for the difference in pore mass. The product mass scales as (m_p,2/m_p,1) = (ξ₂/ξ₁)(d₂/d₁)²; therefore, the product curve for the smaller pore should be scaled by a factor of 1/4 to allow for comparison with the larger pore on a product mass basis. From Fig. 10, the product mass formed by the smaller pore is roughly one-third of that formed by the larger pore after 2 ns of compression time. Considering the slopes of the product curves for t > 1 ns, the mass production rate (⁠$ṁp$⁠) of the smaller pore is found to be ∼0.5 times that of the larger pore.

With reference to the pore closure process illustrated in Fig. 5, the collapse of the smaller pore is very much similar. The main difference is that the jets formed within the smaller pore are slightly slower than those produced by the larger pore. The time-averaged speed of the jet tip is ∼2.74 μm/ns for the larger pore and ∼2.60 μm/ns for the smaller pore. The higher reactivity of the larger pore (on a product mass basis) is therefore attributed to larger amounts of inelastic work and higher jet speeds. Future studies should consider a wider variation in pore size to further scope out its effects on localization and reactivity.

The characteristic length of heat diffusion in a medium with thermal diffusivity λ may be taken as $λt$⁠, where t is the elapsed time. Considering a post-collapse simulation time of 2 ns and the assumed diffusivity, the expected heat diffusion distance is about 20 nm. While this is relatively small compared to the hot spot size in the pore collapse region, it is significant when compared to the shear band width (cf. Fig. 7). The pore collapse simulations have so far neglected conductive heat transfer. To assess this assumption, the crystal model was exercised with Fourier heat conduction. The thermal conductivity of each phase was assumed to be constant and equal to the β-phase ambient value. This would appear to be a conservative estimate, as the measured thermal conductivity of the β phase decreases with temperature.¹¹ Although the pressure dependence has not been measured, analytical efforts suggest non-monotonic behavior wherein the conductivity increases to about 0.75 W/m-K at 6 GPa and then decreases below the ambient value (0.5 W/m-K) when the pressure is increased beyond 8 GPa.⁵² 

The temperature field that was computed with heat conduction was rather similar to that computed using the adiabatic model (Fig. 4). The temperature gradients were, of course, less steep in the conductive case, but the hot spots found in the pore collapse region were of similar size and temperature. There were some subtle differences in the shear band patterning predicted by the conductive model. This is illustrated in Fig. 11(a), which depicts temperature profiles that were extracted from a slice of the crystal. The slice spans the shear bands located in the lower right-hand quadrant of Fig. 4. In the conductive case, there are two fewer shear bands over the same slice and those present exhibit depressed temperature peaks and more diffuse profiles.

The product curves for the adiabatic and conductive cases are plotted in Fig. 11(b). The conductive case is only slightly less reactive on this time scale due to the flux of heat away from the central hot spot. This indicates that the adiabatic assumption is reasonable for the reference case, as the results of interest are not changed in a substantial way when heat conduction is included and the simulation duration is limited to a few nanoseconds. Simulations that seek to address longer time scales (e.g., tens to hundreds of nanoseconds) and hot spot burning should, however, treat heat conduction in a meaningful way.

The sensitivity of the model to the flow strength of solid phase is now considered. The crystal model, which has been used in all simulations up to this point, accounts for the elastic/plastic anisotropy of the β phase and the time-dependent nature of plastic deformation. This is thought to be important because (i) the mechanical responses of low-symmetry crystals tend to be highly dependent on the orientation of the sample relative to the loading axes and (ii) many solids exhibit a sharp increase in flow stress when subjected to a rapid change in strain rate. In regard to the problem of shock-induced pore closure, if the time scale of plastic relaxation is comparable to the time scale of pore collapse, high shear stresses will persist until they are relieved by inelastic deformation. Despite these considerations, a large majority of HE shock simulation efforts view the elastic/plastic responses of crystalline phases as isotropic and insensitive to strain rate. Given the sparsity of strength data and challenges associated with crystal model development and parameterization, it is understandable why simplified models are assumed. It will be shown, however, that conventional strength models can fail to capture the full extent of localization behavior, which is essential to the production of hot spots and the ignition of heterogeneous explosives at moderate stress levels.

To investigate the effects of solid flow strength, a pore collapse simulation was performed using a conventional isotropic/rate-independent strength model for the β phase. Here, the yield strength is taken as Y = Y₀(1 + βϵ^p)ⁿ, where Y₀ is the initial yield strength, ϵ^p is the effective plastic strain, and β and n are hardening parameters. In the literature, it has been common to assume a constant yield strength in the range of 0.060–0.180 GPa.^25,26,65 To consider the effects of somewhat higher strength and weak strain-hardening behavior, the following parameters were selected: Y₀ = 0.300 GPa, β = 0.060, n = 1. All other aspects of the reference-case simulation remained the same.

The temperature fields that are computed when pore collapse is simulated using the isotropic/rate-independent model are shown in Fig. 12. The pore is now collapsed by a single smooth jet that produces a symmetric configuration of hot spots. These temperature fields may be compared to those computed earlier using the crystal model (Figs. 4 and 5). In the crystal calculations, shear banding is a prominent feature of the deformation, whereas in the isotropic/rate-independent calculations, there is no shear banding whatsoever. The differences in localization behavior are largely attributed to the strain-rate-dependence of the strength model. Although rate-dependent flow stresses are often seen as tending to suppress localization (given that higher stresses are required to drive higher strain rates), this notion can break down in the presence of thermal softening and melting. Such is the case in our rate-dependent calculations, where relatively large flow stresses and small plastic strains produce enough dissipation to melt the material located at the tips of the bands. Loads are then shed onto the surrounding material, which allows for the propagation of the localization (melt) bands. In the rate-independent case, the stress state is forced to remain on a strain-rate-independent yield surface and plastic strains are computed according to a consistency condition. For the parameters chosen above, the dissipated work is insufficient to trigger shear banding. It may be possible to induce shear banding in the rate-independent model by prescribing higher yield strengths or a spatial distribution of yield strength in the material around the pore (i.e., seeded weak regions). The former approach would, however, be at odds with quasistatic stress-strain curves⁴² and the observed relaxation behavior under shock loading.⁴⁰ It is possible to induce shear banding in the pore collapse calculations by assigning sufficient strain-rate-dependence to an isotropic model. In this case, the bands form on the planes of maximum shear, whereas, in the crystal model, the details of plastic flow and band structuring are sensitive to the slip plane orientations.

The isotropic/rate-independent model and the crystal model clearly lead to different predictions for the post-collapse temperature field and localization behavior. This difference is quantified by the temperature histograms given in Fig. 13. These histograms were constructed by binning the mass of a sample according to temperature and normalizing those bins by the pore mass. As shown in Fig. 13, the isotropic model is biased toward lower temperatures (notice the large peak in-between 400 and 500 K), whereas the crystal model predicts higher temperatures due to shear localization (600–1200 K) and exothermic reaction (>2000 K). For equivalent energies imparted to each sample, the crystal model predicts higher degrees of localization and, therefore, hot spots that are larger or higher in temperature.

The extent of reaction that is predicted by the isotropic/rate-independent model is effectively zero on this time scale, despite the formation of hot spot temperatures of greater than 1000 K. Therefore, the strength law has a significant effect on the character of the shock-induced plasticity, the peak temperatures that are generated, and the occurrence of decomposition reactions during pore collapse. Given the agreement of the crystal model with experimentally measured precursor wave profiles⁴⁰ and the shear banding that is observed in recovered HE crystals,^58,59 we believe the crystal model provides an improved description of the strength behavior of the solid phase under shock wave loading. The differences in reactivity indicate that crystal flow strength is an important consideration when simulating the behavior of HE crystals at the length scale of individual pores or grains.

The shear bands that grow out of the pore collapse region are filled with liquid-HMX. In the reference case, the liquid viscosity was taken as constant and equal to 5.5 cP. This corresponds to the liquid viscosity that was computed from atomistic simulations at 800 K.⁴⁸ The sensitivity of the model to the assumed liquid viscosity is now considered. As a comparative case, the liquid viscosity was increased to 22.0 cP, which corresponds to the atomistic value at 700 K. All other aspects of the reference case remained the same.

The temperature field of the shock-loaded sample with higher liquid viscosity is depicted in Fig. 14. The shear bands are now wider, fewer in number, and hotter than those of the lower viscosity case (Fig. 4). When the liquid viscosity is increased, the model predicts reactions in both the pore collapse region and some of the shear bands (e.g., the forward-most band in the lower right-hand quadrant of Fig. 14). This is due to higher levels of dissipation in the liquid phase.

The product curves for the higher and lower viscosity cases are given in Fig. 15(a). The initial trajectory of the product curve is controlled by reactions that occur in the liquified pore collapse region. As such, the higher viscosity case is more reactive from the start. At about 2 ns, the higher viscosity case exhibits a sharp increase in overall reaction rate. This corresponds to the onset of reactions in the shear bands. For a post-collapse simulation time of about 2 ns, about 80% of the pore mass is reacted.

The temperature histograms (Fig. 15(b)) confirm that a relatively modest increase in liquid viscosity has a noticeable effect on the amount of mass contained in the unreacted high-temperature tail of the distribution. This is in agreement with previous computations,²² which showed the slope of the high-temperature tail scales with the viscosity of the medium. The model predictions of reactivity are thus quite sensitive to the viscosity of the liquid phase. This is because the liquid contained within a localization band can be subjected to exceptionally high shearing rates. Since the liquid phase experiences temperatures ranging from about 550 K (initial melting at 1 atm) to greater than 1000 K, future modeling efforts should incorporate the temperature-dependent liquid viscosity that was computed in atomistic simulations, i.e., $η=η0 exp(Ta/T)$⁠, where η₀ = 3.1 × 10⁻⁴ cP and T_a = 7800 K.⁴⁸ The effect of pressure on liquid viscosity should also be considered, although we are not aware of any such data. Along these lines, one might also consider the pressure- and temperature-dependence of the product gas phase viscosity.⁶⁶ 

A computational framework has been developed to study the chemical reaction of HMX crystal under shock wave loading. The framework incorporates a material model that accounts for the crystal mechanics and decomposition reactions that are driven by the temperature field. This allows for the direct numerical simulation of coupled thermal/mechanical/chemical responses, including hot spot formation and reaction, under well-defined shock loading conditions. The framework was used to study intragranular pore collapse in β-phase HMX crystal.

A reference case involved a pore diameter of 1 μm, shock loading normal to the $(1¯1¯1)$ plane, a shock stress of 9.4 GPa, and locally adiabatic conditions. In the reference case, shear banding was found to be an important mode of localization. The shear bands grow out of the pore region due to the presence of highly non-equilibrium shear stresses, which drive plastic flow along certain crystallographic slip planes. As a result, the hot spot footprint is noticeably larger than the initial defect size. Crystal melting was found to be responsible for the rapid propagation of shear bands and the asymmetric collapse of the pore, despite its small size. The model predicts a modest extent of reaction (approximately 10% of the pore mass) during a post-collapse compression time of ∼2 ns. The reactions are concentrated in the pore collapse region and are due to elevated temperatures produced by material jet formation and impingement. In this case, the shear bands are not hot enough to react on a nanosecond time scale.

The reference case was then modified in a systematic manner to map out model sensitivities to shock strength, pore size, and various modeling assumptions. The findings of this study are summarized as follows:

• Increasing the shock stress from 6.5 GPa to 9.4 GPa leads to a noticeable increase in the extent of shear localization and reaction on a nanosecond time scale.

• Decreasing the pore size from 1 μm to 0.5 μm does not have a significant effect when reactivity is measured on a relative basis (i.e., when product masses are normalized by the pore mass). When considering an absolute basis, the smaller pore yields a gas production rate (⁠$ṁp$⁠) that is reduced by a factor of two.

• The diffusion of heat away from hot spots in the reference case is an effect that is rather small for nanosecond duration shock simulations. However, modeling efforts that seek to address longer time scales (e.g., tens of nanoseconds) and hot spot burning should account for conductive heat transfer.

• The model predictions of shear localization and reactivity are quite sensitive to the strength behavior of the crystalline phase. Under a shock stress of 9.4 GPa, the viscoplastic crystal model predicts extensive shear banding and a modest extent of reaction; whereas a conventional strength model (i.e., one that is rate-independent and isotropic) predicts no shear banding or reaction whatsoever.

• The model predictions of shear localization and reactivity are quite sensitive to the viscosity of the liquid phase. A modest increase in the assumed viscosity leads to a sharp increase in overall reactivity due to the occurrence of reactions in the shear bands.

Much remains to be done as we progress towards a predictive model of coupled responses in shock-loaded HMX crystal. For example, the material model would be improved through a better treatment of melting and liquid phase behavior. In this regard, one might include kinetics for the β → liquid phase transformation, a more appropriate description of the liquid phase EOS, and explicit kinetics for thermal decomposition from the liquid phase.¹⁶ A refined treatment of certain solid/liquid thermal properties is also needed (cf. Sec. III C). Perhaps the most important (and straight-forward) improvement that should be made is incorporating the pressure- and temperature-dependence of the liquid viscosity, as the simulated reactivity is quite sensitive to this property. To assess hot spot criticality (self-sustained burning), future work should consider conductive heat transfer and shock loading states that persist beyond a few nanoseconds.

The long term goals of this work include the simulation of 3D spherical pore collapse, small arrays of intragranular pores, and HMX grain aggregates. However, we believe these computationally expensive efforts should be deferred until the issues related to material model development and hot spot criticality are more fully addressed.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344 (LLNL-JRNL-664305).