局部高温诱发CL-20/TNT共晶炸药内反应流传播的ReaxFF分子动力学模拟

刘 海, 杨 镇, 何远航

(1. 中国空气动力研究与发展中心超高速空气动力研究所, 四川 绵阳 621000; 2. 北京理工大学爆炸科学与技术国家重点实验室, 北京 100081)

1 引 言

当前共晶技术已广泛应用于生物、医药等领域[1],而在含能材料领域的研究尚处于起步阶段。含能材料的应用水平多依赖于其本身的密度、熔点、分子结构等,且共晶材料的属性可以明显地异于原始单一组分及各组分的物理混合物,因此通过共晶技术可以“取长补短”地修饰当前单一组分含能材料的不足,如感度低,稳定性好但爆压低; 或较好的能量输出,如高爆压、高爆速等,但摩擦、撞击感度高等物性的含能材料。

尽管有很多的潜在替代品,但三硝基甲苯(TNT)仍是非常重要的军用炸药,因为其价廉、对冲击和摩擦均不敏感且易于制备,而相对较低的爆速和爆压则制约了其在高能战斗部中的应用。六硝基六氮杂异伍兹烷(HNIW,CL-20)是具有笼型的高能量密度材料,爆速和爆压均明显优于TNT,但其安定性不够理想,感度较高,因此没有能够得到广泛应用。共晶技术在含能材料领域的应用可以修正这种单一组分表现出的不足。当前已经实现CL-20/TNT[2-3],CL-20/HMX(环四甲撑四硝胺)[4-5],CL-20/BTF(苯并三呋咱氧化物)[6]共晶以及HMX/TATB(三氨基三硝基苯)共晶[7-9],并且实验分析显示,共晶炸药具有较理想的能量输出和安定性。

含能材料晶体高温热分解过程研究的意义,在于通过对含能材料整体施加爆轰柴普曼-柔格(Chapman-Jouguet , CJ)点温度从侧面分析反应区内的化学反应动力学机制[10-12]。而含能材料的热点火过程则使得初始热量集中在某一区域。一旦局部高温诱发的燃烧加速并形成带有冲击波阵面的剧烈化学反应,爆轰则将快速的建立起来。面对含能材料热点火而提出的局部高温引发的化学反应,进而形成的反应流成长和传播问题,可以为预测含能材料非冲击起爆以及爆轰反应区的性质提供帮助[13]。当前利用分子动力学对局部高温引发热传播问题模拟的对象主要集中于金属材料[14-16],而局部高温诱发含能材料内反应流成长和传播的研究还相对较少[17-18]。Zhou通过原子模拟程序GRASP (General Reactive Atomistic Simulation Package)研究了HMX在热冲击和单轴压缩-热冲击耦合作用下的热力学及化学反应问题[17]。Sergeev结合ReaxFF分子动力学方法通过对季戊四醇四硝酸酯(PETN)含能材料局部施加高温模拟了热点(hot spot)成长过程中热流传播及化学反应情况,模拟结果表明PETN中热点温度为2000 K时,热流传播速度最大为0.32 km·s-1[18]。

当前ReaxFF分子动力学在极端条件下含能材料的物理化学响应方面做了很多卓有成效的工作[19-22],其在高能炸药热分解领域已经在初始化学反应路径和化学反应动力学方面提供定量的描述; 冲击起爆的ReaxFF分子动力学模拟也正在宏观爆速、爆压、爆轰波结构,以及微观分子振动模态,能量激发—发展—传播方面为爆轰过程提供信息,并且ReaxFF势函数本身也在不断的完善中,增加长程修正项(Elg)的ReaxFF/lg势函数在平衡密度计算方面具有优越性[23],当前的ReaxFF/lg势函数已经在含能材料领域得到了很好的应用[24-25]。采用大规模原子分子并行模拟器(Large-Scale Atomic/Molecular Massively Parallel Simulator, LAMMPS)同时结合ReaxFF/lg势函数分析了CL-20/TNT共晶炸药内局部高温诱发反应流传播的时空行为以及化学反应情况,利用此结果可为含能材料热点火以及热点成长等问题提供参考。

2 模型及计算细节

