View Article Online View Journal

PCCP Accepted Manuscript

This article can be cited before page numbers have been issued, to do this please use: T. Zhou, H. Song, Y. Liu and F. Huang, Phys. Chem. Chem. Phys., 2014, DOI: 10.1039/C4CP00890A.

This is an Accepted Manuscript, which has been through the Royal Society of Chemistry peer review process and has been accepted for publication. Accepted Manuscripts are published online shortly after acceptance, before technical editing, formatting and proof reading. Using this free service, authors can make their results available to the community, in citable form, before we publish the edited article. We will replace this Accepted Manuscript with the edited and formatted Advance Article as soon as it is available. You can find more information about Accepted Manuscripts in the Information for Authors. Please note that technical editing may introduce minor changes to the text and/or graphics, which may alter content. The journal’s standard Terms & Conditions and the Ethical guidelines still apply. In no event shall the Royal Society of Chemistry be held responsible for any errors or omissions in this Accepted Manuscript or any consequences arising from the use of any information it contains.

www.rsc.org/pccp

Page 1 of 39

Physical Chemistry Chemical Physics View Article Online

Thermal Shock Initiated Thermal and Chemical Responses of HMX

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

Crystal from ReaxFF Molecular Dynamics Simulation Tingting Zhou,a,b Huajie Song,a Yi Liu,c* and Fenglei Huangb Abstract: To gain an atomistic-level understanding of the thermal and chemical responses of condensed energetic materials under thermal shock, we developed a thermal shock reactive dynamics (TS-RD) computational protocol using molecular dynamics simulation coupled with ReaxFF force field. βoctahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocane (HMX) was selected as a target explosive due to its widely usage in military and industry. The results show that thermal shock initiated by large temperature gradient between “hot” region and “cold” region results in thermal expansion of particles and induces a thermal-mechanical wave propagating back and forth in the system with an averaged velocity of 3.32 km·s-1. Heat propagating along the direction of thermal shock leads to temperature increment of the system and thus chemical reaction initiation. Applying a continuum reactive heat conduction model combined with the temperature distribution obtained from the RD simulation, heat conduction coefficient is derived to be 0.80 W·m-1·K-1. Chemical reaction mechanisms during thermal shock were analyzed, showing that reaction is triggered by N-NO2 bond breaking followed by HONO elimination and ring fission. The propagation rates of reaction front and reaction center are obtained to be 0.069 km·s-1 and 0.038 km·s-1 based on the time and spatial distribution of NO2. The pressure effect on thermal shock was also investigated by employing uniaxial compression before thermal shock. We find that compression significantly accelerates thermal-mechanical wave propagation and heat conduction, resulting in higher temperature and more excited molecules and thus earlier initiation and faster propagation of chemical reactions. Keywords: thermal shock; heat conduction, chemical reaction; energetic materials; HMX; molecular 1

Physical Chemistry Chemical Physics Accepted Manuscript

DOI: 10.1039/C4CP00890A

Physical Chemistry Chemical Physics

Page 2 of 39 View Article Online

DOI: 10.1039/C4CP00890A

1. Introduction Condensed energetic materials (EMs), such as HMX, RDX, and TATB, are widely used in military and

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

civilian fields. During the different stages of preparation, manufacture, package operation, transportation, storage, and usage, EMs are subjected to thermal and mechanical loadings, resulting in microstructure changes that can lead to damage, initiation, or even accidental explosion. Understanding the mechanical damage and chemical initiation mechanisms are essential to ascertain the safety of these materials that is critical to reduce the loss of lives and property. As a kind of thermal loading, thermal shock causes thermal expansion and thermal stress in materials that contribute to the formation and growth of micro-crack, even fracture if the temperature changes sharply.1-3 These phenomena are very common especially in brittle materials such as the commonly used high energy compounds. Additionally, temperature gradient between the surface and the inner of explosive induces heat conduction, leading to temperature increment of the material that causes expansion or deformation of the molecular structure of the explosive, weakening the bonds that hold the molecules in formation, which benefits bond breaking.1 The temperature effects on mechanical properties, like modulus, strength, and compressibility, have been widely investigated by experiments and numerical modeling based on finite element method and important conclusions were obtained.1-8 However, the ignition mechanism of EMs under thermal shock is more challenging and is far from being well understood due to the complex and multiscale process involving thermal expansion, endothermic phase change, molecular decomposition, etc. Peterson et al.1 conducted a series of experiments in which a carefully controlled temperature gradient has been applied along a cylinder of PBX 9501 to explore the details of this intricate process. Their results clearly show that thermal damage in PBX 9501 can be classified into two separate temperature regimes—an initial low-temperature regime (155–174°C) dominated by the endothermic β-δ crystalline phase change, thermal expansion, and Ostwald ripening, and a high-temperature regime (175–210°C) dominated by exothermic chemical

2

Physical Chemistry Chemical Physics Accepted Manuscript

dynamics; ReaxFF

Page 3 of 39

Physical Chemistry Chemical Physics View Article Online

DOI: 10.1039/C4CP00890A

morphology and the chemical reactions leading to a potential thermal self-ignition in the explosive. Thermally damaged high explosives can behave much differently than pristine high explosives when

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

subjected to thermal or mechanical stimuli. After thermal shocking, PBX 9501 displays increased shock sensitivity, reduced threshold pressures for violent reaction, and reduced time to self-ignition at elevated temperatures, compared to pristine PBX 9501.9-11 The increased sensitivity may be related to the thermally induced morphological and chemical changes within the explosive.11 To understand the thermally induced microstructure changes, heat conduction, and particularly chemical reactions of EMs at atomistic level, we built a thermal shock model and employed molecular dynamics coupled with reactive force field (ReaxFF) to simulate the thermal shock process. Additionally, we applied uniaxial compression before thermal shock to look into the effect of pressure on these thermal-chemical responses. β-octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocane (HMX) was chosen due to its important and widely usage in powders, propellants, and composite. Molecular dynamics have received increasing acceptance in recent years as a means of obtaining information concerning some of the fundamental properties and processes indicated above. Given an accurate force field and careful consideration of the suitability of the approximation of classical statistical mechanics for predicting some quantity of interest atomistic simulations can yield accurate temperature- and pressure-dependent predictions of EMs.12 The rest of this paper is organized as follows: method and computational procedures are described in Sec. 2, followed by results and discussions in Sec. 3; conclusions are drawn in Sec. 4. 2. Method and Computational Procedures 2.1 ReaxFF reactive force field ReaxFF reactive force field, a bond-order dependent potential parameterized from quantum mechanical (QM) calculations,13 provides accurate descriptions of bond breaking and bond forming. The fundamental difference between ReaxFF and other force fields is that ReaxFF does not need any

3

Physical Chemistry Chemical Physics Accepted Manuscript

decomposition. The results further demonstrate the complex interplay between the evolving sample

Physical Chemistry Chemical Physics

Page 4 of 39 View Article Online

DOI: 10.1039/C4CP00890A

instantaneous interatomic distances which are updated continuously. The electrostatic and van der Waals interactions are also calculated between every atom pairs during the simulation. ReaxFF has been

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

demonstrated to be capable of describing the structures, energy surfaces, and reaction barriers of various materials, including hydrocarbon,14,

15

polymer,16 nitramines,17-28 metals,29-31 and metal oxides.32-34

Combined with molecular dynamics (MD), it provides a direct dynamic approach to study large-scale and complex chemical reactions at the atomistic level with a higher efficiency than QM, considering temperature and pressure effects additionally. It has been successfully applied to study the complicated dynamical processes for EMs including RDX, PETN, TATB, and HMX subjected to various stimuli,1726

providing value information about physical and chemical properties of these organic systems. ReaxFF-lg,35 an extension of ReaxFF, is developed to more accurately account for London dispersion

(van der Waals attraction), which leads to a more accurate description for molecular crystal structure. The improved prediction of density is critical for EMs as density is correlated with detonation performance of these materials. It has been tested for several energetic materials, including RDX, PETN, TATB, nitromethane, as well as HMX,26, 35 showing good agreements with experimental results. The performance of predicting chemical reactions for ReaxFF-lg is as good as that for ReaxFF.35 All simulations in this work were performed using ReaxFF-lg potential. 2.2 Thermal shock reactive dynamics (TS-RD) protocol The thermal shock reactive dynamics (TS-RD) computational protocol is proposed to model the micro thermal shock dynamical process in HMX crystal. The TS-RD procedures are as follows. 1. Structure optimization. A small 5 × 2 × 4 supercell was built and then optimized by running energy minimization. Atom coordinates and cell parameters were then relaxed using isobaric-isothermic molecular dynamics (NPT-MD) simulation for 60 ps to gain equilibrated structure at 300 K and at zero pressure. The equilibrated lattice parameters of β-HMX are: V0 = 527.39 Å3, a = 6.57 Å, b = 11.11 Å, c = 8.76 Å, and β = 124.41°, which are in good agreements with experimental results: V0 = 519.39 Å3, a = 4

Physical Chemistry Chemical Physics Accepted Manuscript

connectivity input for the chemical bonds, but instead calculates the bond order directly from

Page 5 of 39

Physical Chemistry Chemical Physics View Article Online

36

DOI: 10.1039/C4CP00890A

large supercell, containing 560 molecules and 15680 atoms. MD simulations in canonical ensemble (NVT) at T = 300 K was next performed for 2 ps to optimize the large system.

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

