(1.淮阴工学院 化学工程学院,江苏 淮安 223003; 2.湖南云箭集团有限公司 制导武器装备研究院,湖南 长沙 410100; 3.南京理工大学 化学与化工学院,江苏 南京 210094)
量子化学; 含能材料; ReaxFF力场优化; 分子模拟; 理论计算
(1.Faculty of Chemical Engineering, Huaiyin Institute of Technology, Huaian Jiangsu 223003, China; 2.Research Institute of Guidance Weapon Equipment, Hunan Vanguard Group Co. Ltd., Changsha 410100, China; 3.School of Chemistry and Chemical Engineering, Nanjing University of Science and Technology, Nanjing 210094, China)
quantum chemistry; energetic materials; ReaxFF force field optimization; molecular simulation; theoretical calculation
DOI: 10.14077/j.issn.1007-7812.202211009
在高温、高压等极端条件下原位追踪含能材料的超快化学反应极其困难,通过组建炸药体系的特定构型模拟其热解、氧化、爆轰等过程的物理化学特征,预测各种可观测的反应行为,进而获取关键的调控方法尤为重要[1-4]。计算机模拟方法加之合理的分子结构模型和物理原理,以研究含能材料在单一和多因素耦合条件下复杂和快速的微观反应是重要的手段[5]。目前,计算模拟在理解和预测炸药分子/晶体构效关系、反应动力学、微观反应机制等发挥着关键性作用,对解决炸药领域中复杂的科学问题和揭示一般性规律具有指导意义[6-7]。
分子动力学(MD)是探索微观结构和宏观性质的有力工具,它通过一系列近似处理可达到模拟较长时间尺度下的大型系统。经典力场本质上是一组参数化数学方程,旨在捕捉原子相互作用的势函数,例如键、价角、扭转、范德华力和库仑相互作用[8]。这种简化处理极大减低了总体计算成本。然而,经典力场(非反应力场)基于固定键/电荷解决方案,缺乏对化学反应中断键和成键的描述。基于量子力学的分子动力学虽然提供了准确的能量、电荷和反应路径,但由于其计算成本高,在模拟时间和空间尺度方面受到限制,妨碍了对影响材料整体行为的关键因素的理论理解。为了弥合量子力学与经验势方法之间的差距,已开发的反应力场(例如,ReaxFF[9]、COMB[10]、REBO[11]、Tersoff[12])允许化学键在分子动力学模拟中形成与断裂,并且原子可动态地被分配部分电荷,如电负性均衡方法(Electron Equilibration Method,EEM)[13]。这种动态键和动态电荷也使得反应力场的势能表达式比经典的复杂得多。
当前,基于第一性原理的ReaxFF反应力场应用最广泛,能相对可靠地评估原子间配位变化(键级和距离的关系)进而描述化学反应(即键断裂和形成)。ReaxFF反应力场的动力学模拟尺度可达到微米等级,并在多尺度建模中发挥核心作用[14],以相对较低的计算成本描述复杂的多步化学反应、机械和热过程。ReaxFF反应力场是由许多函数形式的经验力场参数组成,主要基于第一性原理的训练集进行优化[9]。ReaxFF-MD模拟在原子水平上描述了化学反应过程中的完整动态,与量子力学计算相比,计算成本显著降低。本文综述了ReaxFF反应动力学在含能材料模拟研究中的方法学和应用领域的进展。
ReaxFF采用结合可极化电荷描述的键级形式来描述原子之间的反应性和非反应性相互作用,键级为使用经验公式直接从原子间距离计算得出。ReaxFF各能量项的总结和其力场训练过程如图1所示。此外,为了弥补ReaxFF力场在有机晶体中对长程vdW相互作用低估的缺点,Liu等[15]引入基于低梯度(Low-Gradient,lg)模型的修正项对色散作用进行修正。
图1 ReaxFF的总势能和力场训练过程
Fig.1 The total potential energy and force field training of ReaxFF
ReaxFF反应力场虽具有强大的化学动力学仿真能力,但毕竟是量子化学方法的一种简化和近似。为了最大程度地精确再现量子力学相互作用下的化学反应过程,ReaxFF反应力场以经典相互作用对其近似只能是特化的,因而具有很差的迁移性。对于不同类型的体系,即使其中某些体系具有相同的物质,也不得不开发出多套力场参数分别使用。例如,有机炸药和生物体系均有大量有机官能团和水、二氧化碳等基础物质,却不能共用一套ReaxFF反应力场参数,因为两种体系分别需要对热分解反应和肽水解反应进行特化,以使精度达到最低要求。而对于研究者而言,一旦获得了目标体系的一套优良的ReaxFF反应力场参数,则代表着其化学动力学研究能力几乎立即从单分子层面达到了介观层次,足以针对界面、孔隙、缺陷、边界、晶格等微观结构效应进行还原,大量的实验现象得以重现,新的机制得以挖掘。
由于对各种特定体系开发ReaxFF反应力场参数的巨大需求,力场参数训练一直都是研究的热点领域。与经典非反应力场不同,ReaxFF反应力场包含了大量的两体、多体相互作用势,力场参数集庞大,训练集包含的化学响应种类多。换言之,训练一套ReaxFF反应力场参数,是一个N维(N=103数量级)M目标(M=102数量级)的最优化问题,因而其训练具有很大的难度。数十年来,一直有新的、更有效的ReaxFF参数优化方法被开发出来,称之为“力场开发的开发(Development of the development of forcefield)”。其发展如图2所示,最初创立的方法为序贯抛物线参数插值法(sequential parabolic parameter interpolation method,SOPPI)[16]。
图2 ReaxFF反应力场优化方法的发展
Fig.2 Development of ReaxFF reaction force field Parameter optimization method
SOPPI对训练集参数需逐一执行单参数搜索,直到满足指定收敛标准为止。SOPPI方法虽简单,但有以下几个缺陷:(1)容易陷入局部最小值;(2)参数量增多时,会使单参数优化量急剧增加;(3)收敛成功与否和初始参数猜想以及各参数优化顺序密切相关。针对SOPPI的缺陷,全局优化方法被使用并展示出高效、自动化的优势,如遗传算法(GA)[17-18]、模拟退火法(SA)[19-20]、进化算法(EA)[21]、粒子群优化(PSO)[22]和机器学习(ML)[23]。全局优化可作为初步优化,在优化接近收敛时,可结合简单的SOPPI方法。
目前,大多数ReaxFF参数是通过执行单一罚函数(penalty function)优化的,罚函数又是通过在指定参数的训练集中评估不同量获得误差的加权组成。在权重的分配上难免需要人为干预来调整优化过程,尤其是涉及到不同维度下组合量的权重,例如,键长、电荷和构型能等。Larentzos和Rice等[24-26]首先引入基于帕累托最优(Pareto dominance)的多目标方法,无需对ReaxFF各参量指定任何权重。在运算中,该方法产生一系列解决方案,原则上帕累托边界可以通过使用随机变化的权重进行重构并使其不断进化,对应于各种任意权重的所有可能值。
多目标进化策略(MOES)框架是一种无需指定目标的后验方法,结合了进化策略(ES)和数据包络分析方法(data envelopment analysis, DEA),被用于训练ReaxFF力场参数[27]。在MOES中,基于帕累托最优提供了一种用于平衡各种不相称量(键长、键角和二面角等)的自动运算策略,从而消除了人为设定权重对优化的影响。力场训练过程包括:(1)准备训练集和初始ReaxFF参数。其中,训练集包含足够多的从头算数据和实验值,如电荷、几何构型(键长、键角和二面角)、生成热、振动频率和能量等,从文献中获取现有力场参数作为初始ReaxFF参数;(2)将初始参数优化的能量和结构特性与训练集进行比对并产生误差,随后输入到MOES中;(3)MOES对误差进行评判并重新创建新一代的力场参数,有相称维度的类似参考数据点被分组在一起,目标函数(误差)的计算公式如下:
式中:j表示目标。
反复执行上述过程直至达到收敛标准。获得的各项ReaxFF收敛参数需经基于帕累托分析的动力学处理得到最优参数。
运用MOES方法对RDX、FOX-7、HMX、TATB、CL-20和PETN的力场参数进行拟合,并获得的最优晶胞参数。优化后的力场可较好地描述这几种炸药的晶格参数、晶体内分子取向和分子结构。这些炸药晶格参数最大误差的绝对值依次为2.28%、3.89%、1.60%、10.30%、5.15%和3.90%。由于TATB中存在强氢键,在长程相互作用下晶格参数的偏移误差导致TATB的密度误差增大,然而ReaxFF-lg模型中的分子更平坦,许多氢键与分子平面呈90°夹角。综上,MOES方法的优势在于:相较于传统的单一目标参数搜索,依托帕累托最优方法分析产生无需人为指定权重的多个有效力场方案,通过对这些有效力场的评估得到最优力场。因此,MOES方法适合用于力场的开发。Larentzos和Rice[28]还利用这种可迁移的拟合力场策略将ReaxFF-lg力场扩展到硝基甲烷的研究。
粒子群优化(particle swarm optimization, PSO)算法也被用于训练ReaxFF力场,是一种用于处理随机群体的演化算法,可避免陷入最小化可用于优化非线性、不可微函数的全局搜索方法。PSO的优化工作开始于多组猜测(或随机生成)的力场参数,将每组待优化参数视为解空间中的一个点粒子,而每次迭代产生的参数变化量数组则被视为该粒子的移动速度,功能主要体现为以下3点:(1)同时接受优等和劣等解决方案,在提高随机搜索效率的同时也避免过早收敛的问题;(2)允许同时存在大量潜在的解决方案,从而提高了在盲搜过程中找到目标的概率;(3)与模拟退火不同,PSO在各搜索个体之间共享信息。搜索个体的位置取决于先前行为的惯性项和个体自身的经验参考。个体的速度和位置表示为:
式中:vik和xik分别为个体i在迭代k时在解空间中的位置和移动速度; pik为个体i的最佳位置; pgk为群体的全局最佳位置; 系数c1和c2为个体和全局项的恒定权重; ω为惯性因子,被用于抑制个体的动量,以防止在优化过程中速度无限增长; 参数ri1k和ri2k表示在[0,1]范围内的均匀分布采样的独立随机标量。
在大多数情况下,可分离函数在旋转下变得不可分离。使用随机标量(而不是向量)使得算法对搜索空间的旋转保持不变(rotation-invariant PSO,RiPSO)[29]。该算法将在可分离和不可分离函数中提供一致的性能。然而,在处理复杂的分子力场时RiPSO算法明显性能不佳[30]。这是由于RiPSO算法缺乏多样性,进而导致群体各成员轨迹塌陷到低维子空间中,即线搜索。Furman等[22]对标准RiPSO算法进行直接的扩展,各向同性的高斯变异(GM)算子被用于增强其全局搜索能力,防止其退化为线搜索,因而被命名为具有高斯突变的旋转不变粒子群优化(RiPSOGM)。变异算子的作用是用群体中最佳个体的变异版本替换当前个体,见下式:
式中:个体的新位置xik+1被设置为以最佳粒子位置pgk为中心的各向同性随机值。步长是从具有零均值和统一的标准偏差的高斯分布G(0,1)中采样的。尺度参数γ与域边界相关,并控制分布的宽度。单位向量t^kk+1通过在d维超球面的表面上均匀拾取点来确保位移的分布是各向同性的。
使用RiPSOGM对色散作用力场参数进行优化,通过DFT-TS计算了HMX、TATB、PETN、FOX-7和NM的状态方程。参数训练中的目标函数是以加权误差平方和的形式表示:
式中:xi,DFT为参考值; xi,ReaxFF为计算值; σi为为训练集中的每个优化目标指定的权重。
Furman等[15]在使用RiPSOGM优化后获得的ReaxFF-lg预测的升华热,并与Liu等拟合的ReaxFF-lg力场作对比。通过对PETN、NM和TATB的测试,发现使用RiPSOGM优化的ReaxFF-lg力场总体上达到了良好效果。对于反应热和活化能等关键反应性质,以TNT、RDX、HMX和CL-20四种典型炸药为例,ReaxFF虽然不能给出准确值,但能获得正确的定性排序,如图3所示。
图3 实验和ReaxFF得到的典型炸药反应热(Er)和活化能(Ea)
Fig.3 Reaction heat(Er)and activation energy(Ea)of typical explosives obtained by experimental and ReaxFF methods
基于神经网络的反应力场是依托了大量的ReaxFF-MD模拟数据,用以拟合原子间势能的方法[31]。这些势能模型与基础物理量无关,将局部原子环境与能量/力建立关系。Yoo等[32]提出了一种数据驱动的力场:神经网络反应力场(NNRF),并较好的描述硝胺类炸药(eg. RDX、HMX、NM)的复杂化学反应。训练集是由密度泛函中的Perdew-Burke-Ernzerhof和处理London色散作用的D2校正计算得到的大量电子结构。双向不同的反馈回路(凝聚相和单分子系统)使得训练集更加丰富。NNRF方法使用二体和三体对称函数与结合局部原子环境的映射。此外,还通过预测测试集中的解离曲线、能量和作用力来验证。
NNRF的输入层是以原子为中心的加权高斯对称函数,公式(6)和(7)中分别为从分子几何结构中提取的径向(wGradi)和角(wGangi)两种高斯对称函数。
式中:Zj是用于区分不同元素的原子量; rs、η、λ和ξ是要指定的超参数; fc(r)是截止函数,其在5Å的截断半径之外平滑衰减为零。
独立的神经网络被用来描述每个元素的能量和力。每个神经网络由一个输入层(42个节点)、两个隐藏层(各50个节点)和一个输出节点(原子能量)组成。其中,将数据随机分为用于优化神经网络的训练集(包含90%的数据)和用于评估收敛性并避免训练期间过度拟合的验证集(10%)。在较宽温度范围内(1500～2500K),NNRF方法被应用到RDX的高温分解的复杂化学性质中,并得到与实验相近的分解产物和活化能。RDX分解会产生大量中间体和产物分子,如CH2O、HCO2、HON、HONO、HOCN、HCN、CO、CO2、NO2和N2等。通过对比PBE-D2、ReaxFF-2014、ReaxFF-2018和ReaxFF-lg以及实验参数计算的形成能[32],发现NNRF的形成能较好地吻合实验值。
NNRF利用神经网络捕获具有径向和角度描述符的短程相互作用(<5Å),并采用ReaxFFs的van der Waals和Coulomb相互作用参数化表达远程相互作用(<10Å)。使用NNRF的迭代训练研究复杂化学具有DFT级别的准确性,并可转移至其他硝胺炸药(HMX、NM、TNT、TATB、CL-20和PETN),为含能材料的大规模热分解和冲击动力学提供了新的解决方案。然而,由于NNRF力场针对体系的几何结构给出力和能量数值的过程是“黑盒”式的,缺乏明确的相互作用形式表达式,难以明确判定化学反应或进行统计,且体系的某些热力学参数还无法从模拟中获得(如熵和自由能),有待进一步解决和发展。
ReaxFF模拟已经广泛应用于探索含能材料的热解和燃烧机制中。研究主要涉及反应物初始分解、中间体/自由基转化和消耗、稳定产物形成的3个阶段[33-35]。其中,炸药初始机理包含N—NO2键断裂、NO2异构化、HONO消去反应、环裂解和二聚体等[36-41]。此外,也对凝聚相炸药中氢转移和碳团簇的形成过程进行了探究[42-43]。图4展示了关于含能材料典型的ReaxFF模拟的内容。
图4 含能材料的ReaxFF模拟
Fig.4 ReaxFF simulations of energetic materials
炸药燃烧传播过程经历从爆燃到爆轰的转变,爆炸物质在周围介质中产生高温高压的化学反应。准确预测和理解炸药在热和压力作用下的基本反应至关重要。ReaxFF方法已经应用到炸药在热和压力耦合环境下的机理研究,其热力学参数见表1。值得注意的是目前加压的模拟主要是通过压缩超晶胞实现的,且炸药反应速率的一般分析方法是将炸药热分解过程进行分段,如初始分解反应(吸热阶段,S1)、次级化学反应(放热阶段,S2)和最终产物生成(S3)。
表1 含能材料ReaxFF模拟的热力学参数
Table 1 Thermodynamic parameters of ReaxFF simulation for energetic materials
单质炸药热解研究对象主要包括TNT、RDX、HMX和CL-20等凝聚态物质。Furman等[5]以TNT为例分别阐释了其在气相和凝聚相中不同的分解机制。计算得出TNT在气相和晶相中的分解活化能分别为约259.5kJ/mol和约146.5kJ/mol,显然晶相中TNT的活化能更低。借助ReaxFF-MD模拟观察到晶相中相邻TNT分子之间率先发生氢转移现象,进一步借助DFT计算证明双分子中氢转移是最低能量路径。Rom等[35]研究了致密的液体TNT热解行为。反应从初始TNT液相开始,先分解为自由基,随后形成碳簇产物,最后形成稳定的气态产物。在机理方面发现,单分子C—NO2键的断裂在较低密度下占主导地位,而二聚体的形成以及分解为TNT衍生物和较小的气体碎片在较高的压缩态占优势。更高的TNT密度则会加剧形成碳堆积,初始气态碎片形成受到抑制。
Sakano等[47]提出了ReaxFF与有限元结合的多尺度模型,并模拟RDX的热传输和化学分解过程。Wang等[34]通过反应动力学模拟和静态势垒计算研究了α-RDX和ε-RDX在高温高压下的分解,发现α-RDX的活化能随体系的不断加压而减小,分解率却随着压力的增加而提高。ε-RDX则表现与α-RDX相反的机制。Joshi等[58]研究RDX晶界处分子能量的转移,发现以弯曲模式为主的高频振动区域主要分布在界面热点处。反应域中的高温气相分子将大量热能转移到界面处的未反应的固体,激发了固体分子内的振动模式并触发局部化学反应,从而导致热点的增长。通过界面RDX分子的功率谱分析表明除了各种弯曲模式外,还激活了N—N、C—H和C—N键的拉伸模式。Sakano等[59]对比在快速均匀加热和热点存在下无定形和晶相RDX的反应性情况。晶相RDX在加热后经历快速吸热过程,晶体有序度降低并造成实际反应温度的降低。在直径为10nm的热点下,非晶相RDX的爆燃速度传播更快。
Wang[57]研究CL-20的热分解机理和特性。CL-20热分解的初始反应路径为N—NO2解离,随后是C—N断裂并导致笼状结构的破坏。此外,温度的升高会导致笼形结构被破坏的更早,并较大影响产物H2O和N2的生成速率。CL-20在起爆过程中的压缩程度对化学反应的影响也被探讨[36]。在高密度下N—NO2键断裂受到抑制,体系容易形成团簇并最终造成气体数量的减少。
Rom等[39]利用ReaxFF-MD研究了不同压缩条件下液体硝基甲烷的热分解机理。在低密度下,CH3NO2→CH3+NO2是主要反应途径。当密度和压力升高至Chapman-Jouget(CJ)爆轰条件(约30%压缩,>2500K)时,主要机制转向为氢转移和N—O键断裂,并形成CH3NO。Chen等[37]研究不同温度下HNS热解机制,揭示HNS主要初始分解机制是C—NO2键解离和硝基的异构化(NO2→ONO)。此外,还总结团簇的形成与温度之间的关系,低温会形成更大的团簇。Song等[38]对NONA的研究也发现相似的初始反应机理,这可能是二者同属于芳香族硝基胺类炸药,具有相似的结构。Zhou等[53]对硅-季戊四醇四硝酸酯(Si-PETN)的热解机理和其自身的高感度原因进行研究,发现低温下Si-PETN分子发生了独特的重排现象(R3Si—CH2—O—R2→R3Si—O—CH2—R2)并导致势垒的显著降低。相比而言,PETN的分解没有此过程。放热阶段,Si—PETN反应速率比PETN高约3倍,活化能降低约29.3kJ/mol。此外,Si—O键的形成释放巨大的能量,进一步推动反应的发生。
混合炸药的ReaxFF模拟包含高分子黏结炸药、液体混合炸药等。Russo等[60]使用ReaxFF模拟了不同组分的DNB和硝酸混合物的燃烧动力学。在某些组成的混合物下,热能会急剧释放并表现出自燃的特征。对关键反应机制的分析发现,在低温下单独模拟这些组分时并没有发生反应。然而一旦这些组分以1:1的质量比混合在一起时,观察到DNB燃烧速率更快。Joshi等[61]研究了HTPB/RDX复合材料的热化学响应。发现HTPB通过延长从点火到爆燃的时间来提供一定的安全缓冲,这种延迟主要是由于HTPB和RDX的热容量的差异。RDX的爆燃进一步导致HTPB的热解。HTPB的热解产物主要包括超支化聚合物、环状分子如环己烯、环己二烯和环戊二烯,以及线性分子如丁二烯和反式丁二烯低聚物。Yan等[62]构建了CL-20和聚合物黏结炸药(PBXs)的热分解模型并预测匀速分解的温度分布。其中,聚合物黏结剂包含聚异丁烯(PIB)、丙烯腈-丁二烯橡胶(NBR)、丁苯橡胶(SBR)、A型氟橡胶和氟黏合剂等。结果表明,在聚合物的作用下,—NO2基团释放O的频率减少,使CL-20更加稳定。与氟聚合物相比,烃类聚合物可以使ε-CL-20在笼状结构破坏前发生更完全的N—NO2断裂,这可用于解释实验中烃类聚合物可以大大缓解ε-CL-20的分解过程,而含氟聚合物对此影响很小的原因。
共晶炸药是协调各现有分子间性能并构建特定功能的结晶策略,被认为是改善现有含能材料性质,并获得低冲击感度或高热稳定性的先进含能材料的方法。共晶炸药的ReaxFF反应动力学模拟以含有CL-20成分的共晶为主。Xue等[55]模拟了CL-20/HMX共晶的热稳定性,发现初始反应占据主导的为N—NO2键的断裂。CL-20分子的热分解速率更快,并加速HMX的衰减。Guo等[56]将CL-20/TNT共晶与TNT和CL-20单晶的热解行为进行比较。结果表明,CL-20/TNT共晶释放能量的速率介于CL-20和TNT单晶之间。相比于CL-20,共晶热解具有更高的势垒,且共晶中形成的碳团簇抑制了活性中间体,阻碍参与二次反应。此外,还发现CL-20/TNT共晶比二者的不定形态混合物释放能量的速率更慢。Ren等[33,63]借助VARxMD捕获CL-20/HMX和CL-20/TNT共晶热解过程中的中间体和化学反应路径,同样也揭示了TNT组分热解产生的中间体和富碳团簇延缓了CL-20/TNT的热解并降低共晶的感度。Sultan等[64]研究了CL-20/DNT共晶体在不同极端温度下的热分解机理,包括热稳定性、热点演变、反应区形成和反应动力学等,发现CL-20/DNT共晶在感度、势能、物种总数、碳团簇数和爆炸能量方面优于纯CL-20。Xiao等[57]研究了CL-20/H2O2溶剂化物中主客体的热分解机制。H2O2的分解产物通过两条重要途径参与了CL-20的分解反应。H2O2分解产生的H和O原子以及OH连接到CL-20的硝基中的O原子上,随后导致N—O键断裂。其次,H2O2分解产生的O原子与C原子相连,导致C—N键断裂,最终破坏了CL-20的笼状结构。
由于尺度效应,纳米炸药在反应动力学和内在反应机制上与块体炸药有所差异。ReaxFF-MD模拟对理解纳米炸药的复杂燃烧机制和表面诱导效应等方面非常重要。Zheng等[65]对比纳米和块状RDX分解过程的动力学和路径。结果表明RDX的热分解不是来源于纳米颗粒表面,而是在整个纳米颗粒内部随机发生。在较小尺寸的纳米颗粒中,从RDX中脱去NO2的反应概率更高。Zhong等[66]使用ReaxFF方法模拟RDX在9种气体环境下的衰减过程。其中,NH3、CO、NO和NO2可以促进RDX分子衰减,尤其是NH3。O2可以抑制衰减,而CO2、H2O、H2和N2不参与RDX的衰变反应,对反应物有稀释作用。She等[67]利用ReaxFF-MD研究了纳米FOX-7的热解。结果表明,温度为3000K时纳米FOX-7的势垒高于2000K时。由温度分布可知,纳米FOX-7分解首先发生在纳米颗粒的中心部位,Han等[68]在研究了粒径3～6nm的CL-20纳米颗粒的热分解过程时也有相似的热解规律。此外,还揭示了纳米FOX-7的3种初始反应途径:氢转移、C—NO2和C—NH2的断裂,其中C—NH2键的断裂是一种新的反应途径。Xiang等[69]研究纳米TATB的热解机制。高温破坏纳米TATB的类石墨烯层π-堆叠结构,并在热解反应之前触发团聚、滑移和分散。纳米TATB的初始分解包括:OH的形成、C—NO2键断裂、NO2→ONO和氢转移。
ReaxFF力场模拟含能材料在冲击加载方案主要有4种:(1)Hugoniostat;(2)多尺度冲击技术(multiscale shock technique,MSST);(3)CS-RD压缩剪切反应动力学;(4)高速活塞撞击。在冲击加载下炸药的ReaxFF模拟研究结果列于表2。
表2 在冲击加载下各种炸药的反应分子动力学模拟结果
Table 2 Reactive molecular dynamics simulation results of various explosives under shock
Hugoniostat分子动力学模拟被用于描述含能材料的冲击载荷。在冲击压力下,Hugoniostat单轴压缩体系并控制总能量满足Hugoniot跳跃条件。由于避开了通过材料的实际冲击波传播,因此计算强度明显较小并且可以使用适度的计算资源进行更长的时间尺度模拟,从而更好地统计化学分解反应。Hugoniostat施加的跳跃条件如下:
式中:ρ为密度; 下标0表示未冲击的初始性质; us和up分别为冲击和粒子速度; PZZ为应力张量在冲击方向上的法向分量。
对于给定期望的冲击压力,Hugoniostat技术使用恒温器和气压调节器来满足跳跃条件。Islam等[72]使用ReaxFF-MD模拟结合Hugoniostat技术研究硝基甲烷和聚乙烯醇硝酸酯(Polyvinyl Nitrate,PVN)对冲击载荷的热力学、化学和光谱响应。通过MD轨迹的后处理得到光谱图,见式(13),并通过时间相关原子速度的傅里叶变换获得[47],最终建立反应机制和光谱特征的关系并与激光驱动冲击实验进行比对。
式中:vj(t)为原子在t时间上的速度; β定义为(kBT)-1; 其中T为温度; τ=nΔt为总采样周期; mj为j原子的质量; Δt为采样率; N为要分析的帧数。
ReaxFF-MD模拟与MSST相结合是基于可压缩流动的Navier-Stokes方程[89]。MSST通过计算盒子中的原子数目和体积的运动方程来模拟稳定的冲击波,分别将冲击传播矢量应力和能量约束到瑞利线(Rayleigh line)和Hugoniot关系。根据穿过冲击波前沿的质量、动量和能量守恒可以得到Hugoniot关系:
E-E0=1/2(P+P0)(V0-V)[90-91](14)
式中:E为能量; P为冲击矢量中应力张量对角分量的负值; V为体积; 下标0表示预激状态,而没有下标的量表示后激状态。
将系统的初始状态连接到其最终Hugoniot状态的热力学路径可以用瑞利线P-P0=U2ρ0(1-ρ0/ρ)来描述(式中:U为冲击速度; ρ为密度)[92]。这两个关系描述了连续介质理论中稳定的平面冲击波,对于稳定的冲击波必须遵守。由于在含有大量氢原子的含能材料中的冲击诱导化学反应具有显著的核量子效应。量子浴多尺度冲击技术(QB-MSST)为冲击加热添加量子力学校正[93]。在应用中,Hamilton等[94]使用量子和经典恒温器研究TATB在冲击载荷下的分解。与经典模拟相比,核量子效应降低了与相关反应的活化势垒。Liu等[76]使用MSST方法对CL-20的冲击起爆进行了反应分子动力学模拟。结果表明,当C—N和C—H键的裂解率分别超过每皮秒3.11%和4.15%时,就会发生爆炸反应。更高的冲击速度会导致这些键的裂解率更高,但也会导致更多的原子形成团簇。不同冲击方向的分子排列存在显著差异,导致不同反应速率和反应路径,Li等[77]预测了TNAZ各方向的冲击感度排序为[100]>[010]>[001]。发现在[010]和[001]方向的分子排列有助于形成“缓冲”效应。MSST方法被应用于探究CL-20/TNT共晶结构和非晶结构在不同速度的冲击波下的响应过程[75]。在等速冲击波的作用下非晶结构比共晶结构更容易压缩,从而实现更高的应力和温度。无定形结构的化学反应更激烈,反应物衰变更快,产物更丰富,中间产品更早完成向稳定产品的转化。
每个冲击方向都会导致几种可能的滑移系。由于空间位阻,每组滑移系都可能导致剪切应力的显著差异,最终导致炸药的各向异性爆轰。最易激活的滑动系预计将具有更大的分切应力(RSS)、更低的临界分切应力和更小的Burgers矢量。CS-RD被设计为一种快速检测方法,用于模拟炸药化学反应与压缩方向和滑动系统的关系,可区分冲击方向的敏感性。Zhou等[83]使用CS-RD研究了δ-HMX的各向异性冲击感度。根据越过剪应力势垒过程中积累的内能,评价7个激波面的敏感程度:(111)/(011)/(110)/(101)>(100)/(010)>(001),其中(111)面最敏感,(001)面最稳定。这种各向异性的敏感度与滑移面上的不同分子间空间排列相关。当HMX剪切形变时,(111)面的分子几何松弛空间更小,而分子碰撞接触更多,进而导致更高的剪切应力势垒、更大的能量和温度增量以及更强的化学反应。Zybin等[81]发现PETN的机械和化学响应因滑动方向而异,这与冲击感度实验一致。特别是CS-RD的剪切应力过冲(由分子的不良空间相互作用引起)与观测的感度相关。这种过冲阻碍剪切形变并导致分子更大的变形,进而提高敏感方向的温度并增加化学键断裂,释放出诱导进一步分解的中间体,最终推动爆炸。相似地,An等[82]也使用CS-RD方法对RDX晶体的低密勒指数面进行单轴压缩,随后以0.5ps-1的恒定剪切速率使超胞变形。预测了(100)和(210)是RDX最敏感的冲击方向,而对(120)、(111)和(110)方向上的冲击是不敏感的。Guo等[84]对比了TNT和CL-20以及二者共晶的冲击敏感性。与纯CL-20晶体相比,CL-20/TNT共晶的冲击敏感性明显降低。在大气压下和高压下(约5GPa),CL-20/TNT比CL-20的剪切敏感度分别降低约70%和约46%。
活塞沿单方向从一端直接撞击超胞,在超胞中形成冲击波并传播到另一端。冲击波传播后,原子加速运动,超胞被压缩。使用活塞直接冲击炸药晶体更接近实际的冲击过程,可以更好地模拟具有各向异性结构的炸药在受到冲击载荷时的化学反应。Budzien等[85]使用大规模ReaxFF-MD模拟研究PETN的冲击诱导反应。激波沿PETN单晶的[100]方向被3或4km/s的活塞冲击,随后反射回来产生了一个稀疏波。强冲击系统在激波前沿后面显示出一个峰值,而弱冲击系统在模拟的后期显示为一个宽峰。在冲击波阵面的后方,除了观察到引发键O—NO2的断裂外,还检测到NO3、CH2NO3和CH2O的片段的形成。Wang等[86]研究CL-20在热和冲击下的化学反应和引发机理。在高温下,NO2是由N—NO2键的直接裂解产生的。相比之下,高压和冲击促进了O原子向与NO2相连的N原子的转移,导致了N—NO2键的断裂。Eason等[95]研究了PETN单晶中定向冲击诱导的非弹性变形机制。平面激波(pRankine-Hugoniot约9GPa)沿敏感性[001]和不敏感性[100]的晶向传播。对于敏感方向,只有弹性压缩发生,导致单个波在材料中传播,而对于不敏感方向,在冲击波阵面处和紧随其后的弹性压缩之后是非弹性变形,导致了一种双波结构。Yang等[87]模拟了在冲击载荷下层状晶体结构炸药ICM-102的各向异性反应机理及敏感度差异。当ICM-102的反应速率由慢变快时,均会观察到临界压力。ICM-102最敏感的方向为:x轴>y轴>z轴。此外,Yang等[88]还将ReaxFF-MD模拟应用于定量预测冲击感度方面。将反应后的炸药分子与初始超级单元中分子总数的比值定义为反应分数,并给出了不同载荷下炸药反应速率的计算方法。然后分析了反应速率和冲击波压力之间的关系,确定冲击下反应阈值的定量方法。最终得到冲击引起的常见硝酸盐炸药(如PETN)、硝酸铵炸药(如CL-20、HMX和RDX)以及硝基苯衍生物炸药(如TNT和TATB)点火和爆炸的临界压力。
ReaxFF力场在带有缺陷的炸药晶体中的应用也非常广泛,本研究内容主要包括点缺陷(空位)、线缺陷(刃型位错和螺型位错)、面缺陷(孪晶)和体缺陷(孔隙)。
(1)点缺陷。Zhou等[96]利用ReaxFF-MD模拟研究了不同温度下分子空位对凝聚相β-HMX分解机理和反应动力学的影响。结果表明,HMX的初始分解机制包括:N—NO2键离解、HONO消除和协同环裂变,这3种机制均受分子空位和温度的影响。分子空位处由于形成热点效应加速了HMX的分解。温度小于2000K时,空位质量分数为12.5%的HMX晶体的反应速率常数几乎是理想晶体的3～4倍,其活化能比理想晶体的活化能小24.3kJ/mol。
(2)线缺陷。Xue等[97]使用MSST方法对比完美的和位错的RDX晶体的冲击响应,还对带有4种位错的RDX进行剪切动力学模拟。发现位错均增强了RDX的冲击感度,尤其对于刃位错。依照剪切应力能垒划分了RDX晶体的迁移率:(010)[001]/screw(s2)>(010)[001]/edge(e2)>(010)1/2[100]/screw(s1)>(010)1/2[100]/edge(e1)。根据温度、压力和反应物衰减率,确认冲击感度遵循 e2>e1>s1≈s2>p的顺序。Deng等[98]研究了“冲击-预热-位错耦合”效应对RDX衰减的影响。结果表明,增加冲击速度和预热温度以及在晶体中出现边缘位错均促进RDX衰减。此外,还发现一些冲击、预热和边缘位错的耦合对RDX衰减具有等效的影响。
(3)面缺陷。Wen等[99]报道HMX孪晶的热解和冲击响应,并与HMX的完美晶体进行比较。结果表明孪晶显著提高了冲击感度,反应中孪晶界面处的温度最高,成为热点区域。在相同的冲击条件下,孪晶表现出较低的冲击起始应力和较高的分解速度。An等[100]将PETN和Si-PETN嵌入HTPB聚合物黏合剂基质中而构建具有锯齿界面的复合炸药模型,并用以研究在激波诱导下热点形成和爆炸行为。当冲击波由炸药传播到聚合物时,在凸形聚合物粗糙处形成热点。而由聚合物传播到炸药时,在凹形聚合物粗糙处先形成一个热点,随后在凸形聚合物粗糙处形成第二个更显著的热点。
(4)体缺陷。Furman等[101]研究了纳米尺度空洞的存在对ETN中的激波传播特性和化学分解的影响,描述孔隙塌缩和反应过程,包括冲击-空隙相互作用、空隙坍塌、反应性加强和ETN分子分解。在直径为5nm的空洞中观察到超音速分子纳米射流的形成,这导致了局部热效应和加快的化学分解速率。纳米空洞塌陷形状是由球形到新月形的形状转变,与瑞利型流体动力气泡坍塌非常吻合。Wang等[102]研究冲击加载下含有纳米空隙的CL-20晶体的分解机理。获得了CL-20分子在空隙周围的初始反应,讨论了冲击速度和空隙尺寸对CL-20初始化学反应的影响。Zhou等[103]以1km/s的速度冲击半径为4nm空隙的β-HMX单晶。冲击强度和空隙尺寸对空隙塌陷和热点形成有显著影响。冲击由弱到强的过程中(Up:0.5→3km/s),在空隙处经历了坍塌-形成流动分子-下方未冲击表面撞击的一系列过程,这种变化完成度受冲击强度的影响。Wood等[104]研究带有圆柱形孔洞的RDX的冲击过程引发的化学反应。发现直径为40nm的孔洞的塌缩会导致爆燃波。坍缩过程中的分子碰撞导致发生的超快和多步化学反应。Zhang等[105]模拟带有立方体空隙的RDX晶体的冲击过程。发现冲击引起的空洞坍塌、热点的形成和生长以及衰落均取决于冲击速度。此外,还揭示了初始热点的形成步骤:一是局部塑性变形引起的温度升高; 二是空隙塌陷过程中上方和下方颗粒碰撞引起的温度升高。
飞秒激光(fs-laser)脉冲在短时间内产生超高能量,此技术具有安全性高、精度高和无污染等优点,且可引发炸药的化学反应。Wu等[106]使用ReaxFF力场研究了HMX在不同飞秒激光能量下的烧蚀机制。结果表明,HMX的飞秒激光烧蚀机制与激光功率密度有关。当激光功率密度足够高时(3.4×1014W/cm2, 1.0mJ/pulse1),HMX在皮秒级(约7.65ps)发生电离或分解反应,并产生高温高压等离子体。随着激光能量的降低,烧蚀产物的电离度降低,其中单原子和离子产物的数量减少,而小分子(如NO、HNO2和HNO)的数量增加。Zhang等[107]研究了电场诱导下CL-20/HMX共晶的分解机制。结果表明,外部电场在加热晶体、拉伸化学键从而导致炸药分解方面起着重要作用。强恒定电场极大地提高了CL-20/HMX的感度,这归因于在电场下CL-20分子的开环反应占据主导。Zhang等[108]还研究了电场诱导CL-20/H2O2和CL-20/N2O的分解。与α-CL-20和ε-CL-20进行了比较,发现在电场下的感度顺序为:α-CL-20/H2O2>ε-CL-20>α-CL-20>α-CL-20/N2O。
铝基炸药燃烧氧化和界面钝化的ReaxFF-MD模拟研究对理解、预测和调控含能材料性能至关重要。Wang等[109]在原始C/H/N/O的ReaxFF参数上引入Al元素,并模拟RDX和Al板的相互作用行为。RDX可在Al(111)面自发反应; 当在(111)面添加不定性氧化铝时,仅稍微增加RDX分子分解的势垒,而RDX分子在α-Al2O3的(001)面的热解却被显著地延缓。在冲击载荷下Al和RDX/Al的动态形变也被研究[110]。当ANP冲击RDX晶面时,RDX与Al的相互作用不影响RDX晶体的各向异性响应顺序,不同厚度的ANP对结构变形的抵抗能力不同。当用3nm厚度的氧化铝涂覆在Al板时,氧化铝层阻碍了RDX分子与活性Al原子的反应。Mei等[111]使用具有低梯度校正的反应力场研究了典型炸药(TNT、RDX、HMX和CL-20)与纳米铝颗粒的反应行为。研究表明,含铝炸药经历了3个热解阶段:吸附期(0～20ps)、扩散期(20～80ps)和形成期(80～210ps),这些阶段依次对应铝与周围炸药分子之间的化学吸附(R-NO2-Al形成)、炸药的分解和O原子扩散到Al纳米颗粒中以及最终产物的形成。Zhong等[112]研究了Al@Al2O3纳米粒子在RDX中开裂行为。发现Al核的原子向外迁移并与RDX反应,从而使氧化铝壳层变厚。增厚的壳层一定程度上保护了铝芯原子,但随着增强的内部应力致使其发生开裂。Hao等[113]研究了纳米铝在TATB中的演变。结果表明,铝添加剂增加整个体系热量释放速度,纳米铝的形态演变取决于其大小和含量,更小的尺寸有助于引发微爆炸。在加热过程中,Al@Al2O3逐渐形成空心铝球。氧原子向铝芯移动,导致最终铝芯的氧原子浓度更大。CL-20和纳米铝在冲击加载下的快速反应机制也被探究[114]。在冲击载荷下,纳米铝先发生形变,然后扩散氧化,最后随着爆轰产物的膨胀而被拉伸和分解,并导致进一步加速了纳米铝的氧化。
由于铝参与到氧化和二次反应中,一些反应机制(如核-壳演化、原子扩散、热传导和微爆炸等)也被揭示。Zeng等[115]研究核-壳结构的Al/Al2O3纳米颗粒的扩散和氧化行为。根据原子均方位移发现,纳米颗粒的初始反应是由氧原子向内扩散引起的。在加热后原子扩散率随着壳厚度的增加而降低。反应后期,扭曲的(AlO)n(n=3、4和5)团簇从系统表面喷出,纳米粒子发生了微爆炸。Liu等[116]研究了铝颗粒从点火到燃烧的热力学和结构特性。研究发现氧化壳阻碍了核心Al原子的热扩散并导致异常熔化。在加热期间,核壳界面处的应力变化最剧烈。燃烧过程由颗粒和外部氧分子中组分的扩散行为控制。壳层较薄的颗粒具有较短的着火延迟和较高的燃烧温度,其直接影响燃烧期间的辐射传热率。Chu等[117]建立了纳米铝和团聚的理论模型用以研究流动时(流速为1～7km/s)颗粒之间的相互作用和热传导行为,模型考虑了流体与颗粒之间的热传导和对流传热、核壳反应、热传导、颗粒弥散和阻力。随着流速的增加,单个纳米铝和团聚体经历了从扩散氧化到各向异性氧化和爆炸氧化的转变。Li等[118]揭示了纳米铝被氧气氧化的细节,并详细阐明微爆炸加速氧化的机制。纳米铝的微爆炸发生在高温和高密度的O2环境中,氧化所释放的热量加速热点和孔隙的形成,同时也促进纳米铝的快速蒸发。此过程具有颗粒分散快、反应物消耗快、热释放和压力增长快、反应和中间体复杂等特点。Zhang等[119]研究了5种不同粒径(3～7nm)纳米铝颗粒的氧化过程,发现粒径越小,铝粉中活性铝含量越低,活性铝含量降低速率越快。Chu等[120]研究了纳米铝高温高压下的传热传质和反应过程,揭示了纳米铝氧化机理。纳米铝经历了预热、熔化、铝成核和壳层氧化4个阶段。较高的初始环境温度和气体压力会缩短预热时间且促进氧化反应。Zhang等[121]利用ReaxFF-MD研究了纳米铝的氧化反应,针对链状氧化物成核和生长等科学性问题进行深入探讨,对调控纳米结构具有现实意义。Song等[122]研究了纳米铝在多种气态氧化物(CO2、CO、NO2和NO)中的氧化行为,阐明了纳米铝形态演化过程中表面氧化、链状产物形成、中空形成和碳沉积等机理。
三氢化铝作为添加剂可以显著提高推进剂比冲和减少喷嘴烧蚀。然而,由于Al—H键较弱导致氢化铝的稳定性差,对环境相非常敏感。近年来,AlH3和炸药的相互作用也被广泛研究。Li等[123]使用ReaxFF-lg研究了RDX和AlH3复合炸药体系的相互作用。首先,AlH3释放H2,随后H2与RDX的分解产物(NO2和CO2)反应,导致H2O、NO和CO的增加。Song等[124]研究了AlH3在氧气中的反应机制。结果表明,在AlH3的氧化中受温度和结构的强烈影响,根据形态可分为球体、缺口球体、大分支和原子链4种典型形状。氧化壳层阻碍内部氢气团的释放,相反温度提供给氢气团足够的动能并使得氧化壳层不稳定,进而导致AlH3在氧化过程中产生丰富的形态。初始快速脱氢会向外排挤氧气进而一定程度上抑制表面的氧化。脱氢反应产物有H2和H原子,高温促进H原子的产生。除了AlH3在氧气中的燃烧,Song等[125-126]还模拟了AlH3与其他含氧气体(H2O、CO2、CO、NO2和NO)的反应。AlH3经历脱氢、铝成核、微爆炸和氧化阶段。相较于AlH3在氧气中的氧化,在这5种氧化物中未表现出形态多样性,主要是AlH3反应初期表面形成的氧化物很薄。
铝在环境相很容易受到腐蚀而失活。Song等[127]用分子动力学研究了碳纳米管封装纳米铝的过程。结果表明,纳米铝通过范德华力和惯性进入碳纳米管中。在管内,纳米铝以2.27Å/ps的速度向碳纳米管帽处移动。由于帽的吸引,其速度增加到3.17Å/ps。封装纳米铝颗粒的氧化过程分为核壳分离、壳损伤和核壳破裂。温度和碳纳米管长度影响内部纳米铝的氧化行为,为纳米铝的自组装和防护提供原子视角。Habibollahi等[128]研究了AlH3被HTPB、癸二酸二辛酯和二乙基醚3种涂层包裹下的燃烧过程,发现温度分别收敛到1007、1018和1012K,对应的势能分别为-502.11、-495.21和-529.17eV。Feng等[129]研究了AlH3和HTPB的复合炸药模型。由于HTPB的包裹,相较于纯AlH3氧化,复合炸药中AlH3的氧化明显延迟。此外,观察到HTPB/HTPB中间体的分解发生在纳米粒子表面,并且一些分解产物与纳米粒子结合。同时,AlH3纳米颗粒通过脱氢或脱氢促进HTPB的引发。Soleymanibrojeni等[130]通过研究环氧树脂、水和纳米铝的相互作用,得到了3种不同交联环氧树脂的玻璃化转变温度。Liu等[131]研究了含氟聚合物PTFE对纳米铝燃烧性能的改善效果。结果表明,在相同的加热速率和燃烧温度下,当温度超过铝的熔点后,裸露的两个纳米铝会烧结在一起。但当PTFE涂层覆盖时,分解产生的Al-F团簇将颗粒推开并增加反应表面积。在加入氧气环境时,覆盖PTFE的纳米铝的烧结明显延迟。
除了上述应用外,在复杂化学反应过程中,ReaxFF-MD模拟结果的分析至关重要且具有一定的挑战性。一方面,ReaxFF-MD模拟的轨迹采用键级的表现形式,相邻原子之间缺少明确的化学键类型; 另一方面,考虑到周期性边界条件,需要将其轨迹坐标转换为正常坐标。最后通过键级截断值判定原子间的成键信息。目前,主要的后处理程序包含反应分子动力学可视化和分析工具(VARxMD)[132]、化学轨迹分析器(ChemTraYzer)[133]和反应网络生成器(ReacNetGenerator)[134]。这些程序能够识别、量化和评估来自单个模拟的大量反应,包括短寿命反应中间体的反应,并得到关于分子碎片、化学反应方程式和反应网络等有价值的信息。此外,ChemTraYzer还可以获得相关反应的动力学分析结果。
ReaxFF反应力场促进了计算化学在凝聚相炸药体系中的应用,拓展了研究体系的规模和时间尺度,其主题涵盖了炸药高温高压下的热解、冲击加载下的响应、激光点火起爆等前沿内容。反应动力学模拟提供探索微观结构反应性、动态演化、初始化学分解路径以及中间和稳定产物分布等途径,尤其是在复杂热化学反应机制方面,如结构官能团、颗粒尺寸、加热速率,为炸药体系非平衡模拟研究提供理论见解,并且在某些情况下提供了半定量预测。然而,从服务于实际的角度,ReaxFF反应力场作为一种计算化学和计算材料学的研究工具,仍存在一些问题:
(1)准确性不足。ReaxFF反应势函数的参数优化目前仍然主要基于已明确的反应过渡态和势能面。应对含能材料分解等复杂化学反应时,ReaxFF对众多未纳入训练集的化学反应表现不佳,因而在整体上限制了其预测活化能和反应热等关键反应参数的准确性。
对此,新兴的机器学习等方法是重要突破方向。在人工智能框架下,由于其能够探索的化学反应种类和数量远远超出人们手动搜索的结果,有望大幅提高训练出的ReaxFF反应力场对于复杂反应的准确性。
(2)尺度迁移性不足。由于ReaxFF的势函数种类和数量远远多于传统非反应力场,因而其模拟的体系规模达到10nm尺度已经是极限,距离含能材料反应的重要影响因素的真实特征尺度0.1～100μm仍有一段距离,如颗粒、晶界、穿晶断裂面和孔隙等。如何与现有的连续体模型如统计热点模型等结合,使其具有含能材料真实建模的性能预测能力,从而服务于实际战斗部和固体发动机等装置的设计和虚拟实验,是一个急需解决的问题。
目前,已经有一些基于实验标定的连续体仿真模型,但其理论基础仍然是唯象的,未能充分考虑微观介观影响因素。若将其中的局域反应度等关键指标的唯象规律替换为基于ReaxFF的反应规律,将是一个包含力学与化学动力学耦合响应的、能够打通含能材料多尺度全流程仿真计算途径的潜在方法。在使用ReaxFF研究时,还应考虑其应用条件,挖掘其内在规律,以及与实验数据的关联性并相互佐证。由于反应力场大量应用,建立可查找、可访问、可操作和可重用数据管理平台,对在线可视化、在线分析以及交叉应用也是必要的。
