冲击作用下CL-20/HMX 共晶力-热-化学耦合响应的动力学模拟

王 宁，苏 晶，关红波，成 杰，程立国，李 军

（1. 湖北航天化学技术研究所，湖北 襄阳 441003；2. 航天化学动力技术重点实验室，湖北 襄阳 441003；3. 湖北航鹏化学动力科技有限责任公司，湖北 襄阳 441003）

1 引言

含能材料的能量和安全性矛盾一直以来受到广泛关注，含能共晶炸药是解决这一矛盾的主要途径之一，其中较为典型代表是CL-20/HMX 共晶炸药［1-4］，该炸药的爆轰性能优于β-HMX，感度与β-HMX 相当［2］，对其感度机理的深入研究有利于新型含能共晶炸药的设计和开发。含能材料分子、晶体结构与其感度的关系，国内外开展了大量研究［4-8］。对于CL-20/HMX 共晶炸药的撞击感度，热力学引发键键级判据和动力学活化能判据［9-10］在单质含能材料中应用较为广泛。由于分子层级混合的含能共晶仍属于不同分子结构的混合，因此这些判据还在进一步发展。对CL-20/HMX 共晶的稳定性，Bolton 等［2］认为与纯组分相比，共晶中的C—H…O 氢键数目明显增多，键长减小。陶俊等［11］根据最大键长和结合能，以范德华为主的作用力“敏化”了CL-20/HMX 共混体系的结构，而根据内聚能密度和径向分布函数，得出CL-20/HMX 共晶体系由于存在长度相对较短的CH…O 氢键而导致其感度较低。若考虑晶体形貌或缺陷影响，用于预示感度的最大引发键键级判据和分子间弱相互作用强弱等将受到限制。随着计算能力提高，从动力学角度模拟晶体结构和感度的关系成为可能。

目前，立足于含能材料分子结构，在晶体结构、颗粒等空间上和介观-微观等时间尺度上，对强冲击下含能材料发生的力-热-化学变化耦合过程的研究已取得较大进展。Jaramillo 等［12］、Cawkwell 等［13-14］和Bedrov 等［15］研究了RDX 不同晶面受到冲击加载时发生的塑性行为，详细分析了位错、剪切带和分子构型的变化规律。Eason 等［16］分析了冲击诱导PETN 时发生的弹塑性行为，而Shan 等［17］则分析了沿PETN 不同晶面冲击加载时化学反应过程。此外，孔洞塌缩形成热点的机制也在不断深入，Wood 等［18］分析了含孔洞缺陷的RDX 发生冲击转爆燃的过程，利用流体动力学方法和反应分子动力学方法（ReaxFF-MD），详细分析了HNS 晶体中孔洞塌缩形成热点的过程［19］，而Shan 等［20］则在更长时间内分析了孔洞塌缩后形成的非对称热点区域。然而含能共晶炸药冲击过程模拟未见报道，因此，本研究着重分析了冲击作用下，有无缺陷CL-20/HMX 共晶发生的力-热结构变化和随后的化学反应过程，为进一步分析冲击波感度机理提供必要基础。

2 计算模型与计算过程

2.1 CL-20/HMX 共晶初始模型及结构优化

含能材料热分解和冲击压缩过程的研究常采用ReaxFF-MD 方法。由于不同时期开发的力场有差异，因此首先需要计算研究对象超胞结构来验证力场的适用 性。本 研 究 依 据ε-CL-20 晶 体［21］、β-HMX［22］和CL-20/HMX 共 晶［2］的X 射 线 单 晶 衍 射 结 果，建 立 了4×2×2 的CL-20 超胞、5×3×4 的HMX 超胞和2×3×2的CL-20/HMX 共晶超胞，再采用ReaxFF-MD 方法［23］进行不同晶体的分子动力学优化，将计算结果和试验结果对比以筛选合适力场。

可能满足CL-20/HMX 共晶计算的力场有ReaxFF-2005［24］、ReaxFF-lg［25］、ReaxFF-2014［18］和ReaxFF-2018［19］等，其中ReaxFF-lg 因添加了低梯度修正项（lg），能够更准确描述晶格参数和密度，在ReaxFF-2014 的基础上，ReaxFF-2018 增加了lg 项，并兼顾了多种含能材料反应途径和燃烧过程。因此，本研究利用ReaxFF-lg 和ReaxFF-2018 分别计算CL-20、HMX 和CL-20/HMX 晶体结构参数，以验证力场适用性。