2. Uniaxial compression. The equilibrated large structure was uniaxially compressed by 0, 15 %, and 30 % of the original volume (ΔV = 0, 15 %, and 30 %) along y direction. The dimensions of these compressed boxes were changed to 32.85 Å × 155.54 Å × 35.04 Å, 32.85 Å × 132.21 Å × 35.04 Å, and 32.85 Å × 108.88 Å × 35.04 Å. These compressed systems were then relaxed with cell parameters fixed using NVT-MD at T = 300 K for 10 ps. The obtained stress tensors were collected in Table S1 of the Electronic Supplementary Information (ESI). 3. Temperature gradient introduction. Two “hot” regions were defined at the two ends of each simulation box in y direction with the dimensions of 32.85 Å × 11.11 Å × 35.04 Å, 32.85 Å × 9.44 Å × 35.04 Å, and 32.85 Å × 7.77 Å × 35.04 Å for systems with 0, 15 %, and 30 % compression, respectively. The formation of “hot” regions is mimicked by fast heating the two ends from 300 K to 3000 K in 1 picosecond while the temperature of the center region (“cold” region) keeps 300 K during this period, forming large temperature gradient between the ends and the middle. The relatively large size and high temperature of “hot” regions used in these simulations is to quickly observe thermal and chemical responses in several picoseconds using affordable computing resources. 4. Thermal shock reactive dynamics. One thermostat was imposed to each “hot” region and NVT-RD simulations were carried out at T = 3000 K. The temperature was scaled using Berendsen method with a coupling constant of 100 fs. Meantime, reactive dynamics simulations in microcanonical ensemble (NVE-RD) without temperature control were performed for the “cold” region to follow the time evolution of the system. The particles near the interface of the two different regions will be subjected to NVE ensemble (or NVT ensemble) once they move to “cold” region (or “hot” region) from “hot” region (or “cold” region), thus, their temperature will be altered. The large temperature gradient drives heat flow from “hot” region to “cold” region as thermal shock during the TS-RD simulation. Fig. 1 shows

5

Physical Chemistry Chemical Physics Accepted Manuscript

6.54 Å, b = 11.05 Å, c = 8.70 Å, and β = 124.30°. The equilibrated system was expanded to 5 × 14 ×4

Physical Chemistry Chemical Physics

Page 6 of 39 View Article Online

DOI: 10.1039/C4CP00890A

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

compression, respectively.

Figure 1. The TS-RD model for HMX single crystal. Different regions are separated by the two vertical black lines

All MD simulations were performed using General Reactive Atomistic Simulation Package (GRASP). A time step of 0.1 fs was used in the integration of the Newtonian equations of motion and threedimension periodic boundary conditions were employed. We performed one-dimensionally binning to analyze the spatial distributions of physical properties (such as density, temperature, and others) during thermal shock. To analyze the numerous reactions that occur during ReaxFF simulations, we carried out fragment analysis. The algorithm of molecule recognition in the fragment analysis uses the connection table and bond orders calculated by the ReaxFF every step22, 25 and records the molecular components and their compositions every 0.1 ps. In this analysis, the bond order cutoff values used to identify molecular fragments in the systems are tabulated in Table S2 of the ESI for various atom pairs. These cutoff values were verified in previous study.26 Thus, any two fragments were considered separate molecules if the bonds between them have bond orders smaller than cutoff values. These data were then used to plot time evolution profiles for reactants, intermediates, and products, providing detailed information about the initiation of chemical reactions. After determining the molecular fragments based on bond order cutoff values, the molecule recognition algorithm assigns a unique identification number to each molecular fragment to trace the reaction pathways.

6

Physical Chemistry Chemical Physics Accepted Manuscript

the TS-RD model. The simulation time was 80 ps, 60 ps, and 30 ps for systems with 0, 15 %, and 30 %

Page 7 of 39

Physical Chemistry Chemical Physics View Article Online

3.1 Thermal-Mechanical Wave Propagation Due to the large temperature gradient between “hot” region and “cold” region, particles near the “hot”

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

region are found to expand from the ends to the center, and induce a mechanical wave. Herein, we denoted such mechanical wave as thermal-mechanical wave (TMW), differing from external mechanical wave. TMWs are initiated from the two ends and move towards each other until they meet at the center of the system at t = 2 ps, resulting in the compression of the center region; the two waves then reflect backwards during the next 2 picoseconds, leading to the expansion of this region. Such compression and expansion occur alternatively and repeat several times until TMWs decay off at about 20 ps. This phenomenon is reflected in the time evolution of microstructure for HMX crystal as shown in Fig. S1 of the ESI. Time-spatial distribution of mass density illustrated in Fig. 2 reveals directly the propagation process of TMWs. We divided the system along y direction into 28 bins equally and the bins are of fixedvolume. Particles flow into or out of the bins during thermal shock, leading to the alteration of density within each bin. The density of “cold” region increases from the two sides to the center with time increasing at t < 2 ps, while the density of “hot” regions decreases gradually during this period (shown in Fig. 2(a)). This verifies that the TMWs initiated by thermal shock propagate from the two ends towards the center of the system within the first 2 ps. During the next 2 ps, the density of “cold” region goes down from ~ 2.15 g·cm-3 to ~ 1.90 g·cm-3, whereas, it goes up gradually for “hot” regions (shown in Fig. 2(b)), demonstrating that the two waves reflect backwards to the ends after they collide at the center during this period. The increase-decrease cycle of mass density repeats several times with a periodicity of 4 ps during the first 20 ps (see Fig. 2(c)), indicating that TMWs oscillate with a periodicity of 4 ps during this period, after which they decay off. Here, we designated the periodicity and time of decay off as tp and td. Note that such periodicity depends on the finite system size in this work, but the observed reflective wave is common in experiments when forward mechanical wave

7

Physical Chemistry Chemical Physics Accepted Manuscript

3 Results and Discussions

DOI: 10.1039/C4CP00890A

Physical Chemistry Chemical Physics

Page 8 of 39 View Article Online

DOI: 10.1039/C4CP00890A

Vs is calculated to be around 3.32 km·s-1, close to the experimental sound velocity of 3.0 km·s-1 for HMX crystal.37 At t = 20 ps, the density of “cold” region is higher than the initial value of 1.90 g·cm-3,

continuously from the two sides to the center, which means that low-density regions expand from the sides to the center (see Fig. 2(d)). These low-density regions are places where chemical reaction occurs, making the microstructure of HMX crystal become amorphous as shown in Fig. S1 of the ESI. (a) 2.4

(c) 2.4

Compression = 0

2.2

Compression = 0 2.2

2.0 -3

Density/gcm

Density/gcm

-3

2.0

1.8 1.6 1.4

t=0 t = 0.5 ps t = 1 ps t = 2 ps

1.2 1.0 0.8

0

20

40

60

t = 0.25 ps t = 0.75 ps t = 1.5 ps

80

100

120

1.8 1.6 1.4 t = 4 ps t = 8 ps t = 12 ps t = 20 ps

1.2 1.0

140

0.8

160

0

20

40

Distance/Å

80

100

120

140

160

(d) 2.4

Compression = 0

Compression = 0 2.2

2.0

2.0 -3

2.2

Density/gcm

-3

60

t = 6 ps t = 10 ps t = 15 ps

Distance/Å

(b) 2.4

Density/gcm

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

suggesting that the total effect of TMWs is compressing. After 20 ps, the density of “cold” region drops

1.8 1.6 1.4 1.2

t = 2 ps t = 3 ps t = 4 ps

1.0 0.8

0

20

40

60

t = 2.5 ps t = 3.5 ps

1.8 1.6 1.4

t = 20 ps t = 40 ps t = 60 ps t = 80 ps

1.2 1.0

80

100

Distance/Å

120

140

160

0.8

0

20

40

60

80

t = 30 ps t = 50 ps t = 70 ps

100

120

140

160

Distance/Å

Figure 2. Time and spatial distribution of mass density under thermal shock for HMX crystal without compression.

The time evolution of particle velocity Vp averaged inside a small segment was calculated to describe the TMW propagation quantitatively (see Fig. 3). This segment is 33.33 Å away from the left “hot” region with the size of 11.11 Å in the y direction as shown in the inset of Fig. 3. The time evolution of Vp shows that the averaged particle velocity fluctuates around zero with a periodicity of about 4 ps,

8

Physical Chemistry Chemical Physics Accepted Manuscript

encounters obstacles or other waves. Based on tp and propagation distance, the averaged TMW velocity

Page 9 of 39

Physical Chemistry Chemical Physics View Article Online

DOI: 10.1039/C4CP00890A

Such behavior continues several times within the first 20 ps, which is further confirmed by the possibility distribution of particle velocities within this bin as depicted in Fig. S2 of the ESI. It is

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

noticed that Vp does not have equal positive and negative amplitudes: the positive velocity is larger and is along the same direction of thermal shock, while the negative one is smaller and is in the opposite direction. Hence, the overall displacements of particles are along the direction of heat propagation leading to the final compression of the system, which is in accordance with the increased mass density of “cold” region. The amplitude of Vp decays rapidly from the maximum 0.5 km·s-1 to 0.1 km·s-1 within the first 20 ps due to energy dissipation, after which particles move irregularly at the velocity within ± 0.1 km·s-1. Such a small particle velocity induced by thermal shock is too low to initiate any chemical reaction. It is previously found that during the shock simulation of RDX, there is little chemistry, and only a few NO2 are formed for impact velocity smaller than 4 km·s-1, some of which are just highly excited N-N bonds that have a chance of appearing as broken in phase space.18

Figure 3. Time evolution of the averaged particle velocity Vp for HMX crystal without compression during the first 40 ps TS-RD simulation. The velocity is averaged over the particles within the bin defined between the two vertical lines.

