RDX/BAMO推进剂结合能、力学性质和能量性能的分子动力学模拟

李苗苗 沈瑞琪 李凤生

(南京理工大学化工学院,南京210094)

RDX/BAMO推进剂结合能、力学性质和能量性能的分子动力学模拟

李苗苗*沈瑞琪 李凤生

(南京理工大学化工学院,南京210094)

利用分子动力学方法研究了著名的含能材料环三亚甲基三硝胺(RDX)、3,3′-双-(叠氮甲基)-氧杂环丁烷(BAMO)和RDX/BAMO推进剂.结果表明,BAMO与RDX(010)面之间分子相互作用最强,其次是(100)和(001)面.以对相关函数g(r)描述了RDX和BAMO之间的相互作用.计算了RDX/BAMO推进剂的弹性系数、模量、柯西压、泊松比等性能.结果表明,BAMO的加入能够改善RDX的弹性力学性能,相对改善效应的顺序为(100)＞

(001)＞(010).RDX/BAMO推进剂的能量性能结果显示,BAMO的加入降低了RDX的比冲,但仍高于著名的双基推进剂的比冲.

分子动力学;RDX;BAMO;结合能;力学性能

1 Introduction

Fig.1 Molecular configurations of RDX(a)and crystal(b)

Fig.2 Molecular configurations of BAMO

Cyclotrimethylene trinitramine(RDX,Fig.1),an important modern molecular explosive,has been widely used as the propellant ingredients in applications such as gun and rocket motors.It offers many advantages for advanced propulsion,including high energy,high specific impulse,low sensitivity,and environment friendliness.In particular,RDX is also attractive because of its low cost.However,RDX has low burning rate (～1.19 cm·s-1at 6.08 MPa)coupled with a relatively high and undesirable burning rate exponent(～0.74).1In the search for new binder materials for the nitramine propellant,a great deal of interest has recently centered on azido compounds,which contain―N＝N≡N bonds in their chemical structure,such as 3,3′-bis-azidomethyl-oxetane(BAMO,Fig.2).This nitramine/ azide combination propellant is characterized by high velocity of combustion,low vulnerability,and good thermal stability. These characteristics are especially useful as rocket propellants,since they provide high specific impulse while generating minimum smoke.Heretofore,a number of experimental measurements and theoretical studies2-12have paid attention to this kind of propellant.However,previous studies mainly focused on thermal decomposition and combustion characteristics of nitramine/azide propellant,and there are few reports on the structure performance for nitramine/azide using the MD simulation.In this paper,we chose RDX/BAMO propellant as an example to investigate the correlations between the structure and performance for propellants.It is hoped that our studies provide some information and guidance for composite propellant formulation design.

2 Modeling and simulation

2.1 Computational model

The initial models were built by using materials studio(MS) package.13The initial RDX structure used in the condensed phase simulations was taken from the experimental result by Choi and Prince,14which crystallized in the orthorhombic space group of Pbca with four independent lattice parameters a=1.3182 nm,b=1.1574 nm,c=1.0709 nm and α=β=γ=90°. There are eight irreducible molecules in the unit cell,see Fig.1 (b).The RDX crystal was cut along three crystalline surfaces

(100),(010),and(001),respectively,with the“cutting”method in MS 3.0.1 and put in the periodic cells with 3.0 nm vacuum layer along the z-axis(c direction),and the periodic MD simulation cells contain 96 RDX molecules,corresponding to (2×2×3)unit cells.BAMO involving six-chain segments was processed by amorphous cell module and simulated by 2.5×106fs using the MD method to get the equilibrium conformation, and the end groups of BAMO polymer were saturated with H and OH groups.According to the practical formulation of nitramine/binder propellant,we embed five equilibrium chains of BAMO into the supercell model in parallel with different RDX crystalline surfaces.A total of three different initial RDX/BAMO configurations with 2631 atoms,and with mass fractions of 80.6%RDX and 19.4%BAMO,respectively,were obtained.

2.2 MD simulations