分子动力学优化过程分为三步，即先后运行20 ps的T=5 K、正则系综NVT 和T=300 K、NVT 进行弛豫，然后选择等温等压系综NPT（T=300 K、p=1.01×105Pa）进行50 ps 分子动力学模拟。采用Nose-Hoover 控温和控压［26］，时间步长为0.1 fs。

2.2 不同晶向CL-20/HMX 共晶模型的建立及冲击加载过程细节

利用Materials Studio 中morphology 模块［27］预测真空状态下CL-20/HMX 共晶形貌，发现（100）晶面面积最大。Ghosh 等［4］根据粉末XRD 和晶体形貌认为（111）晶面是重要晶面。因此，选择（100）和（111）晶面作为冲击加载的主要方向。

考虑模拟体系大小、晶面的尺寸匹配和计算资源等因素，沿CL-20/HMX共晶的［100］晶向构建了10×1×24的模型CL-20/HMX（100），沿［111］晶向构建了5×1×48的模型CL-20/HMX（111）；沿CL-20 的［011］晶向构建了1×10×24 的模型CL-20（011），沿HMX 的［011］晶向构建了1×10×24 的模型HMX（011）。具体尺寸如表1 所示。为减少因周期性边界条件而可能增加的非键弱相互作用，在z 向均增加了2 nm 的真空层。

表1 不同晶向CL-20/HMX 共晶模型详细信息Table 1 Detailed information of various CL-20/HMX models

为分析CL-20/HMX 共晶在冲击加载下发生的力-热结构变化和化学反应，本研究先对构建的CL-20/HMX 共晶模型结构优化，再按照动量镜原理［15］开展了冲击加载模拟。冲击加载过程如下：取［100］晶向的两层和［111］晶向的四层重复单元作为固定壁，壁厚约2～3 nm，其余部分以质点速度2 km·s-1［20，28］撞击固定壁，形成反射冲击波。当反射冲击波传播到尾端时，采用吸收波阵面方法［14］，继续在NVE 系综下开展分子动力学模拟。整个模拟过程的时间步在0.1～0.01 fs之间调整，满足每步电荷收敛平衡需要和捕捉到更全面的化学反应产物。通过设定不同种类原子对的键级阈 值［29］区 分 分子碎片，后采用Sergey Zybin 提 供的mol_fra.c 源代码［30］处理。所有反应分子动力学模拟过程均通过LAMMPS 软件［31］完成。

2.3 含孔洞缺陷的CL-20/HMX 共晶模型的建立

孔洞缺陷是含能材料中最常见的缺陷类型之一。由于孔洞缺陷大小和冲击强度对孔洞塌缩诱导的热点形成有显著影响，因此结合Shan 等［20］的研究结果，本研究建立了含孔洞直径为20 nm、沿［100］晶向的模型CL-20/HMX（100）-void。为便于分析，还建立了含同样孔洞尺寸的沿CL-20［011］晶向的模型CL-20（011）-void，沿CL-20/HMX 共晶［100］晶向的无缺陷模 型CL-20/HMX（100）-large，沿CL-20/HMX 共 晶［111］晶 向 的无缺陷模 型CL-20/HMX（111）-large。具体尺寸如表2 所示。

表2 ε-CL-20、β-HMX 和CL-20/HMX 晶胞参数的实验值和模拟结果Table 2 Comparison of the experimental and calculated lattice parameters of CL-20，HMX，and CL-20/HMX

3 结果与讨论

3.1 CL-20/HMX 力场适用性分析

为筛选适用于CL-20/HMX 共晶冲击过程的力场，需正确描述CL-20/HMX 共晶结构。为此，本研究通过ReaxFF-lg 和ReaxFF-2018 分 别 优 化 了ε-CL-20、β-HMX 和CL-20/HMX 超胞结构，获取CL-20/HMX 共晶结构，各物质的晶胞参数和密度优化结果见表2。由表2 可知，ReaxFF-lg 和ReaxFF-2018 这两套 力 场 在描述含能材料晶体参数和密度时均具有较好的计算精度。