The propagation rate of TMW is markedly accelerated by uniaxial compression, which can be demonstrated by the time and spatial distributions of mass density for compressed systems as shown in

9

Physical Chemistry Chemical Physics Accepted Manuscript

indicating that particles move positively in the first 2 ps and move negatively in the following 2 ps.

Physical Chemistry Chemical Physics

Page 10 of 39 View Article Online

DOI: 10.1039/C4CP00890A

the compressed systems during the propagation of TMW is the same as that for the uncompressed system but the periodicity is remarkably reduced. The periodicity tp and the decay off time td of TMW

ps and 5 ps for the highly compressed system (ΔV = 30 %). The TMW velocities Vs for the 15 % and 30% compressed systems are calculated to be about 5.67 km·s-1 and 9.33 km·s-1, much higher than that for the system without compression. (b1)3.0

(a1) 2.5 Compression = 15 %

2.4

Compression = 30 %

2.9

Density/gcm

Density/gcm

-3

-3

2.3 2.2 2.1 2.0 1.9 t=0 t = 0.4 ps t = 0.8 ps

1.8 1.7 1.6

0

20

40

60

t = 0.2 ps t = 0.6 ps t = 1 ps

80

100

2.8 2.7 2.6 t=0 t = 0.2 ps t = 0.4 ps

2.5

120

140

2.4

0

20

40

Distance/Å

t = 0.1 ps t = 0.3 ps t = 0.5 ps

60

80

100

120

Distance/Å (b2) 3.0

(a2) 2.5 Compression = 15 %

2.4

Compression = 30 %

2.9

Density/gcm

-3

-3

2.3

Density/gcm

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

propagation are estimated to be 2 ps and 14 ps for the system with low compression (ΔV = 15 %) and 1

2.2 2.1 2.0 1.9 t = 1 ps t = 1.4 ps t = 1.8 ps

1.8 1.7 1.6

0

20

40

60

t = 1.2 ps t = 1.6 ps t = 2 ps

80

Distance/Å

100

2.8 2.7 2.6 t = 0.5 ps t = 0.7 ps t = 0.9 ps

2.5

120

140

2.4

0

20

40

60

t = 0.6 ps t = 0.8 ps t = 1 ps

80

100

120

Distance/Å

10

Physical Chemistry Chemical Physics Accepted Manuscript

Fig. 4. Periodic increase-decrease behavior of the density revealing the compression-expansion cycle of

Page 11 of 39

Physical Chemistry Chemical Physics View Article Online

(b3)3.0

(a3)2.5

Compression = 30 %

Compression = 15 %

2.4

2.9

Density/gcm

2.2 2.1 2.0 1.9 t = 2 ps t = 4 ps t = 10 ps

1.8 1.7 1.6

0

20

40

60

t = 3 ps t = 7 ps t = 14 ps

80

100

2.8 2.7 2.6 t = 1 ps t = 3 ps t = 5 ps

2.5

120

2.4

140

0

20

40

60

80

100

120

(b4) 3.0 Compression = 15 %

2.4

Compression = 30 %

2.9

Density/gcm

-3

-3

2.3 2.2 2.1 2.0 1.9 t = 14 ps t = 30 ps t = 50 ps

1.8 1.7 1.6

t = 2 ps t = 4 ps t = 6 ps

Distance/Å

Distance/Å (a4) 2.5

Density/gcm

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

Density/gcm

-3

-3

2.3

0

20

40

60

t = 20 ps t = 40 ps t = 60 ps

80

100

2.8 2.7 2.6 t = 6 ps t = 15 ps t = 25 ps

2.5

120

140

2.4

0

20

40

t = 10 ps t = 20 ps t = 30 ps

60

80

100

120

Distance/Å

Distance/Å

Figure 4. Time and spatial distributions of mass density under thermal shock for compressed HMX crystals: (a1-a4) ΔV = 15 %, (b1-b4) ΔV = 15 %.

As we mentioned above, continuously decrease of densities from the two ends to the center beginning after TMW decaying off indicates the occurrence of chemical reaction. The smaller td for compressed systems suggests earlier reaction initiation compared with uncompressed system. Also, the faster decrease of densities for center region suggests that the decomposition process propagates more rapidly for compressed systems. These are mainly arisen from faster temperature increasing for compressed systems, which will be confirmed in Section 3.2. It is noticed that for the highly compressed system, the densities of the regions near the two ends become higher than that for center region after 5 ps, which is distinct from those for low compressed systems. We propose that high uniaxial compression strongly accelerates heat propagation, leading to rapid temperature increasing and chemical initiation, followed by dramatic exothermic chemical reactions that generate plenty of gas products. These products

11

Physical Chemistry Chemical Physics Accepted Manuscript

DOI: 10.1039/C4CP00890A

Physical Chemistry Chemical Physics

Page 12 of 39 View Article Online

DOI: 10.1039/C4CP00890A

compression of the ends. The behavior of the averaged particle velocity Vp for compressed systems is similar to that for the

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

system without compression but with shorter oscillation periodicity, further confirming the accelerated propagation of TMW. Fig. 5 presents Vp ~ t relationship, illustrating that the oscillation periodicity is about 2 ps for the 15 % compressed system. The velocity amplitude attenuates rapidly from the maximum value 0.26 km·s-1 at t = 0 to 0.05 km·s-1 at t = 14 ps, after which typical oscillation behavior disappears and the velocity changes randomly within the range of ± 0.1 km·s-1. For the highly compressed system, it is hard to observe explicit oscillation behavior because of the very short periodicity. The maximum Vp is about 0.16 km·s-1 at the beginning and particles move intricately at the velocities within the range of ± 0.1 km·s-1 after 5 ps. The oscillation periodicity and the maximum amplitude of the averaged particle velocity for the compressed systems are about less than half of the uncompressed system. As is well known, external compression decreases intermolecular distance, resulting in a smaller free space available for the movements of particles. The reduced free space enhances the probability of collision among particles and thus the movement directions could be altered more frequently, leading to the reduction of oscillation periodicity. This would definitely benefit mass as well as heat conduction, resulting in faster temperature increment.

12

Physical Chemistry Chemical Physics Accepted Manuscript

propagate from center region towards two ends, resulting in the expansion of the center and the

Page 13 of 39

Physical Chemistry Chemical Physics View Article Online

Figure 5. Time evolutions of the averaged particle velocity Vp for compressed HMX crystals during the first 30 ps TSRD simulation: (a) ΔV = 15 %, (b) ΔV = 15 %.

3.2 Heat Conduction The temperature time-spatial distribution during the dynamical process is plotted in Fig. 6 (a), showing that the temperature is high at the ends and low at the center. This large temperature gradient drives heat flow from the ends to the center, heating up “cold” molecules and leading to temperature increment. The temperature of the center increases gradually with the time increasing and reaches 1050 K at t = 80 ps.

13

Physical Chemistry Chemical Physics Accepted Manuscript

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

DOI: 10.1039/C4CP00890A

Physical Chemistry Chemical Physics

Page 14 of 39 View Article Online

Figure 6. Temperature time-spatial distribution diagram for HMX crystal without compression: (a) extracted from TSRD simulation, the two horizontal dashed lines from up to down present the temperatures of reaction center (T = 1780 K) and reaction front (T = 1080 K), (b) calculated from reactive heat conduction model with heat conductivity coefficient of 0.80 W·m-1·K-1. The numbers on curves represent time (ps) elapsed during simulation.

A one-dimensional continuum heat conduction model with heat sources related to chemical reaction was set up to describe heat wave propagation as follows:

T  k  2T  Q t (1) Ei Q   qi Ai exp( )(i  1, 2,3) RT

c

where T is temperature, t time,  density, c specific heat capacity,  heat conduction coefficient, qi heat of reaction, Ai frequency factor, Ei activation energy, and R gas constant. 14

Physical Chemistry Chemical Physics Accepted Manuscript

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

DOI: 10.1039/C4CP00890A

Page 15 of 39

Physical Chemistry Chemical Physics View Article Online

DOI: 10.1039/C4CP00890A

experimental results from Tarver et al.,38 who suggested that the reaction mechanism of solid HMX decomposition is HMX → fragments → intermediate gases → final gases. The following thermal and

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

chemical kinetic parameters were used in the three-step reactive heat conduction (RHC) model:

= 1.9 g·cm-3

c = 2.3012 J·K-1·g-1

 = 0.80 W·m-1·K-1

q1 = 418.4 J·g-1

E1 = 52.7 kcal·mol-1

lnA1 = 48.7 s-1

q2 = -1522.2 J·g-1

E2 = 44.1 kcal·mol-1

lnA2 = 37.3 s-1

q3 = -5020.8 J·g-1

E3 = 34.1 kcal·mol-1

lnA3 = 28.1 s-1

