高温下含缺陷CL-20 初始化学反应的分子动力学模拟

胡靖伟，甘 强，冯长根，李昌霖，朱双飞，程年寿

（北京理工大学爆炸科学与技术国家重点实验室，北京 100081）

1 引言

六硝基六氮杂异伍兹烷（简称HNIW，或CL-20）是已应用能量密度最高的单质炸药，1986 年由美国Nielsen［1］首次合成，其热安定性较好，最大爆速、爆压、密度等均高于奥克托今（HMX），能量输出比HMX高10%～15%［2］。为深入研究CL-20 爆轰机理，近年来许多理论模拟和实验研究CL-20 的初始化学反应，较少考虑缺陷的影响。研究含缺陷CL-20 的初始化学反应，对于了解其复杂的微观反应机理以及安全应用均具有重要意义。

1993 年Patil 等［3］首先采用光谱法研究了CL-20 热分解，计算出CL-20的热分解动力学参数。Tureotte等［4］通过热重分析研究了CL-20 的热分解过程，并对热分解产生的气体进行表征。陈松林等［5］采用布氏压力计法研究了CL-20 在真空条件下的热分解，并估算出70 ℃下的贮存寿命为43.8 年。许丽娟等［6］发现CL-20 热分解后的气体产物主要有NO、CO、CO2和H2O，计算了CL-20 的热力学参数和动力学参数，确定了活化焓、活化熵和活化自由能。目前关于CL-20热分解实验研究已取得较大进展，但只局限于宏观性质如分解温度、活化能等研究，缺少热分解反应中分子结构的变化、化学键的生成和断裂以及反应路径等信息。

近年来理论模拟方法在CL-20 单分子热解、极端高温热解研究中得以广泛应用。2005年Okovytyy等［7］采用量化计算研究了CL-20 单分子热分解反应，认为其初始反应是N—NO2断键形成NO2自由基，之后骨架上C—N 键断裂，导致开环以及H 原子转移。Isayev 等［8］发现CL-20 分解过程主要包括N—N 键断裂并伴有NO2生成、C—C 和C—N 键断裂导致开环反应。张力等［9］采用ReaxFF 力场研究发现，ε-CL-20 超晶胞发生N—N 键断裂，导致硝基官能团脱离，形成NO2分子。对α-CL-20 及其水合物在1000～3000 K 的分解研究发现，含水分子体系反应速率常数增大，水分子促进了CL-20 第二阶段的分解［10］。Wang 等［11］研究了2000～3500 K 下CL-20 热分解初始反应路径，发现温度不影响CL-20 的分解路径，随着温度升高，CL-20的分解速率增加并且笼状结构更早被破坏，温度会极大地提高H2O 和N2的生成速率。Ren 等［12］研究了800～3000 K 温度下CL-20 的热分解过程，同样发现CL-20 初始分解路径主要是N—NO2键断裂反应和C—N 键断裂引起的开环反应，高温下N—NO2键断裂反应显著增多，而C—N 键断裂开环反应随温度升高变化不大。Wang 等［13］采用ReaxFF-lg 力场研究不同压缩密度下ε-CL-20 超晶胞热分解，发现密度越大热分解反应释放的能量越小，反应速率越快，但反应活化能不随密度变化而改变。

含能材料中普遍存在杂质、空位、孔隙等缺陷，会降低黑索今（RDX）、奥克托今（HMX）等含能材料反应活化能，加快发生热分解，同时有助于热点的生成和演化，导致含能材料感度提高［14-16］。2000 年Kuklja等［14-15，17］研究了位错对RDX 热分解反应的影响，发现缺陷会激发其周围的电子，促进RDX 的N—NO2键的解离，从而促进RDX 发生分解。2003 年黄亨建等［18］采用急速重结晶方法研究发现晶体缺陷使HMX 热分解起始温度及相应表观活化能明显降低，热感度增加，热安定性降低。Zhou 等［16］发现分子空位缺陷会大大改变HMX 的N—NO2键断裂、HONO 离解以及分子主环断裂等分子内反应的频次，增大了HMX 的反应速率常数、降低了其反应活化能。含能晶体内部的微观缺陷将在外界能量作用下形成高温“热点”，使含能晶体对外界感度显著增大，从而对热稳定性产生较大影响。当前晶体缺陷对CL-20 热分解反应影响报道较少，有必要展开相关研究。

