Multiscale Modeling of Shock Wave Localization in Porous Energetic Material

The most common approach for studying shock response of heterogeneous materials is through a continuum mechanics representation of material properties. Continuum simulations usually follow one of two approaches in the literature: a hydrodynamic model introduced by Mader mader1965 , or a viscoplastic model introduced by Carroll and Holt carroll1972 , which was later extended by Khasainov et al. khasainov1996 , Butler et al. kang1992 , and Frey frey1972 . In summary of the different models, the type of pore collapse depends on the ratio of inertial to viscous forces, i.e. Reynold’s number. For a cylindrical pore geometry, the Reynold’s number is given as,

where $a_{0}$ is the initial pore radius, $\rho$ and $P$ are the density and pressure following the shock wave, and $\mu$ is the shock viscosity. The shock viscosity is estimated from experimental measurements hambir2001 or theoretical calculations chou1993 , and is usually orders of magnitude lower than normal values. At low Reynold’s number ($Re\ll 1$), the viscoplastic models assume radial symmetry and incompressible plastic flow behind the shock wave, whereas at high Reynold’s number ($Re\gg 1$) the compressible flow equations are solved using a hydrocode, and material strength is often neglected. Lines of constant $Re=1$ are drawn on a shock viscosity-pore radius plot for HNS, assuming different values for the shock pressure and a range of shock viscosities from Chou et al. chou1993 (see Fig. 1). Insets in Figure 1 schematically show these differing pore collapse mechanisms. Unfortunately, the viscoplastic to hydrodynamic transition is ambiguous across a wide range of pore radii (0.1 to 30 $\mu$m) and both mechanisms may be relevant to HNS initiation.

Continuum simulations that naturally capture both viscoplastic and hydrodynamic behavior (including material jetting, viscoplastic heating, and shear banding) require a full stress tensor model. Only a few such models have been developed to model pore collapse in explosives; one of the most relevant to the current work is the $\beta$-HMX crystal plasticity model of Austin et al. austin2015 . One key result from that work is rate-independent strength models are incapable of reproducing shear banding and viscoplastic collapse. Hence, a simple strength model with plastic strain dependence alone cannot give the desired material behavior.

In this work we employ the strain-rate dependent, Steinberg-Guinan-Lund (SGL) viscoplastic strength model steinberg1989 for HNS. For simplicity, a reduced form of the full SGL model was chosen for tuning the yield strength to the strain rate. In summary of this constitutive model for HNS, the stress tensor is decomposed into spherical and deviatoric terms,

where the mean pressure, $\bar{P}$, is given by the thermodynamically complete Mie-Grüneisen equation of state (EOS) from Kittell and Yarrington kittell2016 , and the von Mises yield criteria is assumed to limit the magnitude of the deviatoric stresses,

where $S=|\bm{S}|=\sqrt{S_{ij}S_{ij}}$ (summation implied over repeated indices throughout) and $Y$ is the yield strength given by,

where $Y_{T}$ and $Y_{A}$ are the thermal and athermal components, respectively. The value for $Y_{A}$ is assumed constant, while $Y_{T}$ contains the strain rate-dependence and is given by the implicit relation,

where $C_{1}$, $C_{2}$, $U_{K}$, and $Y_{P}$ are material constants and $\dot{\varepsilon}_{p}$ is the plastic strain rate calculated from the rate of deformation tensor via $\dot{\varepsilon}_{p}=\sqrt{\tfrac{2}{3}{\dot{e}^{p}}_{ij}{\dot{e}^{p}}_{ij}}$. In addition, the total rate of deformation tensor is decomposed into an elastic ($e$) and plastic ($p$) component given by Hooke’s law and the associated flow rule,

and

where $G_{0}$ is the shear modulus calculated from the slope of the Hugoniot elastic limit, $\overset{o}{\bm{S}}$ is the Jaumann corotational derivative, and $\lambda=|{\bm{\dot{e}}^{p}}|/|{\bm{S}}|$ is a scalar used to normalize Eq. 7. Finally, the SGL model assumes a melt curve of the form,

