Hot spot ignition and growth from tandem micro-scale simulations and experiments on plastic-bonded explosives

Heterogeneous energetic materials (EMs) are composites of organic crystals and other materials such as plastic binders, metals, and additives. For example, plastic-bonded explosives (PBXs), used in munitions and propulsion devices, consist of organic CHNO compounds such as octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX)¹ with a viscoelastic binder such as hydroxy-terminated polybutadiene (HTPB)² or a polyurethane (PU) produced by polymerizing HTPB,³ possibly with additives to ensure stable formulations. A PBX is expected to remain unresponsive to unintentional mechanical insult during storage, transportation, and handling, yet maintain effectiveness during designed operational conditions. To balance safety and performance, it is important to be able to predict thermo-mechanical stresses that initiate chemical reactions in a PBX. It is well known that the EM microstructure strongly influences the thermo-mechanical properties that control the initiation and growth of reactions, and, hence, the overall sensitivity and performance.^4,5 Therefore, predicting the structure–property–performance (S–P–P) relationships of EM is a subject of active study in current energetics research.^6–8 Predictive approaches at the system (i.e., macro-scale) level should rely on well-tested and validated mechanical and thermochemical sub-models (micro-/meso-scale) to describe the macroscale response of a PBX to shocks and other transient loads. Therefore, establishing S–P–P linkages for heterogeneous EM with complex stochastic microstructures inherently relies on multi-scale models, where phenomena at smaller scales must be accurately represented at larger scales.^9–13 This work examines the ability of currently available models to provide physically correct descriptions of shock-induced initiation in the HMX-based PBX by conducting head-to-head comparisons of experimental and numerical simulations on single isolated energetic crystals embedded in a binder, with a 20 GPa shock input, examining phenomena at the micrometer spatial and nanosecond time scales. The analysis reveals the extent to which current models can reproduce experimentally measured temperatures and hot spot initiation and growth phenomena and points to advances needed in both the experimental and numerical aspects to establish a model validation framework for the shock initiation of EM.

One of the complexities of predicting critical hot spots and modeling their contributions to the overall performance of an EM is that the overall shock-to-detonation transition (SDT) involves processes separated across disparate length and time scales.¹⁴ The thermo-mechanics of critical hot spot formation and evolution at the microscale is key to the initiation of PBXs. Several mechanisms of hot spot formation at such microstructural features have been hypothesized, including plastic dissipation, shear banding, friction between grains and crack surfaces, and the collapse of pores.^4,15 For strong shock conditions, shock focusing and pore collapse are considered to be the predominant mechanisms for hot spots, with shear bands playing an important role as well, particularly, at moderate shock strengths.^16–18 Critical hot spots,^19–22 local elevated-temperature regions that are large enough and hot enough, can grow and coalesce rapidly, leading to deflagration and subsequent detonation in PBXs. The initial hot spot development and growth associated with shock initiation occurs at length scales of nm–μm and time scales of ns–μs.^21,23,24 Direct observations of such small length scale and short time scale phenomena are sparse in the literature, but necessary to acquire physical insights and to validate computer models. Recent advances in the direct observation of hot spot initiation and growth in a single or small cluster of crystals have been made by Johnson et al.²⁵ Johnson et al. studied both nominally defect-free HMX single crystals and HMX polycrystalline aggregates with a variety of defects. They found that for single crystals, hot spots developed primarily at crystal–binder interfaces, especially along sharp edges or corners, with an initial temperature rise of about 4000 K. These hot spots could grow and consume the entire crystal. For crystal aggregates, hot spots developed at crystal–binder interfaces and also at the locations of heterogeneity, where the initial temperature rise was about 6000 K. However, these 6000 K regions had a small thermal mass and cooled quickly so that subsequently, the defective crystals behaved much like the single crystals. The physics of such hot spot initiation and growth phenomena in single and small collections of crystals needs to be further elucidated to understand, design, and control the effects of microstructural heterogeneities on the macro-scale sensitivity of energetic materials.

Numerical studies of hot spot mechanisms have been extensive and have examined void collapse^18,26–30 and shear band induced localization. Most of those studies assumed idealized geometries for voids such as circular/elliptic/triangular voids in a background homogeneous HMX material.^{18,26–28,31–34} Springer et al.¹⁸ studied the effect of shock pressure and pore morphology on the collapse of voids in HMX and the subsequent growth of hot spots. One of their findings was that for moderate shock pressures (2–10 GPa); the hot spot growth mechanism was mainly through the formation of shear bands. At higher pressures (20–40 GPa), hot spots were formed mostly by hydrodynamic void collapse. Springer et al.¹⁸ also studied the collapse of voids at different locations in an HMX-Kel-F binder system and showed that the hot spot criticality depended on where the voids were located in the crystal–binder system. While those studies are in a 2D setting, 3D void collapse calculations have also been performed^35–37 for voids of various shapes (spherical, cylindrical, etc.), to show the generation of intense baroclinic vortical structures during and following the collapse, correlating the surface area of the voids with the vorticity and subsequent growth rates of reaction fronts. Head-to-head comparisons of void collapse simulations against molecular dynamics (MD) calculations have also been performed recently.³⁸ However, head-to-head comparisons between experimental observations and simulations, particularly for crystals embedded in a polymeric binder, are lacking. A step in this direction was taken in a recent work by Rai et al.³⁹ comparing head-to-head simulations and experimental observation of the collapse of millimeter-sized pores in the inert, transparent, polymeric material poly-methyl methacrylate (PMMA). The 3D simulations of pore collapse in PMMA across a range of shock strengths showed good agreement for pore collapse profiles and times, indicating transitions from shear band dominated to hydrodynamic jetting dominated localization mechanisms as the shock strength increased. Tandem experiments and simulations of reactive hot spots in energetic crystals or crystal–binder mixtures are a much-needed line of work to ascertain whether emerging physics-based (typically, MD-guided³⁸) models of energetic materials accurately capture observed hot spot phenomena, and to guide further improvements in the thermophysical models of EM to make reliable and accurate predictions and S–P–P linkages through numerical simulations.

