Hot Spot Interaction With HTPB Binder in Energetic Composites - The

Subscriber access provided by - Access paid by the | UCSB Libraries

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Hot Spot Interaction With HTPB Binder in Energetic Composites Kaushik L. Joshi, and Santanu Chaudhuri J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b11155 • Publication Date (Web): 02 Apr 2018 Downloaded from http://pubs.acs.org on April 2, 2018

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 41 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

Hot Spot Interaction with HTPB Binder in Energetic Composites Kaushik Joshia, Santanu Chaudhuria,b,* a

Illinois Applied Research Institute, University of Illinois, Urbana-Champaign, Champaign, IL 61820 b

Department of Civil and Materials Engineering, University of Illinois at Chicago, IL 60607

Abstract Binders play a pivotal role in solid propellants and polymer-bonded energetic composites in providing the means for safe transport, storage, and use. The fundamental role of binders is known qualitatively, but the exact role in modification of chemistry and thermal transport is hard to quantify. Thermochemical response of a commonly used binder, hydroxyl-terminated polybutadiene (HTPB), in a model composite comprising of cyclotrimethylene trinitramine (RDX) in presence of a thermal hot spot was studied using Reactive molecular dynamics (RMD) simulations. The results from model interfaces between RDX and HTPB were analyzed to identify the unique role played by the binder phase. It was found that the binder phase provides a safety buffer by extending the time required for the hot spot to undergo ignition to deflagration transition. This delay is mainly due to the disparity in heat capacities of HTPB and RDX. HTPB, which has higher heat capacity than RDX, absorbs energy from hot spot at a rate that is at least four times higher than the corresponding rate for RDX. RMD trajectory was mapped onto 1 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

Eulerian control volumes to capture energy flux across a moving deflagration front. It was found that a deflagration wave created in RDX moves in HTPB at least three times slower than RDX. Due to deflagration in RDX, HTPB is exposed to high temperatures (< 2000K) and high pressures (< 10 GPa). At such high temperatures and pressures, thermal degradation of confined HTPB primarily occurs through chain branching bimolecular reactions. The primary products of thermal degradation include heavy branched polymers (heavier than HTPB), cyclic molecules like cyclohexene, cyclohexadiene, cyclopentadiene and linear molecules like butadiene and trans-butadiene oligomers. The results show that RMD simulations are capable of providing an in-depth understanding and qualitative picture of the events after thermal initiation in PBX. There is still scope for enhancing quantitative accuracy of RMD predictions by further improving the force field parameters. Thus RMD can play significant role in understanding the role of binders in intergranular regions of PBXs.

2 ACS Paragon Plus Environment

Page 2 of 41

Page 3 of 41 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. Introduction Safe handling, storage and transportation of energetic materials (EM) and solid propellants are enabled by polymeric additives called binders1. The energetic composites thus formed, which are also referred as polymer bonded explosives (PBXs), are often subjected to extreme environments or insults like high temperatures, high pressures and shock waves. Understanding the response of a PBX to such insults is important for military, propulsion and mining technology. Depending upon the desired mechanical formability, insensitivity, and ageing properties, binders can be fluoropolymers, elastomers or even energetic polymers. In general, polyurethanes are most widely used binder. Elastomers such as hydroxyl-terminated polybutadiene (HTPB) mixed with plasticizers are preferred for mechanically sensitive explosives like cyclotrimethylene trinitramine (RDX) and Octahydro-1,3,5,7-tetranitro-1,3,5,7tetrazocine (HMX)2-4. The commonly used HTPB is available in three forms with average molecular weight in the range of 2800 g/mol. The plasticizers can also vary in proportion and chemical compositions. Some common forms of plasticizers often mixed with HTPB are dioctyl adipate (DOA), isodecyl pelargonate (IDP), dioctyl secabate and a range of phthalates. The PBX formulations in service today are designed largely by experience in formulation chemistry, careful sourcing of individual components, and expensive standardized tests. Thus, the ability of energetic composites to resist initiations from thermal, mechanical, thermochemical, or thermomechanical insults is key to the safety of PBX from accidental initiation. However, modeling of actual composites in extreme conditions has been a challenge due to the complexity in the chemical composition and coupling between thermo-chemical response and composite microstructure. For example, a common PBX such as PBXN-109 contains approximately 64% 3 ACS Paragon Plus Environment


Page 4 of 41

RDX, 20% Aluminum, 7.5% HTPB, 7.5% DOA, and some other components. As a result, molecular level investigations of the EM and binders treat the ingredients separately due to large molecular weight and large number of possible combinations. Mechanical changes in EMs have been studied extensively using computational methods of different length scales. At high pressures, EM can undergo different types of structural transformations such as phase change, fracture, shear banding and plastic flow. The α-γ phase transition at high pressures in RDX has been studied using Density-Functional-Theory (DFT) based molecular dynamics (MD)5, interatomic potential based molecular dynamics6 and using phase-field model7. Cawkwell et al.8 reported homogenous dislocation nucleation and anomalous plastic hardening in RDX under shock loading using molecular dynamics. Similar structural changes

have

been

reported

in

1,3,5-triamino-2,4,6-trinitrobenzene

(TATB)9

and

Dihydroxylammonium 5,5′-bis(tetrazole)-1,1′-diolate (TKX-50)10. In addition to such structural changes, extreme conditions can also initiate chemical transformations in EM if enough energy gets deposited in a localized region. It is widely believed that chemical reactions in EM are mostly initiated at hot spots. The mechanisms behind formation of such hot spots mostly depend on the crystal structure and orientation of the grains. In single crystal EMs, hot spots can form due to void collapse11-12 and shear band localization13-14. In polycrystalline materials, hot spots are more likely to occur near crystal interfaces due to cracks in microstructure15, grain-grain friction16-17 and due to conversion of plastic work into thermal energy18. A hot spot can trigger chemical reactions and it can either quench or form a self-sustaining deflagration or detonation front19-21. ReaxFF22 reactive molecular dynamics has been used for investigating the chemical transformations due to a hot spot in Ammonium Nitrate23, RDX11,

24-27

, HMX28 and

Pentaerythritol tetranitrate (PETN)29. In mesoscale simulations, different variants of dissipative 4 ACS Paragon Plus Environment

Page 5 of 41 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


particle dynamics (DPD) were used to investigate the shock to detonation transition in liquid nitromethane30 and the thermal decomposition of RDX31. However, these mesoscale methods require a priori knowledge of reaction chemistry and like RMD, both these studies were focused on single crystal EMs. Molecular studies have also been used to investigate mechanical response of different polymeric binders. Hooper et al. performed a series of non-reactive molecular dynamics simulations to obtain high pressure ( < 100 kPa) PVT data of polydimethylsiloxane (PDMS), trans hydroxyl-terminated polybutadiene (HTPB) and Estane-like polymer3. Fröhlich, Sewell, and Thompson32 found that HTPB form glass-like state under shock compression. Recently, Xie et al.33 reported that bending distortions dominate over torsional distortions in polyethylene under shock compression. Chemical response of HTPB has been mostly studied in pyrolysis experiments. It was reported that thermal degradation of HPTB occurs by a two-stage mechanism and is sensitive to heating rate34-36. In the first stage, which is net exothermic, HTPB molecules undergo depolymerization, cyclization and crosslinking. The residues formed in the first stage undergo further thermal decomposition in the second stage which is net endothermic. At low heating rates (< 102 K/min), exothermic cyclization and crosslinking reactions dominate over depolymerization and as a result, thermal decomposition occurs mainly via cyclic intermediates37. At high heating rates, endothermic chain-scission reactions are more prevalent, and trans-butadiene oligomers dominate the overall chemistry36. Tingfa also reported that apparent activation energy increases with increase in heating rate34. However, these experimental studies were done at thermodynamic conditions (temperatures < 900K, pressures < 25 atm) that are orders of magnitude lower than the P-T conditions possible in extreme environments.