The COMPASS force field was used to study the structures and properties of the RDX,BAMO,and RDX/BAMO.Its parameters were debugged and ascertained from the ab initio calculations,optimized according to the experimental values,and parameterized using extensive data for molecules in condensed phase.Its nonbonded parameters were further amended and validated by the thermal physical properties of the molecules in liquid and solid phases obtained using the MD method.Consequently,COMPASS is able to accurately predict the structural,conformational,vibrational,and thermophysical properties for a broad range of compounds in both isolation and condensed phases.15-18Moreover,the nitro(―NO2)functional groups required to model energetic materials have been specifically parameterized and included in the COMPASS force field.17Up to now,this force field has been successfully employed to investigate the nitramine.19-25It is therefore suitable for performing MD simulations on them.

The molecular dynamics simulations of RDX,BAMO,and RDX/BAMO were performed using the COMPASS force field and periodic boundary conditions.Minimizations were initially carried out for 5000 iterations to equilibrate the RDX/BAMO models and then the simulation boxes were compressed slightly(0.3%)along the c direction.Afterward,another 5000 iterations of minimizations were carried out to reach the equilibrium state and the boxes were compressed further along the c direction.The above processes were repeated step by step until the density approached the theoretical maximum value,which can be predicted according to the mass fraction of each component in the propellant.Starting from the above-minimized structures,the MD simulations were conducted at constant volume and constant temperature(NVT)conditions.After an equilibrium run,the module allowed one to collect the results of the dynamics simulation in a trajectory file.Through analyzing trajectory files,the static elastic properties and pair correlation functions were obtained.Considering the condition of equilibrium and spend of CPU time,all the simulation time was added to 4×105fs.A fixed time step size of 1 fs was used in all the cases.In the above-mentioned MD simulations,the Andersen thermostat method26was employed to control the system at the temperature of 298 K.The Coulomb and van der Waals longrange,nonbond interactions were handled by using the standard Ewald and atom-based summation methods,respectively. Nonbonded interactions,spline width,and buffer width were truncated at 0.95,0.1,and 0.05 nm,respectively.All the calculations were implemented on a Pentium IV computer.

3 Results and discussion

3.1 Criteria of system equilibrium

There are two criteria to judge the equilibrium:one is the equilibrium of temperature and the other is the equilibrium of energy.The fluctuations of temperature and energy are in the range of 5%-10%,that is to say,the fluctuation of temperature is within±15 K and the energy is invariable or small fluctuation around the average energy value.For instance,Fig.3 shows the fluctuation curves of temperature and energy in the MD simulation of BAMO on the molecule layers parallel to (100)crystalline surface of RDX.It can be found that both the fluctuation curves of temperature and energy trend to be smooth,and the temperature reaches equilibrium state indeed as it is fluctuating 10 K or so.

3.2 Equilibrium configuration and interactions between constituents

After the MD simulations,one can get the equilibrium configurations of the models.Fig.4 illustrates the equilibrium configurations of RDX/BAMO with BAMO on RDX crystalline surfaces of(100),(010),(001),respectively.As can be seen from Fig.4,BAMO polymer binder is closely contacted with the RDX crystalline surface and consequently extensive interactions exist between BAMO polymer and RDX.Binding energy(Ebinding)can well reflect the capacity of the two components blending with each other,which is defined as the negative value of the intermolecular interaction energy(Einter),Ebinding=-Einter. The intermolecular interaction energy can be evaluated by the total energies of the composite and each component in the equilibrium state.So the Ebindingbetween RDX and BAMO polymer can be determined as follows:

where Etotalis the total energy of RDX/BAMO mixed system. ERDXand EBAMOare the total energies of RDX and BAMO polymer,respectively.

For visualization,Etotal,ERDX,EBAMO,and Ebindingfor different crystalline surfaces are presented in Table 1.It can be found that the RDX(010)surface has a stronger capability to blend with BAMO than(001)and(100)surfaces due to the larger binding energies.In other words,when BAMO polymers are put into the RDX crystal,they tends to concentrate on the RDX(010)surface due to their stronger intermolecular interactions.

Fig.3 Plots of temperature and energy vs simulation time for RDX/BAMO with BAMO on RDX(100)crystalline surface at the temperature of 298 K