CL-20/TNT共晶初始单胞结构来自剑桥结构数据库(Cambridge Structural Database, CSD)[26],单胞内CL-20和TNT的分子数量比为1∶1,初始密度为1.91 g·cm-3。为了充分观察局部高温诱发反应流成长和传播的时空行为,构建20×2×1条形超胞,x/y/z方向的尺寸依次为193.4832 Å,38.738 Å,24.6898 Å。超晶胞内共含有640个分子,其中CL-20和TNT各为320个,体系内总计18240个原子。CL-20/TNT共晶结构以及物理模型如图1所示。

图1晶体结构及物理模型

Fig.1Crystal structure and physical model

首先将初始共晶晶胞通过结构优化,获得最低能量的晶体结构。随后通过NVE系综,同时结合Berendsen温度耦合方法调整系统温度与外界热浴保持一致,这里外界热浴设定为300 K,弛豫10 ps后获得常温下的晶格结构。其次,将CL-20/TNT沿x方向分为三个区域: Hot-1,Hot-2和Cold区,如图1所示。其中0～20 Å为Hot-1区,含有1888个原子,173.4832～193.4832 Å为Hot-2区,此区域内共有1892个原子,这两个区域沿x方向分别位于CL-20/TNT共晶晶胞的两端,并在上述区域施加高温,其空间尺寸均为20 Å×38.738 Å×24.6898 Å,而20～173.4832 Å为Cold区,总计含有14460个原子。区域划分完毕后,对Hot-1和Hot-2区同时施加3000 K热浴,快速的温升及温度控制将引发化学反应流的传播。为了控制温度在目标值附近波动,这里采用Berendsen温度耦合方法进行控温,同时外界热浴和系统温度耦合的特征时间为50 fs。而中间Cold区则采用NVE系综保持能量守恒。为了研究温度差异对含能材料初始热分解速率以及反应流传播规律的影响,另外一组则针对Hot-1和Hot-2区同时施加4000 K热浴,Cold区仍通过NVE系综保持能量守恒,其他参数设置相同。两组模拟的时间步长为0.1 fs,x方向为固定边界条件,而y,z方向均为周期边界条件,模拟时间持续100 ps。

3 共晶炸药内反应流传播的时空响应

3.1 局部高温诱发反应流形成的温度剖面

为了分析局部高温诱发共晶炸药内反应流的传播情况,这里采用分层技术[27]对Cold区域沿x方向(20～173.4832 Å)进行一维空间划分。将Cold区域分为153层,每层厚度约为1 Å,并且随时间输出各剖面的热力学以及动力学参数。

图2是两端热载荷为3000 K时,不同时刻Cold区各剖面的温度分布情况。对整体处于常温条件(300 K)的CL-20/TNT的两端(Hot-1,Hot-2区)施加高温热载荷后,Hot-1和Hot-2区快速升温,在升温的过程中,Hot区内含能材料分子快速发生分解反应(CL-20和TNT热分解演化曲线请见 (图5a,图5b)),热量由两端向中间Cold区传递形成热-机械波[17]。热量的传递使得Cold区域两端形成温度梯度(2800～400 K),但在Cold区两端明显的温度梯度区域内,并不会均发生热分解反应,因为低温区(400 K)内的含能材料分子,在前 10 ps,不会发生分解反应。原因可参照Goddard对实验室密度(1.68 g·cm-3)RDX进行1200 K高温热分解模拟时的情况,结果显示在100 ps后,热分解反应才开始[28]。而CL-20/TNT共晶在2800 K高温下,前10 ps内,分解反应已经开始[29]。所以Cold区内化学反应的剧烈程度由外向内依次递减。

a. 2～10 ps

b. 20～100 ps

图2热载荷为3000 K时不同时刻Cold区的温度剖面

Fig.2Temperature profile of cold region with the hot region of 3000 K at different time

