Correlation between Chemical Bond Cleavage and Detonation of ε-2

Mar 5, 2019 - Researchers have been striving to determine the connection between the microscopic chemical reactions and macroscopic detonation laws of...
0 downloads 0 Views 2MB Size
Subscriber access provided by ECU Libraries

C: Energy Conversion and Storage; Energy and Charge Transport

Correlation Between Chemical Bond Cleavage and Detonation of #-2,4,6,8,10,12-Hexanitrohexaazaisowurtzitane Danyang Liu, Lang Chen, Deshen Geng, Jianying Lu, and Junying Wu J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.9b01975 • Publication Date (Web): 05 Mar 2019 Downloaded from http://pubs.acs.org on March 10, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Correlation Between Chemical Bond Cleavage and Detonation of ɛ-2,4,6,8,10,12-Hexanitrohexaazaisowurtzitane Danyang Liu†, Lang Chen†*, Deshen Geng†, Jianying Lu†, Junying Wu†



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

AUTHOR INFORMATION Corresponding Author *E-mail: [email protected]

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 27

Abstract: Researchers have been striving to determine the connection between the microscopic chemical reactions and macroscopic detonation laws of explosives. In this study, we performed reactive molecular dynamics simulations

of

shock-induced

explosion

of

the

2,4,6,8,10,12-hexanitrohexaazaisowurtzitane (CL-20) explosive. The results show that detonation is mainly determined by rapid irreversible cleavage of the C–N and C–H bonds. Such C–N and C–H bond cleavage determines early formation of N2 and H2O. The detonation reaction occurs when the cleavage rates exceed 3.11%/ps and 4.15%/ps for the C–N and C–H bonds, respectively. Higher shock velocity results in higher cleavage rates of these bonds, but it also leads to more atoms being trapped in clusters. However, the decomposition rate of these clusters is mainly affected by the decrease in the density, not by the shock velocity, indicating that the late detonation reaction is mainly based on the characteristics of the explosive.

ACS Paragon Plus Environment

Page 3 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1. INTRODUCTION The detonation reaction is a unique phenomenon of explosives that rapidly and steadily releases energy, so it can be used in weapons. Many experiments show that the detonation process is sustained by propagation of the detonation wave along the explosive. Classical Zeldovich–Neumann–Döring (ZND) detonation theory1,2 states that the detonation wave consists of a shock wave front and a subsequent reaction zone. In the detonation process, the explosive is first compressed by the shock wave, and rapid chemical reactions then occur accompanied by volume expansion until the explosive completely releases its energy. However, although ZND detonation theory can solve many detonation engineering calculation problems,

the detonation reaction

mechanism and reaction details are still unclear because the extreme reaction environment and rapid reaction process make the chemical reaction details difficult to experimentally observe. Molecular dynamics (MD) simulations based on the reactive force field (ReaxFF)

3−5

offer the possibility of understanding the detonation reaction at

the microscopic level. The ReaxFF is based on bond order theory, in which cleavage and generation of chemical bonds are judged by the distances between the atoms. In 2003, Strachan et al.4 performed MD simulations with the ReaxFF to investigate the cyclotrimethylenetrinitramine (RDX) explosive reaction. This method has subsequently been widely used to analyze the reaction mechanisms of explosives 6



14

. In 2011, based on the ReaxFF

framework, Liu15 revised the description of the dispersive force and proposed a modified ReaxFF (ReaxFF-lg) that better describes the long-range interactions between molecules and can be used to obtain more accurate

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 27

molecular crystal structures16 20. ReaxFF MD can achieve simulation times −

various magnitudes longer than quantum mechanics and provide highly detailed reaction details. However, the larger simulation system also makes it more difficult to analyze the mechanism because the information is buried in millions of reactions. Existing research has mainly focused on the initial reaction details of explosives for the first few hundred picoseconds under shock or heating, and results revealing the connection between the microscopic chemical reactions and the macroscopic detonation laws are limited 21 , 22 . In this study, we investigated the relationship between the chemical

structure

and

the

detonation

characteristics

of

the

2,4,6,8,10,12-hexanitrohexaazaisowurtzitane (CL-20) explosive23,24, which is currently the most energetic explosive, using the MD method with the ReaxFF-lg. We developed a specific processing program to discover the high-frequency reactions and track them by observing snapshots of a unit cell in the system, so the reaction path can be analyzed in a simple and clear way.

2. COMPUTATIONAL METHODS The simulations were performed with the LAMMPS Molecular Dynamics Simulator using the ReaxFF-lg force field. We investigated thermal decomposition of CL-20 at various temperatures and densities16,17. The results show that the ReaxFF-lg force field provides a good description of the CL-20 crystal structure and reaction energy compared with experimental and density functional theory (DFT) results. Another key issue is how to apply a shock load to the explosive crystal. Here, we used the multiscale shock technique (MSST) 25