除优化CL-20/HMX 的共晶结构，为研究共晶的动态加载过程，还需计算相关高压物理特性，包括静水压和体积压缩比p-v/v0和Hugoniot 曲线相关的p-v/v0和T-v/v0及us-up关系等［20］。由 于CL-20/HMX 共晶 的 这些物理特性缺乏试验数据，故本研究以高精度的色散修正的DFT-D［32］计算结果为参考，对比分析了ReaxFF-lg 和ReaxFF-2018 描述CL-20/HMX 共晶高压物理特性的准确性，结果如图1～图3所示。由图1、图2a和图3 可知，与ReaxFF-lg 相比，采用ReaxFF-2018 计算值更接近于DFT 计算值。 由图2b 可知，ReaxFF-2018 计算的高压缩比下CL-20/HMX 共晶温度偏离未反应Hugoniot 曲线，这与部分物质发生反应有关，比如质点速度超过1.787 km·s-1时，类似结果也反 映 在Islam 等［33］的 研 究 工 作 中。 因 此，选 取ReaxFF-2018 开展质点速度为2 km·s-1时CL-20/HMX共晶力-热-化学反应过程研究。

图1 CL-20/HMX 共晶的冷静压下体积压缩比和静水压关系Fig.1 Relationship between compressive volume ratio and pressure for CL-20/HMX cocrystal

图2 CL-20/HMX 共晶受冲击后（a）p-v/v0和（b）T-v/v0关系Fig.2 Curves of p-v/v0（a）and T-v/v0（b）of CL-20/HMX cocrystal under shock loading

图3 计算获得的CL-20/HMX 共晶的us-up关系Fig.3 Calculated us-up relationship of CL-20/HMX cocrystal

3.2 CL-20/HMX 不同晶向响应特性分析

为考察CL-20/HMX 共晶是否具有明显的各向异性，分别沿［100］和［111］晶向冲击加载了CL-20/HMX 共晶。初始质点速度设定为2 km·s-1，其温度和初始物质（包括CL-20 和HMX）分解量随时间的变化关系如图4 所示。为对比分析，同时计算了沿［011］晶向冲击加载的CL-20 和HMX。

图4 沿［100］和［111］冲击加载CL-20/HMX 共晶时（a）温度和（b）未分解分子随时间的变化规律Fig.4 Evolution of temperature（a）and percentage of decomposed molecules（b）vs time along the［100］and［111］shocks in CL-20/HMX cocrystal

由图4a 可知，整个快速冲击响应过程持续了约5 ps，不同模型冲击加载时间与模型加载方向长度有关。快速冲击响应后发生缓慢化学分解过程，表现在不同模型的温度均呈现不同程度的升高。CL-20、HMX 和CL-20/HMX 共晶的温度增长速率大小顺序为：CL-20＞CL-20/HMX（100）＞HMX＞CL-20/HMX（111）。由 图4b 可 知，CL-20、HMX 和CL-20/HMX 共 晶 中CL-20 分子的分解速率大小顺序为CL-20＞CL-20/HMX（100）＞CL-20/HMX（111），HMX 分 子 则 为CL-20/HMX（100）＞HMX＞CL-20/HMX（111）。因此单纯CL-20 的分解速度最快，而CL-20/HMX 和HMX的分解速度相当，CL-20/HMX（111）模型的分解速度甚至比单纯HMX 的还慢，可以初步解释CL-20/HMX共晶的冲击波感度和HMX 更接近的实验现象。与［111］晶向相比，沿［100］晶向加载时CL-20/HMX 共晶分解更快，表明出明显的冲击各向异性。

另外，当快速冲击加载结束时模型会达到最大压缩点，此时若只施加周期性边界条件继续开展分子动力学模拟，将因为末端产生的反向稀疏波而使得体系变得松散，因此在最大压缩点时额外施加一静止活塞［14］，可以有效阻止稀疏波进入波阵面后区域，从而能够进一步分析结构塑性变形和化学反应过程。