where $\gamma_{0}$ is the Grüneisen parameter. When temperatures are found in excess of $T_{m}$, the yield strength is set to zero ($Y=0$).

The SGL model in Eqs. 4 through 8 was implemented in CTH mcglaun1990 , a shock physics hydrocode developed by Sandia National Laboratories. CTH is used to model multidimensional, multimaterial, large deformation shock wave physics, and employs a fixed Eulerian mesh with Lagrangian and remap solution steps.

In a similar fashion to the distinction between hydrodynamic model forms, molecular dynamics (MD) simulations require the selection of a material model in the form of an interatomic potential (IAP) which is a short-ranged description of atomic energies and forces. This is the leading approximation in these simulations. A wide range of different types of IAP have been developed over the last few decadesDawEAM84 ; BaskesMEAM92 ; RappeUFF ; Tersoff88 ; COMB07 . As a general trend the MD community has been moving toward more accurate and more computationally intensive potentialsPlimptonThompson12 ; Plimpton1995 ; Rappe1991 .

A subset of these IAP are known as bond order potentialsREBO02 ; BrennerIAP which dynamically calculate per-atom energies and forces as a function of the bonding environment around each atom, allowing for reactions to progress naturally rather than an ad hoc approach to chemistryBrennerAB ; HerringAB used with classical potentials. The most commonly used of these reactive MD potentials is ReaxFF, which has parameterized force fields for proteins, ceramics, metal-oxides and organic solidsSenftle2016 ; vanDuin2001 . The ReaxFF potential specific for energetic materials has gone through several adaptations since its original implementation by Strachan et al Strachan2003 . Notably, the original nitramines potential was reparameterized by merging training data with the combustion branch of ReaxFF Chenoweth2008 as well as including full reaction paths of common energetics (HMX, RDX, PETN) to form the potential reported by Wood et al Wood2014 . Liu et. al.Liu2011 added a long range, low gradient attractive term to the ReaxFF total energy functional in order to correct for the underestimation of the ambient density. These density corrected reactive potentials are denoted as ReaxFF-lg for the additional low gradient term (Eq. 9 below). This model form is employed for the present study on hexanitrostilbene (HNS) and the force field file needed to run these simulations is included as Supplemental Material.

Using the force field from Shan et. al.RayHNSReax , which is a derivative of the combustion and energetics merger, corrections to the low gradient term were needed in order to accurately predict the ambient density of HNS. This was done by adjusting the constant term, $d$, term in equation 9 from its initial value of 1.0 to 0.9125 while holding the other constant term, $C_{ij}$, at 0.55. This value was optimized by running MD simulations coupled to a thermal reservoir at 300 K and a pressure reservoir at 1 atm and comparing the predicted density to the target experimental density of $1.744\frac{g}{cm^{3}}$.

To confirm that these changes to the reactive potential do not significantly alter the behavior of HNS, we compute the $U_{s}-U_{p}$ shock Hugoniot, a property that is central to results presented here, using a variety of MD methods. In Figure 2, the DFT-MD data from Wixom et. al.WixomMattsonSAND is plotted as the solid black curve along with two equilibrium methods, Multi-Scale Shock Technique (MSST)ReedPRL03 , and Constant Stress Hugoniotstat (NPHUG) RavelPRB2004 . Also included in Figure 2 are direct measurements (NEMD) of the shock velocity from a [010] directed single crystal shock experiment using the simulation details outlined in Shan et. al.ShanPRB16 . Each of these methods results in a good agreement with the DFT reference data in the range of piston velocities ($0.5-1.25\leavevmode\nobreak\ km/s$) of interest here. Details about these two methods are contained in the Supplementary Material.

