Shock-Induced Hot Spot Formation and Spalling in 1,3,5-trinitroperhydro-1,3,5-triazine Containing a Cube Void

The initial reaction mechanism of energetic materials under impact loading and the role of crystal properties in impact initiation and sensitivity are still unclear. In this paper, we report reactive molecular dynamics simulations of shock initiation of 1,3,5-trinitroperhydro-1,3,5-triazine (RDX) crystals containing a cube void. Shock-induced void collapse, hot spots formation and growth, as well as spalling are revealed to be dependent on the shock velocity. The void collapse times are 1.5 and 0.7 ps, for the shock velocity of 2 and 4 km·s–1, respectively. Results indicate that the initial hot spot formation consists of two steps: one is the temperature rise caused by local plastic deformation and the other is the temperature increase resulting from the collision of upstream and downstream particles during the void collapse. Whether hot spots will continue to grow or quench depends on sensitive balance between energy release caused by local physical and chemical reactions and various heat dissipation mechanisms. In our simulations, hot spot would grow for Up = 4 km·s–1; hot spot is weak to some extent for Up = 2 km·s–1. The tensile wave reflected by the shock wave after reaching the free surface causes the spalling, which depends on the initial shock velocity. Typical spalling occurs for the shock velocity 2 km·s–1, while the tensile wave induces the microsplit region in RDX crystals in the case of Up = 4 km·s–1. Chemical reactions are studied for Rankine–Hugoniot shock pressures Ps = 14.4, 57.8 GPa. For the weak shock, there is almost no decomposition reaction of the RDX molecules near the spalling region. On the contrary, there are large number of small molecule products, such as H2O, CO2, NO2, and so forth, around the microsplit regions for the strong shock. The ruptures of N–NO2 bond are the main initial reaction mechanisms for the shocked RDX crystal and are not affected by shock strength, while the microsplit slows down the decomposition rate of RDX. The work in this paper can shed light on a thorough understanding of thermal ignition, hot spot growth, and other physical and chemical phenomena of energetic materials containing voids under impact loading.


INTRODUCTION
Condensed phase energetic materials usually contain defects (such as crystal dislocations, shear bands, defects, voids, etc.), making the energetic material heterogeneous, which is beneficial for hot spots formation and detonation under impact loading. 1−5 Hot spots are localized regions of high temperature and high pressure in materials that serve as nucleation sites for initiating and possibly sustaining rapid chemistry. Various mechanisms have been provided to explain the formation mechanism of "hot spots": void collapse, adiabatic shear bands, fraction, and heating at dislocation pileups. 6−9 Because of the complex physical and chemical coupling effects and the elastoplastic deformation of the energetic material, as well as very short time scale and space scale, the current experimental techniques are difficult to reveal the detailed process of the shock initiation of energetic materials from an atomic and molecular level, while the molecular dynamics (MD) simulation method can be used to reveal the formation process of hot spots and the progress of chemical reactions from the microscopic scale; the latter also provides a new way to analyze the physical and chemical processes of energetic materials under shock loading.
The work of Holian et al. 10 showed that the reason for the formation of hot spots is the evaporation of materials into the pores and subsequent recompression in the process of pore collapse in Lennard-Jones system. Hatano 11 used 3-D Lennard-Jones system simulations containing cuboid voids of varying size and found that increasing cross-sectional area perpendicular to the shock direction increases the number of energetic intermolecular collisions, whereas the increase of the gap length parallel to the impact direction increases the peak temperature resulting from void collapse. MD simulations have also been used to investigate cubane nitrogen with square pore collapse that produces jetting. 12 Herring et al. 13 studied the effect of the arrangement of circular pores with various sizes on the formation process of impact-induced hot spots and found that the wide pressure wave generated by the hot spots seems to be the main mechanism of the detonation transition. Shan and Thompson. 14,15 simulated the hot spots formation and growth in shocked pentaerythritol tetranitrate (PETN) crystals. They found an asymmetrical evolution of hot spots growth and secondary circular shockwave, formed upon void collapse. Long and Chen. 16 simulated the void collapse and hot spots formation in shocked polymer-bonded explosives (PBXs) as well. The results indicated out that the hot spots formation should undergo three processes, namely void collapse, heat transformation, and temperature relaxation. Zhou et al. 17 focused on the hot spots formation of octahydro-1,3,5,7tetranitro-1,3,5,7-tetrazocine (HMX) crystals with various spherical voids size and impact strengths. They found out that increasing void size leads to a slower shock wave velocity and a more curved shock wave front. Jaramillo et al. 18 characterized the atomic-level mechanisms of plastic deformation in HMX crystal under dynamical loading. The results showed that the deformation of the HMX crystal passes through the elastic deformation, the plastic deformation, and the nanoscale shear band formation, and these deformations are strongly dependent on the impact strength. An et al. 19,20 reported on the shock-induced behavior of PBX based on the silapentaerythritol tetranitrate, 1,3,5-trinitroperhydro-1,3,5-triazine (RDX). They found hot spots also occur when a plane shock wave propagates to the heterogeneous cross section of a polymer explosive due to the difference in density of the two materials in the polymer. Eason and Sewell. 21 considered that strong shock produces a fluid-like jet structure during collapse in the RDX crystals. Strachan et al. 22 studied shock propagation via nonequilibrium MD simulations at various collision velocities and found that local overheating caused by impact promotes chemical reaction. Nomura et al. 23 reported million atoms reactive force field MD simulation of shock initiation of RDX crystal with a nanometer-scale void. They indicate that the nanojet formed during void collapse is related to the free volume of the void to enhance the intermolecular collision.
Although quite a few researches have been done by using MDs simulation methods for the initial reaction mechanism of energetic materials under impact loading, the role of crystal properties in impact initiation and sensitivity remains to be understood. In this paper, the void collapse process of RDX with a cube void under various shock velocities is simulated by MDs method with the reactive force field (ReaxFF)/lg 24 force field. The formation of hot spots and shock-induced spallation during the interaction between shock wave and cube void are to be studied in order to explain the detonation mechanism of energetic materials from an atomic level. The work should give insight for thermal ignition and hot spots growth of energetic materials.