Fig.4 Equilibrium structures of RDX/BAMO with BAMO on different crystalline surfaces of RDX at 298 K

Table 1 Binding energies for RDX/BAMO with BAMO on three different crystalline surfaces of RDX at 298 K

The interactions between each constituent can be further analyzed by examining pair correlation function g(r).The pair correlation function gives a measure of the probability of finding another atom at a distance r from a specific atom.It has many applications in structural investigations of both solid and liquid packings(local structure),in studying specific interactions such as hydrogen bonding,and in statistical mechanical theories of liquids and mixtures.

To be compact,only g(r)-r relations with important interaction of the crystalline surfaces(100)are described in Fig.5. The corresponding atoms are named as follows:(1)O,N,and H atoms in RDX are noted as O1,N1,and H1,respectively;(2) O(―OH),N,and H in BAMO are noted as O2H,N2,and H2,respectively.

In general,intermolecular actions include hydrogen bonding action and van der Waals(vdw)force,in which the vdw force is composed of dipole-dipole,induction,and dispersion force. If the distance between atoms is 0.26-0.31 nm,0.31-0.50 nm, or above 0.50 nm,the interaction belongs to hydrogen bonding,strong vdw,or weak vdw force,respectively.Although the hydrogen bonding action is weaker than chemical bond,it is the strongest force among intermolecular actions and can strengthen them.

From Fig.5(a),it is found that probability for N1in RDX and O2Hin BAMO to simultaneously arise in the distance of 0.29-0.39 nm is high to 2.5 or so,predicting the strong van der Waals interaction between them.In the region of 0.29-0.40 nm,a comparative high peak arises in the g(r)describing theatom pair(Fig.5(b)),predicting the strong van der Waals interaction.As seen from Fig.5(c),in the region of 0.25-0.29 nm,the high peak value arises in the g(r)curve ofindicating the strong hydrogen bond interaction between this atom pair;in the region of van der Waals,high peak arises again,predicting that certain van der Waals interaction exists between them.Fig.5(d)shows that mainly hydrogen bond and van der Waals interactions exist in O1-H2.In all,van der Waals and hydrogen bonds are the main interactions between RDX and BAMO,especially the van der Waals and hydrogen bond interactions exist in N1-O2Hand O1

-H2.

3.3 Mechanical properties

The material stress and strain tensor are denoted by σ and ε, respectively.From the statistical mechanics of elasticity,27the generalized Hooke′s law is often written as

[Cij]is symmetric,and hence a maximum of 21 constants is required to fully describe the stress-strain behavior of an arbitrary material.

The effective isotropic compliances in terms of single-crystal compliances averaged over all orientations can be obtained by Reuss average.The effective bulk and shear moduli are given by:

The subscript R denotes the Reuss average.The compliance matrix S is equal to the inverse matrix of elastic coefficient matrix C,i.e.,S=C−1.Note that for the most general crystal struc-ture(all 21 constants are independent)the Reuss modulus depends on only nine of the single-crystal compliances.From the rules of isotropic linear elasticity we have

Fig.5 Pair correlation function of RDX(100)/BAMO propellant

E=2G(1+n)=3K(1-2n) (5) where E is tensile modulus and n is the Poisson′s ratio,so that after the bulk and shear moduli are calculated,the tensile modulus and Poisson′s ratio can be obtained.

Such plastic properties as hardness,tensile strength,fracture strength,and elongation in tension,can be related to the elastic modulus.28Hardness and tensile strength representing the resistance to plastic deformation are proportional to the shear modulus G.Fracture strength is proportional to the bulk modulus K. The quotient K/G indicates the extent of the plastic range(elongation in tension),so that a high value of K/G is associated with ductibility and a low value with brittleness.

The predicted value of elastic constants and effective isotropic mechanical properties(tensile modulus,bulk modulus, shear modulus,and Poisson′s ratio)are summarized in Table 2 and Table 3,respectively.