where , c, qi, Ei, and Ai were derived from Tarver et al.38 Finite difference method was employed to solve Eq. (1) to predict the time evolution of temperature. Requiring that the RHC model leads to temperature profiles close to those observed in the TS-RD simulation, the heat conduction coefficient used in the RHC model can be determined. When  is set to be 0.80 W·m-1·K-1, the predicted temperature of the center rises to ~ 1050 K at 80 ps and the temperature time-spatial distribution is very close to those from TS-RD simulation as shown in Fig. 6(b). The calculated  is 1.59 times as high as the one at ambient condition, 0.502 W·m-1·K-1.23 The heat conduction coefficient derived from the RHC model can be considered as an averaged value of the mixing reacted and unreacted HMX molecules. The decomposed small products facilitate heat and mass conductions via move and collision in a finite space, leading to higher thermal conductivity coefficient. Fig. 7 shows the temperature spatial distributions at various times for compressed systems, revealing that uniaxial compression significantly speeds up the heat flow. The temperature of the center reaches 1550 K at 60 ps for 15 % compressed system and rises to 1780 K at 30 ps for the 30 % compressed one, which are significantly higher than that for the system without compression. The heat conduction coefficients for the compressed systems are derived to be 1.25 W·m-1·K-1 for the 15 % compressed system and 1.82 W·m-1·K-1 for the 30 % compressed system, which are 1.56 times and 2.28 times as large as the uncompressed system. The enhanced heat conduction coefficient further confirms that heat 15

Physical Chemistry Chemical Physics Accepted Manuscript

Three-step reaction for the HMX decomposition process was taken into account according to the

Physical Chemistry Chemical Physics

Page 16 of 39 View Article Online

DOI: 10.1039/C4CP00890A

Figure 7. Temperature time-spatial distribution diagrams for compressed HMX crystals under thermal shock: (a) with 15 % compression, the two horizontal dashed lines from up to down represent the temperatures of reaction center (T = 2100 K) and reaction front (T = 1550 K), (b) with 30 % compression, the two horizontal dashed lines from up to down represent the temperatures of reaction center (T = 2280 K) and reaction front (T = 1780 K). The numbers on curves represent the time (ps) elapsed during simulation.

The accelerated heat propagation and higher temperature under uniaxial compression are partially arisen from the reduced intermolecular distance under compression, leading to higher collision probability among particles and shorter oscillation periodicity that are in favor of heat conduction as discussed in section 3.1. Additionally, external compression would lead to the increment of internal

16

Physical Chemistry Chemical Physics Accepted Manuscript

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

propagation is accelerated by uniaxial compression.

Page 17 of 39

Physical Chemistry Chemical Physics View Article Online

DOI: 10.1039/C4CP00890A

channeled into appropriate molecular vibrational modes, the result may be bond scission and the initiation of chemical reaction followed by exothermic reactions. From the “multi-phonon up-pumping”

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

viewpoint,39 a shock wave produces a bath of highly excited phonons, which flow into the lowest vibrational modes (doorway modes) of the molecules that make up a crystal. Continued flow of energy from the hot phonon bath into the doorway modes leads to higher levels of vibrational excitation. The internal vibrational molecular modes equilibrate very quickly, and the internal temperature of the molecule rises. If the initial phonon temperature is sufficiently high, multi-phonon vibrational uppumping will heat the molecules to temperatures at which chemical reaction initiation. Uniaxial compression enhances the levels of molecular vibration, making these atoms more active. When they supposed to thermal shock that provides initial phonons with high temperature, the multi-phonon vibrational up-pumping will be excited more easily, leading to higher temperature of molecules and earlier bond breaking followed by faster exothermic reactions that in turn contribute to temperature increase. 3.3 Chemical Reaction Initiation and Propagation Heat conduction leads to the temperature increase of “cold” molecules, benefiting chemical reaction initiation. Analyzing chemical species produced during the dynamical process assists in understanding the evolution of chemical reactions. Fig. 8 shows the time evolutions of major chemical products with a large quantity, which suggests that N-NO2 bond breaking, HONO elimination, and concerted ring fission are three predominant initial reaction mechanisms under thermal shock. These are in agreements with recent QM studies on isolated HMX molecule as well as condensed phase40,

41

and previous

experiments reports.42-45

17

Physical Chemistry Chemical Physics Accepted Manuscript

energy via energy conversion from mechanical work done on the system. If enough of this energy is

Physical Chemistry Chemical Physics

Page 18 of 39 View Article Online

DOI: 10.1039/C4CP00890A

(a1)600

Compression = 0

HMX

400 300

NO2 HONO NO

200

HO

HNO

CO2

CO

N2

H2O

100

0

5

10

15

20

25

Number of products

Number of products

500

0

80

CHN 60 CH2N CHNO 40

CH2O

CNO

O2

0

5

10

15

20

25

30

(b2)100 Compression = 15 %

500

Number of products

Number of products

H2

O

Compression = 15 % HMX

400 N2

H2O

300 NO2 200

NO

HONO

HO

HNO CO2 CO

100

0

5

10

15

20

25

80

60

HNO3 H 2

O NO3

CHNO

40

O2

CNO CH O 2

CHN 20

0

30

H

CH2N 0

5

10

15

20

25

30

t / ps

t / ps (c2) 300

(c1) 600

Compression = 30 %

500

N2

HMX

H2O

400 300

Number of products

Compression = 30 %

HO NO2

200

HONO HNO

NO

CO2

100 0

HNO3

t / ps

(b1) 600

0

H

20

0

30

NO3

t / ps

Number of products

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

Compression = 0

CO 0

5

10

15

t / ps

20

25

30

250

H

200 O

150 NO3 100

CH2O H2 HNO3 CHN

CHNO CNO O2

CH2N

50 0

0

5

10

15

20

25

30

t / ps

Figure 8. Time evolutions of major chemical reaction products for HMX crystals with different compressions during the first 30 ps TS-RD simulation. (a1-a2) ΔV = 0; (b1-b2) ΔV = 15 %; (c1-c2) ΔV = 30 %.

N-NO2 bond rupture is the most primary initial reaction mechanism compared with the other two, leading to a larger number of NO2 in comparison with HONO and CH2N2O2, etc. The number of NO2 increases dramatically at the beginning and then decreases gradually after reaching the maxima due to subsequent reactions, leading to the formation of NO3, HNO3, etc. As the products of initial N-NO2 bond breaking from HMX molecule (C4H8O8N8 → C4H8N7O6 + NO2), the time and spatial distributions 18

Physical Chemistry Chemical Physics Accepted Manuscript

(a2)100

Page 19 of 39

Physical Chemistry Chemical Physics View Article Online

DOI: 10.1039/C4CP00890A

equally and count up the number of specific product within each bin. Clearly that the population of NO2 is much larger than that of C4H8N7O6, indicating that multiple successive N-NO2 bonds break after

which may be due to higher energy barrier for HONO elimination (52.3 kcal/mol) than that for NO2 dissociation (47.9 kcal/mol).41 The spatial distribution of HONO at various times shown in Fig. 10 further confirms that the formation of HONO is less favorable than NO2. At 80 ps, NO2 molecules were observed at the regions close to the center, but no HONO appears at these regions. 70

70

t = 5 ps Number of C4H8N7O6

Number of NO2

60 50 40 30 20 10 0

1

2

3

4

5

6

7

8

9

50 40 30 20 10 0

10 11 12 13 14

t = 5 ps

60

1

2

3

4

70

Number of NO2

Reaction center

40 30

Reaction front

20 10

Number of C4H8N7O6

t = 20 ps

50

7

8

9

10 11 12 13 14

t = 20 ps

60 50 40 30 20 10

1

2

3

4

5

6

7

8

9

0

10 11 12 13 14

1

2

3

4

5

6

7

8

9

10 11 12 13 14

Distance (x 11.11 Å)

Distance (x 11.11 Å) 70

70

Number of C4H8N7O6

t = 40 ps

60 50 40 30 20 10 0

6

70

60

0

5

Distance (x 11.11 Å)

Distance (x 11.11 Å)

Number of NO2

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

losing the first NO2. The formation of HONO is slower and the amount is smaller compared with NO2,

1

2

3

4

5

6

7

8

9

10 11 12 13 14

Distance (x 11.11 Å)

t = 40 ps

60 50 40 30 20 10 0

1

2

3

4

5

6

7

8

9

10 11 12 13 14

Distance (x 11.11 Å)

19

Physical Chemistry Chemical Physics Accepted Manuscript

of NO2 and C4H8N7O6 are shown in Fig. 9. We separated the system along y direction into 14 bins

Physical Chemistry Chemical Physics

Page 20 of 39 View Article Online

DOI: 10.1039/C4CP00890A 70

70

Number of C4H8N7O6

Number of NO2

40 30 20 10 0

1

2

3

4

5

6

7

8

9

50 40 30 20 10 0

10 11 12 13 14

1

2

3

4

Distance (x 11.11 Å) t = 80 ps

Number of C4H8N7O6

Number of NO2

6

7

8

9

10 11 12 13 14

70

60 50 40 30 20 10 0

5

Distance (x 11.11 Å)

70

1

2

3

4

5

6

7

8

9

50 40 30 20 10 0

10 11 12 13 14

t = 80 ps

60

1

2

3

4

5

6

7

8

9

10

11

12

13

14

Distance (x 11.11 Å)

Distance (x 11.11 Å)

Figure 9. Spatial distributions of NO2 (left) and C4H8N7O6 (right) at various times under thermal shock for HMX crystal without compression. Reaction fronts are defined at the front boundaries of NO2 distribution and reaction centers at the boundaries with the maximum value of NO2 distribution. 20

20

t = 20 ps Number of HONO

t = 5 ps Number of HONO

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

50

t = 60 ps

60

Physical Chemistry Chemical Physics Accepted Manuscript

t = 60 ps

60

15

10

5

0

1

2

3

4

5

6

7

8

9

10 11 12 13 14

Distance (x 11.11 Å)

15

10

5

0

1

2

3

4

5

6

7

8

9

10 11 12 13 14

Distance (x 11.11 Å)

20

Page 21 of 39

Physical Chemistry Chemical Physics View Article Online

DOI: 10.1039/C4CP00890A

20

Number of HONO

