Shock-induced plasticity and structural changes in energetic molecular crystals are well documented. These processes couple with the leading shock wave and affect its propagation, resulting in long, transient responses that are challenging to capture with all-atom simulations due to their time scale. Hence, the effects of this coupling and the transient shock response on the formation of hotspots and the initiation of chemistry remain unclear. To address these challenges, we investigate the role of shock-induced plastic deformation on shock initiation with a recently developed particle-based, coarse-grain model for 1,3,5-trinitro-1,3,5-triazinane (RDX) that utilizes the generalized dissipative particle dynamics with reactions framework. This model enables reactive simulations at micron length scales, which are required to achieve steady-state shock propagation. The simulations show that the shock Hugoniot response of RDX can involve transient behavior for up to 150 ps before steady-state behavior is achieved for shock strengths above the elastic limit. Pore collapse simulations demonstrate that the intensity of the resulting hotspot will weaken as the shock transitions from transient to steady-state behavior, ultimately affecting the shock-to-deflagration transition. Our results highlight the importance of considering the mesoscopic effects of shock-induced plastic deformation in simulations of shock-to-deflagration transitions of high explosives.
Shock initiation of high explosives (HE) is a complex, multiscale phenomenon involving various physical and chemical processes. A crucial mechanism governing the shock initiation and transition to detonation is the formation of hotspots from the interaction of the propagating shock with the microstructural features of the material. Significant amounts of research efforts have examined various components of such interactions, including the size and shape of defects,1,2 crystal orientation,3,4 and shock strength.5,6 However, an overlooked aspect, especially from simulation studies, is the role of shock-induced plastic deformation on the interplay between shock propagation and hotspot formation.
Shock-induced structural changes occur when the shock stress surpasses the material's strength at the Hugoniot elastic limit (HEL). Experiments7–13 and simulations14–17 have shown that RDX (C3H6N6O6) responds to shock loading with a plethora of processes, such as dislocations and shear-band formation, stacking faults, and an α-to-γ crystal transition, which depend on the shock strength and orientation. These transformations act to relieve non-hydrostatic stresses in the material15,18 and significantly alter the shock Hugoniot response of the HE.14,16,19 Atomistic simulation methods utilizing Hugoniotstat15,18 or shock absorbing boundary conditions2,16,20,21 have been particularly useful for investigating such processes as they provide molecular-resolution data and allow simulation times on the order of nanoseconds. However, these techniques cannot be applied to simulate the interaction between a propagating shock and microstructure, as well as the associated hotspots1–3,6,22 and shock-to-deflagration transitions.5 Non-equilibrium shock simulations can, in principle, describe the interaction between a propagating shock and microstructural features. However, the long transient times required to achieve a steady propagating shock for certain conditions, associated with defect nucleation, suggest that the O(10 ps) duration atomistic simulations carried out to date likely include transient effects.
A viable method for investigating the effect of shock-induced nucleation of defects and other relaxation phenomena and achieving steady states in non-equilibrium shock simulations is particle-based, coarse-grain (CG) modeling. For example, our recent CG model parameterized for RDX23 that utilizes generalized dissipative particle dynamics with reactions (GenDPDE-RX) framework24,25 is a suitable approach as it is designed to capture mechanical, thermal, and chemical processes of atomistic models with orders of magnitude improvement in computational efficiency. Specifically, our model has been developed using molecular force fields26 and chemistry models27 that exhibit good agreement with experiments and atomistic simulations in capturing the crystal properties, liquid structure, phase transition, elastic constants, and reaction kinetics of RDX. In addition, the computational efficiency of the model and other CG models of similar resolution has allowed the investigation of systems up to 25 μm in length.28–31
In this study, we examine the effect of shock-induced relaxation processes on shock-to-deflagration following pore collapse in RDX using our CG model on μm-scale systems. First, the shock Hugoniot behavior as a function of particle velocity and run time/distance is investigated. Our results demonstrate that the material undergoes a transition from elastic Hugoniot response to two-wave regimes, followed by an overdriven wave depending on the shock strength, which is consistent with previous atomistic simulation studies.19 Importantly, in the two-wave regime, we observe significant temporal dependence on the shock speed, with transient behavior lasting ∼150 ps (corresponding to run distances of ∼700 nm) before a steady-state propagation is obtained. Our reactive CG simulations demonstrate that the shock-to-deflagration transition is dependent on the details of the shock, with transient shocks resulting in higher hotspot temperatures compared to steady-state shocks that facilitate the transition to deflagration. These findings highlight the importance of the interaction between mesoscopic relaxation and energy localization mechanisms in shock initiation of HE materials, an effect that has been previously overlooked in atomistic non-equilibrium shock simulations.
II. MODELS AND METHODS
A. CG RDX model
Here, RDX, INT, and PRD represent the reactant RDX, intermediate, and product species, respectively. The INT and PRD species do not correspond to specific species, but represent the cluster of species obtained from dimensionality reduction analysis on reactive atomistic simulation data. The reaction kinetics are governed by an Arrhenius equation , where θ is the internal temperature of the CG particle. The kinetics parameters are as given in Refs. 23 and 27.
Here, , are the position, momentum, velocity, internal energy, and mass of particle i, respectively. are, respectively, the conservative forces, dissipative forces, and random contribution to the momentum due to interactions between particles i and j. The primed variables indicate the state that is one time step after the current state , while the non-primed variables refer to the current state . and are thermal conduction and random contribution to the internal energy. We refer the reader to Ref. 24 for a detailed explanation of GenDPDE-RX and Ref. 23 for the parameters employed.
Here, , and represent the total pressure, excess pressure, and particle local density, respectively. is a unit-less parameter that represents the contribution to stemming from implicit degrees-of-freedom that have been coarse-grained into the CG particle, where in this work, .23, is given by the analytical equation-of-state of the exponential-6 potential33 as a function of , . The Exp-6 parameters ( ) are (12, 0.892 kcal/mol, 9 Å) and (14, 0.595 kcal/mol, 9.5 Å) for INT and PRD, respectively, as parametrized in Ref. 23.
Here, θ refers to the internal temperature representing the temperatures associated with the higher resolution degrees-of-freedom that have been coarse-grained. , , and denote the heat of formation, isobaric heat capacity, and number of molecules in particle i, respectively. indicates the species type (RDX, INT, and PRD for our model), and s is the number of species (three in our chemistry model). and correspond to the reference temperature (298.15 K) and molar fraction of species , respectively. We assume that is independent of θ, i.e., . As parametrized in Ref. 23, and . Lastly, while of Ref. 23 was obtained from DFT calculations, we use a classical expression for the heat capacity: , where 60 corresponds to the degrees-of-freedom of RDX molecule less the translational degrees-of-freedom of the CG particle. We use the classical model to facilitate future comparisons with atomistic molecular dynamics that are described by classical heat capacity models, but the conclusions from this study are qualitatively not affected by the choice of the heat capacity model.
B. Simulation setup and computational details
We used the Large-scale Atomic/Molecular Massively Parallel Simulator34 (LAMMPS, version 3Mar2020) software, with GenDPDE-RX implemented through modification to the USER-DPD package. The CG crystalline RDX samples were prepared by replicating a unit cell with lattice parameters a = 13.4 Å, b = 11.6 Å, and c = 10.6 Å, which is in accordance with the lattice constants for the employed force field of RDX.26 The crystal orientations were aligned such that  was along the x axis and  was along the z axis. The shock simulations were performed by driving the crystalline RDX samples into a momentum reflector.
For non-reactive Hugoniot simulations, the sample dimensions of the homogeneous crystal were 3 μm × 5 nm × 30 nm when shock was applied in the x  direction and 30 nm × 5 nm × 3 μm when shocking in the z  direction. For the reactive simulations, three different samples were prepared that included a 40 nm diameter cylindrical pore with varying distances of 100, 400, and 800 nm from the piston to the leading edge of the pore surfaces (dpore). The dimensions of the three samples were 120 × 5 nm2 in transverse directions and 1.27 μm (dpore = 100 nm), 1.47 μm (dpore of 400 nm), and 2 μm (dpore = 800 nm) for the axis in the shock direction. The purpose of the reactive pore collapse simulations is to compare the formation of hotspots and their subsequent reactivity, following the collapse of a pore driven by shock waves that are at various stages between transient and steady states.
Hugoniot responses were analyzed with one-dimensional Eulerian bins with a grid size of 1 nm in the shock direction. Shock speed (us) at an observation time (τobs) was determined by measuring the difference in the shock front positions at τobs + 10 ps vs τobs – 10 ps. Pore collapse simulations were analyzed with two-dimensional bins with grid sizes of 1 × 1 nm2. Temperature distributions from hotspot formation were analyzed by θ-area plots that depict the cumulative area that is above a certain θ value.
III. RESULTS AND DISCUSSIONS
A. Hugoniot response
The shock Hugoniot curves of the CG RDX model along  and  directions are characterized for impact velocities (up) varying from 0.25 to 2.5 km/s and at an observation time (τobs) up to 300 ps. The transient Hugoniot response of RDX with shock in the  direction is illustrated in Fig. 1. As reported in Ref. 23 and in agreement with previous atomistic simulations,19 we observe three distinct regimes depending on up. When up ≤ 0.5 km/s, RDX is below its HEL and displays an elastic response. In this regime, the steady state of us is reached in a very short time [Fig. 1(a)]. As expected, no shear strain on the shocked material [middle snapshot of each triplet in Fig. 1(e)] is observed and the velocity along the shock direction (vx, left snapshot of each triplet) and the internal temperature (θ, right snapshot of each triplet) are found to be independent of τobs [Fig. 1(e), up = 0.5 km/s].
In the two-wave regime (0.75 km/s ≤ up ≤ 1.5 km/s), we observe a leading elastic precursor followed by a plastic wave. This regime is characterized by a relatively long transient regime required for the plastic wave to fully develop, during which the velocity of the leading elastic wave decreases. Figure 1(a) shows the instantaneous velocity of the leading wave as a function of observation time τobs. For all piston velocities in this regime, steady-state propagation of the leading wave (i.e., constant us) is attained after ∼150 ps. As depicted in Fig. 1(e), the leading elastic precursor region between elastic and plastic shock fronts exhibits homogeneous vz and θ, as well as, by definition, a lack of localized shear strain. This is followed by a plastic region that displays significant localized shear strain and elevated temperature compared to the elastic region. The elastic region expands over time even after steady-state propagation of both elastic and plastic waves is achieved, as the elastic wave is faster than the plastic wave. The temperature, pressure (see Fig. S1), and wave speed of the elastic precursor are nearly independent of up within the two-wave regime, while they increase with up in the plastic regime until the overdriven regime is achieved.
When up > 1.5 km/s, an overdriven regime is observed where the plastic wave is faster than the elastic wave, leading to a single-wave structure as the plastic wave catches up to the elastic wave. In this regime, steady-state propagation of shock is attained after τobs ∼ 50 ps and the change in us is less than those observed in the two-wave regime [see Fig. 2(a)]. In addition, we observe relatively homogeneous vz and θ across the shocked region.
Next, we examine the effect of the crystal orientation on shock propagation by analyzing the Hugoniot along the  direction (Figs. S1 and S2). The Hugoniot response in the  direction is similar to the response observed for shock in the  direction. We observed elastic, two-wave, and overdriven regimes with quantitatively similar values of steady-state temperature and pressure compared to results from the  direction (see Fig. S1). Here, we again refer to the steady-state as τobs in which us of the leading shock wave reaches a constant value. A few distinctions include the HEL at up = 1 km/s and a crisscross pattern of shear-band formation in the two-wave regime. While us of the leading elastic wave remains nearly constant throughout the two-wave regime for the  case as well, we observe a distinct decrease at up = 1 km/s [Fig. S2(b)]. We attribute such a decrease to the significant fluctuation in the plastic wave us that is not fully resolved for this particular case. While our simulation systems are significantly larger than those typically explored with atomistic models, the length scale, especially lateral to the shock direction, may not be sufficient to attain stable shock conditions for this case.
B. Shock initiation and deflagration transition
We evaluate the effect of the transient nature of the shockwave on the shock-to-deflagration transition of the CG RDX model via shock-induced pore collapse simulations. To study shocks at different stages of evolution within the two-wave regime, three μm-scale samples were prepared each with a single 40 nm diameter cylindrical pore located at different distances from the impact surface: 100, 400, and 800 nm from the impact plane (dpore), see Fig. 2(a). The location of the pore at dpore = 100 nm was chosen for two reasons. First, atomistic simulations1,5,6,35–37 generally locate pores of similar size at 40–150 nm from the piston. Second, the shock front reaches the pore in approximately 20 ps, which is insufficient time for the shock to reach the steady-state, for the two-wave and overdriven regime cases. On the other hand, the shock front reaches the pore in approximately 80 ps when dpore is 400 nm, which is sufficient to attain a steady-state shock for the overdriven regime but not the two-wave regime. Lastly, the dpore = 800 nm case allows an investigation of the steady-state shock behavior for all wave regimes.
The cumulative area vs internal temperature (θ) plots for the two systems with varying dpore are depicted in Figs. 2(b)–2(e). These plots depict the area with temperature exceeding a value θ vs θ and are particularly useful for analyzing the shock initiation behavior of hotspots, as their criticality condition is a function of both the temperature and size of the hotspots.38 In these plots, solid lines represent the results corresponding to 5 ps after the total collapse of the pore, while dashed lines depict the results from t = 150 ps after the collapse. As shown in Fig. 2(b), RDX quenches for up = 0.5 km/s in all cases and temperature distributions are nearly independent of the location of the pore. The results for up = 2 km/s, in Fig. 3(e), also show consistent behavior for all cases, in this case, the three systems deflagrate regardless of dpore, with nearly identical temperature distributions. The equivalence of the shock response of the systems in the elastic and overdriven regimes is consistent with our analysis of the transient Hugoniot response, where we observed that the material reaches the steady-state in τobs < 50 ps, and that the changes in shock velocities are not as significant for elastic and overdriven regimes as compared to the two-wave regime cases.
Figures 2(c) and 2(d) depict the temperature distributions for shock strengths within the two-wave regime. Here, we observe a clear distinction in the temperature distributions depending on the location of the pore with respect to the impact plane, dpore. The case with dpore = 100 nm exhibits temperatures 200 ∼ 500 K higher than the systems with dpore = 400 and 800 nm. Furthermore, for up = 1.125 km/s, this difference leads to a divergence in the deflagration behavior; the material deflagrates for the dpore = 100 nm case and quenches for the dpore = 400 and 800 nm cases. Specifically, we observe quenching for all cases at up = 1 km/s [Fig. 2(c)], divergence of response for up = 1.125 and 1.25 km/s [Fig. S3(a)], and deflagration for all cases at up = 1.375 and 1.5 km/s [Figs. S3(b) and S3(c)]. For the critical cases (up ≥ 1.375 km/s), the rates of growth of the hotspots significantly depend on dpore, with the dpore = 100 nm case expanding most rapidly.
The process of hotspot formation in the two-wave regime merits’ additional analysis. Figure 3 shows the temperature profiles at various times during pore collapse for the up = 1 km/s shock along . The upper panels in Figs. 3(a)–3(c) represent local temperature (θ) vs position along the shock direction (x) for the cross section of , which includes the pore (see rightmost sketch). The lower panels depict temperature profiles on a cross section of , which does not include the pore. The plotted simulation times correspond to (1) when the leading shock wave reaches the leading edge of the downstream pore surface (red); (2) the approximate midpoint between the times at which the shock reaches the leading edge of the pore surface and the total collapse of the pore occurs (green); (3) the time at which the total collapse of the pore occurs (blue); (4) 5 ps (purple) after total collapse; (5) 15 ps (yellow) after total collapse; (6) 25 ps (cyan) after total collapse; and (7) 40 ps (black) after total collapse.
The processes of pore collapse and hotspot formation differ significantly depending on the location of the pore. First, the pore collapse rate is faster for dpore = 100 nm (18 ps for total pore collapse) compared to the other cases (∼23 ps) because shock velocity of the leading wave is greater for dpore = 100 nm and decreases with time as the plastic wave develops. This causes the initial hotspot temperature to be approximately 80 K higher for dpore = 100 nm as compared to the dpore = 400 nm and dpore = 800 nm cases, which show similar temperatures. Second, for dpore = 100 nm, the trailing plastic wave reaches the hotspot around the time of the full collapse of the pore, while for greater dpore cases, the plastic wave arrives later and interacts with the hotspot created by the pore collapse. This is because the distance between elastic and plastic wave fronts increases with time. In addition, the trailing plastic wave interacts with the secondary waves generated by the interaction of the leading wave and the pore, most noticeably the rarefaction fan propagating back toward the impact plane. The simulation snapshots in Figs. 3(d)–3(f) show that the trailing plastic wave is weakened in the region upstream from the hotspot by this rarefaction wave. However, the plastic waves surround the hotspot and increase the overall temperature. Quite remarkably, the core temperature of the hotspot does not decrease for several tens of picoseconds after collapse; suggesting that the additional energy delivered by the plastic wave may mitigate the heat dissipation of the hotspots. Future research scrutinizing the interaction between the plastic waves with hotspot formation will contribute to increasing our understanding of HE shock sensitivity.
Lastly, the temperature distributions of the hotspots that result from shock in the  direction are depicted in Fig. S4. Similar to the response shown for shock in the  direction, the temperature distributions for up = 0.5 and 2 km/s [Figs. S4(a) and S4(d)] do not depend significantly on dpore, while the distributions in the two-wave regime are greatly affected by dpore. For this orientation, the hotspots formed with up = 1 km/s [Fig. S4(b)] display divergence of the deflagration behavior compared to up = 1.125 km/s for shock in the  direction.
IV. DISCUSSIONS AND CONCLUSION
This study investigates the role of shock-induced plastic deformation on the formation of hotspots due to pore collapse and criticality. Explicit shock simulations using a particle-based coarse-grained model for RDX demonstrate that the steady-state wave propagation in the two-wave regime can require up to 150 ps and 800 nm of run distances to be established. The transient shockwave regime involves a reduction in the strength of the leading elastic shock as the plastic wave develops. As an important consequence, we observe significant differences in the temperature distributions of the hotspots created from the collapse of pores depending on the transient shock response that can lead to differences in the deflagration behavior of the material.
As with any physics-based model of a material, our model is not without approximations and it is important to assess the robustness of the results. For example, the CG force field has a known deficiency in replicating the morphology of the shear bands observed from atomistic simulations.39 In addition, while atomistic simulations15,18 have predicted that the orientation significantly affects the shock-induced structural transformation and the associated changes in energies, our model does not fully account for such behavior. However, while a direct comparison of our results with atomistic simulations is not possible due to the length scales in this study, many of our results are in qualitative agreement with previous studies employing atomistic simulations. For example, Bedrov et al.15 investigated the time-dependent response of RDX through atomistic, uniaxial Hugoniotstat simulations. Their results demonstrated that RDX above its HEL and compressed in the  direction initially displays elastic response up to 10 ps and then exhibits plastic deformation, with stress relaxation requiring a 100 ∼ 200 ps incubation period. The reported HEL pressure of approximately 8 GPa was also qualitatively similar in magnitude to those from the CG model (6.25 GPa for  and 8.93 GPa for ). These qualitative agreements in the time scale and the HEL support the importance of considering the effect of shock-induced plastic deformation on the shock response of HE. However, as these results are from Hugoniotstat simulations, future studies of μm-scale direct shock simulations with an atomistic model would provide further validation as well as additional insight into the effect of shock-induced plasticity that was not fully resolved by our study.
The findings of this study are especially important to consider when simulating hotspot formation from pore collapse simulations using atomistic models. To our knowledge, systems studied using atomistic models have been limited to approximately 300 nm, while the distance from the piston to the pores is generally in the range of 40–150 nm.1,5,6,35–37 Those length scales are significantly smaller than the required run distance to obtain steady-state shock for both the two-wave and overdriven regimes. This suggests that the temperature of the hotspots formed in those atomistic model simulations may be overestimated in such systems and that the transient Hugoniot response of the HE should be considered in the analysis of those results.
The supplementary material includes pressure and temperature data from Hugoniot simulations, Hugoniot response in the  direction, and additional shock-to-deflagration simulation results for shocks in  and  directions.
We thank Dr. Garrett Tow, Dr. Brenden Hamilton, and Jason Wilkening for their insightful discussion. This research was sponsored by the Army Research Laboratory and was accomplished under Cooperative Agreement No. W911NF-20-2-0189. This work was supported in part by high-performance computer time and resources from the DoD High-Performance Computing Modernization Program. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Laboratory, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.
Conflict of Interest
The authors have no conflicts to disclose.
Brian H. Lee: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Writing – original draft (lead); Writing – review & editing (equal). James P. Larentzos: Funding acquisition (equal); Methodology (supporting); Software (lead); Supervision (equal); Writing – review & editing (equal). John K. Brennan: Funding acquisition (equal); Investigation (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). Alejandro Strachan: Conceptualization (equal); Formal analysis (equal); Funding acquisition (lead); Supervision (lead); Writing – review & editing (equal).
The data that support the findings of this study are available from the corresponding author upon reasonable request.