From the elastic coefficients matrices for RDX in Table 2,it can be found that the diagonal elements Ciiand C12,C13,C23(nine elements totally)are larger but the other coefficients are smaller,which indicate that the crystal has anisotropy to some extent.In addition,the larger C11,C22,and C33imply that,to reach the same strain,RDX needs a larger stress.This characteristic can also be further validated by the differences of elastic coefficients of the RDX/BAMO propellant with BAMO on different crystalline surfaces.Compared with the pure RDX, the smaller diagonal elements C22and C55increase for RDX/ BAMO,while other diagonal elements Ciiall decrease.Moreover,the off-diagonal elements C13and C23increase while C12, C15,C25,C35,and C46decrease or increase to some extent.This evolution tendency of elastic coefficients shows that adding some amounts of BAMO can reduce the anisotropy of the system.

Cauchy pressure(C12-C44)can be used as a criterion to evaluate the ductibility or brittleness of a material.As a rule,the value of(C12

-C44)for a ductile material is positive;on the contrary,that is negative for a brittle material.Meanwhile,the more positive the(C12-C44)value is,the more ductile the material is.According to this,the data in the last column of Table 2 indicate that the pure RDX and the obtained RDX/BAMO are all ductile due to their positive(C12-C44).But,the Cauchy pressures of the RDX/BAMO are all larger than that of the pure RDX crystal.This indicates that the ductibility of RDX is greatly improved by adding some quantities of BAMO polymer. Comparing the values of(C12-C44),we find that the ductibility of RDX/BAMO depends on different surfaces,and it changes in the following order:(010)＞(100)＞(001).

As can be also seen from Table 3,all moduli of the obtained RDX/BAMO decrease in comparison with the pure crystal except K for(010).For example,the tensile modulus E decreases from 8.51 to 6.38 GPa as the system changes from the pure crystal to the RDX/BAMO mixture,the shear modulus G decreases from 3.29 to 2.37 GPa,and the bulk modulus K decreases or increases slightly as compared to the former two moduli.This further indicates that the rigidity and brittleness of RDX/BAMO propellant decrease,while its elasticity and plasticity strengthen.Meanwhile,these variations also suggest that the materials will deform more easily when they are subjected to an external force,because the resistance to plastic deformation is proportional to the elastic shear modulus G and the fracture strength for materials is proportional to the bulk modulus K.29As a whole,the effect of BAMO on improving the mechanical properties of different crystal surfaces is somewhat different and changes in the order of(100)＞(001)＞(010).

In addition,the ratio(bulk modulus/shear modulus,K/G) can be used to evaluate the tenacity of a material.A higher value of K/G is associated with malleability and a lower value with brittleness.According to this,it can be deduced from the K/G values in Table 3 that RDX/BAMO has better malleability than the pure RDX.On the whole,the malleability of RDX/ BAMO propellant with BAMO on three different crystalline surfaces decreases as follows:(100)＞(001)≈(010).

3.4 Energetic properties

Energetic properties are the much important factors to evaluate the propellant performances.The calculations of propellant energetic properties can be approximately divided into foursteps.Firstly,calculate the assumed chemical formula,oxygen balances(OB100)and initial total enthalpy of 1000 g propellant. For1000 g RDX,the assumed chemicalformula is C13.506H27.013O27.013N27.013,and for 1000 g RDX/BAMO(80.6/19.4) propellant,that is C16.659H31.011O22.927N28.701.Next,compute the thermochemical properties of propellant combustion process in the combustion chamber.Propellant isenthalpic combustion decomposes into high temperature working fluid,which is heat balance and chemical equilibrium(energy conservation and mass conservation).According to the principle of minimum free energy function of balance system,the equilibrium composition,chamber temperature(Tc)and other thermodynamic functions in the combustion chamber can be further calculated. Thirdly,calculate the thermochemical properties of combustion products expansion process(isentropic expansion)in the nozzle.Using the entropy difference inserting method,one can get the outlet temperature(Te),outlet species composition,etc. Lastly,calculate specific heat ratio,standard theoretical specific impulse(Isp),etc.

Table 2 Elastic constants(GPa)of pure crystal RDX and RDX/BAMO propellant with BAMO on different crystalline surfaces at 298 K