为研究CL-20/HMX 共晶冲击各向异性产生机理，统计分析了最大压缩点时CL-20/HMX（100）和CL-20/HMX（111）模型中完整的CL-20 和HMX 分子的分布，计算结果如图5 所示。由分子堆积过程可知，每两层CL-20 分子和一层HMX 分子交替堆叠形成共晶，HMX 和CL-20 之间因为氢键作用而较稳定的组装在一起。由图5a 中沿［100］晶向冲击加载时，作用力将垂直于CL-20 和HMX 分子层而逐层传递，而由图5b 可知，沿［111］晶向冲击加载时，CL-20 分子层和HMX 分子层之间产生滑移。与［100］晶向的这种显著差异，直观表现为即使在［111］晶向冲击加载端也有较多保持完整的CL-20 和HMX 分子。结合图4a 可知，与沿［100］晶向冲击加载相比，通过不同活性分子的交替排布和一定量的滑移等能够有效分散冲击加载能量，致使CL-20/HMX（111）模型的响应速度变慢，表现出明显的冲击各向异性。

图5 沿［100］（a）和［111］（b）冲击加载且处于最大压缩状态的CL-20/HMX 共晶结构（红色表示完整的CL-20 分子，绿色表示完整的HMX 分子，其余属于碎片分子）Fig.5 Structures of CL-20/HMX models at the maximum moment along［100］（a）and［111］（b）shock directions（Color ranks：red-intact CL-20 molecules，green- intact HMX molecules，others-fragments）

为进一步分析沿不同晶向冲击CL-20/HMX 共晶模型时能否发生滑移，计算了沿x 向和z 向的局部位移分布（如图6）和z 向的局部速度分布（如图7）。图6a和图6b 中CL-20/HMX（100）模型处于3.2 ps，波阵面位于20 nm，而图6c 和图6d 中CL-20/HMX（100）处于3.8 ps，波阵面位于17.5 nm，两模型的压缩量均为最大压缩量的60%。在CL-20/HMX（111）模型中，z向起始端属于固定壁位置，沿z 向的2.5～17.5 nm 范围属于波后受压区域，17.5～32.5 nm 属于未压缩的CL-20/HMX 共晶区域，32.5～40 nm 范围属于空隙区。由图6a 可知，沿［100］晶向加载时，受压缩区沿z 向位移方向和冲击波传播方向一致，而图6b 中沿x 向位移方向并无明显规律。由图6c 和图6d 可知，沿［111］晶向加载时，受压缩区沿z 向和x 向位移方向均有正有负，而且位移数值也比沿［100］晶向的大。相似规律也反映在图7 中。这些结果表明，沿［111］晶向冲击加载后，受压缩区的确发生了明显的滑移，最终延缓了反应进程，而沿［100］晶向冲击加载时滑移并不显著。

图7 沿不同晶向冲击加载的CL-20/HMX（100）和CL-20/HMX（111）模型的局部速度变化规律Fig.7 Variations of local velocities along different shock directions in CL-20/HMX（100）and CL-20/HMX（111）models

3.3 孔洞缺陷对冲击过程的影响

在分析无缺陷CL-20/HMX 模型的基础上，进一步考虑孔洞缺陷对冲击过程中力-热耦合过程影响。Φ20 nm 孔洞存在下温度、压力以及初始物质分解量随时间的变化关系分别如图8 和图9 所示。由图4a、图8a、图4b、图9 可知，沿［100］和［111］冲击加载无缺陷CL-20/HMX 共晶时，温度和不同分子分解量相对变化规律一致，即沿［111］冲击加载时温度（主要和粒子质点速度相关）降低最多，压力增加缓慢，CL-20 和HMX 分子分解速率最低，表明不同规模模型对计算结果的影响较小。然而，当存在Φ20 nm 孔洞缺陷时，伴随孔洞的压缩，体系温度出现先增后减的“鼓包”现象，体系压力则出现先减后增，对CL-20 和HMX 分子的分解产生一定的影响。下面结合孔洞压缩过程瞬态图做进一步分析。

沿［100］冲击加载含Φ20 nm 孔洞的CL-20/HMX共晶时，孔洞受压缩初期，模型中各物理量分布如图10 所示。沿z 轴正向，含缺陷CL-20/HMX 共晶分别由固定壁、含孔洞的CL-20/HMX 共晶受压区和运动区、空隙区等四部分组成。由图10a 可知，波阵面在碰到孔洞时被隔断，然而以冲击波传播方向为参考，孔洞上游壁面和下游壁面附近区域温度均较高（图10b），这是由于这两部分局部区域对应的粒子质点速度和位移均因孔洞存在而增加了自由度，受周围分子的约束减弱而动能增加的缘故。

图8 沿［100］和［111］冲击加载含Φ20 nm 缺陷CL-20/HMX 共晶模型温度和压力随时间的变化规律Fig.8 Evolution of temperature and pressure along the［100］and［111］shocks in CL-20/HMX cocrystal with 20 nm diameter void

图9 不同CL-20/HMX 模型中未分解CL-20 和HMX 分子数的变化规律Fig.9 Evolution of percentage of undecomposed CL-20 and HMX molecules in different CL-20/HMX models

为详细分析冲击加载CL-20/HMX［100］模型的结构变化，根据孔洞变形阶段，计算了几组典型时刻的密度、温度和压力的瞬态分布，结果如图11 所示。选取的典型时刻已标注于图8 中，即2.86 ps 对应于图10中时刻，5.44 ps 对应于孔洞临近闭合时刻，6.62 ps 对应于近球形孔洞塌缩“热点”产生的应力波正追赶冲击波的时刻，8.30 ps 对应于孔洞塌缩热点产生的应力波和之前冲击波复合形成的加强型冲击波继续向前传播的某一时刻。由图8～图11a 所展示的孔洞塌缩过程和后效分析可知，孔洞塌缩过程导致孔洞上游壁面和下游壁面获得较高的质点速度，两者是相对运动，孔洞塌缩的过程产生较高的温度、降低局部压力、减弱CL-20 和HMX 分子的分解。在图11b 中，孔洞闭合后，以孔洞闭合位置为中心产生新热点，此时CL-20 和HMX 分子的分解速率明显提高，形成了更高的温度和压力区并向周围扩散。在图11c 中，当孔洞塌缩热点产生的应力波和之前冲击波复合时形成了加强型冲击波，继续向前传播。从新热点到两波平齐，含Φ20 nm孔洞的CL-20 模型需要1.61 ps，而CL-20/HMX 模型则需要1.92 ps。

图10 沿［100］方向冲击加载含Φ20 nm 缺陷CL-20/HMX 共晶时孔洞受压缩初期各物理量的瞬态变化Fig.10 Evolution of various parameters of CL-20/HMX with 20 nm diameter void along the［100］shock direction at initial stage，which includes（a）density；（b）temperature；（c）pressure；（d）displacement

由孔洞塌缩形成热点过程可知，孔洞缺陷首先造成受力结构变化，以2 km·s-1质点速度冲击下Φ20 nm 孔洞周围结构没有产生射流现象，而是粘塑性收缩过程，然后孔洞闭合产生的高温高压条件和结构上的粘塑性变形均有效促使CL-20 和HMX 分子发生快速的化学分解，孔洞塌缩新热点的形成反过来进一步增强了冲击加载过程。

4 结论

采用反应分子动力学方法，探讨了有无缺陷CL-20/HMX 共晶在冲击加载下发生的力-热结构变化和随后的初始化学反应过程。结果表明：

（1）ReaxFF-lg 和ReaxFF-2018 均能够较好的描述CL-20、HMX 以及CL-20/HMX 共晶的晶格参数、密度。与ReaxFF-lg 相比，利用ReaxFF-2018 计算的CL-20/HMX 共晶高压物理特性值更接近于DFT-D 计算值。

（2）利用动量镜原理产生平面反射波后，在最大压缩点增加一静止活塞，可有效阻止稀疏波进入波后区域，有利于分析结构塑性变形情况或化学反应过程，不仅可实现高速冲击加载过程模拟，也能够实现低速撞击过程模拟，具有较大应用潜力。

（3）冲击加载无缺陷CL-20、HMX 和CL-20/HMX共晶时，CL-20 比HMX 的分解速度快，CL-20/HMX 共晶的分解速度接近HMX。和［111］晶向相比，沿［100］晶向冲击加载CL-20/HMX 共晶分解更快，这与CL-20 和HMX 分子层的交替排布以及能否产生滑移等因素有关。

（4）以2 km·s-1质点速度冲击下Φ20 nm 孔洞周围结构没有产生射流现象，而是粘塑性收缩过程。孔洞塌缩形成的高温高压条件和结构上的粘塑性变形有效促使CL-20 和HMX 分子发生快速的分解，孔洞塌缩新热点的形成进一步增强了冲击加载过程。