A head-to-head comparison of experiments and simulations is difficult from both the experimental and computational standpoints, but here, we demonstrate such comparisons against experiments conducted using the procedure and setup recently described in Johnson et al.²⁵ The new feature of our experiments is that prior to shock measurements, the structure of the polymer-embedded crystal was characterized with x-ray computed tomography (CT) with sub-micrometer resolution. These high-resolution 3D structures then formed the basis for simulations.

To date, well-resolved numerical studies of microscale crystal–binder combinations have been mostly performed on idealized geometries and in 2D.^{18,26–28,31–34} Furthermore, there remain large uncertainties with regard to several sub-models that are key to physically realistic computations of hot spot initiation and evolution, including material (strength)³⁸ models and reaction chemistry models.⁴⁰ Fundamental, physics-based, MD calculations in quasi-2D systems^41–43 are possible only for small length (10s of nm) and short time (picosecond) scales and cannot approach experimentally observable (micrometer-scale/100 ns) spatial/temporal dimensions. While continuum and atomistic numerical studies aid in analyzing and understanding hot spot mechanisms, they seldom represent the complex ignition phenomena that occur in real PBX microstructures. A PBX can contain aggregates of crystals of various sizes, shapes, and types of defects, as prominently seen in images by Chen et al.⁴⁴ Often, two or more hot spot generation mechanisms may occur simultaneously and at different length scales.^18,45 Mesoscale numerical calculations can capture the variety of shock and chemical reaction localization that occur in real PBX microstructures, provided they adequately resolve the material interfaces to capture both small- and large-scale phenomena. Most meso-scale simulations, however, cannot be performed at sufficient resolution to accurately treat localized physics at interfaces; even for single void calculations, the resolution requirements are quite stringent.³⁶ Some works that have attempted simulating shock initiation of actual microstructures acquired from images are that of Barua et al.^46,47 and Rai et al.^32,48 However, the numerical models used in those works are based on multiple approximations and need to be validated by experiments. Thus, there is a need for well-tested and systematically vetted numerical models for capturing hot spot ignition mechanisms and growth dynamics at the meso-scale in realistic EMs, resolving the intricacies of the microstructure and adequately capturing key interfacial mechanisms at the crystal–void, crystal–crystal, crystal–binder, and binder–void interfaces where energy localization occurs. Such models need to be calibrated and validated against accurate, well-controlled experimental studies of specific crystal structures and binders while directly observing hot spot evolution with high temporal and spatial resolutions. This paper attempts to move toward such a coordinated program of improving both experimental and numerical capabilities through a head-to-head comparison between micro-scale simulation and experiments.

The present paper extends work by Das et al.⁴⁵ on the numerical calculations of high shock strength loading of an isolated HMX crystal embedded in a polyurethane (PU) binder, and experimental work by Johnson et al.²⁵ on the visualization and measurement of hot spots in an HMX–PU system. In the present work, we use x-ray CT characterization of the HMX–PU samples as input to 2D simulations. Visual comparisons of the origin and growth of hot spots in crystal geometries are done both qualitatively and quantitatively. Techniques are developed to compare 2D temperature fields from the simulation with the radiance-based temperature values from the experiment. The setup and methodology used in this work pave the way for future 3D and high-fidelity numerical mesoscale calculations on a PBX that can be validated for accuracy through experiments.

This paper is structured as follows. In Sec. II (Methods), the setup and methodology adopted for the experimental and numerical studies are described, and the data analysis and techniques for the head-to-head comparison are outlined. In Sec. III (Results and discussion), three individual polymer-bonded HMX structures are analyzed, comparing experimental and simulation results both qualitatively and quantitatively, with an eye to revealing insights into specific mechanisms of hot spot initiation and growth. The different types of energy localization mechanisms are discussed. Finally, the key findings of this study and directions for future work, both from a numerical and experimental standpoint, are discussed in Sec. IV (Conclusions).

