Reactive molecular dynamics study on the thermal decomposition reaction of a triple-base solid propellant

The study of the combustion property of newly designed propellant by means of computational simulation is an efficient pathway for assessment and could avoid exposure to hazardous chemicals. An RDX-modified triple-base solid propellant formula was proposed in this study. Reactive molecular dynamics simulations employing ReaxFF-lg force field were performed to explore the thermal decomposition property of the propellant for a variety of temperatures. The reaction kinetics of the system and major ingredients were analyzed, and the apparent decomposition activation energies were calculated. The population of decomposition intermediates and products is thoroughly investigated. H2O is consumed at high temperatures indicating a water–gas reaction that could reduce carbon clusters during the combustion of solid propellant. The water–gas reaction, as well as the population of H2 at high temperature, points out the way of adjusting the formula of the propellant, which is adding fuel and oxidizer to improve combustion temperature and oxygen balance.


Introduction
Solid propellants are a kind of special energy materials extensively used in the military, commercial, and space industries. Compared with liquid propellants, solid rocket motor propellants are widely accepted because of the fact that they are relatively simple to formulate and use, and have excellent storage reliability, and could be transported conveniently. Typical solid propellants are usually formulated with an oxidizing agent, and a binder which may also act as the fuel. Design and assessment of new solid propellants with good performance is the need for the development of munitions and space exploration.
Generally, solid propellants can be classified as being either homogeneous or composite. The homogeneous propellants are further subclassified as being either single base, double base, or triple base. Single-base (SB) propellant consists of nitrocellulose (NC) as the only explosive ingredient. Double-base (DB) propellants comprising NC and an energetic plasticizer such as nitroglycerin (NG) are a kind of important solid propellants. Nitroguanidine (NQ) is a low sensitive energetic compound with a detonation velocity of 7650 m s −1 . When a propellant is made of NC, NG, and NQ, it is termed a triple-base (TB) propellant. TB propellant has the advantage of reduced muzzle flash compared with DB propellant. The introduction of NQ to a propellant composition also results in the reduction of the combustion temperature. When it is used as gun powder, the erosion of the gun will be reduced. Usually, the proportion of NQ is controlled at about 50 − 55% to adjust the gas output, energy, temperature, and burning rate [1,2].
The performance of solid propellant systems can be improved by adding other oxidizers such as ammonium perchlorate and RDX. The addition of oxidizers will result in a higher burning rate, a higher specific impulse, and improved mechanical properties of the propellant composition. RDX is probably the most important high explosive due to its high chemical and thermal stability, high density, and high detonation velocity. Its performance properties are only slightly inferior to those of the HMX. However, RDX is cheap and easily manufactured. Therefore, replacing partial NQ of TB propellant with RDX will result in a new propellant with improved energy.
Understanding the combustion process of solid propellants is important for the design and assessment the performance. However, the combustion of solid-state propellant is an extreme chemical reaction with high pressure and high temperature. The exceedingly complicated reactions take place in the solid, liquid, and gas phases of heterogeneous mixtures. Therefore, it is hard to study the decomposition mechanism of this kind of chemicals experimentally. Benefitting from the development of computational techniques, it is possible to simulate the combustion process of solid-state propellant with large amounts of atoms.
In present technology conditions, quantum chemical methods could handle a system with hundreds of atoms. However, molecular dynamic methods could simulate a system with up to millions of atoms on a supercomputing cluster. The empirical force fields cannot deal with bond breaking and bond formation of chemical reactions. Consequently, the reactive force fields (ReaxFF) are developed to study the chemical reactions that may exist. ReaxFF is based on bond order which has a nearly linear relationship with bond length. A log-log plot of bond order versus dissociation energies is also almost linear. Therefore, ReaxFF correlates with bond order, bond length, and bond energy. The total energy in ReaxFF is as follows: where E bond , E val , E dihe , E vdW , E Coul , E over , E under , E pen , and E conj represent bond energy, valence-angle energy, torsion angle energy, Van der Waals energy, Coulomb interactions energy, over-coordination energy penalty, resonance energies between π − electrons, stability of systems with two double bonds sharing an atom in a valence angle, and conjugation energy. ReaxFF has been successfully applied to simulate the reaction kinetics of chemical systems including combustion of hydrocarbons [3,4], decomposition of energetic materials [5][6][7][8], oxidation of metal nanoparticles [9], and growth of metal nanoclusters [10].
Most of the ReaxFF reactive force fields were developed based on the quantum calculations at B3LYP/6-31G** level. However, these calculations did not account adequately for the London dispersion which is very important in molecular solids. It results in overestimating equilibrium volumes by about 10 to 15%. Therefore, Liu [11] and Larentzos [12] have developed ReaxFF-lg force fields for organic solids containing C, H, N, and O elements, respectively. Later, Ju et al. [13] developed a new ReaxFF force field consisting of C, H, N, O, and Al parameters. These force fields could accurately predict the crystal structures of high-density solids, especially for energetic materials. They are extensively applied for the simulation of high energy density materials [14][15][16].