proposed by Reed et al., which combines MD simulation with macroshock

dynamics calculation theory and thereby greatly reduces the size of the

ACS Paragon Plus Environment

Page 5 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

calculation system. In this technique, a kinetic equation developed by the Hamiltonian function is constructed to constrain the MD simulation system in the macroscopic thermodynamic state by satisfying the Hugoniot shock condition with the Euler equation and the Hugoniot relationship. Thus, shock compression of the system can be achieved by changing the predefined direction volume of the cell according to a specified shock velocity. This method has been applied to detonation calculations of various explosives26 29. −

The ε-CL-20 phase is the most stable CL-20 phase at room temperature, and it is widely used in mixed explosives. Thus, we chose to investigate the ε-CL-20 phase. The single-crystal structure of ε-CL-20 was obtained from X-ray diffraction crystal data 30 and then expanded along the a, b, and c directions to construct a 4 × 3 × 3 CL-20 supercell. A schematic of the CL-20 molecule and supercell are shown in Figure S1 (Supporting Information). Before shock loading, the supercell was relaxed to obtain a stable structure under normal temperature and pressure. The structure was first relaxed by energy minimization, followed by thermalization at 298 K. A MD simulation was first performed for 10 ps with the canonical ensemble (NVT, constant number of particles, volume, and temperature) using the Berendsen thermostat. Pressure relaxation at 0 GPa was then performed for 5 ps using the isobaric–isothermal ensemble (NPT, constant number of particles, pressure, and temperature) with the Nose–Hoover thermostat and barostat. The density predicted

by

the

ReaxFF-lg

is

1.932

g/cm3.

According

to

the

Becker–Kistiakowsky–Wilson thermodynamic calculation program31, when the density is 1.932 g/cm3, the calculated detonation velocity of CL-20 is 8920 m/s, which is close to 9 km/s. Therefore, the relaxed supercell was subjected to

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

shock loading using the MSST at five shock velocities from 7 to 11 km/s along the orthogonal [010] direction. The calculation time was 50 ps. The time step was 0.1 fs for 7 and 8 km/s and 0.05 fs for the shock velocities of 9 km/s and above.

3. RESULTS AND DISCUSSION 3.1 Hugoniots from MSST shock simulations The shock velocity–particle velocity and pressure–volume Hugoniot plots from the ReaxFF-lg MD simulations along with the experimental data of the C-1 mixed explosive (CL-20/binder/96/4) 32 are shown in Figure 1. These calculated CL-20 shock Hugoniot results were obtained by shock loading the CL-20 supercell with different shock velocities from 3.5 to 10 km/s along the [100], [010], and [001] directions. The calculation time was 1 ps and the time step was 0.1 fs. The simulation results are in good agreement with the experimental results, indicating that MD can attain an accurate CL-20 state under shock conditions. At the same pressure, the experimental C-1 explosive (1.92 g/cm3) has a higher degree of compression because it contains a small amount of binder and pores, and thus it is easier to compress under shock.

Figure 1. Shock velocity versus particle velocity (left) and pressure versus volume compression (right) for the CL-20 supercell and C-1 explosive.

ACS Paragon Plus Environment

Page 6 of 27

Page 7 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Assuming that the relationship between the shock velocity Us and particle velocity Up is linear (left of Figure 1) then the Us–Up Hugoniot Us = C0 + sUp. The intercepts C0 and slopes are given in Table 1. The Us–Up Hugoniots for the three different directions are similar, indicating that the mechanical properties of CL-20 are similar in different shock directions. Table 1. Intercepts and slopes of the Us–Up Hugoniots of CL-20 This work

Shock

C-1 experiment

orientation

C0

s

x[100]

3.68

2.03

y[010]

3.3

2.31

z[001]

3.46

2.14

C0

s

2.6

2.24

3.2 Evolution of the pressure and temperature The pressure and temperature profiles of CL-20 under 7–11 km/s shock loading are shown in Figure 2. The results show that the system pressure and temperature instantaneously increase to their shocked state values, relatively slowly increase for a period of time, and then the system tends to be stable. A higher shock velocity results in higher pressure and temperature, as well as a shorter time to reach a steady state. In the ZND model, the detonation wave propagates in the explosive at the detonation velocity and the pressure in the detonation reaction zone decreases from the Neumann compression peak to the Chapman–Jouguet point (end point of the reaction). Figure 2 shows that the system pressure is 58 GPa under 9 km/s shock, which is similar to the experimental Neumann peak pressure of 51 GPa calculated with the parameters in Table 1, indicating that the MSST simulation with the ReaxFF-lg

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

