基于碳骨架氧化放能的设计思路, 高能炸药的边界已经从黑索今(RDX)、奥克托今(HMX)等传统炸药拓展到了八硝基立方烷(ONC)^[1], 六硝基六氮杂异伍兹烷(CL-20)^[2]等新型的高能化合物。这些新型化合物充分利用了笼状结构中蕴藏的张力, 使之分解时释放出更多的能量, 从而提升爆炸威力。由于这些新型化合物普遍具有敏感度高, 合成步骤复杂, 成本高等缺点, 所以人们继而将目光转向了高氮材料, 由于N≡N三键的键能远高于单键或双键, 这些材料在分解形成氮气时能放出大量能量, 又因为产物N₂对环境友好, 故高氮材料被认为是含能材料设计与合成的未来^[3]。其中, 5, 5′-联四唑-1, 1′-二氧二羟铵(TKX-50)以其易于合成, 威力接近CL-20, 感度接近RDX的特性被认为是未来最有前途的含能材料之一^[4]。

自Fischer等^[4]在2012年首次合成TKX-50, 并报道了其优异的物化性能与能量特性以来, 研究者已经在敏感性^[5], 合成方法^[6-7], 热分解机理^[8-13]以及与其他炸药的复合性质^[14]等方向对TKX-50开展了一系列的研究。在针对TKX-50热分解机理的相关研究中, 微观反应机理的研究主要依赖于计算模拟, 研究相对较少, 且受限于后期反应的复杂度以及计算资源的限制, 计算模拟工作主要集中于初始的分解过程, 如Qi^[8]指出质子转移降低唑环分解能垒, Lu^[10]指出压力会抑制质子转移的过程, 均针对起始的质子转移。到目前为止, 对于TKX-50整体反应路径和机理, 特别是后期过程没有特别深入的研究。因此, 本研究利用从头算分子动力学方法对TKX-50在不同温度范围下及不同密度下的热分解全过程进行模拟, 并对产物中各原子的来源进行统计分析, 查明了主要产物N₂的反应路径, 并藉此给出温度和压力对整体反应速率及机理的影响。

采用从头算分子动力学方法, 计算软件为北京应用物理与计算数学研究所的含能材料团队开发的第一性原理计算软件HASEM^[15]。TKX-50晶体单胞由4个(NH₃OH)⁺阳离子与2个氧代偶联四唑环形成的阴离子(C₂N₈O₂)^2-构成(图 1)。构建的模拟体系大小为2×2×1, 四个单胞, 共192个原子。为了排除热源影响, 研究TKX-50分子通过自身分解时能量转移及释放规律, 动力学模拟时采用了微正则系综(NVE), 而非常用的等温等压系综(NVT)。模拟体系的初始温度设定为2000~3500 K, 以保证能在计算模拟可承受的时间范围内观察到整个分解过程。体系在不同压力下的模拟可通过改变对应模型的体积实现, 设计了体积从0.8V₀~1.3V₀的模型(V₀为300 K零压下模拟体系的体积), 与之相对应的压力为17~0 GPa。模拟的动力学时间步长为0.5 fs, 总模拟时间为5 ps。通过测试发现体系最低在2500 K下可在5 ps内达到平衡状态, 因此不同压力下模拟的初始温度统一设定在2500 K。对中间产物的产生及湮灭、原子间成键断键过程以及原子来源追踪的相关分析工作主要由HASEM程序所集成的团簇分析工具所完成, 该工具主要基于键长进行团簇, 原子间成键的识别和判断。

图 1(Fig. 1)

图 1 TKX-50的分子结构和晶体结构示意图 Fig.1 Schematic diagrams of molecular structure and crystal structure of TKX-50

图 2展示了模拟体系的初始温度为3000 K下热分解模拟过程中体系势能、主要中间产物数量随时间的变化。从图 2可知, TKX-50在开始分解之前有一个短暂的吸热过程, 在该过程中发生的是(NH₃OH)⁺基团的振动和旋转, 偶联四唑环的振动以及两个环之间的相对扭转。同时还有质子转移在该阶段发生。而在放热过程中主要为四唑环以及铵盐分解生成N₂及H₂O。最后生成的气态产物为N₂与H₂O以及少量的NH₃, 中间产物还包括了NO以及少量N₂O。通过对中间及最终产物变化的分析发现N₂主要来源于唑环的直接分解; 而H₂O的形成则与NH₃OH⁺的分解相关。据文献[8]报道, 质子转移过程中一方面降低了唑环分解的能垒, 另一方面会将O原子从唑环环裂后留下的C₂N₄O₂₂^-中剥离下来, 生成H₂O, 并帮助氧与体系中的碳结合。在保持体系体积不变的情况下, 化合物体系内约有75%的N转化为N₂。由于体系为负氧平衡的缘故, 其中阳离子中的N约50%形成N₂, 剩下的主要以NH₃或NH₂的形态存在。如果一开始提供的能量不够的话, 即一开始的反应温度不高, 最后生成的N₂中来自阳离子的比例会有所降低, 同时唑环分解生成的残基将很难进行二次反应。

图 2(Fig. 2)

图 2 产物数量和势能随时间的变化曲线(初始温度为3000 K) Fig.2 The change curves of product population and potential energy with time (initial temperature is set at 3000 K)

对于N₂的生成路径, 一般认为是由唑环的环裂直接生成, 即图 1中N(2)和N(3)集体从唑环上脱落, 形成N₂。Qi^[8]采用从头算方法计算了唑环上不同为位置环裂的能量, 同样认为N(2)与N(3)原子脱落形成N₂是能垒最低的环裂方式。需要注意的是虽然该位置处环裂发生所需的能量最低, 但在真实的反应过程中, 高反应能垒的过程仍有一定的几率发生, 特别是在2500 K至4000 K的模拟温度下, 能够观察到唑环环裂生成N₂时所呈现的各种方式和顺序。

由于该过程是一个典型的单分子反应, 故推测其生成速率应与反应物浓度无关, 而从图 2中可以看出其主要产物N₂的演化曲线在反应前、中期均呈线性变化, 只在反应末期其变化速率有略微降低, 似乎印证了之前的推测。不过由于采用了NVE系综, 在整个反应期间, 温度与压力都在实时变化, 因此无法从产物曲线呈线性变化的规律导出N₂的生成过程是零级反应, 即典型单分子反应的结论。

通过追溯N₂分子中N原子的来源, 发现N₂的形成有主要有三条路径。最主要的一条路径对应于四唑环在N(3), N(4)之间以及N(2), N(1)之间断裂, 直接生成N₂, 其分子中的N均来源于N(2)与N(3), 将此路径记为P1。此外还有四唑环在N(2), N(3)之间以及C(1), N(1)之间断裂, 同样直接生成N₂, 分子中的N来自于N(1)与N(2)位置, 将该路径记为P2。除以上两种四唑环裂直接生成N₂的方式之外, 还有N(4)与O(1)形成的NO与NH₃相互作用生成N₂的途径, 将其记为P3, 此外还有氮气分子中两个N均来自N(4), 将其路径记为P4, 最终得到了N₂中N原子来源的统计图(图 3), 图 3显示了P1, P2, P3, P4这几条主要路径及另外三种次要路径(以N原子的来源标记为N(4)/N(3), N(2)/N(5)和N(1)/N(4)过程)以及其他所有可能路径在总生成路径中所占的比例。同时, 由图 3可知, 在2500 K以及3000 K下由P1路径生成的N₂占总量的40%以上, 由P2, P3, P4路径生成的各占10%左右。由于P3过程与N(4)/N(3)过程均需要N(4)原子, 所以两者比例加和在25%左右。3000 K下P3路径占比高达23%, 而N(4)/N(3)过程只有3%, 这可能是模拟当中出现的随机因素导致P3的占比过高。

图 3(Fig. 3)

图 3 N2不同生成路径的配比 Fig.3 Proportion of different paths for N2 formation

图 4描述了不同压力下, 通过P1至P4路径所形成的N₂分子数目随时间的变化关系。由图 4可知P1过程总是最先发生的, 其次为P2, 剩下的过程均在此之后, 由于P1与P2过程是一步的基元反应, 剩下都是间接反应生成, 所以慢于P1与P2过程。而P1与P2过程的先后顺序以及比例, 则反映了唑环上不同键的强度。通过对四唑环以及NH₃OH⁺各个键做COHP(Crystal Orbital Hamilton Population)分析, 可以获得各原子间成键强度的信息。表 1列出了各个键的COHP计算结果及断键过程所对应的N₂生成路径。由表 1可知, N—H键最弱, 对应的过程为质子转移, 其次为(NH₃OH)⁺中的N—键, 是羟基及H₂O生成的主要来源, 然后为唑环上的1号与3号键, 对应于路径P1, 然后是2号键与6号键, 对应于过程P2, 由于6号键键能较高, 不易断裂, 因此相比较P1难以发生。

表 1(Tab. 1)

表 1 TKX-50晶体中各键的COHP及对应的断键方式 Tab. 1 The COHP value and corresponding breaking bond mode in TKX-50 crystal number bond COHP/eV identifier of different reaction path/process 1 N(1)—N(2) -4.38 P1 2 N(2)—N(3) -4.94 P2 3 N(3)—N(4) -4.67 P1 4 N(4)—C(1) -6.12 P3 5 N(1)—C(1) -5.83 P2 6 N(4)—O(1) -4.84 P3/P4 7 N(5)—H(1) -2.84 Proton transfer 8 O(2)—N(5) -3.46 OH- group formation 9 O(2)—H(2) -2.78 Proton transfer

表 1 TKX-50晶体中各键的COHP及对应的断键方式 Tab.1 The COHP value and corresponding breaking bond mode in TKX-50 crystal

图 4(Fig. 4)

图 4 不同体积下从P1-P4路径的N2产物数量随时间变化的曲线 Fig.4 The change curves of N2 product population from different formation path (P1-P4) with time at different volume

通过比较2500, 3000, 3500 K下TKX-50分解过程中体系势能随时间的变化曲线(图 2), 发现温度对TKX-50的分解速率有较大的影响。温度越高, 其势能曲线达到平稳(势能曲线不再变化)所需的时间越短, 意味着反应速率越快。通过拟合不同温度下N₂的生成曲线的斜率, 发现其不同温度下的生成速率近似满足阿仑尼乌斯律^[16]。虽然反应速率有所变化, 然而N₂的总量在不同温度下却没有明显的变化, 均占到了N原子总量的50%左右。而由图 3可知, 唑环开环的位置同样也受到温度的影响, 在2500 K和3000 K下, P1与P2是主要路径。当温度上升至3500 K时, N(1)/N(4)路径以及other路径的占比均由原先的10%不到上升至16%。说明在高温下, 所有键的活性都有所提升, 原先因键能较高不容易断裂的位置也可能发生断裂, 而生成新的反应路径。

模拟结果显示0.8V₀下TKX-50的分解速率高于正常密度以及1.2V₀下的速率(图 4)。由于N₂是主产物且性质稳定, 不容易与其他产物进行二次反应, 因此可以通过N₂分子的数量变化来判断反应的进程。如果将N₂的产率为零, N₂分子数保持不变的时间定义为体系达到平衡的时间, 记为τ, 则有τ_0.8V0=5.0 ps, τ_1.0V0=3.2 ps, τ_1.2V0=2.1 ps。而反应速率k∝1/τ, 由此可知在1.2V₀(低密度)下, 反应速率为正常密度下的1.52倍; 而在0.8V₀(高密度, 压缩)条件下, 反应速率仅为正常压力的0.64倍。

通过对不同路径生成的N₂比例分析, 发现在压力的作用下, 不但不同路径下N₂的生成速率均受到了抑制, 并且P3这条反应路径几乎消失。由3.2.1节的分析可知反应路径P3对应于NO^-与NH$^{+}_{\text{n}}$的相互作用, 由于NO^-与NH$^{+}_{\text{n}}$在空间上相互隔离, 两者的反应需要通过基团扩散实现。而在压力的作用下, 物质之间空隙较小, 扩散速率较常压下低, 因此该反应路径在压力下受到抑制。此外, 利用均方位移(Mean Square Displacement, MSD)方法计算了各压力下不同元素的扩散速率, 结果显示N原子在0.8V₀, 1.0V₀和1.2V₀体积下的扩散速率D分别为2.62×10^-10, 4.09×10^-10 m²/s和6.13×10^-10 m²/s, 其比值为0.64:1.00:1.50。C, O原子的扩散速率与N原子基本一致, 而H原子的扩散速率约为N原子的2倍。

通过统计分析主要产物N₂中N的来源, 发现TKX-50在热解时存在多条N₂的生成路径, 并给出了这些路径在温度和压力变化时所表现出来的行为差异。

(1) TKX-50在热分解生成N₂时有三条主要的反应路径, 当中最主要的一条, 也即目前最为研究者所熟知的为四唑环通过环裂(P1路径)直接生成N₂, 其他两条路径包括四唑环不等价位的环裂(P2路径)以及四唑环与NH₄的相互作用(P3路径)。不同的反应路径在N₂生成过程中所占的比重与其相对应的断键过程中的键的强度相关。

(2) 温度影响TKX-50的反应速率, 通过拟合不同温度下反应达到平衡的时间所得的反应速率模型近似服从阿仑尼乌斯律, 且N₂生成路径P3/P4以及一些次要的反应过程在总反应中的占比随温度的上升而增加。不同温度下N₂的总生成量维持在体系总N原子数目的50%左右, 剩下的N以NH₃, NO以及HCNO的形式存在。

(3) 压力与反应速率呈负相关, 增加体系的压力会降低反应速率。同时, 与扩散相关的反应路径(如P3)会被抑制, 单位时间内通过该路径产生的N₂数目下降, 进而降低总的反应速率。

计算结果还显示不同压力下原子的扩散速率与表观反应速率呈很高的线性相关性, 不过由于本研究样本数目太少, 不足以做出确定的结论。但是认为可以在未来通过对原子及分子扩散速率的进一步研究, 继而将压力与TKX-50的反应速率联系起来。