在0～10 ps内,Hot区在升温的过程中由于热传导作用不断的向Cold区传递热量,而使得Hot区的温度并没有很快上升至目标温度(3000 K),而是在Cold区两端沿x方向40 Å宽的范围内形成2800～400 K的温度梯度,在Cold区域沿x方向60～140 Å的范围内,温度基本维持在400 K左右,并且Cold区两端的温度梯度基本对称。随着热量由两端不断的向中间区域传递,Cold区的温度也逐渐的上升,在20～80 ps内,Cold区左端的温度逐渐上升并维持在3000 K左右,而Cold区右端的温度则在2500～2800 K范围内波动,并且在20～100 ps期间,Cold区的温度梯度等呈不对称分布,至模拟结束时(t=100 ps),整个Cold区内的“低温区”出现在120 Å处,温度约为1800 K。出现这种不对称分布现象的可能原因为: 反应流传播过程中,Cold区两端粒子往复运动致使两端高温区域温度波动进而影响反应流传播的结果。

3.2 粒子密度分布

图3为两种热载荷条件下(3000,4000 K)反应流传播过程中不同时刻Cold区的一维密度分布图。在2 ps时刻,Cold区中间部分存在一个相对较高的密度区,并且热载荷为4000 K时,在Cold中间区域形成的高密度区范围(70～130 Å)大于3000 K时的情况(90～110 Å),但是高密度区的峰值基本相同,均为2.34 g·cm-3。而此时Cold区两端的密度仅为0.25 g·cm-3。由于局部高温使得Cold区两端的粒子向中间运动且形成相向的热-机械波,两波相遇后开始反向运动,反向运动的热-机械波将使得Cold区域中间部分拉伸而降低密度。因此在4 ps时,热载荷为3000 K时中间区域的密度降至约1.70 g·cm-3,热载荷为4000 K时Cold区中间区域的密度则为约2.00 g·cm-3,热-机械波相向碰撞—反向运动引发粒子来回震荡,如此往复,直到机械能完全转化为热能,并且通过Cold区中间区域密度的变化可以看出此过程所需的时间约为20 ps。两端的温度差异致使Cold区域内未分解的含能分子整体向右偏移,而使得密度逐渐出现不对称,并且随着这种效应的累积,至30 ps时,两端密度出现明显的不对称。相比于热在金属内传递的直接传递过程,压缩及膨胀体现了热在含能材料内传播的波动性。

a. 2～30 ps at 3000 K

b. 40～100 ps at 3000 K

c. 2～30 ps at 4000 K

d. 40～100 ps at 4000 K

图3不同热载荷下不同时刻Cold区内的密度剖面

Fig.3Density profiles of Cold region at different time and different hot region

3.3 粒子运动的速度剖面

图4为反应流传播过程中Cold区内粒子平动速度的剖面图。两端热载荷为3000 K时,在2 ps时刻,CL-20/TNT内沿x方向60 Å和140 Å处,粒子瞬时平动速率最大,沿x正方向分别为+0.4 km·s-1和-0.4 km·s-1。当热载荷为4000 K时,通过2 ps和8 ps时的粒子瞬时速度剖面可以看出,Cold区内两端粒子的压缩和膨胀的程度相比3000 K时的情况更为剧烈。同样在2 ps时刻,Cold区内约150 Å处粒子的速度达到-0.5 km·s-1。图4直观反映了局部温度热点向含能材料内传输热量的过程中,各剖面内的粒子首先往复运动,而后动能逐渐衰减至受热后的振动状态而表现出较小速度波动的过程(3000 K: 粒子由0.4 km·s-1衰减至振动模态; 4000 K: 粒子由约0.5 km·s-1衰减至振动模态)。Zhou在对HMX的研究中,局部热点温度为3000 K时,HMX内粒子前期往复震荡过程中,粒子的速度变化为0.5～0.1 km·s-1[17]。

3.4 高温诱发化学反应的产物识别分析

基于键级(Bond Order, BO)理论的产物识别方法为定量分析极端条件下含能材料的初始化学反应路径提供了支撑。反应流传播过程中,晶胞内的分子在高温作用下不断的分解形成各种分子碎片。这里,当任意两个分子碎片中的任意两个原子构成的原子对的键级≥0.3时,则认为新的产物生成[30]。

a. 3000 K

b. 4000 K

图4Cold区内粒子的瞬时速度剖面

Fig.4Particle velocity profiles of Cold region