Computational methods
In this study, the reactive molecular dynamic (RMD) simulations were performed using the large-scale atomic molecular massively parallel simulator (LAMMPS) [17] with ReaxFFlg potential developed by Liu et al. under periodic boundary conditions (PBCs). A cubic periodic box with a size of 4.7 × 4.7 × 4.7 nm was built. An oligomer of NC with six pyran rings is chosen as the molecular model of NC to simplify the modeling. To ensure the mechanical properties of the propellant, the total amount of NC and NG is controlled at around 65%. Half of the NQ in the TB propellant is replaced with the more energetic compound RDX. The mass percentages of NC, NG, NQ, and RDX in the propellant are 27.5%, 26.9%, 22.8%, and 22.8%, respectively. All molecules are filled randomly into the cubic box. The model has a density of 1.70 g cm −3 which is close to that of a real solid propellant as shown in Fig. 1. The whole system contains 10,317 atoms.
The time step of the RMD simulations was set to 0.1 fs for an accurate description of the violent reaction. The structure was relaxed by energy minimization with the conjugation gradient algorithm at 0 K. The molecules were further equilibrated with an NPT ensemble at 300 K for 10 ps. The size of the box changed to 4.846 × 4.846 × 4.846 nm after structural relaxation corresponding to a density of 1.55 g cm −3 . RMD simulations of the relaxed periodic structure were heated from 300 to 1000, 1500, 2000, 2500, 3000, and 3500 K with NVT ensemble for 200 ps, respectively. Nose-Hoover thermostat was used to maintain a constant temperature during simulations, and the dump interval was set to 100 fs. The bond order was calculated using the ReaxFF bond-order utility with a dump interval of 100 fs. The molecular species were analyzed with a script written in C language. The choice of bond-order cut-off may alter the molecular composition of some intermediates and final products. If the default bond-order cutoff value Fig. 1 The initial structure of the periodic box. Color scheme: nitrogen -blue, oxygen -red, carbon -gray, hydrogen -white (0.3 in the LAMMPS suite) is used, chemicals with unreasonable bonds will be observed [18]. Therefore, the bond-order cutoff values are set according to pairs of atoms and listed in Table 1.

Reaction kinetic analysis
The time evolution of potential energy (PE) for temperatures from 1000 to 3500 K is shown in Fig. 1a (black lines). All curves increase rapidly at the earliest stages indicating an endothermic process in the initial period. Therefore, the maximum value of PE increases with temperature. The curves decrease after PE reaches maximum due to the release of heat during the decomposition of energetic materials. PE declines more rapidly when temperature increases. PE curves reach a plateau earlier at higher temperatures. It takes 200, 100, and 45 ps to equilibrate for reaction temperatures at 2500, 3000, and 3500 K. The equilibration times at lower temperatures are longer than 200 ps.
The solid propellant studied in this research contains four ingredients. The expression of the reaction kinetics by the decay of molecules is impracticable. Alejandro Strachan et al. proposed a simple exponential function to describe the evolution of the reaction [19]: where U is the potential energy, U 0 is the asymptotic energy of the products, ΔU is the exothermicity of the reaction, t E−I is the initial equilibration and induction time, and τ is the average characteristic time (i.e., inverse reaction rate, 1/k) of the energy-releasing reactions. The Arrhenius equation (Eq. 2) is employed to describe the relationship between the reaction rate constant and the apparent activation energy  Table 2. Ln(1/τ) versus inverse temperature should be linear relation if the reaction could be described by the Arrhenius equation (Eq. 2). The plot of ln(1/τ) versus inverse temperature shows a turning point at 2000 K indicating two different reaction mechanisms at temperatures below and above 2000 K (Fig. 2b). The red curve in Fig. 2b is the linear fit of the data for temperatures from 2000 to 3500 K. Apparent activation energy (E a ) of the decomposition of the propellant derived from the slope of the fitted curve is 87.69 kcal mol −1 .

Decay kinetics of RDX, NG, and NQ
The propellant studied in this research is composed of NG, NC, NQ, and RDX. The model contains 16 NC molecules which are insufficient for the kinetics analysis. Wherefore, the decomposition of RDX, NG, and NQ is analyzed. The times taken for the decompose of all RDX (t d ) are listed in Table 3. RDX molecules do not decompose completely within 200 ps at 1000 K (Fig. 3a). t d is shortened to 70 ps at 1500 K and to only 0.4 ps at 3500 K. The decay kinetics of RDX could be expressed by the first-order rate Eq. 3, which can be converted to Eq. 4.
where N 0 , N t , t, and k are the initial number of RDX, RDX amount at t (ps), decay time, and decay rate constant. The decay rate constant k could be obtained through the linear fitting of lnN t versus t. As shown in Fig. 3b, the decay of RDX follows a good linear relationship. The decay rate constants and lnN 0 derived from linear fittings at different temperatures are listed in Table 3. The values of lnN 0 are close to each other at different temperatures indicating the results are reliable. The further linear fitting of lnk versus inverse temperature could obtain the apparent activation energy of the decomposition of RDX (114.8 kcal mol −1 ; Fig. 3c).
The decay kinetics of NG and NQ could be fitted through the same procedure. The respective values of E a are 100.4 and 127.9 kcal mol −1 . The relatively low apparent activation energy of the whole system is probable due to the early decomposition of NC.