Number of HONO

15

10

5

0

1

2

3

4

5

6

7

8

9

15

10

5

0

10 11 12 13 14

1

2

Distance (x 11.11 Å)

3

4

5

6

7

8

9

10 11 12 13 14

Distance (x 11.11 Å)

20

t = 80 ps Number of HONO

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

t = 60 ps

15

10

5

0

1

2

3

4

5

6

7

8

9

10 11 12 13 14

Distance (x 11.11 Å)

Figure 10. Spatial distribution of HONO at various times under thermal shock for HMX crystal without compression

CH2N2O2 is one fourth of HMX molecule, which can be derived from the following dissociation mechanisms: 1) C4H8N8O8 → 2C2H4N2O2, which decomposes to two CH2N2O2; 2) C4H8N8O8 → CH2N2O2 + open RDX; 3) C4H8N8O8 → 4CH2N2O2. Sharia et al.41 predicted the activation barriers for the three reactions in gas phase to be 97.23, 70.19, and 82.98 kcal/mol, respectively, all of which are higher than that for N-NO2 bond cleavage and HONO elimination. They did not consider ring breaking reactions for condensed phase since they believe that the unfavorable mechanisms even in the gas phase arisen from high activation barriers and low reaction rates will further be suppressed by crystalline field. Indeed, we see much fewer CH2N2O2 molecules compared with HONO and NO2 during the dynamical process. Furthermore, the majority of CH2N2O2 molecules are founded to be formed directly from the ring fission of C4H8N8-nO8-2n, the product of HMX after losing nNO2, through C4H8N8-nO8-2n → CH2N2O2 + open C3H6N6-nO6-2n (n ≤ 3). Further decomposition of CH2N2O2 and ring fission lead to the productions of various small fragments, such as CH2N, CH2O, CHNO, CHN, and CNO.

21

Physical Chemistry Chemical Physics Accepted Manuscript

20

t = 40 ps

Physical Chemistry Chemical Physics

Page 22 of 39 View Article Online

DOI: 10.1039/C4CP00890A

propagation, confirming again that chemical reactions accompany heat conduction and temperature increase. Since NO2 is the most favorable initial decomposition product, we defined reaction fronts and

were defined at the front boundaries of NO2 distribution, such as 33.33 Å at t = 20 ps, and reaction centers at the boundaries of NO2 distribution with the maximum value, such as 22.22 Å at t = 20 ps. Thus the temperature is about 1080 K for reaction fronts and 1780 K for reaction centers as represented by the two dashed lines in Fig. 6(a). For regions that temperature is lower than 1780 K, less NO2 molecules are generated due to low temperature; for regions that temperature is higher than 1780 K, the number of NO2 molecules is also reduced arisen from subsequent reactions that consume NO2. The propagation rates of reaction front VRF (T = 1080 K) and reaction center VRC (T = 1780 K) are calculated to be 0.069 km·s-1 and 0.038 km·s-1 averaged from 20 ps to 80 ps, 47 times and 86 times slower than TMW propagation velocity (Vs = 3.32 km·s-1), respectively. Reaction front moves 1.8 times as fast as reaction center, leading to the expansion of reacted zone, which is consistent with the expanding lowdensity region and amorphous microstructure after 20 ps. H2O, N2 and CO2 are found to be major final products of HMX decomposition, the time evolutions of whose spatial distributions are plotted in Fig. 11. Distinctly, these products are initiated from “hot” regions and propagate from two sides towards the center, which is in consistence with the propagation direction of heat. This figure also manifests that the amounts and reaction rates for H2O and N2 are comparable, which are much higher than those for CO2.

80 60 40 20 0

1 2 3 4 5 6 7 8 9 10 11 12 13 14

Distance (x 11.11 Å)

100

t = 5 ps

80 60 40 20 0

1 2 3 4 5 6 7 8 9 10 11 12 13 14

Distance (x 11.11 Å)

Number of CO2

t = 5 ps

Number of N2

100

120

120

120

Number of H2O

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

reaction centers based on the time evolution of NO2 distribution. As presented in Fig. 9, reaction fronts

100

t = 5 ps

80 60 40 20 0

1 2 3 4 5 6 7 8 9 10 11 12 13 14

Distance (x 11.11 Å)

22

Physical Chemistry Chemical Physics Accepted Manuscript

These initial reactions occur from the two ends towards the center along the direction of heat

Page 23 of 39

Physical Chemistry Chemical Physics View Article Online

40 20

80 60 40 20 0

1 2 3 4 5 6 7 8 9 10 11 12 13 14

60 40 20

100

60 40 20 0

1 2 3 4 5 6 7 8 9 10 11 12 13 14

Number of N2

Number of H2O

t = 60 ps

60 40 20

100

Number of N2

Number of H2O

80 60 40 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14

Distance (x 11.11 Å)

60 40 20

40 20

100

100

t = 60 ps

80 60 40 20 0

1 2 3 4 5 6 7 8 9 10 11 12 13 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14

Distance (x 11.11 Å) 120

t = 80 ps

80 60 40 20 0

1 2 3 4 5 6 7 8 9 10 11 12 13 14

Distance (x 11.11 Å)

120 t = 80 ps

t = 40 ps

80

Distance (x 11.11 Å)

120

0

t = 60 ps

60

0

1 2 3 4 5 6 7 8 9 10 11 12 13 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14

120

80

Distance (x 11.11 Å) 100

100

0

1 2 3 4 5 6 7 8 9 10 11 12 13 14

120

80

20

Distance (x 11.11 Å)

120

0

t = 40 ps

80

Distance (x 11.11 Å) 100

40

Distance (x 11.11 Å)

Number of CO2

0

60

120

Number of CO2

Number of N2

Number of H2O

t = 40 ps

t = 20 ps

80

0

1 2 3 4 5 6 7 8 9 10 11 12 13 14

120

120

80

100

Distance (x 11.11 Å)

Distance (x 11.11 Å) 100

t = 20 ps

DOI: 10.1039/C4CP00890A

1 2 3 4 5 6 7 8 9 10 11 12 13 14

Distance (x 11.11 Å)

100

t = 80 ps

80 60 40 20 0

1 2 3 4 5 6 7 8 9 10 11 12 13 14

Distance (x 11.11 Å)

Figure 11. Spatial distributions of H2O (left), N2 (middle), and CO2 (right) at various time under thermal shock for HMX crystal without compression

For compressed systems, the time evolutions of major products were plotted in Fig. 8 (b-c). It shows that a great number of isolated H and O atoms appear at the beginning besides NO2 and HONO, which are distinguishable from those for uncompressed system. Additionally, the formations of various intermediates (NO3, HNO3, HO, NO, HNO, CHN, CHNO, CNO, CH2O) and final products (H2O, N2, CO, CO2) are also influenced by uniaxial compression. 23

Physical Chemistry Chemical Physics Accepted Manuscript

60

100

Number of CO2

80

0

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

t = 20 ps

120

Number of CO2

100

120

Number of N2

Number of H2O

120

Physical Chemistry Chemical Physics

Page 24 of 39 View Article Online

DOI: 10.1039/C4CP00890A

30 % compressed system in comparison with the uncompressed one as shown in Fig. 12(a). This is in agreement with the decomposition of RDX during shock simulation, which shows that the amount of

of NO2 is partially arisen from reversible reaction of N-NO2 bond rupture (C4H8N8O8 ⇌ C4H8N7O6 + NO2 or C4H8N7O6 ⇌ C4H8N6O4 + NO2) under high compression and partially ascribed to fast subsequent reactions that consume NO2. The latter factor can be confirmed by the increased quantities of NO3 and HNO3 for compressed systems (see Fig. 12(b-c)) since NO3 is mainly originated from NO2 + O → NO3 or 2NO2 → NO3 + NO and HNO3 forms via NO3 + H → HNO3 or NO2 + HO → HNO3. The accelerated subsequent reactions are mainly resulted from faster temperature increment under compression, making more molecules active. The temperature of the reaction centers reaches 2100 K and 2280 K for 15 % and 30 % compressed systems so that we expect subsequent reactions of these highly excited molecules. (b) 100 Compression = 0 Compression = 15 % Compression = 30 %

160

Compression = 0 Compression = 15 % Compression = 30 %

80

Number of NO3

Number of NO2

(a) 200

120

60

40

20

80

0

20

40

60

0

80

0

20

t / ps

60

80

60

80

(d) 100

Number of HONO

Compression = 0 Compression = 15 % Compression = 30 %

80

60

40

20

0

40

t / ps

(c) 100

Number of HNO3

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

NO2 remarkably decreases when impact velocity is higher than 8 km·s-1.18 The highly reduced number

Compression = 0 Compression = 15 % Compression = 30 %

80

60

40

20

0

20

40

t / ps

60

80

0

0

20

40

t / ps

24

Physical Chemistry Chemical Physics Accepted Manuscript

The number of NO2 is slightly smaller for the 15 % compressed system and highly reduced for the

Page 25 of 39

Physical Chemistry Chemical Physics View Article Online

80

80

60

40 Compression = 0 Compression = 15 % Compression = 30 %

20

0

0

20

40

60

Compression = 0 Compression = 15 % Compression = 30 %

60

40

20

80

0

0

20

t / ps

40

60

80

t / ps

Figure 12. Comparisons of the time evolutions of NO2 (a), NO3 (b), HNO3 (c), HONO (d), NO (e), and HNO (f) under thermal shock for HMX crystals with different compressions