Table 3 Elastic properties of RDX crystal and RDX/BAMO mixture system with BAMO on different crystalline surfaces at 298 K

Table 4 Energetic properties of RDX/BAMO propellant

Among above parameters,standard theoretical specific impulse(Isp),equal to units of thrust per unit mass of propellant consumed per unit time,is the most important parameter to evaluate energy characteristics of propellants.Isphas great influence on the rocket range,the higher Isp,the farther rocket range.The calculation standard of Ispis as follows:(1)the combustion chamber pressure equals to 6.86 MPa;(2)environment pressure equals to 0.1 MPa;(3)nozzle expansion ratio is optimal,that exit pressure equals to environment pressure;(4)nozzle exit spread angle equals to 0°.The predicted energetic properties are listed in Table 4.

From Table 4,it can be seen that Ispof RDX is relatively high(2608 N·s·kg-1)due to its high density(ρ),oxygen balances(OB100),and chamber temperature.Compared with RDX, BAMO has positive heat of formation(ΔHf,2460.0 kJ·kg-1), yet it also has large negative oxygen balances(OB100=-123.8), which will cause the incomplete combustion and reduce the chamber temperature,and so BAMO has lower Isp(2137 N·s· kg-1).When adding BAMO polymer to RDX crystal,the OB100, Tc,and Ispof RDX/BAMO propellant are smaller than those of the pure RDX crystal.But,compared with the famous double base(DB)propellant(Isp=2100-2270 N·s·kg-1)30,it is obvious that RDX/BAMO propellant has better energetic performances than DB propellant.

4 Conclusions

In this study,we have performed MD simulations to study the binding energies,mechanical properties,and energetic properties of RDX/BAMO propellant.The effects of BAMO on different crystalline surfaces of RDX were investigated,the major findings can be summarized as follows:

(1)The binding energies between RDX and BAMO are obtained,and the order of binding energies is(010)＞(100)＞(001).

(2)The pair correlation function g(r)analysis shows that hydrogen bond and van der Waals interactions exist between RDX and BAMO.

(3)The mechanical properties of the pure RDX can be effectively improved by adding BAMO ploymer.On three different cyrstalline surfaces,the effect of BAMO on improving the mechanical properties is approximately(100)＞(001)＞(010). Whereas the improvement in the ductibility(C12-C44)made by BAMO polymers changes in the sequence of(010)＞(100)＞

(001).

(4)The calculations on energetic properties for RDX/ BAMO propellant show that its energetic properties can be lowered by blending some amount of BAMO polymer,but it is still superior to the double base propellant and can be used as advanced propellant.

In a word,the MD simulations on the RDX and RDX/ BAMO propellant provide us much information about its mechanical properties,binding energies,and energetic properties. These may be useful for guiding composite propellant design.

(1) Luca,L.D.;Cozzi,F.;Germiniasi,G.Combust.Flame 1999, 118,248.

(2) Oyumi,Y.;Brill,T.B.Combust.Flame 1986,65,127.

(3) Chen,J.K.;Brill,T.B.Combust.Flame 1991,87,157.

(4) Miyazaki,T.;Kubota,N.Propell.Explos.Pyrotech.1992,17,5.

(5)Oyumi,Y.;Inokami,K.;Yamazaki,K.;Matsumoto,K.Propell. Explos.Pyrotech.1993,18,62.

(6)Shen,S.M.;Chiu,Y.S.;Wang,S.W.;Chen,S.I.Thermochim. Acta 1993,221,275.

(7) Kimura,E.;Oyumi,Y.Propell.Explos.Pyrotech.1995,20,322.

(8) Kubota,N.J.Propul.Power 1995,11,677.

(9) Liu,Y.L.;Hsiue,G.H.;Chiu,Y.S.J.Appl.Polym.Sci.1995, 58,579.

(10) Oyumi,Y.;Kimura,E.;Nagayama,K.Propell.Explos. Pyrotech.1998,23,123.

(11) Pisharath,S.;Ang,H.G.Polym.Degrad.Stabil.2007,92,1365.

(12) Zhai,J.;Yang,R.;Li,J.Combust.Flame 2008,154,473.