SIMULATION DETAILS
2.1. Details of ReaxFF. The first principles based ReaxFF reaction force field can be used to the formation and breakage of chemical bonds in energetic materials under extreme conditions by using the relationship among bond length, bond order and bond energy, and the dynamic transfer of charges above atoms. 25 The ReaxFF is a bond-order potential based from quantum mechanical (QM) calculations, which provides accurate descriptions of complex chemical reactive processes (see bond-order cutoff in Supporting Information Table S1). As an extension, the ReaxFF/lg 24 added a longrange correction term using a low-gradient model Elg based on ReaxFF (see more details about ReaxFF/ g in Supporting Information S2). The simulation results of the crystal structure and state equation of energetic materials using ReaxFF/lg potential are in good agreement with experimental results. 17,26−28 2.2. Details of MD Simulation. The initial RDX crystal unit cell contains 8 molecules. 29 In order to obtain calculation accuracy and study the microscopic response of local shock wave propagating in RDX materials, 30 a supercell containing 400 unit cells was built by extending the initial RDX unit cell in the x, y, and z directions (40 × 10 × 1). The dimensions in the x, y, and z directions are 530.235, 115.740, and 10.709 Å respectively, and the initial density is around 1.80 g·cm −3 . The unit cell contains a cube void with a size of 60 × 60 × 10.709 Å 3 to study the shock-induced hot spots formation of RDX crystals as shown in Figure 1. A nonperiodic boundary condition is applied in the x direction, while the periodic boundary conditions are applied in the y and z directions. The structure of the supercell is first optimized by energy minimization, and then the relaxation process is carried out at 300 K using the NVT ensemble for 2 ps. The particles in the system are driving against the reflection wall at the y−z plane at velocities of 2 and 4 km·s −1 in the negative x direction to generate a shock wave propagating along the x-direction. Reflection wall is applied to the left boundary of the system, which means that if an atom moves outside the wall on a timestep by a distance delta, then it is wrapped back inside the other face by the same delta, and the sign of the corresponding component of its velocity is flipped. 31 The shock simulations were performed using ReaxFF/lg potential in large-scale atomic/molecular massively parallel simulator 31,32 with a 0.1 fs time step under the NVE ensemble.

RESULTS AND DISCUSSION
3.1. Void Collapse and Hot Spot Formation. Figure 2 shows the spatial−temporal temperature distribution in process of void collapse for the case of U p = 2 km·s −1 (a,c,e) and U p = 4 km·s −1 (b,d,f). The temperature in each simulation is calculated after subtracting the center-of-mass (COM) velocity of all the atoms using the previous method. 18,33 More details about the method are shown in the Supporting Information S3. Once the shock wave propagates to the upstream surface of the void, the molecules on the upstream surface begin to enter the void, causing local plastic deformation. Figure 2a,b shows the temperature distributions when the shock wave reaches the upstream surface of the void under the impact velocity of 2 and 4 km·s −1 . The temperature rise at this time is mainly caused by shock compression and local plastic deformation. Comparing Figure 2e,f, it can be seen that at low impact velocity (2 km·s −1 ), the upstream free surface starts to collide with the downstream interface at around 3.5 ps, until the void is closed at around 5.0 ps; at high impact velocity (4 km·s −1 ), the void is fully closed at around 3.2 ps, with the corresponding void collapse time of 1.5 and 0.7