The procedure for our head-to-head comparison is illustrated in Fig. 1, where starting from x-ray CT images of a deconstructed PBX consisting of HMX crystals and PU binders [Fig. 1(a)], the right branch illustrates the experimental route [Figs. 1(e)–1(g)] and the left branch [Figs. 1(b), 1(c), 1(d) and 1(g)], the simulation route. Experiments [Fig. 1(e)] used 25 μm thick, 3.5 km/s laser-launched Al flyer plates to produce a planar impact at the polymer surface. The input shock pressure to the polymer was 20 GPa, and the pressure in HMX was estimated at 25 GPa. Figure 1(f) shows hot spots observed via their intense thermal emission. The CT images were the basis for numerical simulations performed using an in-house code SCIMITAR3D.^49,50 The simulations were 2D due to computer limitations, so they were performed on selected cross-sectional CT images from each PBX sample [Fig. 1(b)]. The simulations used physical parameters for HMX and for Estane® 5703 (Lubrizol) polymer, because it has similar chemical, shock, and mechanical properties to PU used in experiment, and there are reliable material properties for Estane in the open literature^51,52 The simulations predicted the experimental observables, the locations where hot spots were produced, the radiance of the emitted light and the time-dependent temperature profiles.

SCIMITAR3D, an in-house multi-material flow solver was used for all simulations. It is a well-tested and validated solver for computing the reactive shock dynamics of energetic materials^34,36,37,48 and employs a Cartesian grid-based sharp-interface framework for compressible multi-material flows.^49,50 The governing equations for the motion and deformation of the crystal and binder materials are cast in the Eulerian form⁵³ and are solved explicitly using high-order shock-capturing schemes. Chemical reactions are modeled using Arrhenius equations and an operator splitting technique is employed to implicitly integrate the stiff reactive source terms. Unlike reactive burn models,⁵⁴ the heat of reaction is not incorporated via detonation Hugoniots. Instead, the shock Hugoniot (resulting from the same equation of state) is used for the products and the reactants, and the heat of reaction is added as a source term in the global conservation law for specific total energy. The numerical framework used in this work has been described in detail in previous works^{32,39,50,55–57} and validated against experiments^39,50,51 and molecular dynamics simulations³⁸ for high-speed multi-material shock and impact problems; details of the numerical treatment can be obtained from the above-cited previous works. To solve the reactive shock dynamics of the HMX-Estane system in the present work, the CT images of embedded HMX were input to the simulation using image processing approaches described in a previous work.⁵⁸ The interfaces between the crystal, binder, and void were embedded in the Cartesian grid and tracked using the level-set method,⁵⁹ and the ghost-fluid method (GFM)⁶⁰ was used to impose appropriate boundary conditions at the material–material or material–void interfaces.⁵⁰ An adaptive local mesh refinement (LMR) scheme was employed to adequately resolve shocks, reaction fronts, and interfaces.⁶¹ 

An MD-guided rate-dependent elastoplastic material model of the Johnson-Cook form developed recently³⁸ was employed for HMX; Estane was modeled as a viscoelastic material.⁵¹ Both HMX and Estane material models employed are described in the supplementary material. HMX was modeled as a reactive material, employing a 3-equation reduced-order reaction model of Tarver;²¹ the Estane binder was modeled as an inert material. Note that the material models for HMX are the most up-to-date ones currently available in the literature. The rate-dependent J2 plasticity (Johnson-Cook) model employed here for HMX was shown in a recent work to produce shear bands and hot spot behavior that closely reproduce those obtained in the MD simulations of pore collapse in HMX.³⁸ Among the available reaction chemistry models for HMX,⁴⁰ the Tarver 3-equation model was also shown to most closely reproduce experimental data on criticality curves for neatly pressed HMX microstructures in recent works.^7,8 Thus, the set of properties used for HMX in the present work can be considered the best available model for that material, although further work is necessary to refine these models so that they are robust and applicable to a wide range of temperature and pressure regimes.³⁸ The present work exercises this recently assessed model for HMX to determine its ability to capture the experimentally observed broad features of hot spot initiation and growth.

The process of creating levelsets from images is illustrated in Figs. 1(b) and 1(c). Cross sections of the CT scans of crystals were examined and the sections with the most representative defect structures and populations, determined by visual inspection, were employed. The cross-sectional images of the crystal [Fig. 1(b)] were first pre-processed prior to being supplied to the flow solver. This involved rotation, cropping, and denoising. These preprocessing steps were conducted using MATLAB™ to orient the desired section of the image in the plane of the incident shock. The processed output is a grayscale intensity field. An image-processing framework⁵⁸ in SCIMITAR3D was used to create levelsets from this grayscale image field. First, the processed image was denoised using a speckle reducing anisotropic diffusion (SRAD) method.⁶² Then, the crystal shape was extracted using active contouring⁶³ to obtain levelsets delineating the edges of all material interfaces. Once the initial image was placed in the computational domain [Fig. 1(c)] and described in the form of a levelset contour, the subsequent evolution of the interfaces under shock loading was tracked using the evolution equation⁵⁹ 

where $ϕ$ represents the levelset field and $un$ is the velocity of the corresponding material interface. Details of the implementation of the narrow-band levelset tracking in the current framework can be found in Refs. 49, 50, 57, and 64. In this work, the levelset approach for tracking isolated interfaces was modified to accommodate the crystal–binder interface, which involves two coincident interfaces that can locally separate by debonding, viz., the crystal–void and binder–void interface. For this purpose, separate levelset fields were used to describe the crystal–void and binder–void interface. When the binder and crystal are in cohesion, the crystal–void and binder–void interfaces are exactly coincident, with no intervening void space. When the binder and crystal are not in perfect contact, a void space may be formed in the debonded region. Further details on handling the crystal–binder, crystal–void, and binder–void interfaces using the levelset approach can be found in a previous work.⁴⁵ 