(13) Material Studio 3.0 discover;Accelrys Inc.:San Diego,2004.

(14) Choi,C.S.;Prince,E.Acta Crystallogr.B 1972,28,2857.

(15) Sun,H.;Ren,P.;Fried,J.R.Comput.Theor.Polym.Sci.1998, 8,229.

(16) Bunte,S.W.;Sun,H.J.Phys.Chem.B 2000,104,2477.

(17)Yang,J.;Ren,Y.;Tian,A.M.;Sun,H.J.Phys.Chem.B 2000, 104,4951.

(18)Mcquaid,M.J.;Sun,H.;Rigby,D.J.Comput.Chem.2004,25, 61.

(19) Sun,H.J.Phys.Chem.B 1998,102,7338.

(20) Zhu,W.;Xiao,J.;Zhu,W.;Xiao,H.J.Hazard.Mater.2009, 164,1082.

(21) Xu,X.J.;Xiao,H.M.;Xiao,J.J.;Zhu,W.;Huang,H.;Li,J.S. J.Phys.Chem.B 2006,110,7203.

(22)Qiu,L.;Zhu,W.H.;Xiao,J.J.;Zhu,W.;Xiao,H.M.;Huang, H.;Li,J.S.J.Phys.Chem.B 2007,111,1559.

(23)Zhu,W.;Wang,X.;Xiao,J.;Zhu,W.;Sun,H.;Xiao,H. J.Hazard.Mater.2009,167,810.

(24) Xiao,J.;Huang,H.;Li,J.;Zhang,H.;Zhu,W.;Xiao,H.J.Mol. Struct.-Theothem 2008,851,242.

(25) Qiu,L.;Xiao,H.J.Hazard.Mater.2009,164,329.

(26)Andersen,H.C.J.Chem.Phys.1980,72,2384.

(27) Weiner,J.H.Statistical Mechanics of Elasticity;John Wiley: New York,2002.

(28) Pugh,S.F.Philos.Mag.Series 7 1954,45,823.

(29) Weiner,J.H.Statistical Mechanics of Elasticity;John Wiley: New York,1983.

(30) Tian,D.Y.;Liu,J.H.Energetics Calculation of Chemical Propellants;Henan Scientific and Technical Publishers: Zhengzhou,1999.[田德余,刘剑洪.化学推进剂计算能量学.郑州:河南科学技术出版社,1999.]

January 21,2011;Revised:March 28,2011;Published on Web:April 15,2011.

Molecular Dynamics Simulation of Binding Energies,Mechanical Properties and Energetic Performance of the RDX/BAMO Propellant

LI Miao-Miao*SHEN Rui-Qi LI Feng-Sheng
(School of Chemical Engineering Nanjing University of Science and Technology,Nanjing 210094,P.R.China)

Molecular dynamics(MD)simulations were performed to investigate the well-known energetic material cyclotrimethylene trinitramine(RDX)crystal,3,3′-bis-azidomethyl-oxetane(BAMO)and the RDX/ BAMO propellant.The results show that the binding energies of RDX with BAMO on different crystalline surfaces change as follows:(010)＞(100)＞(001).The interactions between RDX and BAMO were analyzed by pair correlation functions g(r).The mechanical properties of the RDX/BAMO propellant,such as the elastic coefficients,modulus,Cauchy pressure,and Poissonʹs ratio,were obtained.We find that the mechanical properties are effectively improved by adding some BAMO polymer and the overall effect of BAMO on the three crystalline surfaces of RDX changes as follows:(100)＞(001)＞(010).The energetic performance of the RDX/BAMO propellant was also calculated and the results show that compared with the pure RDX crystal,the standard theoretical specific impulse(Isp)of the RDX/BAMO propellant decreases but it is still superior to that of the double base propellant.

Molecular dynamics;RDX;BAMO;Binding energy;Mechanical property

O641

∗Corresponding author.Email:lmm406@126.com;Tel/Fax:+86-25-84315942.

The project was supported by the Jiangsu Postdoctoral Sustentation Fund,China(0902018C).江苏省博士后基金(0902018C)资助项目