With a properly trained reactive force field in place, we now turn our attention toward the dynamics of shock induced pore collapse. It is worth noting that the extra cost associated with running the MD simulations with the reactive versus a non-reactive IAP is necessary in order to capture shock induced chemistry and the natural evolution of the hot-spot in the energetic materialWood2015 .

In order to properly calibrate the SGL strength model, both codes must share an observable quantity that represents the shock response of the porous material. This quantity should be sensitive enough to capture the characteristics of a hydrodynamic versus viscoplastic shock response within the limited time and length scales accessible to MD. Furthermore, the training metric should directly or indirectly exercise many, if not all, of the free parameters in Eqs. 4 through 8 in order to ensure a uniqueness in the fitted solution. To satisfy these constraints, we have designed the scale independent simulation geometry shown as an inset to Figure 3, and will use the pore collapse rate as the shared metric between either simulation code. In this simulation setup, the material is impacted at the left surface with a fully supported piston of variable velocity. For the MD simulation of this geometry, the shock is always directed along the [010] crystallographic direction in the (100) plane. More details about the initial conditions for either code can be found in the Supplemental Material. To simplify the analysis for every pore diameter and piston velocity pairing, we define a characteristic time ($\tau=D/U_{s}$), where $D$ is the original pore diameter and $U_{s}$ is the shock velocity for the given piston velocity, which represents the time for the shock to travel a single pore diameter.

The slopes of normalized pore area ($\hat{A}=\frac{A(t)}{A(t=0)}$) versus $t/\tau$ are collected from the NEMD runs as the set of points in Figure 3. While the pore size does have a noticeable effect on the collapse rate, the main driver behind the transition from purely viscoplastic $(\lim_{U_{p}\to 0}\frac{\partial\hat{A}}{\partial\tau})$ to hydrodynamic $(\lim_{U_{p}\to\infty}\frac{\partial\hat{A}}{\partial\tau})$ is the piston velocity that drives the shock wave. This can be seen by comparing the range of pore collapse rates at a given piston velocity to the range of values for a given pore size.

To generate these data from the continuum code, a mesh resolution was achieved using 400 zones across the pore diameter, and the domain was held constant at 1600 by 1200 cells; symmetry conditions were imposed on the impact wall and periodic boundary conditions in the lateral direction to mimic the simulation geometry used in the atomistic code.

The fitted SGL model parameters are summarized in Table 1, and were found by matching the training metric shown in Figure 3. The use of a Latin hypercube sampling (LHS) algorithm mckay1979 and $>$5,000 simulation runs were required to fine-tune the model fit; however, further improvement could be made to strength correct the EOS. Despite this, the fitted parameter values are physically realistic; for example, the value of the drag coefficient $C_{2}$ is similar to the value obtained from experimental measurements of void collapse in polymethylmethacrylate (PMMA) hambir2001 , and the Peierls stress for HNS is greater than the yield stress of 140 MPa, but less than the shear stress of 5686 MPa.

One of the main advantages of training the strain rate dependent (SRD) SGL model to reproduce the MD shock response of HNS is the larger range of defect sizes and shock strengths that can be accurately simulated with the continuum shock physics approach. To demonstrate this, Figure 4 collects several thousand individual CTH simulations that span a wide range of shock pressures and pore sizes of the same simulation cell geometry used during training. The color axis in this figure is the scaled pore collapse rate that was defined in the discussion of Figure 3. Additionally, the measured pore size distribution of HNSDammHNS is shown on a common horizontal axis having a mean and standard deviation of $0.9\pm 7\mu m$ that matches experimentally observed microstructures. As shown in Figure 4, the primary factor controlling the transition from viscoplastic (blue, region A) to hydrodynamic collapse (red, region B) is the input shock pressure, i.e. moving from region A to B. However, there does exist a size effect that is most apparent at lower shock pressures, i.e., moving from region A to region C.