Appropriate boundary conditions were applied at the embedded material interfaces⁵⁰ and the domain boundaries to accord with the experimental loading conditions. Frictionless contact was assumed at the crystal–binder interfaces, i.e., the continuity of normal velocity and normal traction was enforced at the interface between the two materials;⁵¹ the tangential components of velocity and traction were assumed to be discontinuous at the interface to allow frictionless sliding.⁵¹ A free-surface condition as described in Sambasivan et al.⁵⁰ was used at the crystal–void interfaces. Interfacial conditions were enforced using a modified ghost-fluid method.⁶⁰ The interfacial conditions and their implementation in the current framework are described in detail in previous works.^50,51

The domain boundary conditions were prescribed to mimic the experimental conditions. Only a sub-region of the overall experimental region was captured in the simulations to economize on the computational effort. The full extent of the experimental setup²⁵ (the 3 mm diameter well, which housed the PBX samples) went well beyond the domain represented in the simulations $(300×300μm2)$⁠. Because the computational domain occupies a smaller lateral extent than the experimental one, the east, west, and south boundaries of the computational domain were imposed with outlet type boundary conditions to allow for the free exit of waves at those boundaries. The north boundary, on the other hand, was given a Dirichlet velocity boundary condition to approximate the pressure pulse induced by the flyer impact. To compare to the experiments where the particle velocity of the PU polymer at the impact surface [Fig. 1(e)] was estimated through an impedance calculation using the flyer velocity, a particle velocity $(Up)$ of 2.3 km/s was imposed on the binder for the initial 4 ns of simulation, after which, the velocity was set to 0.

It is to be noted that the 1D pulse velocity boundary condition employed at the north boundary in the simulations presented here is an idealization of the actual behavior of the binder–air interface impacted by the flyer in the experiments. The assumption here is that the thickness of the flyer $(25μm)$ is large enough that the release wave reflected from the free end of the flyer does not catch up to, and weaken, the main shock in the PBX before the shock traverses ∼1/2 the length of the crystal, i.e., before t ≈ 12 ns. For the cases presented here, the locations of prospective hot spots (upper surface and internal voids) lie roughly within this distance. We also know from previous numerical tests⁴⁵ that the leading edge (with respect to the shock) of the crystal is much more susceptible to reactions at the shock-facing crystal–binder interface as compared to the lower surface. Thus, in the interest of computational tractability, it is assumed that the main shock from the impact is able to reach the regions of interest in the crystal, creating hot spots and triggering a reaction, before being weakened by the ensuing rarefaction. While a boundary condition incorporating the release wave in the binder would be physically more accurate, an application of such a condition will necessitate the simulation of flyer impact on the crystal–binder system. In the experiments, it is observed that flyer impact leads to very large deformation of the flyer including shattering, and hence, modeling the impact, fragmentation, and separation of the flyer is quite complicated in our current numerical framework. The portrayal of interaction between the flyer and binder at the top surface is a subject of current work to extend the framework.

The computational domain for simulations of the shocked HMX crystals was 300 μm² with a uniform grid spacing of 0.25 μm. The level 4 tree-based cell division of the LMR employed for the current cases resulted in the finest grid cell of size $0.0325μm$⁠. The grid sizes and levels of LMR were based on experience from previous works.³⁶ During the simulation runtime, the LMR procedure⁶¹ dynamically and adaptively refined the mesh. A typical estimate of the number of grid points with LMR was ∼8 M cells for the 2D cases. The high-fidelity simulations were run using 30 computing nodes. Each node had a 2.8-GHz Intel Xeon Broadwell 44 core and 128 GB of DDR4 memory. The simulations took about 50–70 wall clock hours or about 66 000–92 000 CPU hours to yield ∼50 ns simulation runs.

The shock experiments on deconstructed PBX consisting of a polymer-encased single HMX grain, either a single crystal or a defective polycrystalline aggregate (collectively referred to as “crystals” in this paper), have been described in detail previously,^25,65 and will be summarized briefly here.

HMX crystals were grown by evaporation on a crystallization dish from a saturated solution of HMX in spectroscopic grade acetone. Warning: HMX is a powerful explosive that should be handled only by trained personnel, and HMX solutions are flammable. Only HMX crystals ranging from ∼100–350 μm in size were selected. The crystals were embedded in fluid HTPB which was subsequently polymerized into PU using an isophorone diisocyanate crosslinker and ferric acetylacetonate catalyst. During polymerization, the HMX was allowed to drift toward the PU surface at the bottom so that all the HMX/PU samples (after being inverted) had the HMX about 10 μm below the PU surface.

CT scans of the embedded crystals were carried out on a high-resolution x-ray computed tomography instrument (Rigaku Nano3DX). Samples were imaged using an 8 keV x-ray emission from a Cu anode. The HMX crystal with its embedding polymer was mounted on a 1 mm diameter rod. The system's L0270 lens was used to achieve a 0.55 μm voxel size. A total of 800 projections was acquired over 180 degrees of rotation, and each projection took 8 s exposure time. The resulting stack of 2D CT images was processed with a 3D image segmentation software (Simpleware ScanIP) to generate a 3D view of the crystal such as what is shown in Fig. 1(a).