In addition to modeling EM and binder separately, some previous studies have also studied multi-component PBX systems in molecular simulations. Those studies reported that hot spots are more likely to form near non-uniform interfaces between fuel and binder. For example, Shi and Brenner studied the thermomechanical response of a PBX system containing nitrogen cubane crystal and inert binder38. They found that shock focusing and local compression at the fuel-binder facets can form a hot spot. An and co-workers39 used RMD to study thermochemistry behind the formation of a hot spot in PBX under shock compression. Their PBX model consisted of RDX, HTPB and di-octyl adipate (DOA). They found that, under shock compression, a hot spot is formed at the non-uniform fuel/binder interface due to shear localization. They also reported chemical initiation of RDX due to large temperature rise caused by the hot spot. A similar mechanism of hot spot formation and subsequent chemical initiation was reported in PETN and silicon pentaerythritol tetranitrate (Si-PETN) by the same group of researchers40. For the PETN systems, the binder component of PBX was represented by HPTB. Instead of studying thermochemistry due to a hot spot, Yan et al.41 used RMD to simulate direct thermal decomposition of PBX systems containing CL-20 and different types of polymer binders. However, these RMD studies were performed for very small simulation durations (< 10 ps) and were restricted only up to the formation of a hot spot near fuel/binder interface. Hence, there is a lack of understanding about the interactions between a hot spot and a polymeric binder. Lack of in-situ experimental probes have created a gap in the understanding of the role of polymeric binder in this complex thermochemical process. The energy up-pumping at a fuel/binder interface in a highly heterogeneous thermochemical environment and the effects of such energy up-pumping on the kinetics of hot spot growth are important links to connect the exact steps involved in the chemical degradation of the binder. Since a hot spot can lead to a self-sustaining 6 ACS Paragon Plus Environment

Page 6 of 41

Page 7 of 41 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


reaction front, a deeper understanding of the thermochemistry of a fuel/binder interface will greatly benefit design of PBXs. The objective of this work is to investigate how a thermal hot spot will interact with polymeric binder in RDX/HTPB energetic composite using RMD. The computational models used for studying RDX/HTPB composite will be discussed in the second section. Capturing thermomechanical and thermochemical responses of both RDX and HTPB with the same set of force field parameters is a challenging task. In our previous studies, we have investigated energy up-pumping42 and thermochemistry26 near a hot spot interface in condensed phase RDX using ReaxFF force field parameters by Wood et al.11. As a result, we have used those force field parameters in this work for the sake of consistency. However, to the best of our knowledge, those force field parameters have never been tested for HTPB binder. Hence, we will first present RMD simulations which were performed to benchmark the force field for mechanical and thermophysical properties of HTPB. Although the aim of this work is not to re-parameterize the force field for the binder, we believe that such benchmarking is essential for an enhanced understanding of the observed trends in RMD simulations of composite systems. Next, we will discuss the non-equilibrium dynamics between a thermal hot spot in RDX and HTPB binder. Our calculations on non-equilibrium dynamics indicate that HTPB binder is subjected to high temperatures (> 1000K) and high pressures (> 10 GPa) if a hot spot undergoes ignition to deflagration transition. Hence, in the third section, we will present RMD-based chemistry of the thermal degradation of HTPB binder at high temperatures and high pressures. The details about the RMD simulation protocols for each simulation type are provided in the corresponding section.



2. Computational Models 2.1 PBX Models Developing a realistic PBX model in molecular scale is challenging. Hence, a simplified model of the interfaces is attempted in this work for systematically assessing the role of interfaces. The binder in our models is represented by trans 1,4 HTPB molecule with 32 backbone carbon atoms (C32H50O2). We assume that the lower molecular weight model of the binder used in the study is a reasonable representative for the local thermchemical initiation and energy transfer dynamics. However, some differences in mechanical properties are expected. It is well known that building a perfect polymer structure is a challenging task due to very long relaxation times and chain entanglementes and hence is beyond the scope of this work. The glass transition and melting temperatures of HTPB are 190K and 284K respectively43. Hence, at temperatures above room temperature (> 298K), HPTB is most likely to occur in amorphous melt state. To represent such HTPB melt, 198 HTPB polymeric molecules (each polymer with 32 carbon atoms) were initially randomly placed in an orthogonal box using PACKMOL44. The HTPB in our model is not cross-linked and can cause a lower estimate for bulk modulus and yield strength. This initial structure was equilibrated at 300K. After thermal equilibration, the box was compressed for 500 ps to obtain an amorphous melt structure with an equilibrium density of 0.94 gm/cm3 at 1 atm. In the PBX system, the explosive is represented by the γ-phase of RDX with starting density of 2.26 gm/cm3. The choice of the -phase is based on recent 8 ACS Paragon Plus Environment

Page 8 of 41

Page 9 of 41 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


observations which indicated that the chemical decomposition is mostly initiated in the γ-phase over the α-phase under moderate shock compression

5, 45

. The unit cell of -phase was energy

minimized using ReaxFF to relax atomic positions. We then constructed a RDX supercell by replicating the unit cell 3 times in the x and y directions and 40 times in the z-direction. The resulting dimensions of the RDX supercell with 2862 RDX molecules were 3.768 nm×2.844 nm ×43.72 nm. The longitudinal length of 43.72 nm was chosen to ensure that a steady deflagration front can be developed during the RMD simulations as observed in our previous work46.

Figure 1: Snapshot of PBX system consisting of -phase RDX and HTPB with chain length nc = 32. The PBX system is 58.2 nm long in the z-direction and has the cross-sectional area of 3.768 nm ×2.844 nm. Color scheme: cyan- carbon, blue – nitrogen, redoxygen, white – hydrogen Then we constructed the composite system between the RDX and HTPB by creating a planar interface along (001) orientation of -RDX. The resulting composite system is 58.2 nm long in the z-direction and consists of approximately 79% RDX and 21% HTPB by atoms (by weight, 87.32% RDX and 12.68% HTPB). The composite system was first energy minimized for 9 ACS Paragon Plus Environment


local relaxation. After energy minimization, the composite was slowly heated to 500K to form the stable interface. The equilibrium pressure of the composite system is approximately 7.25 GPa (see Figure S1 and S2 of SI) and is higher than the critical pressure of 4 GPa above which RDX undergoes α to phase transition5. This observed pressure is probably over-estimated and the addition of van der Waals correction to the force field may be required for the improved description of the P-T state47. Additional discussion on the choice of the force field is provided in the next section. 2.2 Reactive Molecular Dynamics Simulations The RMD simulations were performed using the LAMMPS MD48 suite with the ReaxFF force field. ReaxFF is a bond order dependent reactive force field which allows bond formation and bond breaking during a MD simulation. In ReaxFF, energy of a system is described using both bonded and non-bonded interatomic interactions. The bonded interactions include bond energy, valence angle energy, torsion angle energy, under and over coordination energies and lone pair energy. Every component of bonded interactions is made bond order dependent to facilitate smooth transitions between bond formation and bond breaking. The non-bonded interactions include Coulomb and van der Waals interactions. More details about each energy term can be found in Chenoweth et al.49. The ReaxFF parameterization for nitramines was originally developed by Strachan et al.

24

. The revised version of that parameterization was

released in Wood et al.50 The revised force field parameters were obtained by combining the original ReaxFF nitramine (C/H/N/O) training set with C/H/O combustion training set49. The nitramine training set was predominantly parameterized for the unimolecular chemical reactions. In the condensed phase, bimolecular reactions can play an important role in the overall chemical 10 ACS Paragon Plus Environment

Page 10 of 41

Page 11 of 41 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


evolution. In our earlier work, we have shown that the revised force field describes such biomolecular reactions within a reasonable accuracy expected from the underlying DFT methods26. The CHO combustion parameters have been widely used to study combustion chemistry of different hydrocarbons like n-dodecane51, n-heptene52 and coal53. Besides this original ReaxFF formulation, Liu et al.47 developed a new formulation, called ReaxFF-lg, by adding a third component to non-bonded interactions to incorporate London dispersion effects. They reported that ReaxFF-lg is better than the original version in describing crystal structure and P-T state of energetic materials. However, authors did not report the performance of ReaxFF-lg for the combustion chemistry of hydrocarbons and polymeric binders like HTPB. Hence, in this work, RMD simulations were performed using the force field parameters by Wood et al.50 (non-lg formulation) because it provides an enhanced description of combustion chemistry of both nitramines and hydrocarbons. 3. Results 3.1 Mechanical and thermophysical properties of HTPB The mechanical and thermophysical response of HTPB can be described by various factors such as Hugoniot Locus, equation of state, specific heat and thermal conductivity. The comprehensive estimation of such mechanical and thermophysical properties is beyond the scope of this work. However, for benchmarking the force field parameters, we performed RMD simulations on a single HPTB configuration to evaluate bulk modulus ( 0), radial distribution function (RDF), constant volume specific heat (Cv) and thermal expansion coefficient (). The starting structure for these RMD simulations corresponds to the configuration shown in Figure 1 (198 HTPB molecules, density = 0.94 gm/cm3). For the bulk modulus, NPT simulations were 11 ACS Paragon Plus Environment