describes the stable detonation reaction of the CL-20 explosive. At this time, the system temperature is 4245 K.

Figure 2. Pressure (left) and temperature (right) evolution under different shock velocities. 3.3 Reaction path To explore the high frequency reactions and analyze the reaction paths, a data processing program was developed using the Python language. The program first analyzes the bond information of each atom with other atoms in each time step. By comparing the time steps before and after, bond cleavage and formation can be determined. Finally, the chemical reactions occurring in the interval can be inferred. However, in the late stage of the reaction, the types of reactions in the system are very complicated, which makes the analysis very difficult. Therefore, in this study, the reaction of CL-20 was mainly analyzed in the first 6 ps. The main reaction paths, reaction frequencies, and reaction times under different shock velocities are given in Table 2. The analysis interval was 0.05 ps.

ACS Paragon Plus Environment

Page 8 of 27

Page 9 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 2. Main reaction paths, reaction frequencies, and reaction times of CL-20 under different shock velocities a Shock Frequen velocity

Time (ps)

Elementary reaction

121

0.1-6

C6H6N12O12+C6H6N12O12 →C12H12N24O24

168

0.1-6

Illustration

-cies (km/s)

C6H6N12O12+C12H12N24O24→

Oligomerization

C18H18N36O36

of CL-20

7 79

0.2-6

C12H12N24O24 →2C6H6N12O12

165

0.25-6

C18H18N36O36→C6H6N12O12+C12H12N24O24

2

2.1-6

C12H12N24O24→C12H12N23O22+NO2

50

0.05-4.7

C6H6N12O12+C6H6N12O12 →C12H12N24O24

29

0.5-5.5

C6nH6nN12nO12n →C6nH6nN12n-1O12n-2+NO2

21

1.35-4.7

C6nH6nN12nO12n+NO2→C6nH6nN12n+1O12n+2

20

3.15-5.9

C6nH6nN12n-1O12n-2+NO2→C6nH6nN12nO12n

36

0.05-0.25

C6H6N12O12+C6H6N12O12 →C12H12N24O24

b

20

0.45-0.9

C6nH6nN12nO12n→C6nH6nN12n-1O12n-2+NO2

c

11

0.5-1.95

C6nH6nN12nO12n+NO2→C6nH6nN12n+1O12n+2

11

0.7-3.1

C6nH6nN12n-1O12n-2+NO2→C6nH6nN12nO12n

3

5.1-5.9

N2+O2→N2O2

Generate NO2 b

Oligomerization

c

8 Generate NO2

Oligomerization

Generate NO2

9

Generate N2 3

5.15-5.95

N2O2→N2+O2

19

0.05-0.1

C6H6N12O12+C6H6N12O12→C12H12N24O24

3

0.2-0.9

C6nH6nN12nO12n→C6nH6nN12n-1O12n-2+NO2

2

0.6-0.85

C6nH6nN12nO12n+NO2→C6nH6nN12n+1O12n+2

10

0.3-2.2

NO2+NO2→N2O4

9

1.3-3.75

N2O4→NO2+NO2

12

1.75-4.4

N2+O2→N2O2

4

3.1-4.0

N2O2→N2+O2

b

Oligomerization

c

Generate NO2 10

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 27

Shock Frequen velocity

Time (ps)

Elementary reaction

17

2.05-5.95

HN2O→N2+OH

21

2.45-5.9

N2+OH→HN2O

12

2.2-5.85

N2+NO2→N3O2

15

3.0-5.85

N3O2→N2+NO2

4

4.9-5.95

CN2O2→N2+CO2

2

5.55-5.95

N2+CO2→CN2O2

11

0-0.05

C6H6N12O12+C6H6N12O12→C12H12N24O24

6

0.3-1.3

NO2+NO2→N2O4

3

0.8-1.9

N2O2→N2+O2

2

0.85-1.95

N2+O2→N2O2

87

1.1-5.95

N2+OH→HN2O

86

1.35-5.95

HN2O→N2+OH

3

1.25-1.5

N2+NO2→N3O2

44

2.15-5.95

N2+CO2→CN2O2

45

2.4-5.95

CN2O2→N2+CO2

19

2.6-5.9

N2+H2O→H2N2O

24

2.8-5.9

H2N2O→N2+H2O

14

2.75-5.75

N2+H2N2O→H2N4O

8

4.15-5.65

H2N4O→N2+H2N2O

Illustration

-cies (km/s)

10

