Effect of shock-induced plastic deformation on mesoscale criticality of 1,3,5-trinitro-1,3,5-triazinane (RDX)

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₃H₆N₆O₆) 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.⁵ 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²³ 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²⁶ and chemistry models²⁷ 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.¹⁹ 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.

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⁡ ( − E R x k B θ )$⁠, 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δ 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 $( Δt+t)$⁠, while the non-primed variables refer to the current state $(t)$⁠. $q ˙ i j$ and $δ 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.

Here, $π i, π i e x$⁠, and $n i$ represent the total pressure, excess pressure, and particle local density, respectively. $χ$ is a unit-less parameter that represents the contribution to $π i$ stemming from implicit degrees-of-freedom that have been coarse-grained into the CG particle, where in this work, $χ=13.9878$⁠.^23, $π i$ is given by the analytical equation-of-state of the exponential-6 potential³³ as a function of $n i$⁠, $θ i$⁠. The Exp-6 parameters (⁠ $α,ε, 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.

Here, θ refers to the internal temperature representing the temperatures associated with the higher resolution degrees-of-freedom that have been coarse-grained. $Δ 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. $ξ$ 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 ξ$ correspond to the reference temperature (298.15 K) and molar fraction of species $ξ$⁠, respectively. We assume that $Δ 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, $Δ H r x n , R 1=69.415 kcal / mol$ and $Δ H r x n , R 2=−545.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.

We used the Large-scale Atomic/Molecular Massively Parallel Simulator³⁴ (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.²⁶ 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² 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². Temperature distributions from hotspot formation were analyzed by θ-area plots that depict the cumulative area that is above a certain θ value.

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,¹⁹ 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.

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.³⁸ 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≤z≤80 nm$⁠, which includes the pore (see rightmost sketch). The lower panels depict temperature profiles on a cross section of $10 nm≤z≤30 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.

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.³⁹ 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.¹⁵ 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.

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.

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.

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.