ACS Omega
Article ps, respectively. It also shows that the time intensity of void collapse depends on the velocity of collision. The shock wave velocity U s averaged before the void collapse is calculated to be 4.0, 8.5 km·s −1 for the case of U p = 2 km·s −1 , U p = 4 km·s −1 , respectively.
In the process of void collapse, the upstream and downstream particles collide continuously inside the void, chemical reactions occur, and heat is released. The internal temperature of the void is higher than 2500 K, and the hot spot where the shape remains radially symmetrical is formed at the end of void collapse. The initial hot spot formation mainly consists of two stages: the temperature rise caused by local plastic deformation and the temperature increase resulting from the collision of upstream and downstream particles during the void collapse. For U p = 2 km·s −1 , the maximum temperature of the hot spot region during the void collapse at t = 5.0 ps is 3300 K, while the peak temperature of the hot spot region after the void closure is above 6000 K at the case of U p = 4 km·s −1 . The upstream and downstream particles collide with each other to induce chemical reactions, heat release, as well as secondary compression waves generated by upstream particles colliding with downstream surfaces, resulting in temperatures up to 8000 K at individual locations in void, as shown in Figure 2f.
As the void collapse process ends, the shock wave propagates forward in the form of convex plane wave at both collision velocities, which is different from the concave plane wave formed by spherical void collapse. 17 Figure 3 shows the averaged atom velocities V x , V y during the void collapse for the case of the U p = 4 km·s −1 . As shown in the figure, the points A and B, C and D used to observe the convex shock wave are located at the upper and lower ends on the left side of the void,

ACS Omega
Article the right side of the void, respectively. As the shock wave arrives at the left side of void, the velocities in the x and y directions of point A is 2.19 and 1.77 km·s −1 , respectively, as shown in Figure 3a,d. According to the theorem for composition of velocities, the direction of resultant velocity of point A is with the angle of approximately 39°to the x positive direction, while the resultant velocity of point B is at about −55°to the x positive direction. It can be seen that the shock wave propagates in a convex plane. The hydrodynamic fluid jets 21,23 were not observed in the process of the void collapse for the two reasons: void size was too small and the particle velocity along the x direction was always greater than the particle velocity along the y direction.
3.2. Hot Spot Growth. Nucleation and growth of hot spots formed after void collapse is an important process to induce chemical reactions in the system. 15,16 It can be seen from Figure 2e,f that the various shock velocities lead to the different sizes and temperature of the hot spot when the voids are closed. Therefore, it is necessary to study the growth process of hot spots after void collapse. Figure 4 presents color maps of temperature for both shock velocities at three times of interest: initial void collapse (a,b), at the time of maximum compression (c,d), and the time when spalling occurs (e,f). As can be seen from Figure 4a,b, the size of the hot spot increases after the void collapse, while the temperature in the hot spot region drops a little, which it is still higher than 2000 K. Once the system reaches the maximum compression, the temperature in the central region of the hot spot decreases slightly because of the heat released by the chemical reaction in the collision of molecules in the void. At this time, the size of the hot spot still grows. The shape of hot spot is close to hemispherical and remains radially symmetrical for the case of U p = 4 km·s −1 , as shown in Figure 4d. When the shock wave arrives the free surface of the right end, it emits a reverse stretching wave and interacts with the shock wave along the x direction in the crystal, which makes the hot spot grow. When spalling occurs, the growth of hot spots slows down and the temperature gradually decreases for U p = 2 km·s −1 . The wavefront temperature is about 1500 K, as shown in Figure 4e; however, the hot spots growth is not observed. The strong interaction between the shock wave and the stretching wave results in more kinetic transformation into heat in the shock region, and the shape of the hot spot changes to asymmetry at t = 9.0 ps. The hot spot becomes larger and hotter as shown in Figure 4f.
Whether hot spots will continue to grow or quench depends on the balance between energy release caused by local physical and chemical reactions and various heat dissipation mechanisms. 8 Combining with the time evolution of the size and peak temperature of hot spot, it can be seen that hot spot would grow for U p = 4 km·s −1 , while quench to some extent for U p = 2 km·s −1 in our simulation.