Generate N2, OH

Generate CO2

11

b

Oligomerization Generate NO2

Generate N2, OH

Generate CO2

Generate H2O

Generate N2 a

The frequencies of the reactions excluded from the table are less than 10.

b

Two or more (less than four) CL-20 molecules can oligomerize to generate CL-20 oligomers,

such as C42H42N84O84. Only C6H6N12O12 + C6H6N12O12 → C12H12N24O24 is listed in the table to represent the oligomerization reaction of CL-20. c

In addition to single CL-20 molecules, CL-20 oligomers can also decompose and generate

NO2. This type of reaction path is described by C6nH6nN12nO12n → C6nH6nN12n−1O12n−2 + NO2, where n is an integer (n ≥ 1).

From Table 2, under different shock velocities, the reaction sequence of

ACS Paragon Plus Environment

Page 11 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

CL-20

is

almost

the

same,

which

involves

oligomerization

and

redecomposition of CL-20, generation of NO2, generation of N2, generation of OH, and generation of CO2. As the shock velocity increases, the time when each reaction occurs is advanced and the reaction frequency increases. However, most of the reactions in Table 2 simultaneously occur in both reaction directions, which makes it difficult to determine the accurate reaction trend. Therefore, it is necessary to combine molecular snapshots for further analysis. Nevertheless, if the number of atoms is too large, many reactions will occur at the same time in the whole system. To easily distinguish the key reaction, only the information concerning the positions and bonding in a unit cell of the system was tracked. Because there are only four CL-20 molecules in the unit cell, the difficulties in observation are greatly reduced. Snapshots of a CL-20 unit cell at different times under 9 and 11 km/s shock are shown in Figure 3. The CL-20 molecule deforms and decomposes under shock. Higher shock velocity results in more severe deformation and destruction. At 5 ps, the system under 9 km/s shock generates NO2 and releases some free oxygen atoms. More small molecule products are produced in the system at 11 km/s, including a large amount of N2, and the N atoms in N2 may originate from the same CL-20 molecule. In conjunction with Table 2, the generation paths of five typical small molecule products and the previous reactions determined from these snapshots are shown below the snapshots. We found that successive generation of NO2, N2, OH (H2O), and CO2 is caused by the following bond cleavage sequence of CL-20 molecules: N–N, N–O, C–N, C–H, and C–C bonds. The details of these paths are as follows:

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1) Formation of NO2 and free O atoms results from cleavage of the N–N and N–O bonds. The CL-20 molecule has a cage structure and the N atoms are each connected to a nitro group by a N–N bond (a–c in Figure 3). When the molecule is shocked, the N–NO2 bond is broken, and thus NO2 forms. In addition, cleavage of the N–O bond of the nitro group generates more free O atoms. 2) After formation of NO2, two C–N bonds on the five-membered ring are successively broken, the cage structure of CL-20 is destroyed, and N2O is released (d–g). N2O then combines with a free O atom to form N2O2, which decomposes to N2 and O2, indicating that N2 can be generated from one CL-20 molecule. 3) The C atoms in the CL-20 cage structure are connected to H atoms by C–H bonds. After collapse of the cage structure, the C–H bonds begin to break (h–k). As a result, H dissociates and then combines with the nearest free O atom to form OH. Subsequently, OH reacts with N2 in the system to form HN2O. HN2O can further decompose to N2 and OH. 4) Next, formation of CO2 occurs (i–o). Interestingly, after the H atom detaches, an O atom replaces the H atom and begins to bind to the C atom. This results in formation of a C–O bond, and a CO molecule then forms after the C–C bond breaks. The CO molecule then reacts with a nearby N2O molecule by removing an O atom to form CO2. Therefore, the reaction path is CN2O2 → N2 + CO2. 5) While generating CO2, as the reaction progresses, the numbers of free H atoms and N2 in the system gradually increase (p–s). Therefore, the previously formed OH can combine with a free H atom to form H2O. H2O can

ACS Paragon Plus Environment

Page 12 of 27

Page 13 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

react with N2 and generate unstable H2N2O. For this reason, N2 + H2O → H2N2O and H2N2O → N2 + H2O are present in Table 2.

Figure 3. Snapshots of the CL-20 unit cell of the system under shock (top) and the generation paths of five typical small molecule products (bottom). C, H, O, and N atoms are colored gray, white, red, and blue, respectively. These

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