It is also possible to use the size transition between region A and C in Figure 4 as a criterion to estimate the shock viscosity. Solid black lines of constant $Re=1$ for Newtonian fluids with viscosity, $\mu$, are plotted in shock pressure-pore radius space from a manipulation of Eq. 1. In contrast to the wide range of viscosities shown in Figure 1, we predict a much smaller range of shock viscosities to define the viscoplastic-hydrodynamic transition, on the order of $1Pa\cdot s$. In turn, these new estimates of the shock viscosity, in conjunction with the data in Figure 1, provide a critical pore size the separates viscoplastic or hydrodynamic style of pore collapse. From the calibrated SGL model, we find that pore sizes in the range 0.1–0.5 $\mu$m define this transition region. In the next section, we will explore these limiting cases of shock response with both codes by comparing qualitative and quantitative measures that were not used as training points for the SGL strength model.

To properly test the accuracy of the strain rate dependent SGL model, another set of metrics common to both codes is needed that are not directly used as training. In the previous sections we focused on defining the characteristics of viscoplasic and hydrodynamic styles of pore collapse simply through the rate of collapse. In this section we detail these differing mechanisms by analyzing the temperature and strain fields around the collapsed pores. The aim here is to gauge the ability of CTH to match the large-scale MD prediction of the same simulation geometry. To exemplify the improvements to CTH, we show the MD results for experimentally relevant pore sizes (0.1 $\mu$m) alongside the strain rate independent (Hydro) and dependent (SGL) forms of the HNS strength model. Of course, we cannot expect that the agreement between MD and CTH be exact, but rather we aim to capture the main features of the strain and temperature fields.

Where appropriate, the use of CTH to model the shock response of HNS is strongly advantageous because the time to solution in CTH is many orders of magnitude faster than MD. For example, using the geometry shown in Figure 3, a MD simulation cell with a 0.1 $\mu$m pore contains 20.6 million atoms and, for the slowest piston velocities, requires 200ps elapsed time before the shock reaches the free surface. A single one of these detailed MD simulations requires a significant allocation of computing resources often unavailable or unfeasible given the amount of compute time needed. Utilizing the Intel Knights Landing hardware partition on the Trinity machine at Los Alamos National Lab, these MD simulations required approximately $25\cdot 10^{6}$ cpu$\cdot$hrs to complete. Further computational details and timing data for these large ReaxFF MD runs can be found in the Supplemental Material. In contrast, the same simulation can be run with the calibrated SGL strength model through CTH in only 20 cpu$\cdot$hrs, a reduction by a factor of $1.25\cdot 10^{6}$ in time to solution.

Figure 5 collects the local temperature and plastic strain for a shock that was generated with a piston moving at 0.75 km/s. Each panel in this figure shows the simulation frame just before ejecta impact. For both quantities, the results are shown for both CTH strength models and the large-scale MD simulation. Inspecting the shape of the ejecta that is formed, it is clear that the strain rate independent strength model, Fig. 5 A and B, predicts a strong fluid-like jet that originates from the centerline of the pore. In contrast, the strain rate dependent CTH simulation, Fig. 5 C and D, and the MD prediction, Fig. 5 E and F, show that the ejecta is formed equally at two axially offset locations of the pore surface that flow toward the centerline of the pore. This phenomenon is best exemplified by the regions of highest plastic strain which in turn have the highest local temperatures in the MD simulation. A detailed look at the material flow around the pore surface is given in the Supplemental Figure S2. This is an important distinction between the strength models in CTH because the purely hydrodynamic pore collapse is seen for all pore sizes and shock strengths for the strain-rate independent form. Meanwhile the SGL model is capable of capturing both the strong jet formation in the strong shock limit and smoothly transition toward a radially symmetric pore collapse in the weak shock limit, giving it a much larger range of applicability of the space shown in Figure 1.

