Modeling pore collapse and chemical reactions in shock-loaded HMX crystals

The localization of deformation in shock-loaded crystals of high explosive material leads to the formation of hot spots, which, if hot enough, initiate chemical reactions. The collapse of microscopic pores contained within a crystal is one such process that localizes energy and generates hot spots. Given the difficulty of resolving the details of pore collapse in shock compression experiments, it is useful to study the problem using direct numerical simulation. In this work, we focus on simulating the shock-induced closure of a single pore in crystalline β-HMX using a multiphysics finite element code. To address coupled thermal-mechanical-chemical responses, the model incorporates a crystal-mechanics-based description of thermoelasto-viscoplasticity, the crystal melting behavior, and transformation kinetics for a single-step decomposition reaction. The model is applied to stress wave amplitudes of up to 11 GPa to study the details of pore collapse, energy localization, and the early stages of reaction initiation.


Introduction
The shock initiation of high explosive (HE) crystals, and aggregates thereof, is sensitive to localized regions of elevated temperature (hot spots) that are generated during shock deformation [1]. Therefore, it is desirable to develop a fundamental understanding of processes that localize energy in HE crystals under shock wave loading. Localization in single crystals of HE material stems from both the interaction of the shock wave with microscopic defects contained within the crystal (material heterogeneity) and the intrinsic heterogeneity of inelastic deformation mechanisms (e.g., crystallographic slip). Since this information is exceedingly difficult to extract from shock compression experiments, we set out to study processes involved in the shock initiation of HE crystals using direct numerical simulation.
A continuum model is developed to simulate the thermal, mechanical, and chemical responses of β-HMX (C 4 H 8 N 8 O 8 , monoclinic lattice) under shock wave loading. As a model problem, we consider the shock-induced closure of a single pore in β-HMX crystal. The problem of pore collapse has been considered in previous work [2][3][4][5]; however, these efforts have focused on HE crystals that are modeled as isotropic and/or inert. The progress of this work is the incorporation of a crystal-mechanics-based description of anisotropic, thermoelastic-viscoplastic material response as well as chemical reactions that are initiated by the deformation and heating of the material. The continuum model is exercised for stress wave amplitudes of up to 11 GPa to study the early stages of reaction initiation. In the following, crystallographic directions and planes are given in P 2 1 /c notation.

Reactive crystal model
The material model is an extension of an inert thermoelastic-viscoplastic crystal model [4], which now makes use of a reactive equation-of-state (EOS). The crystal model handles the coupling between the distortional and volumetric aspects of the deformation, utilizing elements of traditional crystal models, as described in previous work [4]. Briefly, the deformation gradient within the solid crystalline phase is written in terms of a multiplicative decomposition for finite deformations (F = V RF p ), wherein F p describes the plastic shearing of the lattice, R is a rigid body rotation, and V is a thermoelastic stretch. In this framework, plastic deformation occurs by the motion of dislocations on slip systems (i.e., on specific crystallographic planes and in specific directions). The crystal model has been calibrated to single-crystal shock wave profiles that were measured in gas gun experiments [6].
The volumetric responses and chemistry are handled by the thermochemical code, Cheetah [7]. In these calculations, the material points are considered as multi-phase volumes with properties calculated and dynamically tabulated according to the species concentrations. Given the composition, density, and internal energy of a material point, the pressure and temperature are calculated according to the equations-of-state of the individual species. The species concentrations are updated according to the prescribed reaction path and kinetics. The following single-step reaction is considered, The global decomposition rate is modeled using an Arrhenius rate law that was determined by Henson et al. [8]. The solid phase (β) is described using a Murnaghan EOS, whereas each of the product phase gas species are described using a Buckingham exponential-6 potential model. The product gases are collectively maintained in chemical equilibrium with trace amounts of other species, such as CO 2 and H 2 O, which may arise during reaction.
To account for loss of strength due to melting, a compression-dependent melt curve has been generated to calculate internal energy at the onset of melting, as predicted by the Lindemann law. As the material progresses through the melting phase transformation, a mixture rule is used to estimate the deviatoric responses of the solid/liquid mixture, wherein the liquid phase is modeled as a viscous Newtonian fluid. Once the internal energy exceeds the melt energy by an amount equal to the latent heat of melting, the material no longer possesses any strength. The strength is reduced in a similar manner for material points that have undergone reaction and contain solid/gas mixtures.

Simulation of pore collapse under shock wave loading
The computational domain is rendered by orienting the lattice within a rectangular block of crystal and inserting a single air-filled pore near the center of the block. To reduce the cost and complexity of these exploratory calculations, a two-dimensional problem is considered by invoking the plane strain condition. Therefore, these simulations address the closure of cylindrical pore. Planar shock waves are generated by prescribing the velocity of the left-hand boundary of the sample. The top and bottom surfaces of the crystal are periodic and the right-hand boundary is restrained by a symmetry plane. The shock simulations are performed using the arbitrary Lagrangian-Eulerian finite element code, ALE3D. In this code, the material and the mesh are permitted to undergo independent motions, with advection algorithms accounting for the flux of material among the computational zones (elements). This allows for simulation of the large deformations involved in pore collapse, while avoiding the excessive mesh distortions that cripple pure Lagrangian calculations. The coupling among thermal, mechanical, and chemical responses is achieved using operator splitting. The calculations have been performed under adiabatic conditions to obtain a baseline response. As a reference case, we considered a pore diameter (d) of 1 µm, a peak particle velocity of 1 km/s (which generates a stress of 9.4 GPa), and a crystal orientation that is obtained as follows: the normal of the (111) plane is aligned with the shock direction (x-axis); the crystal is then rotated about the shock direction until the [110] direction forms an angle of θ with the y-axis. To allow for comparison with previous work [4], we selected θ = −38.8 • . For a pore diameter of 1 µm, convergence studies indicate that a mesh length of 8 nm was suitable for our purposes. To mesh the entire crystal sample, about 10 million elements were required. These calculations were run in parallel on 512 processors for a wall clock time of about 72 hours.