ACS Omega
Article 3.3. Spalling. The shock wave reflects from the free surface at the right end of the RDX crystal to form a tensile wave and thus results in local tensile stress on the free surface at the right end. If the tensile stress satisfies the dynamic fracture criterion, the material will fail and the fragments would fly out with a certain momentum. Figure 5 illustrates the snapshots of the spalling process of RDX crystals under impact velocity of 2 km· s −1 . At 12.5 ps, a plastic deformation region is formed at the right end of the crystal under the action of tension wave (marked in Figure 5a). At 13.5 and 14.2 ps, the typical spallation phenomena, that is the fracture of materials initiated by the tensile stress of RDX crystals satisfying the dynamic fracture criteria, can be observed. This is due to the fact that the hot spot strengthens the leading shock wave and reflects it to form a strong tensile zone. At 15.6 ps, the secondary shock wave generated by the hot spot catches up with the initial shock wave, and then the spallation phenomenon disappears.
(Also see the Supporting Movie SM1).
As shown in Figure 5, the spalling occurs at 330−400 Å in the X direction for U p = 2 km·s −1 . Figure 6 shows the spatial and temporal distribution of density in the above regions. There are no atoms in the spalling region, and the density is close to zero, as indicated by the dark blue in Figure 6. Once the shock wave propagates to the free surface right end (t = 8.0 ps), the reflected tensile wave propagates along the X negative direction. Under the action of tensile wave and shock wave, the spalling occurs at a specific location. In this paper, the initial position of spallation which can be calculated by the average sound velocity inside the RDX crystal is at 330 Å at 12.3 ps (see in Figure 6a). The calculated sound velocity is 2.25 km· s −1 , which is close to the 2.8 km·s −1 (for ρ = 1.80 g·cm −3 ) in the handbook 34 The difference is due to the inhomogeneity of density within the crystal. The spalling region expands continually under the dual action of shock waves and tensile waves and reaches a maximum at 14.0 ps (Figure 6d). After that, the secondary shock wave drives the atoms in the crystal to propagate along the X positive direction, and the spalling gradually disappears (Figure 6f).
To further characterize the expand and disappearance of spallation regions, the spalling region analysis was implemented using a series of PYTHON scripts. Figure 7 shows the time-dependent evolution of the projected area of region A (in Figure 6a) on the XY plane. It has been seen that the expand and decay times of spalling are almost the same in Figure 7. The relationship between the projection area in the XY plane and time is similar parabolic by the polynomial fitting method (see more details about the parabolic formula in Supporting Information S4). The rate of growth of the spalling region can be calculated about 70 Å 2 /ps, and the decay rate is 82 Å 2 /ps. Figure 8 shows the snapshots of microsplit and secondary reflection tension of RDX crystals under the impact velocity of 4 km·s −1 . For U p = 4 km·s −1 , there are many obvious microsplit regions in the inner cleavage zone of RDX crystals, which shows a loose porous state (marked with circles in Figure 8b, see the density distribution for U p = 4 km·s −1 in Supporting Information S5). Comparing with Figures 5 and 8, it can be seen that the structure of spallation zone strongly depends on