NO2 spatial distributions at various times for compressed systems shown in Fig. 13 elucidate that the propagation rate of N-NO2 bond breaking is highly accelerated by uniaxial compression. NO2 molecules are detected in the center of the system at 60 ps and 30 ps for the 15 % and 30 % compressed systems, which are much earlier than that for the uncompressed system. The accelerated propagation rate of N-NO2 bond rupture can explain well the sharp increase of NO2 population after 50 ps for 15 % compressed system and after 20 ps for 30 % compressed system as shown in Fig. 12(a), where secondary reactions consuming NO2 do not start. The propagation velocity of reaction front VRF and reaction center VRC are 0.094 km·s-1 and 0.063 km·s-1 averaged from 20 ps to 60 ps for ΔV = 15 % and are 0.155 km·s-1 and 0.108 km·s-1 averaged from 5 ps to 30 ps for ΔV = 30 %. Thus, for the 15 % and 30 % compressed systems, VRF is 1.36 times and 2.25 times as fast as the uncompressed system, and VRC is 1.66 times and 2.71 times as rapid as the uncompressed system. The temperatures of reaction fronts TRF and reaction centers TRC are 1550 K and 2100 K for ΔV = 15 % and are 1780 K and 2280 K for ΔV = 30 % based on the time-spatial distributions of NO2 and temperature.

25

Physical Chemistry Chemical Physics Accepted Manuscript

(f) 100

Number of HNO

(e) 100

Number of NO Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

DOI: 10.1039/C4CP00890A

Physical Chemistry Chemical Physics

Page 26 of 39 View Article Online

DOI: 10.1039/C4CP00890A

70

t = 5 ps Number of NO2

50 40 30 20

0

t = 5 ps

60 50 40 30 20 10

10

0

1

2

3

4

5

6

7

8

9

10 11 12 13 14

1

2

3

4

t = 20 ps

50 40 30 20

9

10 11 12 13 14

t = 10 ps

50 40 30 20

0

1

2

3

4

5

6

7

8

9

10 11 12 13 14

1

2

3

4

5

6

7

8

9

10 11 12 13 14

Distance (x 7.77 Å)

Distance (x 9.44 Å) 70

70

t = 40 ps

60 50 40 30 20

t = 20 ps

60

Number of NO2

Number of NO2

8

10

10

10

50 40 30 20 10

1

2

3

4

5

6

7

8

9

0

10 11 12 13 14

1

2

3

4

Distance (x 9.44 Å)

5

6

7

8

9

10 11 12 13 14

Distance (x 7.77 Å) 70

70

t = 60 ps

60

40 30 20

t = 30 ps

60

Number of NO2

50

50 40 30 20 10

10 0

7

60

Number of NO2

Number of NO2

60

0

6

70

70

0

5

Distance (x 7.77 Å)

Distance (x 9.44 Å)

Number of NO2

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

Number of NO2

60

1

2

3

4

5

6

7

8

9

10 11 12 13 14

Distance (x 9.44 Å)

0

1

2

3

4

5

6

7

8

9

10 11 12 13 14

Distance (x 7.77 Å)

Figure 13. Spatial distributions of NO2 at various times under thermal shock for HMX crystals with 15 % compression (left) and with 30 % compression (right)

26

Physical Chemistry Chemical Physics Accepted Manuscript

70

Page 27 of 39

Physical Chemistry Chemical Physics View Article Online

DOI: 10.1039/C4CP00890A

high compression (ΔV = 30 %). The number of HONO for the low compressed system is close to that for the uncompressed system, but it is smaller for the high compressed system as shown in Fig. 12(d).

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

The origin of HONO is suggested to be nitro group abstracting hydrogen atom from methyl group or isolated NO2 directly attracting hydrogen atom from other radicals. The reduced amount of HONO is partially arisen from the opposite bending of nitro group and methyl group under compression. It is found that half of H-C-N-N dihedrals vary from 32.7ºfor the uncompressed system to 35.0ºfor ΔV = 15 % and to 37.2ºfor ΔV = 30 %, leading to longer non-bonded distance of H···O, making it more difficult to form HONO for compressed systems. The decomposition of HONO leads to the formations of NO and HNO via HONO → HO + NO and NO + H → HNO, the amounts of which reduce remarkably for the high compressed system as shown in Fig. 12(e-f), corroborating that HONO elimination and subsequent decomposition are impeded by high compression. Although the quantity of HONO is reduced, the propagation rate of HONO elimination is accelerated by uniaxial compression (see Fig. S3 of the ESI). The populations of products generated from ring fission (CH2N2O2, CH2O, CH2N, CHN, CHNO, CNO) decrease exceedingly under compression as depicted in Fig. 8 (b-c). This is mainly ascribed to accelerated subsequent reactions that quickly deplete these unstable initial products or intermediates, leading to the increase of final products, such as N2 and CO2. The amounts of N2 and CO2 go up monotonously with the increase of compression as shown in Fig. 14(a-b), which are the same as the decomposition of RDX under shock loadings showing that the quantities of N2 and CO2 rise with the increase of impact velocity.18 It is found that all eight C-N bonds as well as six C-N-C angles reduce under compression that lead to stronger intramolecular interactions, indicating that the changes of molecular structure in bonds and angles arisen from uniaxial compression do not favor C-N bond breaking. Thus, the promoted C-N bond rupture is mainly ascribed to higher temperature for compressed systems as discussed above. Furthermore, N-N bond breaking also contributes to the 27

Physical Chemistry Chemical Physics Accepted Manuscript

HONO elimination is hardly influenced by low compression (ΔV = 15 %), however, it is hindered by

Physical Chemistry Chemical Physics

Page 28 of 39 View Article Online

DOI: 10.1039/C4CP00890A

reactions under compression leads to more N2 formed in shorter time. The increased number of isolated oxygen atom under compression (Fig. 14(c)) benefits the formation of CO2 via oxidation of CO,

(a) 800

(b) 400 Compression = 0 Compression = 15 % Compression = 30 %

Number of CO2

Number of N2

600

400

200

0

0

20

40

60

Compression = 0 Compression = 15 % Compression = 30 %

300

200

100

0

80

0

20

t / ps

80

60

40

Compression = 0 Compression = 15 % Compression = 30 %

80

Number of CO

Number of O

60

(d)100 Compression = 0 Compression = 15 % Compression = 30 %

80

20

0

40

t / ps

(c) 100

60

40

20

0

20

40

60

0

80

0

20

t / ps

60

80

(f) 400 Compression = 0 Compression = 15 % Compression = 30 %

Compression = 0 Compression = 15 % Compression = 30 %

300

Number of H

600

400

200

100

200

0

40

t / ps

(e) 800

Number H2O

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

resulting in the decrease of CO (Fig. 14(d)) while increase of CO2.

0

20

40

t / ps

60

80

0

0

20

40

60

80

t / ps

28

Physical Chemistry Chemical Physics Accepted Manuscript

formation of N2, therefore, the enhanced reaction rate of N-N bond breaking as well as subsequent

Page 29 of 39

Physical Chemistry Chemical Physics View Article Online

Compression = 0 Compression = 15 % Compression = 30 %

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

Number of O2

Number of HO

300

200

Compression = 0 Compression = 15 % Compression = 30 %

80

100

60

40

20

0

20

40

t / ps

60

80

0

0

20

40

60

80

t / ps

Figure 14. Comparisons of the time evolutions of N2 (a), CO2 (b), O (c), CO (d), H2O (e), H (f), HO (g), and O2 (h) under thermal shock for HMX crystals with different compressions

The number of another primary final product H2O rises for HMX crystal with low compression, but declines for high compressed system as shown in Fig. 14(e). The amounts of H and HO fragments go up dramatically under compression (Fig. 14(f-g)), leading to the increase of H2O (H + HO → H2O). However, this reaction is found to be reversible under high compression, resulting in the reduction of H2O. The reaction propagation rates for the three final products CO2, N2, and H2O are accelerated by compression, the time evolutions of whose distributions under compression are presented in Fig. S4 ~ Fig. S6 of the ESI. It should be noticed that uniaxial compression promotes the formations of many small fragments, such as H, O, HO, O2, etc, which are in coincidence with those observed during the decomposition of RDX supposed to shock loading at high impact velocities.18 Small free space under compression leads to large deformation of molecules and strong intermolecular contact and extrusion, making it easier for atoms out of the ring to be crowded out that results in the formations of isolated H and O. These small species can move freely in small free space and easily collide with each other, forming HO, O2, etc. The origin of HO is mainly from the decomposition of HONO (HONO → HO + NO) and partially from the bond formation of H and O (H + O → HO) at atmospheric pressure, but these reaction mechanisms are varied by compression. The reduced number of HONO under compression leads to the reduction of HO and NO via HONO decomposition, thus, the increased HO is mainly arisen from the combination of H 29

Physical Chemistry Chemical Physics Accepted Manuscript

(g) 400

0

DOI: 10.1039/C4CP00890A

(h) 100

Physical Chemistry Chemical Physics

Page 30 of 39 View Article Online

DOI: 10.1039/C4CP00890A

compression contributes to the increment of HO. Although the quantities of some products are reduced by uniaxial compression, the propagation rates

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