为探讨晶体缺陷对CL-20 热分解机理和反应动力学的影响，采用分子动力学方法和ReaxFF-lg 力场，研究在不同温度（1500～3500 K）下含空位缺陷CL-20 的热分解反应过程，获得了CL-20 势能的演变规律和反应动力学参数，分析不同高温下初始反应路径和产物演变，讨论了空位浓度对CL-20 热分解反应的影响，期望对CL-20 热分解研究提供参考。

2 计算方法

2.1 ReaxFF-lg 力场势函数

ReaxFF-lg 反应力场通过键级判定原子间距离，并用于描述化学键断裂和生成。当化学键断裂时，与价键相关的能量和力变为零，采用该反应力场计算得到系统的总能量为［19］：

式中，Elp、Eover、Eunder、Ebond为与化学键伸缩有关的势能项，Eval、Epen、Ecoa为与键角弯曲有关的势能项；Etors、Econj为与二面角扭转有关的势能项，EH-bond为氢键能、EvdW、ECoulomb为非价键作用力，最后是描述分子间长程作用的修正项，各项单位均为kJ·mol-1。

2.2 模型的构建

从剑桥数据库中获得ε-CL-20 单晶数据［20］，将单胞沿a、b和c方向分别扩展4倍、3倍和3倍，构建4×3×3的ε-CL-20 超晶胞模型，包含36 个单胞，144 个CL-20 分子，共5184 个原子，记为CL-20 完美晶体。参考已报道含空位缺陷的HMX 模型［16，21-22］，从CL-20 超晶胞中分别剪切移除2 行、4 行、6 行CL-20 分子，即8、16、24个分子，构建空位浓度分别为5.56%、11.2%、16.7%的CL-20 超晶胞模型，分别含有4896、4608、4320 个原子，如图1 所示，任意两组空位的距离较远，以消除空位间的相互作用。

2.3 分子动力学模拟

分子动力学计算均使用大规模原子分子并行模拟程序（Large-scale Atomic/Molecular Massively Parallel Simulator，LAMMPS）［23］和ReaxFF-lg 反应性力场。采用共轭梯度算法［24］，根据麦克斯韦-玻尔兹曼分布［25］给原子分配初始速度。采用等温等压系综（isothermal-isobaric ensemble，NPT），在300 K 和一个大气压条件下驰豫10 ps 得到超晶胞的平衡结构。采用正则系综（canonical ensemble，NVT）和Berendsen控温方法［26］，将弛豫后的体系分别在1500～3500 K 下进行分子动力学计算。整个计算过程采用周期性边界条件，时间步长为0.1 fs。采用C++语言自编程进行产物分析。

2.4 反应速率分析

采用Arrhenius 方程［27］描述分析含缺陷CL-20 晶体热分解动力学参数：

式中，t 为时间，ps；T 为温度，K；α 为反应进度；f（α）为反应模型；k（t）为与温度相关的反应速率常数，ps-1；A为指前因子，s-1；E 为反应活化能，kJ·mol-1；R 为理想气体常数，8.314 J·mol-1·K-1。

图1 CL-20 完美晶体和不同空位浓度CL-20 超晶胞模型Fig.1 CL-20 perfect crystal and its supercell models with different vacancy concentrations

根据Arrhenius 定律将反应速率常数与反应活化能联系起来：

根据Rom 等［21］提出的反应速率的分析方法，将不同温度下CL-20 的热分解过程分为吸热反应、放热反应和终态产物生成三个阶段进行分析。在吸热反应阶段，将CL-20 分子数量随时间的变化按照一阶衰减指数式进行拟合：

式中，t 为反应时间，ps；α（t）为t 时刻的反应进度；N（0）和N（t）分别是为初始时刻和t 时刻CL-20 的分子数量；k 为吸热阶段的反应速率，ps-1。将式（4）得到的反应速率常数k 代入式（3）并进行线性拟合，得到反应活化能E 和指前因子A。

3 结果与分析

3.1 空位缺陷对CL-20 热分解反应动力学的影响

考察1500～3500 K 下空位缺陷对CL-20 热分解反应速率常数、活化能和指前因子的影响，CL-20 完美晶体和缺陷模型反应速率常数见表1，反应速率常数的对数与温度的倒数关系见图2，线性拟合得到反应活化能E 和指前因子A 见表2。

已有研究表明CL-20 初始化学反应为N—NO2键的断裂，该反应遵循一级化学反应速率方程［10-11，28］。由图2 可见，在反应初始阶段，完美和空位缺陷模型ln（k）与1000/T 的关系符合一阶衰减指数规律。从表1 可见，在1500～3500 K 温度下，反应速率常数k 均随着空位浓度的增加而增大，这说明空位缺陷加速了CL-20 的初始化学反应速率。在相对较低温度（1500、2000 K）下，缺陷对反应速率的促进作用更加明显，其中1500 K 时含16.7%空位缺陷CL-20 反应速率常数是完美晶体的2 倍以上。

表1 不同温度下含空位缺陷CL-20 反应速率常数Table 1 Reaction rate constant of CL-20 with different vacancy concentrations at different temperatures

图2 含空位缺陷CL-20 热分解反应速率常数的对数与温度倒数的对数拟合关系Fig.2 The fitting curve between the logarithm of the thermal decomposition rate constant of CL-20 with vacancy defectsand the logarithm of the reciprocal temperature

表2 含空位缺陷CL-20 晶体热分解的反应动力学参数Table 2 Reaction kinetic parameters of thermal decomposition of CL-20 crystal with vacancy defects

由表2 可见，空位浓度越大，CL-20 热分解反应的活化能越低，其中含16.7%空位CL-20 的热分解活化能比完美晶体低11.80 kJ·mol-1，这表明空位降低了CL-20 热分解反应能垒，促进CL-20 发生分解。本研究得出的CL-20 完美晶体热分解的活化能，与Wang 等［11］的ReaxFF-lg 分子动力学计算和Isayev 等［8］的量化计算结果接近。许丽娟等［6］测定了CL-20 热分解的活化能，由于实验通常低于1000 ℃，且与外界有温度和压力交换，导致与模拟结果存在一定差异。空位缺陷的引入降低了CL-20 的热分解反应活化能，使CL-20 热分解速率显著增大，与含位错、空位RDX、HMX 热分解规律［14，16］类似。

3.2 含空位CL-20 势能和产物种类演化规律

1500～3500 K 热分解过程中CL-20 分子平均势能随时间变化曲线见图3，相对产物种类随时间的变化曲线见图4。为了消除CL-20 初始体系中分子数不同的影响，分别对各体系CL-20 分子数进行归一化处理。

由图3 可见，在各温度下，不同空位浓度下CL-20分子平均势能迅速升高后缓慢下降，表明初始分解反应以吸热反应开始，势能达到最大值后初始分解产物发生次级反应，体系的势能逐渐降低，与已报道CL-20势能变化规律一致［11-12］。温度显著影响CL-20 平均分子势能达到平衡的时间，在3000～3500 K 温度下势能曲线在40～60 ps 后平衡，1500～2500 K 时势能100 ps内分子势能未完全达到平衡，可见在100 ps 时间尺度内，反应温度为3000 K 及以上时含空位CL-20 充分反应，而3000 K 以下温度梯度未完全反应。选取势能变化趋势差异较大的两个温度（2000 K 和3000 K），进行热分解深入分析。对比CL-20 完美晶体可见，空位缺陷浓度越高，CL-20 势能降低速度越快，这说明空位缺陷降低了分解反应的能垒，促进了CL-20 的热分解，温度较低时势能差异更显著。

由图4a 可见，反应温度为2000 K 时，在20 ps 产物种类达到最大值，之后略有下降，说明反应还停留在中间产物阶段，没有大规模生成终态产物。100 ps 时含16.7%空位CL-20 分子平均产物种类约为1.5 个，而完美晶体只有1.2 个，这说明空位缺陷使CL-20 热分解产物种类增加，增大了晶体的破碎程度，促进了CL-20 热分解。由图4b 可见，反应温度为3000K 时，产物种类到达峰值在4.7 ps，相比于2000 K 产物种类增大；100 ps 时产物种类趋向平衡，CL-20 分子平均分解产物为0.8 个左右，这比2000 K 时产物种类减少。是因为3000 K 时CL-20 完美晶体和缺陷体系产物更接近于终态产物，反应比2000 K 时更加完全，因而物质种类更少，与Wang［11］对CL-20 完美晶体热分解研究结果一致。

图3 1500～3500 K 时含空位缺陷CL-20 平均分子势能随时间变化曲线Fig.3 Curve of average molecular potential energy of vacancy defect CL-20 with time at 1500-3500 K

图4 含空位缺陷CL-20 分子产物种类随时间变化曲线Fig.4 Curves of products species of vacancy defect CL-20 with time

3.3 空位缺陷对CL-20 初始反应路径的影响

对含缺陷CL-20 初始化学反应进行统计分析，2000 K 和3000 K 时CL-20 发生初始化学反应时高频基元反应式、持续时间和反应频次见附表S1、S2，表4、图5～图7 是根据附表统计的结果。其中，表4 为统计出的CL-20 主要初始化学反应类型及其反应式，图5是CL-20 两种脱硝基反应示意图，图6 和图7 分别是2000 K 和3000 K 下CL-20 热分解反应类型分布图（扣除逆反应）。

由表4 可见，CL-20 初始化学反应主要有脱硝基反应（包括步步脱硝基和连续脱硝基两种类型反应）、CL-20 与小分子反应、CL-20 与大基团的反应。CL-20热分解初始化学反应主要是CL-20 分子环上发生N—NO2断裂、脱去一个硝基，其次发生其余硝基脱落或者主环的裂解反应，少部分CL-20 发生多重脱硝基反应。 随着硝基脱离CL-20 分子，生成NO2和C6H6N10O8、C6H6N9O6等中间体，会发生CL-20 与NO2、CL-20 以及中间体的反应，但后者发生频次不高。表S1 和表S2 的统计结果也表明，脱硝基反应是频次最高的热分解反应，CL-20 与小分子或大基团的反应频次很低，随着空位浓度增大，反应频次有一定提高。

由图5 可以看出，2000 K 时CL-20 初始分解反应主要为步步脱硝基反应，CL-20 依次脱去一个硝基分别生成化合物Ⅰ（C6H6N11O10）、化合物Ⅱ（C6H6N10O8）和化合物Ⅲ（C6H6N9O6）。这与Isayev 等［8］得到的单分子CL-20 分解路径相同。3000 K 时出现一定频次连续脱硝基反应，即CL-20 连续脱落两个硝基生成化合物Ⅱ，连续脱落三个硝基生成化合物Ⅲ。

表4 CL-20 初始化学反应类型及其主要反应Table 4 Initial chemical reaction types and main reactions of CL-20

图5 2000 K（a）和3000 K（b）时CL-20 分子脱硝基反应Fig 5 Denitration reaction of CL-20 at 2000 K（a）and 3000 K（b）

图6 含空位缺陷CL-20 热分解反应类型分布Fig.6 Thermal decomposition reaction distributions of CL-20 with vacancy defects

由图6 可以看出，空位缺陷对CL-20 各类反应频次产生显著影响。2000 K 和3000 K 下主要发生脱硝基反应，随着空位缺陷浓度增大，反应频次略有下降。另一种重要反应是CL-20 及较大中间产物的开环反应，反应频次随着空位浓度增大逐渐增加。说明空位缺陷加速了CL-20 体系的裂解，使CL-20 开环反应频次增加，促进了CL-20 分解反应。相比于2000 K，3000 K 下各类分解反应及其逆反应频次减少，CL-20热分解进程提前。

Zhou 等［16］研究表明，在高温下HMX 主要分解路径为N—NO2键断裂、HONO 离解以及分子主环断裂，空位缺陷促进N—N 键裂解和开环反应，但阻碍了HONO 的形成。本文结果表明空位的引入并不改变CL-20 主要初始反应路径，含空位CL-20 热分解仍以N—NO2键断为主，其他反应比例有一定增加。凝聚态含能材料热分解时，外界压力也会产生重大影响。文献［29-30］显示，硝基甲烷和HMX 在低压下的初始分解以分子内反应为主，而高压下以分子间反应为主。Wang 等［13］指出高压下CL-20 热分解以分子内反应为主，本研究表明常压下CL-20 完美晶体和缺陷模型仍主要发生分子内反应，高压下反应则需进一步研究。

3.4 空位缺陷对CL-20 小分子产物的影响

图7 和图8 分别是2000 K 和3000 K 时，完美和含空位CL-20 主要小分子产物含量随时间变化曲线。从分解产物变化可见，完美晶体和含空位CL-20 初始分解反应是N—NO2键的断裂生成NO2，这与文献报道一致［7-8］。

空位缺陷不会改变CL-20 主要初始热分解反应类型，但会对产物分布产生一定影响。如图7 所示，2000 K 时，完美CL-20 最主要的初始分解产物是NO2，次级反应消耗掉大量的NO2，生成N2、NO3等中间产物。其中NO3是仅次于NO2的产物，进一步分解生成N2、H2O 等终态产物，与Wang 等［11］、Ren 等［12］以及Isayev 等［8］的模拟结果基本一致。CL-20 分子骨架上C—N 键断裂，发生开环反应，主要小分子产物为CO2，该反应在CL-20 的热分解初期不占主导地位，20 ps 以后CO2数量剧增，说明骨架上C 原子在快速脱离，这与Xiang 等［31］的研究结果一致。随着CL-20 中缺陷浓度的增大，NO2峰值逐渐增大，其他小分子产物无显著变化。相比之下，如图8 所示，3000 K 时初始产物NO2峰值较2000 K 显著增大，由2 上升到2.5，N2、H2O、CO2、H2等终态产物含量增加，其中N2最大值由3 增加到4.5 左右、CO2最大值由0.5～1 增加到2.5 左右、H2O 最大值由1 增加到1.5 左右且提前达到平衡，中间产物NO3峰值急剧减小。相比于2000 K 时NO、N2O 等含量极少。说明在较高的温度下CL-20的分解程度更高，副产物也更少。该温度下空位缺陷未见明显影响，不同空位浓度CL-20 热分解小分子产物差异不显著。

图7 2000 K 时含空位缺陷CL-20 分子主要分解产物相对含量随时间变化曲线Fig.7 Curve of relative content of the main decomposition products of vacancy defect CL-20 with time at 2000 K

图8 3000 K 时含空位缺陷CL-20 分子主要分解产物相对含量随时间变化曲线Fig.8 The curve of relative content of the main decomposition products of vacancy defect CL-20 with time at 3000 K

4 结论

采用ReaxFF-lg 力场和分子动力学方法，研究了1500～3500 K 下空位缺陷对CL-20 初始热分解反应、反应动力学和热分解产物的影响规律，得出以下结论：

（1）高温下CL-20 的初始化学反应以硝基脱落为主，引入空位缺陷后，脱硝基反应比例略有下降，开环反应比例增大，反应首先发生在空位缺陷周围，在CL-20 中可能形成热点。

（2）空位缺陷降低了CL-20 热分解活化能，增大了反应速率常数，温度升高和空位浓度的增加，都会使热分解反应加快。空位浓度越大对热分解促进作用越显著，且随着温度升高影响逐渐减小，2000 K 和3000 K时含16.7% 空位CL-20 反应速率常数分别是完美CL-20 的1.7 倍和1.4 倍，说明较低温度下空位缺陷的促进作用更为明显。

（3）高温下CL-20 主要发生分子内反应，当温度从2000 K 上升到3000 K 时，小分子产物NO2、N2、H2O、CO2等产物峰值显著增加，但中间产物NO3等的含量降幅较大。2000 K 时空位提升了NO2峰值，其他小分子产物无显著变化，3000 K 时空位对CL-20 热分解小分子产物影响不显著。