图5为局部高温引发CL-20/TNT初始化学反应随时间的演化情况。其中图5a和图5b为两种热载荷条件下CL-20和TNT以及总物种数量随时间的变化。在两端施加热浴后,两端Hot区的CL-20和TNT随着温度的快速上升而迅速分解,并且温度上升至高温的时间尺度约为1 ps。3000 K热载荷条件下,1 ps后,体系内的CL-20和TNT由初始的320个分解至约200个; 而升至4000K时,CL-20和TNT分解至200个和180个,并且此过程中系统的新物种数量呈直线上升。随后,随着反应流的传播,CL-20和TNT逐渐分解。由于CL-20的感度高于TNT[5],因此,两种热载荷条件下,反应流传播过程中CL-20的分解速率一直高于TNT。当热载荷为3000 K时,CL-20完全分解的时间尺度约为92 ps,而至模拟结束时,体系内仍有约45个TNT分子未发生分解,此时,Cold区域内的最低温度约为1800 K。当热载荷温度为4000 K时,由于反应流传播速度明显高于3000 K时的情况,因此CL-20和TNT的完全分解的时间尺度显著变短,且分别为65 ps和70 ps。对比图5a与图5b可以看出当CL-20和TNT分解完毕时,系统内的新物种数量均达到峰值,并且两种温度载荷下,新物种的峰值数量几乎相同(约380种)。当热载荷为4000 K时,物种的数量达到峰值后,随着系统内的平均温度进一步升高,初级分解反应形成的小分子碎片逐渐聚合,进而形成最终产物而导致物种的数量开始衰减。

a. evolution curves of CL-20 and TNT molecule and total number of species at 3000 K

b. evolution curves of CL-20 and TNT molecule and total number of species at 4000 K

c. evolution of NO2 and NO and final products at 3000 K

d. evolution of NO2 and NO and final products at 4000 K

e. evolution of intermediate at 3000 K

f. evolution of intermediate at 4000 K

图5高温诱发CL-20/TNT共晶初始化学反应的产物分布

Fig.5Initial chemical products of CL-20/TNT co-crystal decomposition induced by hot spots

CO2的形成主要来自于芳香环的断裂,产生的CHON和CH2O反应并形成CO,CO与O原子结合形成稳定的CO2分子。N2主要来源于(1)CL-20和TNT上—NO2键断裂后,发生NO2+H→NO+OH,而NO上的O和OH或H原子结合进一步形成H2O,或和CO结合形成CO2,而N,N原子结合形成N2; (2)CL-20笼状主环断裂形成部分N—NO2基团,随后两个基团结合形成N2和NO2,而NO2进一步参与到NO2+H→NO+OH反应中。对比图5c与图5d可看出，局部高温诱发CL-20和TNT共晶化学反应最终的产物为N2,CO2和H2O,这与CL-20/TNT共晶高温热分解得到的结果一致[29]。

4 结 论

对CL-20/TNT晶体两端施加高温热载荷,模拟了局部高温在含能材料内引发反应流成长和传播的时空响应,并定量的揭示了初始化学反应情况。模拟结果表明:

(1) 相比于高温在金属内传递,此过程主要分为两个阶段: ①热-机械波传播使得Cold区内两端粒子向中心往复运动而出现密度梯度,并且前期粒子运动速度为百米量级。随着热-波机械波的反复碰撞、衰减,粒子动能逐渐转化为热能,而最终处于振动状态而呈现较小的速度波动。②含能分子在热传递的过程中而逐渐发生分解反应,并且上述两个过程均紧密的依赖于两端的热载荷大小。

(2) 温度的差异作用表现在反应流的传播速度方面,这将直接影响粒子往复运动的速度以及CL-20和TNT的分解速率,另外NO2和N2,H2O,CO2分别是局部高温诱发共晶炸药热分解的初始产物和最终产物。

含能材料的爆轰过程涉及复杂的物理、化学以及流体动力学问题。当前已经出现利用ReaxFF反应力场模拟极端条件下数百万原子体系内化学反应的研究成果。未来应当在反应模型、计算能力、势函数、化学反应动力学和力学相结合的交叉学科理论等方面共同努力,详细的了解含能材料能量释放的时空分布及其与外界条件,材料参数,力学运动耦合关系的细节,并期待从原子水平上揭示极端条件下真实含能材料热解、燃烧以及爆轰的本质。