performed at 300K at different pressures that range from 1 atm to 5 GPa. For each target pressure, the system was compressed to the desired pressure in 250 ps and then isothermalisobaric simulation at the desired pressure was performed for 750ps. For the constant volume specific heat ( ), the MD simulation was performed using NVT ensemble at 300K for 1 ns. For calculating the thermal expansion coefficient, the system was heated from 300K to 600K in 1 ns using NPT ensemble to measure the change in the system volume as a function of temperature. These RMD simulations were performed with a time step of 0.25 fs. Table 1 lists the calculated mechanical and thermophysical properties of HTPB from RMD simulations. The equilibrium density was obtained from the NPT-MD simulation at 1 atm. The equilibrium density predicted by ReaxFF force field parameters is in good agreement with the previously reported values in the literature. The pressure-compression curve was used for calculating the bulk modulus. Figure 2a compares the ReaxFF pressure-volume curve with the

(a) Pressure-compression curve

(b) Radial Distribution Function


Page 12 of 41

Page 13 of 41 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


Figure 2: (a) Pressure-compression behavior of HTPB. Hooper et al. data is based on nonreactive MD simulations3. Millett et al. data is based on single-stage gas gun plate impact experiments4. (b) RDF as a function of distance of HTPB. The distribution function includes all atom types. previous experimental4 and molecular simulations3. Following Hooper et al.3 and Fröhlich et al.32, we fitted our data to empirical Tait equation of state (EOS) 54 which is given by following equation:

() = 1 − × 1 +

(1)

where is equilibrium volume and and are material-dependent parameters which were evaluated using least square fitting technique. The zero-pressure bulk modulus was calculated using; κ0 = / and is listed in Table 1. The bulk modulus obtained from the RMD simulations is lower than the modulus reported by the non-reactive MD3, 32 and experimental4 studies in the literature. The under-prediction of bulk modulus can be due to the lower molecular weight and absence of any cross-links in HTPB model. The discrepancy also indicates that further improvement in the force field is needed since it was primarily trained only for CHO combustion chemistry. Since force field under predicts bulk modulus, it will allow higher compression of the binder at higher pressures and will under predict the speed of sound in HTPB. Table 1: Summary of calculated mechanical and thermophysical properties of HTPB using RMD Property

ReaxFF


Literature


Page 14 of 41

Density, 0 (gm/cm3)

0.94

0.9153, 0.83032, 0.95055

Bulk Modulus, 0 (GPa)

0.9691

1.4143, 1.36832

Constant volume specific heat, (J/Kg.K)

4210

286055

Thermal expansion coefficient, (10-4/K)

10.1

7.556, 8.057

The radial distribution function is calculated using equation 2.

() =