Evolution of the decomposition intermediates and products
The time evolutions of total species for temperatures from 1000 to 3500 K are plotted in Fig. 3d. The reaction propagates quickly in the first 5 ps resulting in the fast increase of the population of intermediates and products. The number of decomposed fragments increases slowly at 1000 K, and the  total amount is less than 1000 at 200 ps. It increases rapidly at temperatures higher than 1500 K. The maximum population at 200 ps is about 3400 for temperatures higher than 2500 K. It takes about 200, 100, and 50 ps to reach maximum species for temperatures at 2500, 3000, and 3500 K indicating the reaction reaches equilibrium. This is in accordance with the evolution of PE. Figure 4 shows the evolution of intermediates and stable gaseous products for temperatures from 1000 to 3500 K. At low temperatures, the formation of abundant unstable intermediates is observed. Note that NO 2 is the major intermediate at all simulated temperatures which comes from the decomposition of NG, NC, and NQ. It is well-known that the cleavage of N − NO 2 and O − NO 2 bonds is normally the trigger of the decomposition of this kind of energetic materials due to the low bond dissociation energies [5,20,21]. NO 2 captures proton forming HNO 2 which could be found at temperatures from 1000 to 2500 K. NO 3 and HNO 3 are also products of NO 2 observed at low temperatures.
NO is another important intermediate at high temperatures which appears shortly after NO 2 (Fig. 4c-f). It is reported that NO 2 will rearrange to nitrite (− ONO) which undergoes bond cleavage to form NO [22]. This means the rearrangement is preferred at high temperatures. HNO appears after NO indicating it comes from NO who captures a proton from an adjacent fragment.
With the consumption of unstable intermediates, the population of stable gaseous molecules increases. The stable reaction products include H 2  To figure out the population of reaction products, the time evolution of four major gases H 2 O, CO 2 , H 2 , and N 2 for temperatures from 1000 to 3500 K is shown in Fig. 5. The evolution of H 2 O shows a one-way increasing trend for temperatures from 1000 to 2000 K within 200 ps. At higher temperatures, it increases first and then decreases. The maximum reaches earlier at higher temperatures. The populations of CO 2 and H 2 are similar which increase at the temperature range of 1000 to 3000 K. At 3500 K, the amounts of these two gases grow and then reach a plateau. We note that the amount of H 2 and CO 2 increase during the consumption of H 2 O, especially at 3000 K in the time range of 50 to 150 ps. After the decay curves of H 2 O flatten out, the amount of H 2 and CO 2 do not increase (temperature and time regions: 3000 K, 150 − 200 ps, and 3500 K, 75 − 200 ps). It indicates a water-gas reaction between H 2 O and carbon clusters derived from the decomposition of the propellant (Eq. 5) [23]. Carbon clusters are observed during the decomposition of some energetic compounds [7,19,22]. In this study, clusters with carbon atoms of 24, 28, 31, 33, 38, and 44 are observed but the number of total clusters is small. These clusters are consumed during the reaction proceeds. Therefore, a higher reaction temperature is beneficial to reduce carbon clusters.
In contrast, the population of inert gas N 2 does not show dependency on the evolution of H 2 O. An RMD study of the decomposition of RDX showed that the population of CO and CO 2 is density-dependent. Essentially, all C atoms form CO molecules for the low-density crystal model while it forms almost no CO for a density of 2.11 g cm −3 [19]. In this study, CO is not a major reaction product at any reaction temperature. Therefore, the water-gas reaction is regarded as Eq. 5. This reaction points out the way to reduce carbon clusters during the combustion of solid propellant. However, the massive presence of fuel gas H 2 in the reaction products indicates the oxygen balance should be further improved.

Conclusions
In this study, we proposed a new triple-base solid propellant formula consisting of NC, NG, NQ, and RDX. The decomposition reactions of the propellant at different temperatures were investigated with reactive molecular dynamics simulations. The reaction kinetic analysis and decay kinetics of RDX, NG, and NQ were investigated. The results reveal that the apparent activation energy of (5) 2H 2 O + C = CO 2 + 2H 2 Fig. 5 Time evolution of a H 2 O, b CO 2 , c H 2 , and d N 2 for temperatures from 1000 to 3500 K the whole system is close to that of three major ingredients. The evolution of decomposition fragments shows the unstable intermediates propagate at low temperatures. They degrade to stable gaseous products H 2 O, N 2 , CO 2 , and H 2 . CO is very few to be found in the products. The time evolution curves of major decomposition products suggest a water-gas reaction at high temperatures between H 2 O and carbon clusters. We infer that adding fuel such as powder aluminum to improve combustion temperature may be a way to reduce carbon clusters and promote complete combustion. The population of H 2 at high temperatures is massive. It indicates the oxygen balance should be further improved through the addition of an oxidizer such as ammonium perchlorate to further improve the energy output of the propellant. This research has shown that reactive molecular dynamics simulation could be a convenient and useful method for the design of solid propellants. The information obtained from reactive molecular dynamics simulations could guide the rational regulation of solid propellant formulations.