Results
When the shock wave reaches the pore, a release wave is emitted from the crystal-air interface. This release wave is followed by a secondary shock that is generated when the pore is fully closed. The simulations are run up until the secondary shock begins to interact with the boundaries. This allows for a post-collapse simulation time of approximately 2 ns. The post-collapse temperature field that is calculated for the reference case is shown in figure 1. The temperature field illustrates energy localization within the crystal. Near the center of the sample, where the pore has collapsed, is an intense hot spot that has been brought about by mechanical deformation and chemical reaction (discussed further below). The narrow regions of high temperature that emanate from the pore collapse region (and span across nearly the entire sample) are shear bands. The shear banding occurs due to localized deformation on specific crystallographic planes. When the shock wave is propagated normal to the (111) plane, the (011) and (102) planes are well-oriented for slip, as shown in figure 1. The localization of plastic deformation on these planes drives up the temperature and softens the material. This deformation eventually becomes unstable and material in the bands is heated beyond its melting point. The shear bands, which are filled with molten viscous fluid, are akin to "melt cracks" that propagate across the crystal sample. These features are also seen in full 3D calculations [4], which indicates the melt cracks are not merely an artifact of the plane strain condition.
The propensity of low symmetry HE crystals to shear band under shock wave loading has been documented in previous experimental work [9,10]. For example, crystals of RDX (C 3 H 6 N 6 O 6 , orthorhombic lattice) that were shocked up to 13 GPa exhibited shear bands on the surfaces of recovered crystals. Furthermore, beaded-up volumes on the surfaces are thought to be molten material that was squeezed out of the shear bands and re-solidified. These experimental observations are encouraging to us, as this is precisely what is predicted by the model: shear bands that are filled with molten material. Shear banding is also found in molecular dynamics simulations of shock wave loading of α-HMX [11], whereas shear cracking is thought to be necessary to explain observed wave profiles in sterically hindered orientations of PETN crystal [12].
An enlarged view of the pore closure process is shown in figure 2 over a period of about 0.5 ns. In figure 2: (a) small regions of localized deformation develop (due to stress concentrations) as the wave front reaches the pore; (b) melt-cracking in the vicinity of the pore contributes to the formation of two fluidized jets; (c) the jets impinge on the right-hand side of the pore wall and are re-directed towards each other; (d) the jets impact one another and create a backward-moving jet (from right-to-left), which collides with a small tertiary jet that forms at the left-hand side; (e) the collision of these jets generates a pressure of about 30 GPa, which drives up the temperature enough to react a small amount of material; (f) the pore is fully closed and the reaction zone is beginning to grow in size. Hence, the shock-induced closure of a micron-sized pore in this low symmetry crystal involves a complicated sequence of events, with fluidized jets traveling in multiple directions.   The only region of the crystal that undergoes chemical reaction on this time scale is the small hot spot located near the center of the collapsed pore. The reactivity of the sample is quantified by calculating the relative mass of the product phase, φ p , which is defined as the total mass of products divided by the mass of crystal that would have fit inside the initial pore. Time histories of φ p are plotted in figure 3 for three shock loading intensities. For a peak longitudinal stress of 6.5 GPa, there is very little reaction on this time scale. At 9.4 and 10.7 GPa, nearly 10 percent of the "pore mass" is reacted. It is tempting to speculate that the initiation threshold of this configuration is around 9.4 GPa; however, longer simulation times and the inclusion of heat transfer are required to predict the initiation threshold with confidence. In experiments, the HMX-based composite PBX 9501 requires roughly 3.5 µs at a shock stress of 3 GPa to achieve ignition [13]. This indicates that reactions can be initiated at stresses less than 9.4 GPa for time scales longer than a few ns, or through mechanisms distinct from pore collapse, such as interactions between grains.

Summary and conclusions
A continuum model has been developed to study pore collapse and energy localization in shockloaded crystals of β-HMX. The continuum model makes use of a constitutive model that accounts for the anisotropic, thermoelastic-viscoplastic response of the crystal as well as temperaturedriven decomposition reactions. The responses of a single configuration (d = 1 µm, (111) impact plane) were simulated for stress wave amplitudes of up to 11 GPa and a post-collapse run time of about 2 ns.
The simulations indicate that shear banding is an important mode of deformation when these (defective) crystals are subjected to shock wave loading. The shear bands nucleate during the initial stages of pore collapse and continue to grow away from the pore after it is fully collapsed. Although the shear bands reside at elevated temperature (approximately 800-1200 K), only very small degrees of reaction are observed within the shear bands on this time scale. Full chemical reaction does occur at the hot spot located in the pore-collapsed region. The model suggests a shock initiation threshold of about 9.4 GPa for the configuration under consideration, however, longer simulation times are necessary to further study the initiation and propagation of reactions in these crystals.