of initial reactions and subsequent reactions are well accelerated. This is consistent with previous experimental results showing that the reaction propagation rates of HMX increase fairly monotonically with pressure.46 The enhanced propagation rates of chemical reactions are arisen from the accelerated heat conduction under uniaxial compression, leading to higher temperature of the system and thus more molecules are excited to be active. This promotes the initial bond breaking and chemical reaction initiation followed by subsequent exothermic reactions, benefiting the formation of intermediates and final products that in turn promote the increase of temperature. 4. Conclusions We developed a thermal shock reactive dynamics (TS-RD) computational protocol to investigate the thermal expansion, heat conduction, and chemical initiation and propagation of HMX crystal under thermal shock. These nonequilibrium MD simulations provide some insight into the microscopic thermal-chemical response during thermal shock process in condensed explosives. Thermal shock initiated mechanical wave —TMW propagates back and forth with a velocity of 3.32 km·s-1, leading to the compression-expansion cycle of the system. It decays off after several periods, at which mass density of “cold” region begins to decrease and chemical reactions initiate. Uniaxial compression accelerates the propagation and attenuation of TMW and benefits heat and mass conductions, leading to earlier reduction of mass density and initiation of chemical reaction. “Hot” regions located at the two ends serving as heat resources heat up “cold” molecules, initiating heat propagation from the ends to the center that leads to temperature increment of the system. Combining the simulated temperature distribution with the continuum reactive heat conduction model, the heat conduction coefficient is derived to be 0.80 W·m-1·K-1. Uniaxial compression remarkably speeds up heat propagation, resulting in higher heat conduction coefficient and temperature. 30

Physical Chemistry Chemical Physics Accepted Manuscript

and O by means of multimolecular process. In addition, the decomposition of H2O promoted by high

Page 31 of 39

Physical Chemistry Chemical Physics View Article Online

DOI: 10.1039/C4CP00890A

fission, among which N-NO2 bond rupture is most favorable. The amounts of NO2 and the initial products originated from ring fission decrease under uniaxial compression ascribed to accelerated

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

subsequent reactions. HONO elimination is hindered for highly compressed system, resulting in the reductions of HONO and subsequent reaction products. Small molecules (like H, O, and HO) generated from multimolecular reactions become dominant products as compression increases. Chemical reactions propagate from “hot” region to “cold” region along the propagation direction of heat, and such propagation rates are significantly accelerated by uniaxial compression, in agreement with previous report. Additionally, the propagation rates as well as temperatures of reaction fronts and reaction centers are obtained. The whole thermal shock process can be divided into two stages. The first stage is mainly characterized by physical responses, specifically the thermal expansion of particles and periodic mechanical wave propagation. The second stage is featured by chemical responses represented by reaction initiation and propagation followed by subsequent reactions leading to the formations of a variety of intermediates and final products that release energies which in turn contribute to temperature increment and sustain decomposition and combustion. Heat propagation proceeds through the whole process, resulting in temperature increasing of the system and chemical reaction initiation. Acknowledgments This work was supported by the National Natural Science Foundation of China (Grant No. 11172044 and Grant No. 11372053) and by the Open Grant of State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, P. R. China (Grant No. KFJJ12-6M). Y. Liu. acknowledges the supports from “Shanghai Pujiang Talent” program (Grant No. 12PJ1406500), Science and Technology Commission of Shanghai Municipality; “Key Program of Innovative Scientific Research” (Grant No. 14ZZ130), Education Commission of Shanghai Municipality; “Recruit Program

31

Physical Chemistry Chemical Physics Accepted Manuscript

Chemical reaction is initiated by N-NO2 bond breaking, followed by HONO elimination and ring

Physical Chemistry Chemical Physics

Page 32 of 39 View Article Online

DOI: 10.1039/C4CP00890A

performance computing facilities at the University of Shanghai for Science and Technology (USST), Shanghai Supercomputer Center, and National Supercomputing Center in Shenzhen, P. R. China.

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

Notes and references a

Institute of Applied Physics and Computational Mathematics, Beijing, 100094, P. R. CHINA.

b

State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, Beijing,

100081, P. R. CHINA. c

School of Materials Science and Engineering, University of Shanghai for Science and Technology,

Shanghai, 200093, P. R. CHINA. E-mail: [email protected] † Electronic Supplementary Information (ESI) available: Microstructure profile for HMX crystal without compression at various times under thermal shock is plotted in Fig. S1. Probability distribution of particle velocity under thermal shock for HMX crystal without compression is shown in Fig. S2. Spatial distributions of HONO, CO2, N2, and H2O at various times under thermal shock for compressed HMX crystals are presented in Fig. S3 ~ Fig. S6. Stress tensors for the systems after uniaxial compression obtained from NVT-MD are collected in Table S1. Bond order cutoff values for various atom pairs used in the algorithm of molecule recognition are tabulated in Table S2. Several important physical parameters to present a direct and clear cognition about the effects of compression are summarized in Table S3. 1. P. D. Peterson, J. T. Mang, and B. W. Asay, Quantitative analysis of damage in an octahydro1,3,5,7-tetranitro-1,3,5,7-tetrazonic-based composite explosive subjected to a linear thermal gradient, J. Appl. Phys., 2005, 97, 093507–1–7. 2. W. B. Zhang, Y. Tian, M. P. Wen, and Y. Hao, Experimental Study on the Thermal Shock Damage

32

Physical Chemistry Chemical Physics Accepted Manuscript

of Global Expert” (“Thousands plan”) in Shanghai. Some computations were carried out using the high

Page 33 of 39

Physical Chemistry Chemical Physics View Article Online

DOI: 10.1039/C4CP00890A

3. Y. Tian, S. H. Luo, and W. B. Zhang, Water-Bathed Thermal Shock Damage of PBX JOB-9003 and

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

its Ultrasonic Characteristics, Chinese Journal of Explosives and Propellants, 2002, 3, 17–19. 4. T. M.Willey,L. Lauderbach, F. Gagliard,B. Cunningham, K. T. Lorenz, J. R. I. Lee, T. van Buuren, R. Call, L. Landt, and G. Overturf, Comprehensive Characterization of Voids and Microstructure in TATB-based Explosives from 10nm to 1cm: Effect of Temperature Cycling and Compressive Creep, The 14th International Detonation Symposium, Coeur d'Alene, ID, United States, April 11–16, 2010, p530–538. 5. Q. Lan, J. Y. Lu, M. Zhang, and L. Yong, Properties for PBX Cylinder during Temperature Rising, Chinese Journal of Energetic Materials, 2008, 16, 693–697. 6. G. T. Gray Ⅲ, D. J. Idar, W. R. Blumenthal, C. M. Candy, and P. D. Peterson, High and low-strain rate compression properties of several energetic material composites as a function of strain rate and temperature, The 11th International Detonation Symposium, Snowmass VILLAGE, CO, United States, Aug 31–Sep 4, 1998, p229–231. 7. D. A. Wiegand, Mechanical properties and mechanical failure of composite plastic bonded explosives and other energetic materials, The 11th International Detonation Symposium, Snowmass VILLAGE, CO, United States, Aug 31–Sep 4, 1998, p85–88. 8. Y. Tian, W. B Zhang, M. P. Wen, Z. F. Yang, Y. Hao, and J. M. Li, Research on correlation of thermal shock damage of PBX JOB-9003, Chinese Journal of Energetic Material, 2004, 12, 174–177. 9. H. L. Berghout, S. F. Son, C. B. Skidmore, D. J. Idar, and B. W. Asay, Combustion of damaged PBX 9501 explosive, Thermochim. Acta, 2002, 384, 261–277. 10. B. W. Asay, P. M. Dickson, B. F. Henson, L. Smilowitz, and L. Tellier, Effect of Temperature 33

Physical Chemistry Chemical Physics Accepted Manuscript

of Explosive by Ultrasonic Testing, Chinese Journal of Energetic Materials, 2004, 12, 85–88.

Physical Chemistry Chemical Physics

Page 34 of 39 View Article Online

DOI: 10.1039/C4CP00890A

Matter-2001: Atlanta, GA, 2001, American Institute of Physics, New York, 2002, p1065–1068. 11. B. W. Asay, B. F. Henson, L. B. Smilowitz, and P. M. Dickson, On the Difference in Impact

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

Sensitivity of Beta and Delta HMX. J. Energ. Mater., 2003, 21, 223–235. 12. D. Bedrov, J. B. Hooper, G. D. Smith, and T. D. Sewell, Shock-induced transformations in crystalline RDX: A uniaxial constant-stress Hugoniostat molecular dynamics simulation study, J. Chem. Phys., 2009, 131, 034712–1–12. 13. A. C. T. van Duin, S. Dasgupta, F. Lorant, and W. A. Goddard III, ReaxFF: A Reactive Force Field for Hydrocarbons, J. Phys. Chem. A, 2001, 105, 9396-9409. 14. K. Chenoweth, A. C. T. van Duin, and W. A. Goddard III, ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A, 2008, 112, 1040-1053. 15. J. E. Mueller, A. C. T. van Duin, and W. A. Goddard III, Application of the ReaxFF Reactive Force Field to Reactive Dynamics of Hydrocarbon Chemisorption and Decomposition. J. Phys. Chem. C, 2010, 114, 5675-5685. 16. K. Chenoweth, S. Cheung, A. C. T. van Duin, W. A. Goddard III, and E. M. Kober, Simulations on the Thermal Decomposition of a Poly(dimethylsiloxane) Polymer Using the ReaxFF Reactive Force Field. J. Am. Chem. Soc., 2005, 127, 7192-7202. 17. K. Nomura, R. K. Kalia, A. Nakano, and P. Vashishta, Reactive Nanojets: Nanostructure Enhanced Chemical Reactions in a Defected Energetic Crystal, Appl. Phys. Lett., 2007, 91, 183109-1-3. 18. A. Strachan, A. C. T. van Duin, D. Chakraborty, S. Dasgupta, and W. A. Goddard III, Shock Waves in High-Energy Materials: The Initial Chemical Events in Nitramine RDX, Phys. Rev. Lett., 2003, 91, 098301-1-4. 34

Physical Chemistry Chemical Physics Accepted Manuscript

Profile on Reaction Violence in Heated and Self-ignited PBX 9501, Shock Compression of Condensed

Page 35 of 39

Physical Chemistry Chemical Physics View Article Online

DOI: 10.1039/C4CP00890A

decomposition of RDX from reactive molecular dynamics, J. Chem. Phys., 2005, 122, 054502-1-10. 20 A. C. T. van Duin, Y. Zeiri, F. Dubnikova, R. Kosloff, and W. A. Goddard III, Atomistic-Scale

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

Simulations of the Initial Chemical Events in the Thermal Initiation of Triacetonetriperoxide. J. Am. Chem. Soc., 2005, 127, 11053-11062. 21. S. V. Zybin, W. A. Goddard III, P. Xu, A. C. T. van Duin, and A. P. Thompson, Physical mechanism of anisotropic sensitivity in pentaerythritol tetranitrate from compressive-shear reaction dynamics simulations, Appl. Phys. Lett., 2010, 96, 081918-1-3. 22. L. Zhang, S. V. Zybin, A. C. T. van Duin, S. Dasgupta, W. A. Goddard III, and E. M. Kober, Carbon Cluster Formation during Thermal Decomposition of Octohydro-1,3,5,7-tetranitro-1,3,5,7tetrazocine and 1,3,5-Triamino-2,4,6-trinitrobenzene High Explosives from ReaxFF Reactive Molecular Dynamics Simulations, J. Phys. Chem. A, 2009, 113, 10619-10640. 23 T. T. Zhou, and F. L. Huang, Effects of Defects on Thermal Decomposition of HMX via ReaxFF Molecular Dynamics Simulations. J. Phys. Chem. B, 2011, 115, 278–287. 24. J. Budzien, A. P. Thompson, and S. V. Zybin, Reactive Molecular Dynamics Simulations of Shock Through a Single Crystal of Pentaerythritol Tetranitrate. J. Phys. Chem. B, 2009, 113, 13142–13151. 25. Q. An, Y. Liu, S. V. Zybin, H. Kim, and W. A. Goddard III, Anisotropic Shock Sensitivity of Cyclotrimethylene Trinitramine (RDX) from Compress-and-Shear Reactive Dynamics. J. Phys. Chem. C, 2012, 116, 10198–10206. 26. T. T. Zhou, S. V. Zybin, Y. Liu, F. L. Huang, and W. A. Goddard III, Anisotropic Shock Sensitivity for β-HMX Energetic Material under Compressive-Shear Loading from ReaxFF-lg Reactive Dynamics Simulations. J. Appl. Phys., 2012, 111, 124904–1–11. 27. Q. An, S. V. Zybin, W. A. Goddard III, A. J. Botero, M. Blanco, and S. N. Luo, Elucidation of the dynamics for hot-spot initiation at nonuniform interfaces of highly shocked materials, Phys. Rev. B, 35

Physical Chemistry Chemical Physics Accepted Manuscript

19. A. Strachan, E. Kober, A. C. T. van Duin, J. Oxgaard, and W. A. Goddard III, Thermal

Physical Chemistry Chemical Physics

Page 36 of 39 View Article Online

DOI: 10.1039/C4CP00890A

28. Q. An, W. A. Goddard III, S. V. Zybin, A. J. Botero, and T. T. Zhou, Highly Shocked Polymer Bonded Explosives at a Nonplanar Interface: Hot spot Formation leading to Detonation. J. Phys. Chem.

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

C, 2013, 117, 26551–26561. 29. T. T. Jarvi, A. Kuronen, M. Hakala, K. Nordlund, A. C. T. van Duin, W. A. Goddard III, and T. Jacob, Development of a ReaxFF description for gold. Eur. Phys. J. B, 2008, 66, 75–79. 30. M. R. LaBrosse, J. K. Johnson, and A. C. T. van Duin, Development of a Transferable Reactive Force Field for Cobalt. J. Phys. Chem. A, 2010, 114, 5855–5861. 31. J. E. Mueller, A. C. T. van Duin, and W. A. Goddard III, Development and Validation of ReaxFF Reactive Force Field for Hydrocarbon Chemistry Catalyzed by Nickel. J. Phys. Chem. C, 2010, 114, 4939–4949. 32. K. Joshi, A. C. T. van Duin, and T. Jacob, Development of a ReaxFF description of gold oxides and initial application to cold welding of partially oxidized gold surfaces. J. Mater. Chem., 2010, 20, 10431– 10437. 33. A. C. T. van Duin, V. S. Bryantsev, M. S. Diallo, W. A. Goddard III, O. Rahaman, D. J. Doren, D. Raymand, and K. Hermansson, Development and Validation of ReaxFF Reactive Force Field for Hydrocarbon Chemistry Catalyzed by Nickel. J. Phys. Chem. A, 2010, 114, 9507–9514. 34. W. A. Goddard III, J. E. Mueller, K. Chenoweth, and A. C. T. van Duin, ReaxFF Monte Carlo reactive dynamics: Application to resolving the partial occupations of the M1 phase of the MoVNbTeO catalyst. Catal. Today, 2010, 157, 71–76. 35. L. C. Liu, Y. Liu, S. V. Zybin, and W. A. Goddard III, ReaxFF-lg: Correction of the ReaxFF Reactive Force Field for London Dispersion, with Applications to the Equations of State for Energetic Materials, J. Phys. Chem. A, 2011, 115, 11016-11022.

36

Physical Chemistry Chemical Physics Accepted Manuscript

2011, 84, 220101-1-5.

Page 37 of 39

Physical Chemistry Chemical Physics View Article Online

DOI: 10.1039/C4CP00890A

by neutron diffraction, Acta Crystallogr., 1970, B26, 1235-1240. 37. K. A. Gonthier, Modeling and analysis of reactive compaction for granular energetic solids,

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

Combust. Sci. and Tech., 2003, 175, 1679-1709. 38. C. M. Tarver, S. K. Chidester, and A. L. Nichols III, Critical Conditions for Impact- and ShockInduced Hot Spots in Solid Explosives, J. Phys. Chem., 1996, 100, 5794-5799. 39. A. Tokmakoff, M. D. Fayer, and D. D. Dlott, Chemical reaction initiation and hot-spot formation in shocked energetic molecular materials, J. Phys. Chem., 1993, 97, 1901-1913. 40. D. R. Chakraborty, P. Muller, S. Dasgupta, and W. A. Goddard III, Mechanism for Unimolecular Decomposition of HMX (1,3,5,7-tetranitro-1,3,5,7-tetrazocine), an ab Initio Study, J. Phys. Chem. A, 2001, 105, 1302-1314. 41. O. Sharia, and M. M. Kuklja, Modeling Thermal Decomposition Mechanisms in Gaseous and Crystalline Molecular Materials: Application to β-HMX, J. Phys. Chem. B, 2011, 115, 12677-12686. 42. R. Behrens, Thermal decomposition of energetic materials: temporal behaviors of the rates of formation of the gaseous pyrolysis products from condensed-phase decomposition of octahydro-1,3,5,7tetranitro-1,3,5,7-tetrazocine, J. Phys. Chem., 1990, 94, 6706-6718. 43. R. Behrens, Identification of octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX) pyrolysis products by simultaneous thermogravimetric modulated beam mass spectrometry and time-of-flight velocity-spectra measurements, Int. J. Chem. Kinet., 1990, 22, 135-157. 44. R. Behrens, and S. Bulusu, Thermal decomposition of energetic materials. 2. Deuterium isotope effects and isotope scrambling in condensed-phase decomposition of octahydro-1,3,5,7-tetranitro-

37

Physical Chemistry Chemical Physics Accepted Manuscript

36. C. S. Choi, and H. P. Boutin, A study of the crystal structure of β-cyclotetramethylene tetranitramine

Physical Chemistry Chemical Physics

Page 38 of 39 View Article Online

DOI: 10.1039/C4CP00890A

45. T. B. Brill, P. E. Gongwer, and G. K. Williams, Thermal Decomposition of Energetic Materials. 66.

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

Kinetic Compensation Effects in HMX, RDX, and NTO, J. Phys. Chem., 1994, 98, 12242–12247. 46. A. P. Esposito, D. L. Farber, J. E. Reaugh, and J. M. Zaug, Reaction Propagation Rates in HMX at High Pressure, Propellants, Explosives, and Pyrotechnics., 2003, 28, 83-88.

38

Physical Chemistry Chemical Physics Accepted Manuscript

1,3,5,7-tetrazocine, J. Phys. Chem., 1991, 95, 5838-5845.

Page 39 of 39

Physical Chemistry Chemical Physics View Article Online

DOI: 10.1039/C4CP00890A

Thermal Shock Initiated Thermal and Chemical Responses of HMX Crystal from ReaxFF Molecular Dynamics Simulation

Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.

(Tingting Zhou, Huajie Song, Yi Liu, Fenglei Huang)

Heat propagating along the direction of thermal shock leads to temperature increase and chemical reaction initiation. Propagation rates of reaction front and reaction center are obtained based on NO2 distribution.

1

Physical Chemistry Chemical Physics Accepted Manuscript

Table of Contents

Shock initiated thermal and chemical responses of HMX crystal from ReaxFF molecular dynamics simulation.

To gain an atomistic-level understanding of the thermal and chemical responses of condensed energetic materials under thermal shock, we developed a th...
4MB Sizes 1 Downloads 3 Views