参考文献:

[1] Thomas J M. Crystal engineering: origins,early adventures and some current trends[J].CrystEngComm, 2011, 13 (13): 4304-4306.

[2] Bolton O,Matzger A J. Improved stability and smart-material functionality realized in an energetic cocrystal[J].AngewChemIntEd, 2011, 50(38): 8960-8963.

[3] 杨宗伟, 张艳丽, 李洪珍, 等. CL-20/TNT 共晶炸药的制备、结构与性能[J]. 含能材料, 2012, 20(6): 674-679.

YANG Zong-wei, ZHANG Yan-li, LI Hong-zhen, et al. Preparation,structure and properties of CL-20/TNT cocrystal[J].ChineseJournalofEnergeticMaterials(HannengCailiao), 2012, 20(6): 674-679.

[4] Bolton O, Simke L R, Pagoria P F, et al. High power explosive with good sensitivity: A 2:1 cocrystal of CL-20: HMX[J].CrystalGrowth&Design, 2012, 12(9): 4311-4314.

[5] 孙婷, 刘强, 肖继军, 等. CL-20/HMX共晶及其为基PBX界面作用和力学性能的MD模拟研究[J]. 化学学报, 2014, 72 (9): 1036-1042.

SUN Ting, LIU Qiang, XIAO Ji-jun, et al. Molecular dynamics simulation of interface interactions and mechanical properties of CL-20/HMX cocrystal and its based PBXs[J].ActaChimSinica, 2014, 72(9): 1036-1042.

[6] YANG Z W, LI H Z, ZHOU X Q, et al. Characterization and properties of a novel energetic-energetic cocrystal explosive composed of HNIW and BTF[J].CrystalGrowth&Design, 2012, 12(11): 5155-5158.

[7] WEI C X, HUANG H, DUAN X H, et al. Structures and properties prediction of HMX/TATB co-crystal[J].Propellants,Explosives,Pyrotechnics, 2011, 36(5): 416-423.

[8] SHEN J P, DUAN X H, LUO Q P, et al. Preparation and characterization of a novel cocrystal explosive[J].CrystalGrowth&Design, 2011, 11(5): 1759-1765.

[9] 卫春雪, 段晓惠, 刘成建, 等. 环四甲撑四硝胺/1,3,5-三氨基-2,4,6-三硝基苯共晶炸药的分子模拟研究[J]. 化学学报, 2009, 67(24): 2822-2826.

WEI Chun-xue, DUAN Xiao-hui, LIU Cheng-jian, et al. Molecular simulation on co-crystal structure of HMX/TATB[J].ActaChimSinica, 2009, 67(24): 2822-2826.

[10] ZHOU T T, HUANG F L. Effects of defects on thermal decomposition of HMX via ReaxFF molecular dynamics simulations[J].TheJournalofPhysicalChemistryB, 2010, 115: 278-287.

[11] Rom N, Hirshberg B, Zeiri Y, et al. First-principles-based reaction kinetics for decomposition of hot, dense liquid TNT from ReaxFF multiscale reactive dynamics simulations[J].TheJournalofPhysicalChemistryC, 2013, 117(41): 21043-21054.

[12] Furman D, Kosloff R, Dubnikova F, et al. Decomposition of condensed phase energetic materials: interplay between uni-and bimolecular mechanisms[J].JAmChemSoc, 2014, 136(11): 4192-4200.

[13] Levesque G, Vitello P, Howard W. Hot spot contributions in shocked high explosives from mesoscale ignition[J].JApplPhys, 2013, 113: 233513.

[14] XIA R H, TIAN X G, Shen Y P. MD simulation of a copper rod under thermal shock[J].ActaMechSin, 2010, 26: 599-604.

[15] FU B Q,LAI W S,YUAN Y,et al. Molecular dynamics study of thermal stress and heat propagation in tungsten under thermal shock[J].ChinPhysB, 2013, 22(12): 126601.

[16] XIA R H, TIAN X G, SHEN Y P, et al. Molecular dynamics simulation of thermoelastic coupling characteristic for low-dimensional nano-rod under thermal shock[J].AdvancedMaterialsResearch, 2012, 581-582: 444-447.

