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/gcm
Density/gcm
-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/gcm
-3
60
t = 6 ps t = 10 ps t = 15 ps
Distance/Å
(b) 2.4
Density/gcm
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/gcm
Density/gcm
-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/gcm
-3
-3
2.3
Density/gcm
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/gcm
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/gcm
-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/gcm
Published on 24 April 2014. Downloaded by University of California - Berkeley on 05/05/2014 19:04:23.
Density/gcm
-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