Three diagnostics were used: a 4-frame 20 ns camera with 2 μm spatial resolution, a 2 ns optical pyrometer, and an 8 GHz photon Doppler velocimeter with 2 ns time resolution. Thermal emission images from the camera were processed using an in-house MATLAB™ code, which converts 16-bit grayscale images to a false color palette. In some images, such as Fig. 1(f), the false-colored thermal images were superimposed on the unshocked crystal image in order to indicate where the hot spots were produced. Time-dependent superimposed images are intended for reference only, since during the shock measurements, the crystal does not retain its original structure. The optical pyrometer measured the spectral radiance, which was converted to a time-dependent history of the spatial-averaged temperature using graybody approximation. The details of the procedure can be found in previous works.^66,67

It is important to note the orientation differences between experiments and simulations. As illustrated in Fig. 1(e), the camera acquires images in the plane perpendicular to shock propagation, whereas the 2D simulations are computed on planes parallel to shock propagation. These different viewpoints necessitate processing the images as well as numerical simulation data to enable quantitative comparison of the measured and calculated quantities of interest. Since the experimentally acquired images are available at discrete time intervals, constrained by the shutter speed of the camera, optical access, spatiotemporal averaging, etc., and the numerically calculated and visualized 2D fields impose limits as well, care has to be exercised in making meaningful head-to-head comparisons between the experiments and simulations. These post-processing measures are described next.

From a library of CT scans, three crystals denoted as Crystal 1, etc., were chosen for comparison to the simulations. Figure 2 shows 3D renderings of the CT scans of the polymer-encased crystals, along with the cross sections used as simulation inputs, that were each chosen as representative of the entire crystal. Crystal 1 is highly defective, with internal voids and cracks. Crystal 2 does not have visual cracks, but it does have clusters of internal voids. Crystal 3 is nominally a defect-free single crystal.

The first line of visual and qualitative comparison with experimental data is done through contour plots of the temperature field. Figure 1(d) shows an example of a simulation-derived contour plot of the temperature field. The range of temperatures shown covers the extent reached in the crystal and binder, and in cases where initiation occurs, the hot spots and final combustion products. Temperature plots are used here primarily to track the hot spots in space and are displayed together with thermal emission images as they evolve in time, as in Figs. 1(d) and 1(f). Other contour plots used in this work are numerical Schlieren images and the final reaction product species mass fraction (Y₄) distribution. The numerical Schlieren is computed from the material density gradient, appropriately scaled to highlight density jumps at shocks and other flow features.⁶⁸ This provides a visualization of the shock dynamics and reaction front progression obtained from the simulations. Incident and reflected shockwaves and their relative strength can be discerned from Schlieren plots. The Y₄ field is useful in tracking reaction fronts and the rates at which the front progresses to consume the unreacted HMX.

The head-to-head comparison focused on two quantities measured in the experiment, the spatially-averaged temperature [Fig. 1(g)], and the spatially-dependent radiance. The spatially-averaged temperature was computed directly from the simulation. To compute the spatially-dependent radiance from simulation, the 2D temperature field at a particular time instant was first interpolated onto a uniform Cartesian mesh of size 1000 × 1000 (from the quadtree structure of the adaptively refined mesh). This 2D spatial distribution of temperature was then processed into a histogram, giving the statistical distribution of the number of grid cells (i.e., area of the domain) occupied by different temperature levels or bins. An example histogram is shown in Fig. 3, corresponding to Crystal 2 in Figs. 2(c) and 2(d). For each bin, we computed the spectral radiance using Planck's blackbody law. We then summed up the spectral radiances and integrated them over wavelength to obtain the radiance at each location. The temperatures obtained through this method are biased toward higher temperatures in the domain (both computational and experimental). This is because the measured thermal emission intensity is a strong function of temperature (⁠$radiance∝T$⁶) and is dominated by the peak hot spot temperatures in the material.

Simulation results for Crystal 1 are summarized in Fig. 4, which were performed on the 2D slice shown in Fig. 2(b). Figure 5 presents the head-to-head comparison of simulation and experiment. The binder–air interface shown in Fig. 2(b) is the site of the flyer impact. Figure 4(a) shows the Schlieren rendering of the simulation at different times. Figure 4(b) shows the computed temperature field. Note how at later times the original shape of the HMX crystal has changed due to compression and plastic deformation from the shock. The outlines in Figs. 4(b) and 4(c) indicate the dynamic boundary between HMX and polymer. Figure 4(c) shows the spatial dependence of the computed mass fraction of reaction products, denoted Y₄.

As seen in Fig. 4(a), the shock wave from the top boundary passes through the crystal, compressing the material in its wake. The voids in the crystal collapse, generating blast wave patterns seen via the Schlieren rendering. The collapse leads to high-temperature zones in the vicinity of the voids, with temperatures exceeding 6000 K in highly localized regions for short durations, resulting in chemical reactions and product formation [Fig. 4(c)]. The large excursions in temperatures are observed only in the instants following a void collapse, for example, between the snapshots shown at 9.60 and 13.71 ns in Fig. 4(b). Note that these high temperatures exist only in small regions and on short time scales during a void collapse. Figure 4(b) also shows small regions of high temperature in the upper left corner of the crystal as well as other local regions on the top surface of the crystal leading to chemical reactions. Such interfacial-initiated reactions are seen only on the surface that faces the incident shock and are not observed on any of the other surfaces of the crystal. These areas of high temperature are caused by shock energy pileup created by the HMX–binder shock impedance mismatch, as shown in a study of mechanisms of energy localization in idealized crystals by Das et al.⁴⁵ Note that the temperatures of the interfacial hot spots are lower than at the void collapse sites, as observed in simulations⁴⁵ as well as in the experimental observations of Johnson et al.²⁵ 