[17] ZHOU T T, SONG H J, LIU Y, et al. Shock initiated thermal and chemical responses of HMX crystal from ReaxFF molecular dynamics simulation[J].PhysicalChemistryChemicalPhysics, 2014, 16(27): 13914-13931.

[18] Yanilkin A V, Sergeev O V. Molecular dynamics simulation of combustion front propagation in a PETN single crystal[J].Combustion,ExplosionandShockWaves, 2014, 50(3): 323-332.

[19] Nomura K I, Kalia R K, Nakano A, et al. Dynamic transition in the structure of an energetic crystal during chemical reactions at shock front prior to detonation[J].PhysicalReviewLetters, 2007, 99(14): 148303.

[20] ZHANG L Z, Zybin S V, Duin A C T, et al. Modeling high rate impact sensitivity of perfect RDX and HMX crystals by ReaxFF reactive dynamics[J].JournalofEnergeticMaterials, 2010, 28: 92-127.

[21] WEN Y S, XUE X G, ZHOU X Q, et al. Twin Induced sensitivity enhancement of HMX versus shock: a molecular reactive force field simulation[J].TheJournalofPhysicalChemistryC, 2013, 117(46): 24368-24374.

[22] 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].JPhysChemA, 2014, 118(5): 885-895.

[23] LIU L C, LIU Y, Zybin S V, et al. ReaxFF-lg: Correction of the ReaxFF reactive force field for london dispersion,with applications to the equations of state for energetic materials[J].TheJournalofPhysicalChemistryA, 2011, 115(40): 11016-11022.

[24] 宋华杰, 周婷婷, 黄风雷, 等. 低压长脉冲载荷下β-HMX单晶滑移系的微观物理化学响应[J]. 物理化学学报, 2014, 30(11): 2024-2034.

SONG Hua-jie, ZHOU Ting-ting, HUANG Feng-lei, et al. Microscopic physical and chemical responses of slip systems in the β-HMX single crystal under low pressure and long pulse loading[J].ActaPhysChimSin, 2014, 30(11): 2024-2034.

[25] ZHOU T T, Zybin S V, LIU Y, et al. Anisotropic shock sensitivity forβ-octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine energetic material under compressive-shear loading from ReaxFF-lg reactive dynamics simulations[J].JournalofAppliedPhysics, 2012, 111(12): 124904.

[26] http://www.ccdc.cam.ac.uk (accessed Feb 24, 2013).

[27] Budzien J, Thompson A P, Zybin S V. Reactive molecular dynamics simulations of shock through a single crystal of pentaerythritol tetranitrate[J].TheJournalofPhysicalChemistryB, 2009, 113(40): 13142-13151.

[28] Strachan A, Kober E M, Van Duin A C T, et al. Thermal decomposition of RDX from reactive molecular dynamics[J].TheJournalofChemistryPhysicals, 2005, 122(5): 054502.

[29] 刘海, 李启楷, 何远航. CL-20-TNT 共晶高温热解的ReaxFF/lg反应力场分子动力学模拟[J]. 物理学报, 2013, 62(20): 208202.

LIU Hai, LI Qi-kai, HE Yuang-hang. Pyrolysis of CL-20/TNT cocrystal from ReaxFF/lg reactive molecular dynamics simulations[J].ActaPhysSin, 2013, 62(20): 208202.

[30] Van Duin A C T, Dasgupta S, Lorant F, et al. ReaxFF: a reactive force field for hydrocarbons[J].TheJournalofPhysicalChemistryA, 2001, 105(41): 9396-9409.

[31] Isayev O, Gorb L, Qasim M, et al. Ab Initio molecular dynamics study on the initial chemical events in nitramines:thermal decomposition of CL-20[J].TheJournalofPhysicalChemistryB, 2008, 112: 11005-11013.

[32] Cohen R, Zeiri Y, Wurzberg E, et al. Mechanism of thermal unimolecular decomposition of TNT (2,4,6-trinitrotoluene): A DFT Study[J].JPhysChemA, 2007, 111: 11074-11083.

ReactiveFlowPropagationinCL-20/TNTCo-crystalExplosiveInducedbyLocalHighTemperatureZones: