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.

## I. INTRODUCTION

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). Experiments^{7–13} and simulations^{14–17} have shown that RDX (C_{3}H_{6}N_{6}O_{6}) 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 material^{15,18} and significantly alter the shock Hugoniot response of the HE.^{14,16,19} Atomistic simulation methods utilizing Hugoniotstat^{15,18} or shock absorbing boundary conditions^{2,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 hotspots^{1–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 RDX^{23} that utilizes generalized dissipative particle dynamics with reactions (GenDPDE-RX) framework^{24,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 fields^{26} and chemistry models^{27} 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

^{23}that utilizes the GenDPDE-RX method. A detailed description and assessment of the accuracy of the model is provided in Ref. 23. Briefly, this model describes each RDX molecule as a single CG particle that acts as a continuously stirred chemical reactor, where within each CG particle, the molecule decomposes according to a prescribed chemistry model. To describe RDX decomposition chemistry, we utilize a reduced order chemical kinetics model

^{27}that represents the species observed from reactive atomistic simulation of RDX decomposition as three components via a dimensionality reduction method, as follows:

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 $ k R x= A R xexp\u2061 ( \u2212 E R x k B \theta )$, where *θ* is the internal temperature of the CG particle. The kinetics parameters are $ A R 1=32.104/ ps, A R 2=11.81/ ps, E R 1=24.4 kcal/ mol, E R 2=23.5 kcal/ mol$ as given in Refs. 23 and 27.

Here, $ r i, p i$, $ v i,$ $ u i, and m i$ are the position, momentum, velocity, internal energy, and mass of particle *i*, respectively*.* $ f i j C, f i j D, and\delta p i j R$ 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 $ ( \Delta t+t)$, while the non-primed variables refer to the current state $(t)$. $ q \u02d9 i j$ and $\delta u i j R$ 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.

^{26}For the reacted RDX, the forces are described by a many-body potential based on an exponential-6 equation-of-state (MB-Exp-6-EOS)

^{32}given by

Here, $ \pi i, \pi i e x$, and $ n i$ represent the total pressure, excess pressure, and particle local density, respectively. $\chi $ is a unit-less parameter that represents the contribution to $ \pi i$ stemming from implicit degrees-of-freedom that have been coarse-grained into the CG particle, where in this work, $\chi =13.9878$.^{23}^{,} $ \pi i$ is given by the analytical equation-of-state of the exponential-6 potential^{33} as a function of $ n i$, $ \theta i$. The Exp-6 parameters ( $\alpha ,\epsilon , r m$) 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.

*θ*) and the internal energy (

*u*) of the CG particle by the following expression:

Here, *θ* refers to the internal temperature representing the temperatures associated with the higher resolution degrees-of-freedom that have been coarse-grained. $ \Delta H f , i$, $ C P , i 0$, and $ N m o l$ denote the heat of formation, isobaric heat capacity, and number of molecules in particle *i*, respectively. $\xi $ indicates the species type (RDX, INT, and PRD for our model), and *s* is the number of species (three in our chemistry model). $ T r$ and $ f i \xi $ correspond to the reference temperature (298.15 K) and molar fraction of species $\xi $, respectively. We assume that $ \Delta H r x n$ is independent of *θ*, i.e., $ C P , R D X 0= C P , I N T 0= C P , P R D 0$. As parametrized in Ref. 23, $ \Delta H r x n , R 1=69.415 kcal / mol$ and $ \Delta H r x n , R 2=\u2212545.852 kcal / mol$. Lastly, while $ C P 0$ of Ref. 23 was obtained from DFT calculations, we use a classical expression for the heat capacity: $ C p 0=60 k B$, 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 Simulator^{34} (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 [100] was along the *x* axis and [001] 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* [100] direction and 30 nm × 5 nm × 3 *μ*m when shocking in the *z* [001] 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 (*d*_{pore}). The dimensions of the three samples were 120 × 5 nm^{2} in transverse directions and 1.27 *μ*m (*d*_{pore }= 100 nm), 1.47 *μ*m (*d*_{pore} of 400 nm), and 2 *μ*m (*d*_{pore }= 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 (*u*_{s}) 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 nm^{2}. 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 [100] and [001] directions are characterized for impact velocities (*u*_{p}) 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 [100] 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 *u*_{p}. When *u*_{p} ≤ 0.5 km/s, RDX is below its HEL and displays an elastic response. In this regime, the steady state of *u*_{s} 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 (*v*_{x}, left snapshot of each triplet) and the internal temperature (*θ,* right snapshot of each triplet) are found to be independent of τ_{obs} [Fig. 1(e), *u*_{p }= 0.5 km/s].

In the two-wave regime (0.75 km/s ≤ *u*_{p} ≤ 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 *u*_{s}) is attained after ∼150 ps. As depicted in Fig. 1(e), the leading elastic precursor region between elastic and plastic shock fronts exhibits homogeneous *v*_{z} 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 *u*_{p} within the two-wave regime, while they increase with *u*_{p} in the plastic regime until the overdriven regime is achieved.

When *u*_{p} > 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 *u*_{s} is less than those observed in the two-wave regime [see Fig. 2(a)]. In addition, we observe relatively homogeneous *v*_{z} and *θ* across the shocked region.

Next, we examine the effect of the crystal orientation on shock propagation by analyzing the Hugoniot along the [001] direction (Figs. S1 and S2). The Hugoniot response in the [001] direction is similar to the response observed for shock in the [100] direction. We observed elastic, two-wave, and overdriven regimes with quantitatively similar values of steady-state temperature and pressure compared to results from the [100] direction (see Fig. S1). Here, we again refer to the steady-state as τ_{obs} in which *u*_{s} of the leading shock wave reaches a constant value. A few distinctions include the HEL at *u*_{p }= 1 km/s and a crisscross pattern of shear-band formation in the two-wave regime. While *u*_{s} of the leading elastic wave remains nearly constant throughout the two-wave regime for the [001] case as well, we observe a distinct decrease at *u*_{p }= 1 km/s [Fig. S2(b)]. We attribute such a decrease to the significant fluctuation in the plastic wave *u*_{s} 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 (*d*_{pore}), see Fig. 2(a). The location of the pore at *d*_{pore }= 100 nm was chosen for two reasons. First, atomistic simulations^{1,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 *d*_{pore} is 400 nm, which is sufficient to attain a steady-state shock for the overdriven regime but not the two-wave regime. Lastly, the *d*_{pore} = 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 *d*_{pore} 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 *u*_{p }= 0.5 km/s in all cases and temperature distributions are nearly independent of the location of the pore. The results for *u*_{p} = 2 km/s, in Fig. 3(e), also show consistent behavior for all cases, in this case, the three systems deflagrate regardless of *d*_{pore}, 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, *d*_{pore}. The case with *d*_{pore }= 100 nm exhibits temperatures 200 ∼ 500 K higher than the systems with *d*_{pore} = 400 and 800 nm. Furthermore, for *u*_{p} = 1.125 km/s, this difference leads to a divergence in the deflagration behavior; the material deflagrates for the *d*_{pore }= 100 nm case and quenches for the *d*_{pore} = 400 and 800 nm cases. Specifically, we observe quenching for all cases at *u*_{p }= 1 km/s [Fig. 2(c)], divergence of response for *u*_{p} = 1.125 and 1.25 km/s [Fig. S3(a)], and deflagration for all cases at *u*_{p} = 1.375 and 1.5 km/s [Figs. S3(b) and S3(c)]. For the critical cases (*u*_{p} ≥ 1.375 km/s), the rates of growth of the hotspots significantly depend on *d*_{pore}, with the *d*_{pore} = 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 *u*_{p }= 1 km/s shock along [100]. The upper panels in Figs. 3(a)–3(c) represent local temperature (*θ*) vs position along the shock direction (*x*) for the cross section of $40 nm\u2264z\u226480 nm$, which includes the pore (see rightmost sketch). The lower panels depict temperature profiles on a cross section of $10 nm\u2264z\u226430 nm$, 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 *d*_{pore }= 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 *d*_{pore }= 100 nm and decreases with time as the plastic wave develops. This causes the initial hotspot temperature to be approximately 80 K higher for *d*_{pore }= 100 nm as compared to the *d*_{pore }= 400 nm and *d*_{pore }= 800 nm cases, which show similar temperatures. Second, for *d*_{pore }= 100 nm, the trailing plastic wave reaches the hotspot around the time of the full collapse of the pore, while for greater *d*_{pore} 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 [001] direction are depicted in Fig. S4. Similar to the response shown for shock in the [100] direction, the temperature distributions for *u*_{p} = 0.5 and 2 km/s [Figs. S4(a) and S4(d)] do not depend significantly on *d*_{pore}, while the distributions in the two-wave regime are greatly affected by *d*_{pore}. For this orientation, the hotspots formed with *u*_{p} = 1 km/s [Fig. S4(b)] display divergence of the deflagration behavior compared to *u*_{p} = 1.125 km/s for shock in the [100] 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 simulations^{15,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 [100] 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 [100] and 8.93 GPa for [001]). 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.

## SUPPLEMENTARY MATERIAL

The supplementary material includes pressure and temperature data from Hugoniot simulations, Hugoniot response in the [001] direction, and additional shock-to-deflagration simulation results for shocks in [100] and [001] directions.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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