The second outside metric that we will compare between the atomistic and continuum methods is aimed at characterizing hot-spots that result from viscoplastic pore collapse. Specifically, we aim to provide a physical understanding of how much heat is generated from the ejecta impact versus the shear banding and other plastic material flow around the pore. In the limiting case of a purely viscoplastic pore closure, the ejecta formation will be suppressed due to irreversible plastic deformation around the pore. We have confirmed this behavior for the weakest piston velocity impacts with the MD data we have generated here.

Figure 6 shows the distribution of temperatures that are present in the MD simulation of a $U_{p}=0.75$ km/s shock at various stages during its compression. Each data series has local temperatures averaged in $1nm^{2}$ square pixels in the viewing plane shown in Figure 5 which is then collected as a histogram with bin width of 8.5 K. Only the material that has been compressed by the leading shock wave is included here. Therefore, there is no peak in these histograms corresponding to $T=300K$ of the un-shocked material. The initial heating of the sample is caused by the shock front compressing and adiabatically heating the material. This temperature distribution is shown in purple and labeled as ’Pore Collapse Begins’ in Figure 6. For reference, the red region in 6 labeled ’Pore Collapse Complete’ corresponds to panel E) in Figure 5. At this time during the simulation, there is extra heat generated from the viscous flow of HNS which comprises a very large area of the sample with temperatures ranging from 500-1500K. Each data series after the pore has fully collapsed incorporates multiple heating mechanisms, namely the collisions of ejected molecules and energy release due to chemical reactions. To show the importance of these additional heat sources, the temperature distribution for 5ps and 15ps after the pore collapse has completed is shown in blue and green, respectively, in Figure 6. Relative to the viscous heating regions, these additional sources represent a small mass fraction of simulated material, but are placed at much higher temperatures 1500K-4000K. However, this small volume of rapidly reacting material contributes strongly to the growth of the hot-spot as exemplified by the bolstered peak near 3500K after 15ps post pore collapse. For comparison, we have included the same data for a much more viscoplastic case ($U_{p}=0.50km/s$) in the supplemental material where it can be seen that the ejecta impact shoulder on the viscous heating is much smaller and the exothermic reactions occur at a much slower pace than Figure 6.

Less the adiabatic heating from the shock front, the shaded area indicating viscous heating has a mean temperature of 860K, and given the integrand of this segment, would equate to an area of 226 $nm^{2}$ at this temperature. The ejecta impact contributes significantly less to the overall hot-spot that is formed, its mean temperature being 1240K, which can be equated to an area of just 27 $nm^{2}$ at this mean temperature. Comparing these two mechanisms, the ejecta impact at this (moderately low) piston velocity contributes a factor of six less $Area\cdot Temperature$ than that of the viscous heating mechanism. In addition, the temperature histograms in Supplemental Figure S2 show that this ratio of hot-spot potency favors viscous heating by a factor of twenty-five. This trend continues up to piston velocities of 1.25km/s where the ratio decreases down to 1.78, but still in favor of viscoplastic heating being the dominant mechanism. Of course the time delay used to capture the ejecta impact is somewhat arbitrary and needs to be carefully chosen since HNS begins to react promptly after the impact occurssidecomment1 .

Up to this point, the discussion of hot spot mechanisms has been focused on idealized pore geometries, but with the SGL strength model in CTH we are able to address much more realistic shock responses of HNS by providing mesoscopic porous microstructures as input geometries. The insets to Figure 7 show a $15\mu m$ by $25\mu m$ slab of HNS with a pore size distribution that is artificially generated such that it matches experimentally observed microstructures. The sample is shocked at the bottom surface with a piston moving at 0.60km/s for both the Hydro and SGL strength models, each simulation snapshot and temperature histogram corresponds to when the shock wave reaches the free (top) surface. While Figure 7 is the culmination of heat generated from many collapsed pores, the SGL model clearly shows additional heat generated in the 400-600K range that is due to viscoplastic heating, see also Figure 5 A),C). Accurately predicting hot spot temperatures is a crucial step in predicting the detonation performance of HNS and many other energetic materials, and the SGL model presented here is a significant advancement for these simulation methods.