representations are considered in the following figures. 3.4 Bond cleavage and the main products The relative bond change, which is the ratio of the number of changing bonds to the number of initial bonds, was then investigated. A positive value indicates bond formation and a negative value indicates bond cleavage. Time evolution of the relative bond changes of the C–H, C–N, N–N, and N–O bonds and the numbers of the main product molecules (NO2, N2, H2O, and CO2) under 9 and 11 km/s shock are shown in Figure 4. The results for 8 and 10 km/s shock are shown in Figure S2. In Figures 4 and S2, the chemical bonds and numbers of product molecules reach equilibrium in shorter time at higher shock velocity (8 km/s, >50 ps; 9 km/s, 31 ps; 10 km/s, 16 ps; 11 km/s, 9 ps). There is a high correlation between the bond changes and the numbers of product molecules. Initial oligomerization of CL-20 results in rapid formation of N–O bonds. The N–N bonds then cleave, which generates a large number of NO2 molecules. However, the number of NO2 molecules then decreases and finally reaches zero, indicating that NO2 is an intermediate product. More details about the intermediate products found in the reaction are given in Figure S3. NO2 is the first intermediate product formed, followed by NO3, O2, and HNO3. Moreover, NO2 has the longest existence time. The N–O bond begins to break after a short time of aggregation and the degree of cleavage is approximately 60%, indicating that there are many N–O compounds in the system. The relative bond changes in the CL-20 unit cell with time under 9 km/s shock are shown in Figure 5. In this figure, only the bond relationship of the atoms contained in the unit cell was calculated. The results show that the N–O bonds are completely cleaved, which indicates that after the N–O bonds

ACS Paragon Plus Environment

Page 14 of 27

Page 15 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

are broken, the related atoms in the unit cell combine with the N or O atoms of surrounding unit cells. Because cleavage of the N–O bond mainly produces the O atom, the above results indicate that the O atom has high mobility during the reaction. The cleavage degree of the N–N bonds in the unit cell is 46%. As shown in Figure 3, a single CL-20 molecule can generate N2 without destroying the N–N bond, so the damaging effect of the shock condition on the N–N bond of CL-20 is not serious. In contrast to the N–O and N–N bonds, the C–N and C–H bonds are irreversibly and essentially completely broken. Higher shock velocity results in higher cleavage rate. The variation in the N2 and H2O molecules is consistent with their breaking rules. As the C–N and C–H bonds break, the main final products N2 and H2O successively form, followed by generation of CO2. Generation of these products also results in an increase in the temperature, as shown in Figure 2.



Figure 4. Time evolution of the relative bond changes of the C–H, C–N, N–N, and N–O bonds and the numbers of NO2, N2, H2O, and CO2 molecules.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 27

Figure 5. Relative bond changes in the CL-20 unit cell under 9 km/s shock.

From Figure 4, as the shock velocity increases, the cleavage rates of C–N and C–H bonds increase, so the formation rates of the detonation products also increase. The main difference between the detonation reactions and other reactions, such as combustion, is that the explosive can rapidly react and release energy in a very short time. Therefore, we conclude that the rapid breaking of the chemical bonds induced by shock is the key factor for explosive detonation. The cleavage rates of the C–N and C–H bonds with different shock velocities are given in Table 3. The results indicate that the detonation reaction occurs in CL-20 when the cleavage rates of the C–N and C–H bonds are greater than 3.11%/ps and 4.15%/ps (results of 9 km/s shock), respectively. Table 3. Cleavage rates of the C–N and C–H bonds Shock velocity (km/s) 8

C-N bond

C-H bond

Relative bond change(%)

Cleavage time (ps)

54.3

49.95

a

Cleavage rate (%/ps)

Relative bond change (%)

Cleavage time (ps)

1.09

51.3

49.95

ACS Paragon Plus Environment

a

Cleavage rate (%/ps) 1.03

Page 17 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Shock velocity (km/s)

C-N bond

C-H bond

Relative bond change(%)

Cleavage time (ps)

Cleavage rate (%/ps)

Relative bond change (%)

Cleavage time (ps)

Cleavage rate (%/ps)

9

92.5

29.75

3.11

94.0

22.65

4.15

10

95.0

12.35

7.69

97.1

8.35

11.63

11

92.6

6.5

14.24

95.4

3.75

25.43

a

The C–N and C–H bonds in the system under 8 km/s shock are still breaking at 50 ps, as

shown in Figure S2, so we only give the ratios of the relative bond changes and calculation time.

Interestingly, the increase in the reaction rate resulting from the increasing shock strength is not sustainable. After the C–N and C–H bonds are completely broken, the product formation rate of the system decreases and the number of molecules tends to be stable. In addition, the total number of products decreases with increasing shock velocity. However, the numbers of the final products in Figure 4 are much lower than those of the thermal equilibrium state (Table S1), indicating that the reactions are not complete. By analyzing the species in the system, we found that in addition to the small molecule products, there are many clusters, which may reduce the reaction rate of the system. The numbers of C, H, N, and O atoms contained in the clusters at 50 ps and the ratios of these numbers to the total number of atoms are given in Table 4. These clusters bind a large number of atoms with the order C > O > N > H. In particular, more than half of the C and O atoms are trapped. A larger number of atoms remain in the clusters with higher impact velocity. This behavior may result from the space constraints at high pressure, as shown in Figure 2, which facilitates bonding of the atoms with nearby atoms to aggregate and form clusters. However, this does not represent the termination state of the explosive reaction because in the real detonation

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 27

reaction, the reaction pressure will gradually decrease and it will not be maintained at a high-pressure state of several tens of gigapascals. Table 4. Numbers of C, H, N, and O atoms contained in the clusters and the ratios of these numbers to the total number of corresponding atoms ε a Shock velocity 0ps

50ps

ε (%)

C

864

550

63.7

H

864

273

31.6

N

1728

609

35.2

O

1728

974

56.4

C

864

655

75.8

H

864

377

43.7

N

1728

860

49.8

O

1728

1162

67.2

C

864

708

81.9

H

864

424

49.1

N

1728

941

54.5

O

1728

1250

72.4

(km/s)

9

10

11

a

A molecule with a greater molecular weight than CL-20 is considered to be a cluster.

3.5 Cluster stability and product variation after expansion In our previous experimental research33, we found that the reaction time of CL-20 (high CL-20 content mixed explosive) is about 48 ns, which is too long for MD simulation. To reduce the calculation time, the cell length in the y direction was increased at a rate of 0.02 km/s to simulate the subsequent reaction of the system during rapid expansion. This calculation is only a qualitative analysis of the cluster stability and product variation. The calculation was stopped when the density of the calculation system decreased

ACS Paragon Plus Environment

Page 19 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

to 1 g/cm3. The maximum molecular weight of the molecules in the system as a function of the density during expansion is shown in Figure 6. As the density decreases, the maximum molecular weight in the system rapidly decreases, indicating rapid decomposition of large clusters. Most of the N and H atoms escape from the clusters, so the clusters mainly consist of C and O atoms in the last stage. When the density increases to 1.48 g/cm3 (9 km/s), 1.73 g/cm3 (10 km/s), and 1.67 g/cm3 (11 km/s), there are no clusters larger than CL-20 in the system.

Figure 6. Maximum molecular weight of the molecules as a function of the density.

The change in the number of typical small molecule products versus the density is shown in Figure 7. With rapid expansion of the system, the numbers of small molecular products, such as N2, H2O, and CO2, rapidly increase. CO begins to appear when the density decreases to 2 g/cm3, as shown in Figure

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

S4, indicating that CO is a late stage product in the CL-20 detonation reaction. At the end of the reaction, the final stable product generated from the C atoms is CO2 or CO, so evolution of the number of stable C atoms (the sum of the number of CO2 and CO molecules) is also given in Figure 7. Interestingly, the change rules of these late products at different shock velocities are similar. The formation rate is not significantly affected by the shock velocity, but it is affected by the decrease in the density. This indicates that after complete destruction of the molecule structure, the late reaction rate is limited and mainly based on the characteristics of the explosive. Therefore, a shock higher than the detonation velocity may not be supported.

Figure 7. Numbers of small molecule products as a function of the density under different shock velocities. 4. CONCLUSIONS

ACS Paragon Plus Environment

Page 20 of 27

Page 21 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

We have attempted to determine the connections between the changes in the chemical bonds and the macroscopic detonation laws of CL-20 explosive by ReaxFF MD simulations. We found that there is a strong correlation between shock-induced bond cleavage and the detonation reaction. CL-20 has the same breaking order under shock along the [010] orientation in the range 8–11 km/s: the N–NO2 bond is cleaved first followed by the C–N, C–H, and C–C bonds. Among these bonds, rapid and irreversible cleavage of the C–N and C–H bonds induced by shock is the key factor for detonation. It significantly affects formation of N2 and H2O, which produces a large amount of heat and increases the temperature of the system. The detonation reaction of CL-20 occurs when the cleavage rates of the C–N and C–H bonds are greater than 3.11%/ps and 4.15%/ps, respectively. Higher shock velocity results in higher cleavage rates of these bonds and higher generation rates of the relative products. However, after the molecule structure is completely destroyed, the high pressure makes the atoms aggregate to form clusters and therefore slows down the reaction. More atoms are trapped in clusters at earlier time with higher shock velocity. Although the clusters will decompose as the volume expands, the reaction rate is mainly related to the density and it is not significantly affected by the shock velocity. This means that the late detonation reaction is mainly based on the characteristics of the explosive. This study reveals the detonation mechanism of the CL-20 explosive and provides a new way to understand the classical detonation model at the microscopic level.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ASSOCIATED CONTENT Supporting Information Calculation details, Figures S1 – S4, and Table S1 can be found in the Supporting Information.

AUTHOR INFORMATION ORCID

Danyang Liu: 0000-0002-8990-8511 Lang Chen: 0000-0003-2623-1061 Notes

The authors declare no competing financial interest.

ACKNOWLEDGMENTS This work was supported by National Natural Science Foundation of China (Grant No. 11832006). L.D.Y. appreciates the calculation guide by Dr. Wang Fu Ping and Wang He Qi.

ACS Paragon Plus Environment

Page 22 of 27

Page 23 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

REFERENCES

(1) Von Neuman J. Theory of Detonation Waves; Institute for Advanced Study: Princeton, New Jersey, 1942. (2) Shepherd, J.E.; Thibault, P. Review of "The Detonation Phenomenon". AIAA J. 2009, 47, 1310−1311. (3) Van Duin, A. C.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A Reactive Force Field for hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (4) Strachan, A.; Van Duin, A. C.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A. Shock Waves in High-energy Materials: The Initial Chemical Events in Nitramine RDX. Phys. Rev. Lett. 2003, 91, 098301. (5) Senftle, T. P.; Hong, S.; Islam, M. M.; Kylasa, S. B.; Zheng, Y.; Shin, Y. K.; Junkermeier, C.; Engel Herbert, R.; Janik, M.J.; Aktulga, H. M.; et al. The ReaxFF Reactive Force-Field: Development, Applications and Future Directions. NPJ Comput. Mater. 2016, 2, 15011. (6) Chaban, V. V.; Fileti, E. E.; Prezhdo, O.V. Buckybomb: Reactive Molecular Dynamics Simulation. J. Phys. Chem. Lett. 2015, 6, 913−917. (7) Fileti, E.E.; Chaban, V. V.; Prezhdo, O. V. Exploding Nitromethane in Silico, in Real Time. J. Phys. Chem. Lett. 2014, 5, 3415−3420. ( 8 ) Wood, M. A.; Van Duin, A. C.; Strachan, A. Coupled Thermal and Electromagnetic Induced Decomposition in the Molecular Explosive αHMX; A Reactive Molecular Dynamics study. J. Phys. Chem. A 2014, 118, 885−895. (9) Rom, N.; Hirshberg, B.; Zeiri, Y.; Furman, D.; Zybin, S. V.; Goddard, W. A.; Kosloff, R. First-Principles-Based Reaction Inetics for Decomposition of Hot,

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Dense Liquid TNT from ReaxFF Multiscale Reactive Dynamics Simulations. J. Phys. Chem. C 2013, 117, 21043−21054. (10) Zhang, L.; Chen, L.; Wang, C.; Wu, J. Y. Molecular Dynamics Study of the Effect of H2O on the Thermal Decomposition of α Phase CL-20. Acta Phys. Chim. Sin. 2013, 29, 1145−1153. (11) Zhou, T. T.; Huang, F. L. Effects of Defects on Thermal Decomposition of HMX via ReaxFF Molecular Dynamics Simulations. J. Phys. Chem. B 2011, 115, 278−287. (12) Guo, F.; Zhang, H.; Cheng, X. Molecular Dynamic Simulations of Solid Nitromethane under High Pressure. J. Theor. Comput. Chem. 2010, 9, 315−325. (13) Budzien, J.; Thompson, A. P.; Zybin, S. V. Reactive Molecular Dynamics Simulations of Shock Through a Single Crystal of Pentaerythritol Tetranitrate. J. Phys. Chem. B 2009, 113, 13142−13151. (14) Strachan, A.; Kober, E. M.; Van Duin, A. C.; Oxgaard, J.; Goddard, W. A. Thermal Decomposition of RDX from Reactive Molecular Dynamics. J. Chem. Phys. 2005, 122, 54502. (15) Liu, L.; Liu, Y. S.; Zybin, V.; Sun, H.; Goddard, W. A. ReaxFF-lg: Correction of the ReaxFF for London Dispersion, With Applications to the Equations of State for Energetic Materials. J. Phys. Chem. A. 2011, 115, 11016−11022. ( 16 ) Wang, F.; Chen, L.; Geng, D.; Wu, J.; Lu, J.; Wang, C. Thermal Decomposition Mechanism of CL-20 at Different Temperatures by ReaxFF Reactive Molecular Dynamics Simulations. J. Phys. Chem. A 2018, 122, ACS Paragon Plus Environment

Page 24 of 27

Page 25 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry



3971−3979. (17) Wang, F.; Chen, L.; Geng, D.; Lu, J.; Wu, J. Effect of Density on Thermal Decomposition Mechanism of ε-CL-20: A ReaxFF Reactive Molecular Dynamics Simulation Study. Phys. Chem. Chem. Phys. 2018, 20, 22600−22609. (18) Yang, Z.; He, Y. H. Pyrolysis of Octanitrocubane via Molecular Dynamics Simulations. Acta Phys. Chim. Sin. 2016, 32, 921−928. (19) Guo, D.; An, Q.; Zybin, S. V.; Goddard III, W. A.; Huang, F.; Tang, B. The Co-Crystal of TNT/CL-20 Leads to Decreased Sensitivity Toward Thermal Decomposition from First Principles Based Reactive Molecular Dynamics. J. Mater. Chem. A 2015, 3, 5409−5419. (20) Liu, H.; Dong, X.; He, Y. H. Reactive Molecular Dynamics Simulations of Carbon-Containing Clusters Formation during Pyrolysis of TNT. Acta Phys. Chim. Sin. 2014, 30, 232−240. (21) Guo, D.; Zybin, S.V.; An, Q.; Goddard III, W.A.; Huang, F. Prediction of the Chapman–Jouguet Chemical Equilibrium State in a Detonation Wave from First Principles Based Reactive Molecular Dynamics. Phys. Chem. Chem. Phys. 2016, 18, 2015−2022. (22) Guo, D.; Guo, D.; Huang, F; An, Q. Influence of Silicon on the Detonation Performance of Energetic Materials from First-Principles Molecular Dynamics Simulations. J. Phys. Chem. C 2018, 122, 24481−24487. (23) Olah, G. A.; Squire, D. R. Chemistry of Energetic Materials; Academic Press: San Diego, CA, 2012. (24) Xue, X.; Wen, Y.; Zhang, C. Early Decay Mechanism of Shocked ε-CL-20:

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



A Molecular Dynamics Simulation Study. J. Phys. Chem. C 2016, 120, 21169−21177. (25) Reed, E. J.; Fried, L. E.; Joannopoulos, J. D. A Method for Tractable Dynamical Studies of Single and Double Shock Compression. Phys. Rev. Lett. 2003, 90, 235503. (26) Manaa, M. R.; Reed, E. J.; Fried, L. E.; Goldman, N. Nitrogen-Rich Heterocycles as Reactivity Retardants in Shocked Insensitive Explosives. J. Am. Chem. Soc. 2009, 131, 5483−5487. (27) Xue, X.; Wen, Y.; Long, X.; Li, J.; Zhang, C. Influence of Dislocations on the Shock Sensitivity of RDX: Molecular Dynamics Simulations by Reactive Force Field. J. Phys. Chem. C 2015, 119, 13735−13742. (28) Wen, Y.; Xue, X.; Zhou, X.; Guo, F.; Long, X.; Zhou, Y.; Li, H.; Zhang, C. Twin Induced Sensitivity Enhancement of HMX versus Shock: A Molecular Reactive Force Field Simulation. J. Phys. Chem. C 2013, 117, 24368−24374. (29) Shan, T. R.; Wixom, R.; Mattsson, A.; Thompson, A. Atomistic Simulation of Orientation Dependence in Shock-Induced Initiation of Pentaerythritol Tetranitrate (PETN). J. Phys. Chem. B 2013, 117, 928−936. (30) Woińska, M.; Grabowsky, S.; Dominiak, P. M.; Woźniak, K.; Jayatilaka, D. Hydrogen atoms can be located accurately and precisely by x-ray crystallography. Sci. Adv. 2016, 2, 1600192. (31) Mader, C. L. Numerical Modeling of Explosives and Propellants; CRC Press: Boca Raton, FL, 2008.

ACS Paragon Plus Environment

Page 26 of 27

Page 27 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry



(32) Chen, L.; Pi, Z. D.; Liu, D. Y.; Yang, K.; Wu, J. Y. Shock Initiation of the CL-20-Based Explosive C-1 Measured with Embedded Electromagnetic Particle Velocity Gauges. Propellants Explos. Pyrotech. 2016, 41, 1060−1069. (33) Liu, D.; Chen, L.; Wang, C.; Wu, J. Detonation Reaction Characteristics for CL-20 and CL-20-based Aluminized Mixed Explosives. Cent. Eur. J. Energ. Mater. 2017, 14, 573−588.

TOC GRAPHIC

ACS Paragon Plus Environment