ACS Omega
Article the velocity and formation of shock wave: typical spallation occurs at low shock velocity (Figure 5b,c), while microsplit occurs at high shock velocity (Figure 8b). The interaction between the shock wave and the closure void would again reflect the formation of the unloading wave propagating along the x negative direction. The shock wave propagating along the x positive direction is continuously strengthened under the action of the hot spot, and it then reflects as a strong unloading wave after arriving at the right end free surface. The compression wave interacts with the strong unloading wave in the hot spot region, and it continuously stretches the sparse high-density region. The shock wave drives the left side of the crack to collide with the right side to catch up with the right spalling zone and make the crack disappear (Figure 8c,d).
(Also see the Supporting Movie SM2).
Because of the instability of energetic materials, the spalling phenomenon of energetic materials can hardly be observed. In recent years, research studies on spallation are mostly focused on metals. 35,36 For example, Luo et al. 35 indicate that spallation occurs without tensile melting due to the weak solid shocks of Cu; in addition, spallation may be accompanied or preceded by partial melting at the case of stronger shocks. Xiang et al. 36 revealed that classical spallation occurs where materials keep in solid state during shock compression and tensile stages, while microspallation is a spallation where melting occurs during shock compression and tensile stages, as in the case of Pb. Unlike metals, there are many chemical reactions in energetic materials by shock loading, which complicates the mechanism of spallation.
3.4. Chemical Reaction. The threshold of impact initiation of RDX crystal is 34.7 GPa. 34 The shock pressures P s , calculated using the Rankine−Hugoniot relation P s = ρ 0 U s U p , where ρ 0 is the initial density of the crystal, were 14.4 and 57.8 GPa for the impact velocities of U p = 2 and 4 km·s −1 , respectively. Obviously, the shock pressure is less than the impact initiation threshold for U p = 2 km·s −1 . Only a little RDX molecules react and concentrate around the hot spot region (quantity is too small to be listed). As seen in Figure 9, there is almost no decomposition reaction of the RDX molecules near the spalling region. On the contrary, there are large number of small molecule products for U p = 4 km·s −1 . These gaseous products, such as H 2 O, CO 2 , NO 2 , NO 3 , O 2 , and so forth, which were around the microsplit region ( Figure  10a,b) may be responsible for the occurrence of microsplit.
The strong shock wave induces the decomposition of RDX into gas products, which move in the crystal under the action of shock wave, and makes the RDX crystal appear a loose porous state.
The ruptures of N−NO 2 bond are the main initial reaction mechanisms for the shocked RDX crystal and are not affected by shock strength. NO 2 and C 3 H 6 O 4 N 5 appear earlier, and their amounts are much higher than the other products. Note that there are more NO 2 molecules than C 3 H 6 O 4 N 5 , suggesting that some RDX molecules lost two NO 2 groups, which is confirmed by the presence of C 3 H 6 O 2 N 4 ( Figure  11b). It can be seen from Figure 11a that the decomposition rate of RDX molecules can be divided into three stages. First, RDX molecules decompose rapidly under the action of shock wave before the maximum compression time. Second, the decomposition rate of RDX molecules is slow at the time of maximum compression, while the formation rate of small molecular products such as NO 2 , CO 2 , HO increases. This also confirms that the increase and mobility of small gas products are the main reasons for the formation of microcracks in the

CONCLUSIONS
In this paper, the physical and chemical phenomena of RDX crystal with a cube void are simulated under different shock velocities by using ReaxFF/lg force field. The void collapse process, hot spot formation and growth, and spallation phenomena are analyzed.
Once the initial shock wave arrives upstream surface of the void, the particles upstream of the void begin to enter the void, which lead to void collapse. The time intensity of void collapse depends on the velocity of collision. The corresponding void collapse times are 1.5 and 0.7 ps, in the case of U p = 2 km·s −1 and U p = 4 km·s −1 , respectively. The shock wave velocity U s averaged before the void collapse is calculated to be around 4.0 and 8.5 km·s −1 . The internal temperature of the void is higher than 2500 K, and the hot spot where the shape remains radially symmetrical is formed after void collapse. The shock wave initiated by hot spot propagates forward in a convex plane way. Nucleation and growth of hot spots are also important processes that induce chemical reactions in the system. In a limited simulation time in case of U p = 2 km·s −1 , the hot spot features certain amounts of attenuation with the shock wave propagation, and thus the peak temperature of the hot spot gradually decreases from the initial 3300−1900 K; while in the case of U p = 4 km·s −1 , the hot spot area is constantly expanding, and the temperature maintains at a higher temperature. The tensile wave reflected by the shock wave after reaching the free surface results in stretching of the material and spalling that depends on the initial shock velocity. Typical spalling occurs for the shock velocity 2 km·s −1 , while the tensile wave induces the microsplit region in RDX crystals in the case of U p = 4 km·s −1 . The secondary reflection of shock wave in the material further expands the material and catches up with the spallation zone, which makes the crack disappear. Shock pressures can be calculated using the Rankine− Hugoniot relation ware 14.4, 57.8 GPa for U p = 2 and 4 km· s −1 , respectively. For the weak shock, there is almost no decomposition reaction of the RDX molecules and around the spalling region are entire RDX molecules. On the contrary, there are large number of small molecule products, such as H 2 O, CO 2 , NO 2 , and so forth, around the microsplit regions for the strong shock. The formation and movement of gaseous small molecules may be the main reason for the formation of microcracks. The initial reaction mechanism of the RDX crystal which is rupture of N−NO 2 bond is independent of the impact strength, while the microsplit slows down the decomposition rate of RDX.
The results in this paper can shed light on thermal ignition, hot spot growth, and other physical and chemical phenomena of energetic materials under impact loading. Future researches will focus on the effects of void shape, size, and void interactions on hot spots and chemical reactions to better understand the relationship between crystal properties and the impact response of energetic materials.

ACS Omega
Article