1 〈 ( − !" )〉 4 !#"

(2)

In equation 2, is the total number of atoms, is the ambient number density and !" is the interatomic distance. The RDF was calculated for all atom types. Figure 2b shows the RDF of the equilibrium HTPB configuration at 1 atm and 300K. The non-zero value of () near 0.96 Å is due the O-H bond of the hydroxyl termination. The next three peaks at 1.35 Å, 1.46 Å and 1.52 Å corresponds to C=C, C-O and C-C bonds respectively. The RDF obtained from the ReaxFF force field agrees very well with the previously reported RDF from the non-reactive MD simulations32. The constant volume specific heat is calculated using energy fluctuations58 as shown in equation 3.

=

〈 % 〉&'( )* +

(3)

,ℎ.., 〈 % 〉 = 〈% 〉 − 〈%〉


Page 15 of 41 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 thermal expansion coefficient is calculated using following equation:

=

1 0 0+ &1(

(4)

Table 1 compares the ReaxFF predictions of and with the existing data in the literature. Since ReaxFF predicts lower bulk modulus (softer binder), it over-predicts the thermal expansion coefficient of HTPB. The current force field parameters also over-predict the constant volume specific heat( ). As a result, during energy up pumping, HTPB will absorb more energy for thermal equilibration. This will be more evident in hot spot/HTPB interaction calculations which are discussed in the next section. 3.2 Non-equilibrium dynamics at RDX/HTPB interface due to a moving deflagration front In order to develop an understanding of the non-equilibrium dynamics at RDX/HTPB interface, it is important to recreate initiation conditions that might exist in a PBX and analyze the role of the binder in the time evolution of a deflagration wave from the energetic solid through the binder phase. A thermal hot spot in the RMD simulations was created in the RDX phase. In order to model a hot spot, 144 RDX molecules (5% of the total energetic material) were subjected to a temperature pulse of 2000K using an NVT ensemble. During the application of the thermal pulse, rest of the system was allowed to evolve using an NVE ensemble. After the thermal pulse duration is over, the entire system was simulated using NVE ensemble to mimic adiabatic conditions. As discussed in the introduction section, hot spots can be formed near or away from the interface region. To take into account both of these possibilities, we performed RMD simulations on two systems (Figure 3). In the first system, the thermal hot spot is located nearly 20 nm away from the interface (central hot spot) and in the second system, the thermal hot 15 ACS Paragon Plus Environment


spot was positioned at the RDX/HTPB interface (interface hot spot). The RMD simulations were performed with periodic boundary conditions in all three directions and with a conservative time step of 0.1 fs.

Figure 3: Snapshots of starting structures of Central hot spot, Interface hot spot and Only RDX configurations. For each configuration, the hot spot is shown by larger red-colored spheres. The size of hot spot is same for all configurations. Color scheme: cyan- carbon, blue – nitrogen, redoxygen, white – hydrogen 3.2.1 Effect of HTPB binder on ignition to deflagration transition The key role of a binder in safe handling of a PBX is in preventing accidental initiation from mechanical and thermal initiation events. It is important to understand if a hot spot created in RDX can be severely weakened before it travels to another RDX grain. The temperature evolution during ignition is followed as a means to identify the onset of ignition and thermal runaway. Figure 4 shows respective temperature profiles of hot spot, non-ignited RDX and 16 ACS Paragon Plus Environment

Page 16 of 41

Page 17 of 41 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


HTPB regions for both central and interface configurations. For comparison, we have also included temperature profiles of the “only RDX” system which does not contain any binder molecules. The

(a) Hot Spot

(b) Non-ignited RDX

(c) HTPB

Figure 4: Temperature profiles of hot spot, non-ignited RDX and HTPB molecules as a function of temperature during ignition to deflagration transition. The observed temperature profiles indicate that ignition to deflagration transition is delayed for the interface configuration due to the energy loss to the binder. “only RDX” configuration is explained in detail in our earlier work46 and is relevant for void collapse inside a single grain of EM which has been studied by other research groups11-12. For the central configuration, a thermal pulse of 15 ps was enough for forming a steady deflagration wave. For the interface configuration, a longer thermal pulse of 30 ps was provided to ensure that a steady deflagration front will be observed within the accessible time scale of RMD. This longer thermal pulse is required due to the heat loss from the hot spot to HTPB at the interface. More details about the energy transfer at the interface are provided later in this section. The temperature profiles indicate that the ignition to deflagration transition occurs for both central 17 ACS Paragon Plus Environment


and interface configurations via ignition delay phase. The ignition delay phase is characterized by energy transport from the hot spot to the surrounding regions via heat diffusion. This is evident from the linear temperature profiles of the non-ignited RDX and HTPB molecules. As the temperature of the RDX molecules surrounding the hot spot increases, it triggers local exothermic chemical reactions. Due to such series of local exothermic explosions of increasing intensity, the energy transport by mass diffusion eventually exceeds the heat diffusion thereby forming a steady deflagration front. The temperature profile of the central hot spot configuration closely resembles the temperature profile of “only RDX” configuration (Figures 4a and 4b). This trend indicates that the presence of binder has nominal effect on the ignition to deflagration transition if a hot spot is formed away from the interface. The effect of binder on system temperature is more evident in the interface configuration which shows significant delay in ignition to deflagration transition. For example, during the ignition delay phase, the hot spot region equilibrates to an intermediate temperature of 1450K after 200 ps (Figure 4a). The corresponding equilibrium temperatures of the non-ignited RDX and HTPB regions are 750K (Figure 4b) and 1025K (Figure 4c) respectively. Figure 4b indicates that the ignition delay of interface configuration is approximately 2.25 times the ignition delay of central hot spot configuration. This increase in the ignition delay is mainly due to the heat loss from the hot spot to the HTPB binder and is discussed in more detail in the next section. 3.2.2 Intermolecular energy transfer at the hot spot/HTPB interface To obtain an estimate on how much energy is lost to the binder at the interface, we calculated the change in specific total energy (potential + kinetic) of the hot spot and the binder for the time period between 15ps to 250ps. The respective energies of the hot spot and the 18 ACS Paragon Plus Environment

Page 18 of 41

Page 19 of 41 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


binder at the end of the thermal pulse (15ps) are assumed as reference energies. Figure 5a shows the change in specific energies up to 250 ps while the Figure 5b shows how the up pumped energy in HTPB is further partitioned into the potential and kinetic energy components. It should be noted that the molar mass of the binder is nearly three times the molar mass of the hot spot. Figure 5a shows that the energy of the hot spot reaches its intermediate equilibrium state in 200 ps and considerably more energy is up-pumped into HTBP binder than the surrounding RDX molecules.

(a) Specific Energy profiles of the hot

(b) Partitioning of adsorbed specific energy in

spot and HTPB

HTPB

Figure 5: Energy up-pumping at RDX/HTPB interface for interface hot spot configuration. (a) Specific energy profiles of the hot spot and binder. The specific energy includes both the kinetic and potential energy components. It should be noted that the molar mass of binder is nearly three times the molar mass of hot spot. (b) Partitioning of the up-pumped energy into kinetic and potential energy components. 19 ACS Paragon Plus Environment


Page 20 of 41

Such an uneven distribution of energy between the adjacent RDX and HTPB regions can be due to the different factors such as interface thermal conductivity, interface geometry and vibrational excitations during energy up-pumping. Near room temperature, thermal conductivities of HTPB and RDX are 0.22 (W/m-K)59 and 0.20 (W/m-K)60 respectively. If thermal conductivities of the two materials remain comparable at higher temperatures, then excitation of vibrational modes under thermal impulse will play significant role in energy uppumping at interfaces. To a first order approximation, vibrational excitations in a material are described by heat capacity. The experimental specific heat capacities of HTPB and solid RDX near 300K are 2860 (J/Kg-K)55 and 1078 (J/Kg-K)61 respectively. Hence, in the initial stages of energy up-pumping, more energy will be up-pumped into HTPB than in RDX due to its higher capacity. Since ReaxFF parameters overestimate heat capacity of HTPB, we anticipate that more energy is transferred to HTPB in RMD simulations than reality. The heat capacity-based reasoning may not be fully valid at different temperatures because heat capacity of a material changes with temperature. To circumvent this problem, we further probed the role of vibrational excitations during energy up-pumping using intermolecular energy transfer rate at the interface between the hot spot and HTPB. The intermolecular energy transfer rate is calculated using following expression:

κ =

∆% ∆4 × ∆+

(5)

where ∇E is change in potential energy, ∇T is the corresponding change in temperature and ∇t is time. For HTPB, we obtained intermolecular energy transfer rate of 1.50×1015 J/(K-mol-s). In our earlier work, we have reported that the energy transfer rate from hot spot to RDX crystal is 20 ACS Paragon Plus Environment

Page 21 of 41 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


approximately 3.40×1014 J/(K-mol-s)42. The observed rates show that, in RMD simulations, the rate of energy transfer from a hot spot to HTPB is at least 4 times higher than the corresponding energy transfer rate between a hot spot and RDX. These calculated energy transfer rates further support earlier heat capacity based conclusion that more energy will be up-pumped into HTPB than in RDX. It should be noted that RMD simulations do not include excitation of electronic modes and quantum effects on vibrations of light atoms like hydrogen. Inclusion of these quantum effects might alter the ratio of energy transfer rates. In our earlier work, we have shown that at higher temperatures, inclusion of quantum effects does not have significant impact on excitation of different types of vibrational states42. Hence, we anticipate that influence of such quantum effects will be nominal on energy transfer ratio at the interface at simulated temperature range. As a result, non-ignited portion of RDX of the interface configuration equilibrates to a lower intermediate temperature than the central hot spot configuration which delays ignition to deflagration transition. Since ReaxFF over predicts energy up-pumped into HTPB, we anticipate that it will also over predict ignition delay. Hence, actual ignition delay observed in experimental systems might be shorter than the one observed in the interface RMD simulation. The results from the previous two sections show that an interface hot spot requires higher energy deposition and longer time to undergo ignition to deflagration transition. Therefore, our results suggest that a void collapse inside an RDX grain can be more consequential in triggering ignition than a hot spot at the binder interface with the same temperature. Hence chemical stability of PBXs can be improved by selecting optimal combination of fuel and binder. The dynamic measure is related to the widely different thermal up-pumping observed at the model interface for HTPB.



3.2.3 Effect of HTPB binder on a moving deflagration front Once the ignition to deflagration transition occurs, a steady deflagration front propagates through rest of the RDX phase. The most challenging part of tracking a moving deflagration front is to predict mass and energy flux due to local chemical transformations. To calculate flux values across a deflagration front, we mapped RMD trajectory onto Eulerian Control Volumes (CVs). Such mapping provides a reasonably good estimate of the chemical source energy of reactive


Page 22 of 41

Page 23 of 41 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


Figure 6: Temporal evolution of specific potential and kinetic energies (source) in microscopic control volumes for the central hot spot configuration. The system was divided in 30 CVs along z-direction. The legend of each frame shows the location of microscopic control volume along zdirection. combustion fronts. More details about the RMD to CV methodology can be found in our earlier work46. For both configurations, the entire system was divided into 30 CVs along the z-direction. Each CV has a width of 1.94 nm. Time averaging over 5 ps intervals was used to reduce the magnitude of fluctuations that are inherent in the MD simulation data. The integration time interval for the CV analysis was 5ps. Figure 6 shows the CV analysis for the central hot spot configuration and Figure 7 shows the corresponding analysis for the interface hot spot configuration. The approximate location of RDX/HTPB interface is indicated in each figure. However, the exact location of the interface might shift during the propagation of deflagration front due to the compression of HTPB at high pressures. Figure 6 shows that, in central configuration, the deflagration wave is formed after 320 ps and it travels further downstream for approximately 10 nm before hitting the HTPB interface. The energy profiles of the CV located at 37.8 nm show that the deflagration wave is characterized by the sharp drop in the potential energy with the simultaneous sharp rise in the local temperature. The source KE profile near deflagration front is broader than source PE profile for all CVs. This indicates the rate of release of chemical energy is higher than the rate of increase in temperature. However, the total chemical energy released due to deflagration (area under the curve of source PE) matches closely with total increase in source KE indicating near complete interconversion of potential energy into temperature. Due to the periodic boundary conditions, the deflagration waves hit the binder from both sides. By the time deflagration waves 23 ACS Paragon Plus Environment


reaches the binder, most of the RDX is converted into final gas products. This hot mixture of gas products then heats the HTPB molecules at an approximate rate of 1.28×1013 K/s (Figure 4). However, unlike RDX, the deflagration wave does not incite localized thermal explosions in HPTB. This is evident from the rise in both potential and kinetic energy components of the HTPB CVs during the propagation of the diminishing deflagration wave (Figure 6). As a result, deflagration wave slowly loses its strength as it moves into the HTPB region. The deflagration wave travels at an approximate speed of 184 m/s in RDX region before reaching the binder region. This observed speed for the central configuration is lower than “only RDX” system wherein the similar wave moves at an approximate speed of 214 m/s. This can be due to the pressure difference between the two systems. The system pressure increases from 9 GPa to 25 GPa during the propagation of the deflagration wave in central configuration. For “only RDX” system, the pressure increases from 16 GPa to 35 GPa in the equivalent time period. In multigrain PBX system with multiple RDX and HTPB regions, such a deflagration wave has to pass through a HTPB region before igniting the next grain. To investigate this scenario further, we performed an additional simulation on a multigrain system that consisted of two grains of RDX separated by HTPB binder (see figure S3 in SI). In this system, a hot spot was assumed to be away from the interface. We found that, energy up-pumping and thermochemistry of deflagration front in multigrain system is similar to the single grain central hot spot configuration. In contrast to single grain configuration, a multigrain configuration can provide an estimate about how much time a wave will spend in HTPB before igniting the next RDX grain. The observed temperature profiles show that minimum of 300 ps were required before the RDX molecules of the next grain attain ignition temperature of 2000K (Figure S4). Hence, the speed of the wave is at least 3x slower than the corresponding wave speed in RDX. 24 ACS Paragon Plus Environment

Page 24 of 41

Page 25 of 41 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


As a result, HTPB impedes the rate of growth of hot spot in the multigrain system and can lead to extinction of hot spots providing a safety mechanisms against accidental initiation. Figure 7 shows the results of the CV analysis for the interface hot spot configuration. For this configuration, the deflagration front appears after 600 ps at a location that is approximately 10 nm from the interface. The deflagration wave then further travels inwards into the RDX region (towards left as per Figure 2) for 35 nm before reaching the interface. The main similarity between central and interface configurations is in the thermochemistry of the propagating deflagration front



Figure 7: Temporal evolution of specific potential and kinetic energies (source) in microscopic control volumes for the interface hot spot configuration. The system was divided in 30 CVs along z-direction. The legend of each frame shows the location of microscopic control volume along z-direction. in RDX region. The PE and KE profiles of the propagating deflagration front in the interface configuration closely match with the corresponding profiles of the central configuration. This trend shows that the location of a hot spot has minimal effect on the energy and mass flux across a steady deflagration wave in RDX. Figures 6 and 7 also show some key differences between the two configurations. For instance, unlike the central hot spot configuration, only one deflagration wave is formed in the system and as a result, the deflagration wave hits the binder only from one side (from right side as per Figure 2). Figure 7 shows that this deflagration wave barely reaches the binder region during the simulated period of 800 ps. As a result, unlike central configuration, we did not observe any direct wave-binder interactions in the interface configuration for the simulated duration. Instead, HTPB molecules are heated by the energy transport from the hot spot region. This mode of energy up pumping is more gradual compared to energy up pumping due to the deflagration wave. E.g. For interface configuration, HTPB molecules are heated at rates that vary between 5.2×1011 - 46×1012 K/s (see Figure 4). This heating rate is one order of magnitude lower than the central configuration. As a result, for the interface configuration, we did not observe any sudden changes in the source PE and KE energies of HTPB CVs (51.3 nm and 53.2 nm panels in Figure 7). For the interface configuration, the observed propagation speed of the deflagration wave is approximately 160 m/s. Such a lower wave speed is expected because the rise in pressure in the interface configuration (from 7 GPa to 20 GPa) is lower than the central and “only RDX” configurations. 26 ACS Paragon Plus Environment

Page 26 of 41

Page 27 of 41 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


3.3 Chemistry of thermal decomposition of HTPB There is a lack of understanding about how the chemistry progresses and how a deflagration wave can decompose the binder before reaching the next grain of RDX. Most of the known post-detonation dynamics is related to combustion of a binder. However, the design of insensitive PBX depend on knowing the exact chemistry in the binder phase undergoing thermal decomposition in the confined intergranular mechanism.

Figure 8a shows the rate of

consumption of RDX and HTPB molecules as a function of time for two hot spot configurations. The BO cut-off 0.5 was used for the molecular population analysis. As discussed earlier, the combustion chemistry is mostly a multi-step process and the disappearance of RDX and HTPB molecules from the population profile does not imply sudden conversion into the final product state. Figure 8 shows that, for both configurations, the rate of consumption of RDX due to a moving deflagration front is significantly higher than the pre-deflagration phase. As expected, for the central configuration, HTPB molecules begin to undergo chemical transformations only after the deflagration front reaches the RDX/HTPB interface (after 400 ps). Figure 8a shows that only 36% of HTPB has reacted in 500 ps. For the interface configuration, the HTPB molecules undergo chemical transformations for the entire simulation duration due to the continuous heating from the hot spot region. However, the rate of chemical transformations of HTPB increases sharply in post-deflagration regime. For example, Figure 8b shows that only 24% HTPB has reacted before deflagration (in 585 ps). The remaining 70% reacts in the post deflagration duration of 250 ps.



(a) Central Hot Spot

(b) Interface Hot Spot

Figure 8: Rate of consumption of RDX and HTPB during the RMD simulations of the central and interface configurations. The molecular populations were obtained using BO cut-off of 0.5. The disappearance of RDX and HTPB molecules from the population profiles do not imply sudden conversion into the final products. We can obtain more chemical insights about the chemical transformations of HTPB under extreme conditions from these two RMD trajectories. However, near the interface region, isolating HTPB molecules from the hot gaseous mixture from the deflagration is non-trivial. Also, in these two RMD simulation, temperature and pressure of HPTB vary continuously during NVE evolution. As a result, the direct analysis of the RMD trajectories of these two configurations cannot reveal the temperature and pressure dependence of the thermal degradation chemistry. Hence, we performed separate pyrolysis simulations at different temperatures and pressures. For each temperature and pressure, the starting structure corresponds to the equilibrium configuration obtained from the corresponding equation of state calculation. The pyrolysis simulations were performed using NPT ensemble. For studying temperature dependence, the pressure was set to 1 atm and temperature was varied from 1400K to 2000K in 28 ACS Paragon Plus Environment

Page 28 of 41

Page 29 of 41 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


steps of 200K. For pressure dependence, the temperature was set to 2000K and pressure was varied from 1 atm to 10 GPa. For each simulation, starting with the desired pressure, the system was heated from 300K to the desired temperature in 125 ps. After reaching the desired temperature, isothermal-isobaric simulation was continued for 200 ps or more depending on the chemical composition of the system. We analyzed RMD trajectories of the temperature dependent simulations to identify the dominant pyrolysis mechanism. Figure 9 shows the schematic of the dominant reaction mechanism observed in the RMD simulation at 1 atm and 2000K. The BO cut-off of 0.5 was used for identifying the molecular species. In Figure 9, the “- C4H6” on the top of an arrow indicates that the product was formed by eliminating butadiene from the reactant molecule. The number below each arrow indicates how frequently that reaction was observed during the simulation. It can be seen that unimolecular depolymerization reactions dominate HTPB pyrolysis chemistry at 2000K and 1 atm. At 2000K, all HTPB molecules were activated by cleavage of the aliphatic C-C bond.



Figure 9: Schematic of the dominant high temperature pyrolysis mechanism of HTPB for pressures up to 1 GPa. The “- C4H6” indicates unimolecular butadiene elimination reaction. Each number shows the frequency of that reaction. The cleavage of the central C-C bond is less preferred over other aliphatic bonds. The larger allylic radicals then further undergo successive β-scission reactions until they form the most abundant intermediate radical C16H25O. This radical then undergoes further β-scission reactions until all molecules are converted into C4H7O and 1,3 butadiene. We found that nearly 50% of C4H7O further eliminates OH to form 1,3-butadiene. The observed trend towards butadiene formation via chain fragmentation is mainly due to high temperatures which prefer entropically favored reactions52. Previous experimental studies have reported that, above 750K, extent of depolymerization and C-C cleavage increases with increase in temperature34, 36. For other three temperatures (1400K–1800K), similar reaction mechanism dominated the RMD evolution. We also calculated apparent activation energy using temporal evolution of HTPB population (see figure S5 of SI). The current force field parameters predict the apparent activation energy to be 42.7 kcal/mol. This calculated value is lower than the previously reported experimental activation energy of 60.1 kcal/mol34. This can be due to the absence of cross-linked versions of HTPB used in PBX and points to the need for improving the force field which is trained only for CHO combustion environment. To investigate the pressure dependence on thermal degradation mechanism, we calculated population of HTPB and butadiene molecules from RMD trajectories. Figure 10 shows the population of HTPB and butadiene as a function of time at various pressures. The HTPB population profiles do not show much dependence on pressure. For the simulated pressure 30 ACS Paragon Plus Environment

Page 30 of 41

Page 31 of 41 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


range, most of the HTPB molecules get converted into intermediate radicals in 200 ps. However, the butadiene formation gets progressively suppressed at higher pressures. Figure 10b shows that the effect of pressure on butadiene suppression is nominal up to 1 GPa. However, population of butadiene decreases by 50% at 5 GPa and by 85% at 10 GPa. This trend clearly indicates a shift in the pyrolysis mechanism at pressures above 1 GPa.

(a) HTPB Population

(b) Butadiene (C4H6) Population

Figure 10: Temporal evolution of molecular populations of HTPB and butadiene during pyrolysis simulations at different pressures. The molecular populations were obtained using BO cut-off of 0.5. The population analysis was performed only for constant temperature duration. At higher pressures (5 GPa and above), we found that nearly 51% of HTPB is activated by the four initiation reactions shown in Figure 9 (depolymerization). However, these activated molecules predominantly undergo bimolecular reactions instead of butadiene elimination. The observed bimolecular reactions include hydrogen transfer between the radicals and with adjacent HTPB molecules and chain cross-linking which leads to branched polymers. To quantify the extent of chain branching, we calculated number of carbon atoms with three carbon neighbors 31 ACS Paragon Plus Environment


(denoted as C3 here after). We performed this analysis on the pyrolysis system at 5 GPa and the interface hot spot system. Figure 11 shows the temporal history of number of C3 atoms for both systems. For the pyrolysis system, we extended the simulation time up to 800 ps to ensure that such branched fragments are stable for longer durations.

(a) Pyrolysis, 5 GPa

(b) Interface Hot Spot

Figure 11: The number of carbon atoms with three carbon neighbors as function of time. The BO cut-off of 0.5 was used for identifying chemical fragments. (a) Pyrolysis simulation at 5 GPa and 2000K, (b) HTPB molecules of the interface hot spot configuration. The comparison between Figures 10a and 11a shows that the branching begins after more than 80% of the HPTB molecules are chemically activated. The temporal population profile of C3 atoms shows an exponential-like behavior with sharp increase in the population up to 400ps. This trend further supports the earlier observation that, at high pressures, the intermediate allylic radicals prefer bimolecular chain branching reactions over β-scission reactions. The population profile of the interface hot spot configuration shows a sudden rise in the population of C3 atoms 32 ACS Paragon Plus Environment

Page 32 of 41

Page 33 of 41 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


when the deflagration wave is developed in the system. This sudden rise occurs in the same temporal interval in which the temperature of the binder increases beyond 2000K and the pressure rises above 15 GPa. It clearly shows that the chemical evolution of HTPB in PBX system closely resembles the pyrolysis chemistry at high pressures. In addition to branched fragments, other major pyrolysis products at higher pressures include 4-vinyl cyclohexene (C8H12), cyclohexadiene (C6H8), cyclopentadiene (C7H10) and ethylene (C2H4). Most of the experimental pyrolysis studies are done at much lower temperatures (T < 600K) and pressures (less than 15 atm). E.g. Arisawa and Brill36 found that trans-butadiene oligomers (t-BDO) are more dominant than butadiene and 4-vinyl cyclohexene for the temperature range of 723-882K. Ganesh et al.37 proposed that, near 773K, HTPB degradation occurs primarily via cyclic compounds and linear oligomers. Replicating experimental temperatures is beyond the accessible time-scale of RMD. However, it is worthwhile to mention that most of the linear and cyclic compounds that are observed in high pressure RMD simulations have been reported in those experimental studies. Besides these organic compounds, we also observed formation of water at high pressures. The water molecules are mostly formed by hydrogen transfer reactions. These hydrogen transfer bimolecular reactions are possible at high pressures (above 15 GPa) due to closer proximity of the molecules. These temperature and pressure dependent pyrolysis simulations clearly indicate that, above 1400K, thermal degradation mechanism of HTPB is dependent upon pressure. At low pressures, thermal degradation is dominated by depolymization reactions. However, at high pressures (5 GPa and above), cyclization and chain cross-linking paths are as important as depolymerization. This high pressure mechanism dominates the overall thermal degradation chemistry of HTPB in PBXs which are exposed to deflagration. 33 ACS Paragon Plus Environment


4. Conclusions Understanding thermochemistry of binders and additives in extreme environments is important for designing insensitive and mechanically formable PBX. The thermochemistry is especially challenging when the PBX undergoes a transition from thermal ignition to an ultrafast kinetically controlled deflagration regime. However, the safety and performance balance is closely linked to the non-equilibrium dynamics at the interface of brittle energetic crystallites and the rubbery binders.

As a way to provide systematic understanding of the role of the

binders, RMD simulations were performed to investigate the thermochemistry in RDX/HTPB composite exposed to a moving deflagration front – a condition that determined if the coalescence of many such fronts will lead to a controlled condensed phase flame in a propellant versus an uncontrolled detonation. In order to understand such events using RMD simulations, same CHNO force field needs to describe the RDX and HTPB with a reasonable degree of accuracy. We have previously used this force field to understand deflagration chemistry of RDX. The same force field was used to perform RMD simulations of the binder phase and benchmark thermo-mechanical properties of HTPB. We found that ReaxFF reproduced the equilibrium density and the equilibrium structure of HTPB with good agreement with the literature. The same force field under predicts the bulk modulus and over predicts the constant volume specific heat capacity. For modeling non-equilibrium dynamics due to a moving deflagration front, we performed RMD simulations on a model interface representing a simplified PBX system consisting of 79% RDX and 21% HTPB. Undoudtedly, this is highly simplified compared to faceting of crystallites and presence of different defects in an in-service PBX. The goal of this work is to understand and quantify the role of a commonly used binder such as HTPB exposed to hot spots. The simulations were performed on two configurations, viz, “central hot spot” with hot 34 ACS Paragon Plus Environment

Page 34 of 41

Page 35 of 41 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


spot located away from the RDX/HTPB interface and the “interface hot spot” with the hot spot located at the RDX/HTPB interface. The RMD simulations show that ignition to deflagration transition occurs in both configurations. However, if a hot spot is formed near the interface, then longer thermal pulse with higher energy is required to observe ignition to deflagration transition. In addition, the transition is significantly delayed for the interface hot spot configuration. This is mainly because HTPB binder takes in more energy from the hot spot than surrounding RDX. This is key to improving the safety of PBX using more complex formulations of binders. We found that intermolecular energy transfer rate from a hot spot to HTPB is at least four times more than the corresponding rate from a hot spot to RDX. This is mainly because heat capacity of HTPB is higher than RDX but have deeper molecular level linked to thermal up-pumping mechanisms available in HTPB vs. RDX. The energy loss from the hot spot to the binder more than doubles the ignition delay in the interface configuration. The simulations show that thermochemistry of an accidental initiation caused inside a PBX crystal due to void collapse can be significantly different than an initiation caused at a non-homogenous fuel/binder interface. We have also developed quantitative methods for tracking the deflagration transition and the deflagration front. To assess the energy flux associated with a moving deflagration front, we mapped the RMD trajectories onto Eulerian CVs. We found that a deflagration wave is characterized by rapid conversion of local chemical potential energy into kinetic energy (high temperature gas products). The observed energy profiles show that the location of a thermal hot spot has nominal effect on the energy and mass flux across a moving deflagration front. In central configuration, HTPB binder is heated due to direct interaction with the moving deflagration front and hence temperature (approximate rate of 1.28×1013 K/s) and pressure of the binder rises sharply due non-equilibrium energy transfer. The CV analysis of the binder region 35 ACS Paragon Plus Environment


also shows considerable weakening of the deflagration wave during its passage through the HTPB region. For the interface configuration, the thermal energy from the hot spot starts heating the HTPB binder before deflagration front is developed and this diffusion based mode of energy up-pumping is an order of magnitude lower than the heating due to deflagration front. For both the configurations, the final temperature of HTPB is above 2000K and the final pressure is above 15 GPa. To identify the chemistry of thermal degradation of HTPB under such extreme conditions, we performed separate pyrolysis simulations at different temperatures and different pressures. We found that thermal degradation mechanism of HTPB is sensitive to pressure at high temperatures. At pressures lower than 1 GPa, thermal degradation primarily occurs via unimolecular depolymerization β-scission reactions which leads to formation of butadiene. However, at pressures higher than 5 GPa, activated HTPB molecules undergo cyclization and chain cross-linking reactions. Due to this mechanistic shift, butadiene formation is suppressed at high pressures. Instead of butadiene, thermal degradation products include heavy (heavier than HTPB) branched polymers and cyclic compounds like 4-vinyl cyclohexene, cyclohexadiene and cyclopentadiene. Like any molecular study, we expect some alterations in our predictions with a change in force field parameter set. For example, we anticipate that the actual energy transfer rate at a fuel/binder interface will be lower than the one predicted by RMD simulations because the force field parameters over-predict the HTPB heat capacity and system pressure. Similarly, we anticipate that our RMD simulations over predict the deflagration speed due to higher pressure. However, the overall thermochemical trends observed from the RMD simulations will remain unchanged with different set of force field parameters. The study shows that non-homogenous


Page 36 of 41

Page 37 of 41 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


fuel/binder interfaces have considerable influence on the hot spot growth and chemical stability of PBX. Acknowledgment The authors acknowledge discussions with colleagues from different institutions during this project. Among many valuable inputs, discussions on the fundamental physics and chemistry of energetics materials with Thomas Sewell, D. Scott Stewart, Joseph Zaug, Dana Dlott, Nick Glumac, Betsy Rice, and Igor Schweigert have enlightened our views. Funding from DTRA grant number HDTRA1-15-1-0034 is used to perform the research. DOD supercomputing hours and help from Ernie Satterwhite in obtaining the necessary allocations on time are also gratefully acknowledged. Supporting Information Further details about interface formation, multigrain PBX simulation, Arrhenius plot of pyrolysis of HTPB References 1. Baer, M. R., Modeling Heterogeneous Energetic Materials at the Mesoscale. Thermochim Acta 2002, 384, 351-367. 2. Chidester, S. K.; Tarver, C. M.; Green, L. G.; Urtiew, P. A., On the Violence of Thermal Explosion in Solid Explosives. Combust Flame 1997, 110, 264-280. 3. Hooper, J. B.; Bedrov, D.; Smith, G. D.; Hanson, B.; Borodin, O.; Dattelbaum, D. M.; Kober, E. M., A Molecular Dynamics Simulation Study of the Pressure-Volume-Temperature Behavior of Polymers under High Pressure. J Chem Phys 2009, 130, 144904. 4. Millett, J. C. F.; Bourne, N. K.; Akhavan, J., The Response of Hydroxy-Terminated Polybutadiene to One-Dimensional Shock Loading. J Appl Phys 2004, 95, 4722-4727. 5. Sorescu, D. C.; Rice, B. M., Rdx Compression, Α→ Γ Phase Transi>on, and Shock Hugoniot Calculations from Density-Functional-Theory-Based Molecular Dynamics Simulations. J Phys Chem C 2016, 120, 19547-19557. 6. Munday, L. B.; Chung, P. W.; Rice, B. M.; Solares, S. D., Simulations of High-Pressure Phases in Rdx. J Phys Chem B 2011, 115, 4378-86. 7. Rahul; De, S., A Phase-Field Model for Shock-Induced Α-Γ Phase Transition of Rdx. Int J Plast 2017, 88, 140-158. 8. Cawkwell, M. J.; Ramos, K. J.; Hooks, D. E.; Sewell, T. D., Homogeneous Dislocation Nucleation in Cyclotrimethylene Trinitramine under Shock Loading. J Appl Phys 2010, 107, 063512.



9. Mathew, N.; Sewell, T. D., Generalized Stacking Fault Energies in the Basal Plane of Triclinic Molecular Crystal 1,3,5-Triamino-2,4,6-Trinitrobenzene (Tatb). Philos Mag 2015, 95, 424-440. 10. An, Q.; Cheng, T.; Goddard, W. A.; Zybin, S. V., Anisotropic Impact Sensitivity and Shock Induced Plasticity of Tkx-50 (Dihydroxylammonium 5,5′-Bis(Tetrazole)-1,1′-Diolate) Single Crystals: From LargeScale Molecular Dynamics Simulations. J Phys Chem C 2015, 119, 2196-2207. 11. Wood, M. A.; Cherukara, M. J.; Kober, E. M.; Strachan, A., Ultrafast Chemistry under Nonequilibrium Conditions and the Shock to Deflagration Transition at the Nanoscale. J Phys Chem C 2015, 119, 22008-22015. 12. Shan, T.-R.; Wixom, R. R.; Thompson, A. P., Extended Asymmetric Hot Region Formation Due to Shockwave Interactions Following Void Collapse in Shocked High Explosive. Phy Rev B 2016, 94, 054308. 13. Robertson, D. H.; Brenner, D. W.; White, C. T., Split Shock Waves from Molecular Dynamics. Phys Rev Lett 1991, 67, 3132-3135. 14. Cawkwell, M. J.; Sewell, T. D.; Zheng, L.; Thompson, D. L., Shock-Induced Shear Bands in an Energetic Molecular Crystal: Application of Shock-Front Absorbing Boundary Conditions to Molecular Dynamics Simulations. Phy Rev B 2008, 78. 15. Rai, N. K.; Udaykumar, H. S., Mesoscale Simulation of Reactive Pressed Energetic Materials under Shock Loading. J Appl Phys 2015, 118, 245905. 16. Barua, A.; Kim, S.; Horie, Y.; Zhou, M., Ignition Criterion for Heterogeneous Energetic Materials Based on Hotspot Size-Temperature Threshold. J Appl Phys 2013, 113, 064906. 17. Barua, A.; Kim, S. P.; Horie, Y.; Zhou, M., Computational Analysis of Ignition in Heterogeneous Energetic Materials. Mater Sci Forum 2013, 767, 13-21. 18. Long, Y.; Chen, J., A Molecular Dynamics Study of the Early-Time Mechanical Heating in ShockLoaded Octahydro-1,3,5,7-Tetranitro-1,3,5,7-Tetrazocine-Based Explosives. J Appl Phys 2014, 116, 033516. 19. Zeldovich, Y. B., Regime Classification of an Exothermic Reaction with Nonuniform Initial Conditions. Combust Flame 1980, 39, 211-214. 20. Tarver, C. M.; Chidester, S. K.; Nichols, A. L., Critical Conditions for Impact- and Shock-Induced Hot Spots in Solid Explosives. J Phys Chem 1996, 100, 5794-5799. 21. Zaug, J. M.; Young, C. E.; Long, G. T.; Maienschein, J. L.; Glascoe, E. A.; Hansen, D. W.; Wardell, J. F.; Black, C. K.; Sykora, G. B., Deflagration Rates of Secondary Explosives under Static Mpa-Gpa Pressure. AIP Conference Proceedings 2009, 1195, 420-423. 22. Senftle, T. P., et al., The Reaxff Reactive Force-Field: Development, Applications and Future Directions. Npt Comp Mat 2016, 2, 15011. 23. Shan, T.; van Duin, A. C. T.; Thompson, A. P., Development of a Reaxff Reactive Force Field for Ammonium Nitrate and Application to Shock Compression and Thermal Decomposition. J Phys Chem A 2014, 118, 1469-1478. 24. Strachan, A.; van Duin, A. C.; Chakraborty, D.; Dasgupta, S.; Goddard III, W. A., Shock Waves in High-Energy Materials: The Initial Chemical Events in Nitramine Rdx. Phys Rev Lett 2003, 91, 098301. 25. Li, Y.; Kalia, R. K.; Nakano, A.; Nomura, K.; Vashishta, P., Multistage Reaction Pathways in Detonating High Explosives. Appl Phys Lett 2014, 105, 204103. 26. Joshi, K. L.; Chaudhuri, S., Reactive Simulation of the Chemistry Behind the Condensed-Phase Ignition of Rdx from Hot Spots. Phys Chem Chem Phys 2015, 17, 18790-801. 27. Sergeev, O. V.; Yanilkin, A. V., Hydrogen Transfer in Energetic Materials from Reaxff and Dft Calculations. J Phys Chem A 2017, 121, 3019-3027. 28. Zhou, T.; Song, H.; Liu, Y.; Huang, F., Shock Initiated Thermal and Chemical Responses of Hmx Crystal from Reaxff Molecular Dynamics Simulation. Phys Chem Chem Phys 2014, 16, 13914-13931.


Page 38 of 41

Page 39 of 41 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


29. Sergeev, O. V.; Yanilkin, A. V., Molecular Dynamics Simulation of Combustion Front Propagation in a Petn Single Crystal. Combust. Explos. Shock Waves 2014, 50, 323-332. 30. Maillet, J. B.; Bourasseau, E.; Desbiens, N.; Vallverdu, G.; Stoltz, G., Mesoscopic Simulations of Shock-to-Detonation Transition in Reactive Liquid High Explosive. Europhys Lett 2011, 96, 68007. 31. Brennan, J. K.; Lisal, M.; Moore, J. D.; Izvekov, S.; Schweigert, I. V.; Larentzos, J. P., Coarse-Grain Model Simulations of Nonequilibrium Dynamics in Heterogeneous Materials. J Phys Chem Lett 2014, 5, 2144-9. 32. Fröhlich, M. G.; Sewell, T. D.; Thompson, D. L., Molecular Dynamics Simulations of Shock Waves in Hydroxyl-Terminated Polybutadiene Melts: Mechanical and Structural Responses. J Chem Phys 2014, 140, 024902. 33. Xie, F.; Lu, Z.; Yang, Z.; Hu, W.; Yuan, Z., Mechanical Behaviors and Molecular Deformation Mechanisms of Polymers under High Speed Shock Compression: A Molecular Dynamics Simulation Study. Polymer 2016, 98, 294-304. 34. Tingfa, D., Thermal Decomposition Studies of Solid Propellant Binder Htpb. Thermochim Acta 1989, 138, 189-197. 35. Chiaverini, M. J.; Harting, G. C.; Lu, Y.-C.; Kuo, K. K.; Peretz, A.; Jones, H. S.; Wygle, B. S.; Arves, J. P., Pyrolysis Behavior of Hybrid-Rocket Solid Fuels under Rapid Heating Conditions. J Propul Power 1999, 15, 888-895. 36. Arisawa, H.; Brill, T. B., Flash Pyrolysis of Hydroxyl-Terminated Polybutadiene (Htpb) I: Analysis and Implications of the Gaseous Products. Combust Flame 1996, 106, 131-143. 37. Ganesh, K.; Sundarrajan, S.; Kishore, K.; Ninan, K. N.; George, B.; Surianarayanan, M., Primary Pyrolysis Products of Hydroxy-Terminated Polybutadiene. Macromolecules 2000, 33, 326-330. 38. Shi, Y.; Brenner, D. W., Molecular Simulation of the Influence of Interface Faceting on the Shock Sensitivity of a Model Plastic Bonded Explosive. J Phys Chem B 2008, 112, 14898-14904. 39. An, Q.; Zybin, S. V.; Goddard III, W. A.; Jaramillo-Botero, A.; Blanco, M.; Luo, S.-N., Elucidation of the Dynamics for Hot-Spot Initiation at Nonuniform Interfaces of Highly Shocked Materials. Phy Rev B 2011, 84. 40. An, Q.; Goddard III, W. A.; Zybin, S. V.; Jaramillo-Botero, A.; Zhou, T., Highly Shocked Polymer Bonded Explosives at a Nonplanar Interface: Hot-Spot Formation Leading to Detonation. J Phys Chem C 2013, 117, 26551-26561. 41. Yan, Q.-L.; Zeman, S.; Sánchez Jiménez, P. E.; Zhang, T.-L.; Pérez-Maqueda, L. A.; Elbeih, A., The Mitigation Effect of Synthetic Polymers on Initiation Reactivity of Cl-20: Physical Models and Chemical Pathways of Thermolysis. J Phys Chem C 2014, 118, 22881-22895. 42. Joshi, K.; Losada, M.; Chaudhuri, S., Intermolecular Energy Transfer Dynamics at a Hot-Spot Interface in Rdx Crystals. J Phys Chem A 2016, 120, 477-489. 43. Fröhlich, M. G.; Sewell, T. D., Pivot Algorithm and Push-Off Method for Efficient System Generation of All-Atom Polymer Melts: Application to Hydroxyl-Terminated Polybutadiene. Macromol Theory Simul 2013, 22, 344-353. 44. Martinez, L.; Andrade, R.; Birgin, E. G.; Martinez, J. M., Packmol: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J Comput Chem 2009, 30, 2157-64. 45. Dreger, Z. A.; Gupta, Y. M., Phase Diagram of Hexahydro-1,3,5-Trinitro-1,3,5-Triazine Crystals at High Pressures and Temperatures. J Phys Chem A 2010, 114, 8099-8105. 46. Joshi, K.; Chaudhuri, S., Observation of Deflagration Wave in Energetic Materials Using Reactive Molecular Dynamics. Combust Flame 2017, 184, 20-29. 47. Liu, L.; Liu, Y.; Zybin, S. V.; Sun, H.; Goddard III, W. A., 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.



48. Plimpton, S., Fast Parallel Algorithms for Short-Range Molecular Dynamics. J Comput Phys 1995, 117, 1-19. 49. Chenoweth, K.; van Duin, A. C. T.; Goddard III, W. A., Reaxff Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J Phys Chem A 2008, 112, 1040-1053. 50. Wood, M. A.; van Duin, A. C. T.; 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. 51. Wang, Q.-D.; Wang, J.-B.; Li, J.-Q.; Tan, N.-X.; Li, X.-Y., Reactive Molecular Dynamics Simulation and Chemical Kinetic Modeling of Pyrolysis and Combustion of N-Dodecane. Combust Flame 2011, 158, 217-226. 52. Joshi, K. L.; Raman, S.; van Duin, A. C. T., Connectivity-Based Parallel Replica Dynamics for Chemically Reactive Systems: From Femtoseconds to Microseconds. J Phys Chem Lett 2013, 4, 37923797. 53. Castro-Marcano, F.; Kamat, A. M.; Russo, M. F.; van Duin, A. C. T.; Mathews, J. P., Combustion of an Illinois No. 6 Coal Char Simulated Using an Atomistic Char Representation and the Reaxff Reactive Force Field. Combust Flame 2012, 159, 1272-1285. 54. Nanda, V. S.; Simha, R., Theoretical Interpretation of Tait Equation Parameters. J Chem Phys 1964, 41, 1884-1885. 55. Cai, W.; Thakre, P.; Yang, V., A Model of Ap/Htpb Composite Propellant Combustion in RocketMotor Environments. Combust Sci Technol 2008, 180, 2143-2169. 56. Tsolou, G.; Harmandaris, V. A.; Mavrantzas, V. G., Temperature and Pressure Effects on Local Structure and Chain Packing Incis-1,4-Polybutadiene from Detailed Molecular Dynamics Simulations. Macromol Theory Simul 2006, 15, 381-393. 57. Sharma, P.; Roy, S.; Karimi-Varzaneh, H. A., Validation of Force Fields of Rubber through GlassTransition Temperature Calculation by Microsecond Atomic-Scale Molecular Dynamics Simulation. J Phys Chem B 2016, 120, 1367-79. 58. Allen, M. P.; Tildesley, D. J., Computer Simulation of Liquids; Clarendon Press, 1989, p 385. 59. Thomas, R.; Boudenne, A.; Ibos, L.; Candau, Y.; Thomas, S., Thermophysical Properties of Ctbn and Htpb Liquid Rubber Modified Epoxy Blends. J Appl Polym Sci 2010, NA-NA. 60. Miller, M. S., Thermophysical Properties of Rdx. Lab, A. R., Ed. Aberdeen Proving Ground, MD, 1997. 61. Beckstead, M. W.; Puduppakkam, K.; Thakre, P.; Yang, V., Modeling of Combustion and Ignition of Solid-Propellant Ingredients. Prog Energy Combust Sci 2007, 33, 497-551.

TOC Figure


Page 40 of 41

Page 41 of 41 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



Hot Spot Interaction With HTPB Binder in Energetic Composites - The

Recommend Documents