Figure 4(c) shows that the reactions initiated by both the void collapse and the surface hot spots lead to hot spot growth, that is, they are “critical” hot spots in the sense described by Tarver et al.²¹ However, because of the higher initial temperatures at the void collapse sites, the collapse-generated hot spots grow at a faster rate and occupy a larger area than the surface-generated hot spots. Note that the strain-rate dependent Johnson-Cook model for HMX⁶⁹ leads to the presence of shear bands which offer preferred paths for the hot spot evolution, as seen from the branched structure of the hot spots. For example, in the third column (27.87 ns) of Fig. 4, the reaction front can be seen to develop a branch; the branching is due to a high-temperature region or melt crack that forms due to shear localization. The reaction front progresses faster along the melt crack leading to the observed branch. Thereafter, the reaction front spreads along this branch as observed in the hot spot at 61.02 ns. Interestingly, the crack in the crystal [see Fig. 2(b)] does not appear to generate a hot spot. The existence of material (binder) inside the crack arguably changes its behavior compared to internal voids, which are not filled with any material for the simulations. In addition, the orientation of the crack is such that the shock impinges on it in an unfavorable direction for the formation of a hot spot. As shown in Rai et al.³² elongated voids (or cracks) that are oriented with their long side parallel to the shock front (i.e., perpendicular to the direction of the shock propagation) fail to produce hot spots; this is due to the generation of a rarefaction zone when the shock impinges on the free surface of an elongated void.

The existence of hot spots in the experiments on Crystal 1 is clearly demonstrated in the emission images in Figs. 5(a) and 5(b). The cross section used for the simulation is located toward the middle of the crystal denoted by the yellow line in Fig. 5(a), which indicates the plane of the 2D crystal slice used in the simulations. This plane contains major cracks/voids in the crystal and is also the region with bright illumination in the experimentally obtained hot spot pictures shown at 17 ns in Fig. 5. It is, however, difficult to gauge from the images if the hot spots are internal or at the surface. It can be noted that the brightness of the hot spot is highest near the center of the image at the 17 ns time instant. Furthermore, the temperatures at 107 ns [Fig. 5(b)] are low enough that no localized high-temperature hot spot is visible at that time. It is surmised, based on the simulated temperature contours in Fig. 5(c) at 115.8 ns, that most of the HMX crystal has combusted, and the temperature is seen to have fallen to values in the vicinity of 2000 K, which is the noise floor of the pyrometer.

Figure 6 shows how we compared the experimental measurements of hot spot locations in Crystal 1 to simulation. Figure 6(a) shows the regions of low density (i.e., voids or cracks) within the HMX crystal. Such regions are clearly evident in the x-ray CT scans due to large differences in x-ray scattering cross sections between HMX and void. Figure 6(a) is based on the plane used in the 2D simulation, and it shows the projection of these void regions along the line shown in Fig. 5(a). Figure 6(b) shows the predicted radiances from simulation along the same line. Figure 6(c) shows the experimentally measured image radiance along this line. It is to be noted that Fig. 6(a) indicates the defect locations in the undisturbed crystal, whereas Figs. 6(b) and 6(c) are extracted from the simulated and imaged plots, in this case, Figs. 5(c) and 5(b), respectively. Hence, in general, hot spots corresponding to defect locations may not be seen in the simulated and imaged radiance plots if the shock has not reached the void/s. Alternatively, if the data are collected at a later time when the reaction zone has branched and spread, then the radiance plots may show peaks additional to the defect locations.

Figure 6 indicates a general, but imperfect correspondence between voids in the crystal and the simulated and measured radiance from hot spots. All the plots show peaks at the middle and right dashed lines. These radiance peaks result from the collapse of the voids indicated in the defect-location plot [Fig. 6(a)]. There are some disagreements among the plots as well. In particular, the simulation in Fig. 6(b) predicts one large hot spot not predicted by the measured void distribution in Fig. 6(a), and the measurement in Fig. 6(c) does not show either of the two large hot spots predicted in Fig. 6(b). The discrepancy along the left line may indicate that the hot spots are present for only very short time durations, which makes them difficult to be captured, given the 20 ns time resolution of the imaging camera. The high radiance on the right-hand side of the image plot may also overshadow the peaks on the left. The isolated peak on the left of the simulated radiance plot likely points to the shear-driven branched reaction front discussed in Fig. 4. This mismatch of the reaction zone between the computed and imaged radiance may be due to the 3D effects (a higher gradient of temperature elsewhere in the 3D crystal causing the reaction to branch into a different plane), or shortcomings of the Johnson-Cook strength model used in the calculations. Further work is needed to examine the causes for such discrepancies between the experiments and simulations.

Figure 7 compares the predicted and measured time-dependent temperature profiles for Crystal 1. Again, there is a general but imperfect correspondence between the experiment and simulation. The experimental temperatures are hotter, and the temperature decay is slower. The lower temperatures in simulation might indicate that the model for HMX mechanics is imperfect, or it might be due to comparing a 2D simulation to a 3D experiment. The experiment might be seeing hot spots not present in the 2D slice used in the simulation. Alternatively, the small regions and short time durations occupied by the hot spots in the simulations, may not register strongly enough in the numerically computed radiance from the simulated hot spot fields in the case of Crystal 1. The highly transient nature of the localized high-temperature regions in the simulations can be seen in Fig. 4(d), where the histogram of the temperatures in the domain is shown at each instant of time featured in Fig. 4(a). The highest temperatures in the simulations (approaching 6000 K) are seen to be present for a very short duration (between t = 9.6 ns and t = 13.71 ns) and occupy small areas of the crystal confined to the vicinity of the collapsed voids.

The analysis in this subsection parallels the analysis for Crystal 1 above. Crystal 2 and its 2D rendering are shown in Figs. 2(c) and 2(d). Figure 8 shows simulation results in the form of Schlieren (a), temperature profile (b), and reaction products (c). The most significant difference between Crystal 1 and Crystal 2 is that Crystal 2 has multiple voids clustered together and no large crack. Comparing the simulated reaction temperatures and product distributions in Figs. 4 and 8 shows that Crystal 2 produces higher temperatures and faster reactions.

Experiment and simulation for Crystal 2 are compared in Figs. 9–11. Figure 9 shows Crystal 2 and its thermal images during shock, along with the corresponding temperature prediction. The yellow line in Fig. 9(a) is the line upon which we projected the void locations and the predicted and measured radiances. Figure 10 shows the defect locations, and the predicted and measured radiances along the yellow line in Fig. 9(a). As with Crystal 1, there is again an overall qualitative but imperfect correspondence between the experiment and simulation. In this case, we find more instances of agreement between the locations of defects and hot spots as compared to Crystal 1. This applies not only to the experimental and simulated radiances in Fig. 10, but also to the temperature vs time data shown in Fig. 11. In Fig. 10, as seen on the 1st, 3rd, and 4th dotted lines from the left, experiments and simulations indicate hot spots that correspond to the defect locations in the section shown in Fig. 9(a). The simulated cluster of peaks around the first dotted line (centered around x = 100 μm) appears to be due to the combination of void collapse and surface reactions that can be observed to have been initiated at the top surface [see Fig. 9(c)]. The surface reaction produces lower temperatures and may contribute to the relatively shorter peaks in that cluster. The 3rd and 4th peaks in the numerical results [Fig. 10(b)] are evidently due to void collapse and show strong single peaks. The main difference is that the numerical predictions show strong concentrated peaks at those locations, while the experimental data shows a broader, less strong distribution of radiance. As pointed out earlier, this difference may be due to the inability of the 20 ns camera snapshots to catch short-lived intense bursts of radiance and the spatiotemporal averaging inherent in the experimental data processing of the emitted radiance. While, overall, the computed radiance signatures at the 1st, 3rd, and 4th peaks are in good agreement with the experimental observations in Fig. 10(c), there is a notable discrepancy with regard to the location identified by the second dashed line (centered around x = 125 μm) in that the simulation in Fig. 10 fails to predict a radiance peak at that location while a strong radiance peak is observed in the experiment. This discrepancy can be attributed to the different ways in which the numerical and experimental data on thermal emission intensity is collected. At the instant for which the radiance data in Fig. 10 is plotted, Fig. 9(c) shows that the computed shock clearly has not reached and collapsed the void located at x ≈ 125 μm in the crystal. Since the experimental data are collected at intervals of 20 ns and spatiotemporally averaged, a strong peak due to the collapse of the void located at x ≈ 125 μm can be seen in the experimental radiance data in Fig. 10(c). These challenges in coordinating and unifying the computational and experimental data collection and processing techniques will need to be overcome to improve the precision of head-to-head comparisons of the shock response and hot spot evolution.

Figure 11 compares measured and computed time-dependent temperature profiles for the case of Crystal 2. The computed temperatures are in better agreement with those obtained experimentally for Crystal 2 as compared to Crystal 1. This better agreement of the peak temperatures in the hot spots is aligned with that of the radiance information for Crystal 2 when contrasted with Crystal 1 (compare Figs. 6 and 10). However, the decay of high temperatures is still slower in the experiment than in the simulation. This slower decay of temperature may once again indicate that thermal transport in 3D needs to be captured in simulations to better match with experimental hot spot growth behavior. Overall, Case 2 has shown that some features of the hot spot phenomena in the crystal are captured to a fair degree by the 2D simulations performed in the paper. This case also highlights the aspects of experiments (data gathering at higher temporal and spatial resolutions) and computations (3D simulation and improved material models) that need to be advanced to obtain a closer agreement in future work.

The nominally defect-free Crystal 3 is shown in Figs. 2(e) and 2(f). The defect-hot spot correlation for this crystal is not of significance unlike those with the defect-laden Crystal 1 and Crystal 2. In fact, the simulation of Crystal 3 shows minimal hot spots and chemical reactivity [see Fig. 13(a)]. In the experiment as well, Crystal 3 did not show the existence of hot spots in imaging or pyrometry, in the duration of ∼100 ns from impact as recorded for this study.

Note that Crystals 1 and 2 did show surface hot spots that continued to grow into critical hot spots [Figs. 4(b) and 8(b)], indicating that such a mechanism can operate even for nominally defect-free crystals. Das et al.⁴⁵ also showed that inter-/intra-crystal voids are not necessary for shock initiation; shock focusing at the crystal–binder interface can also lead to initiation. To look more closely at the behavior of “clean” (i.e., defect-free) crystals in comparison with defect-laden ones, the CT images of Crystals 1 and 2 were image processed in MATLAB™ to obtain the images of crystal sections without any internal defects (voids, cracks), i.e., by “cleaning up” the crystals, to obtain internal defect-free crystals (called Crystals 1* and 2*, respectively) as shown in Figs. 12(b) and 12(c).

Simulations of shock loading at the identical particle velocity of 2.3 km/s were performed on Crystals 1* and 2* as well. The temperature contours from the shock simulations are shown in Fig. 13, where the top row corresponds to Crystal 3, the middle row to Crystal 1*, and the bottom row to Crystal 2*. Crystals 1* and 2* show prominent surface reactions. Between the first two snapshots in the middle and bottom rows in Fig. 13, localized hot spots are seen to develop on the surface, with temperatures momentarily reaching values approaching 4000 K. These critical hot spots produce sustained reactions, with growing hot spots displaying temperatures ∼2000 K. The hot spots sustain the surface reaction and consume appreciable areas of the crystal in the timeframe of 70 ns, as shown in Figs. 13(b) and 13(c).

The hot spot growth in Crystals 1* and 2* is compared with their original counterparts Crystals 1 and 2, respectively, in Fig. 14 for longer simulation times. Figures 14(a) and 14(b) correspond to Crystals 1 and 1*, respectively, with the progression of hot spots shown from left to right at time instants t = 14, 37, 77, and 115 ns. Figures 14(c) and 14(d) correspond to Crystal 2 and 2*, respectively, with the progression of hot spots shown from left to right at time instants t = 9, 19, 37, and 61 ns. It is seen that the hot spots from the void collapse events within the crystal initiate at higher temperatures, and the burn fronts due to such events proceed to consume the crystals faster than the surface-initiated hot spots. The two types of reaction fronts also spread into and merge across regions of the crystals. However, the occurrence and sustenance of surface hot spots are not altered by the void collapse events in the crystals, i.e., the surface-generated hot spot mechanism is independent of other (void collapse or reaction wave) hot spot mechanisms that occur within the crystal. This indicates that the surface reaction can occur even in nominally “clean” crystals under suitably strong shocks. In addition, by comparing the evolution of hot spots between Crystal 1 and 2, it is seen that the reaction wave mechanism that arises in the latter crystal due to the cluster of voids produces much faster growing hot spots resulting in the rapid completion of reactions in the crystal. Therefore, the three different mechanisms of hot spot formation in the cases studied produce distinctly different reaction progress rates, as pictured in Fig. 14.

Head-to-head comparisons were presented for the deconstructed PBX consisting of HMX crystals embedded in polyurethane (experiment) and Estane binder (simulation). The temperature fields obtained from the simulations and experiments as well as the points of initiation and the subsequent evolution of hot spots were compared for three selected crystals, bringing to light different types of hot spot initiation mechanisms: (1) strong, high-temperature hot spots formed due to pore collapse; (2) lower temperature hot spots initiated at the surface of the crystal due to corners and asperities; and (3) rapidly propagating high-temperature reaction waves that propagate faster than the speed of the incident shock, within the crystal, leading to a fast consumption of the energetic material.

The overall agreement between experiment and simulation was qualitatively correct but imperfect. The predicted hot spots appeared where the CT scans indicated defects were present, but the predictions also showed hot spots not observed by the experiment. This is likely not so much a failure of the model but a limitation on the part of the experiments to capture short-lived hot spots. The time-dependent temperature predictions were similar to or a bit lower than those measured by the experiment. The prediction that the greatest reactivity would occur in the order Crystal 2 > Crystal 1 > Crystal 3 was supported by the experiment as well.

The differences between experiments and simulations are attributed to the comparison of 2D simulations with experiments on real 3D crystals, and possibly, the simplified homogenized model for HMX mechanics used in the simulations. These are clearly two areas that require improvements, i.e., further refinement of the mechanical models for HMX and the use of 3D simulations. In addition, approximations were made in applying the boundary conditions in the numerical simulations, most notably in the way in which the flyer impact on the top surface of the binder was represented. More realistic boundary conditions need to be implemented to simulate the impact of the flyer on the crystal–binder system. The tandem simulations and experimental results presented in this paper point to areas of improvement on both fronts, which will make it possible to distinguish between different models for the initiation of HMX-based PBXs, leading to the creation of reliable validated models for shock initiation.

See the supplementary material for the material and reaction models for HMX and Estane used for the simulations. Simulation movies of the evolution of Schlieren, temperature, and reaction products, as well as temperature histograms for the crystals, are also included.

The research described in this study is based on the work at the University of Illinois Urbana-Champaign and the University of Iowa supported by AFOSR-MURI Grant No. FA9550-19-1-0318, program manager: Dr. Martin Schmidt. The nano-CT characterization was carried out in Beckman Institute at the University of Illinois at Urbana-Champaign.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding authors upon reasonable request.