PYX基高聚物粘结炸药界面相互作用及力学性能的分子动力学模拟

2022-05-11任圆圆刘思佳

原子与分子物理学报 2022年5期

陈芳, 任圆圆, 何磊, 刘思佳

(中北大学化学工程与技术学院, 太原 030051)

1 引言

高聚物粘结炸药(Polymer Bonded Explosives, 简称PBX), 以新型单质炸药为主体, 以高分子材料为粘合剂并添加助剂制成的, 具有耐热和易加工成型、安定和爆轰性能好等优势, 广泛用于国防和民用领域[1, 2]. 要研制出性能优异的配方就要选择与主体炸药结合性能好, 粘合强度大的高聚物粘结剂. 对于筛选优质粘结剂, 计算机模拟是一个不错的选择. 计算机模拟虽然不能完全代替实验, 但可以指导实验和验证某些理论假设. 分子动力学(Molecular Dynamics, MD)模拟方法是基于对原子或分子的排列和运动进行统计的原理, 可以在原子尺度上模拟高聚物粘结炸药的性质[3, 4]. 目前, 在高聚物粘结炸药领域研究较多的是以黑索金(RDX)[5], 奥克托今(HMX)[6, 7], 六硝基六氮杂异伍兹烷(CL-20)[8-11]等为主体炸药与粘结剂之间的相互作用以及PBX体系的结合能, 内聚能密度(CEDs)、力学性能等性质.

耐热炸药是指在高温环境中经过较长时间仍能可靠使用的一类炸药, 其熔点高, 蒸汽压低, 在高温条件下具有良好的安定性, 已经广泛应用于宇宙航行、空间探索、深井钻探、石油开采及其他特殊用途.为了满足更高的耐热性要求, 需要寻找具有更加耐高温性的炸药, 例如六硝基菧(HNS), 三氨基三硝基苯(TATB)[12-14]和2,6-双(苦氨基)-3,5-二硝基吡啶(PYX)等. PYX是目前世界上耐热性较好的单质炸药[15]. 但是, 在PYX单独使用时, 其机械感度和摩擦感度较高. 因此通过以PYX为主体, 添加少量高分子材料的粘结剂构成高聚物粘结炸药(PBX), 能够显著改善PYX的安全性能和使用性能. 粘结剂既保证了PBX复合材料的完整性, 使炸药具有塑性和钝感, 又影响了炸药的能量水平. 本文采用MD方法, 计算3种高聚物粘结剂(Estane 5703, EPDM, F2311)与PYX(011)晶面的相互作用以及PYX基PBXs复合体系的力学性能, 为基于PYX的高聚物粘结炸药的设计提供理论参考.

2 理论和计算方法

2.1 模型构建

从剑桥晶体数据中心(CCDC, Cambridge Crystallographic Data Centre)获得PYX晶胞结构(CCDC号为1429070)[16], 晶胞结构如图1所示. 选用COMPASS力场[17], 利用Forcite模块对初始的PYX晶体进行结构优化, 优化前后的晶胞参数及其误差见表1. 从表1可知, PYX优化后的结果与实验结果的最大误差为0.06(b值), 在合理的误差范围之内, 后续研究都选择COMPASS力场进行模拟研究. 利用Materials Studio软件构建Estane 5703, EPDM, F2311的分子模型, 如图2所示, 并对聚合物链进行几何优化, 获得能量最小的稳定构型.

图 1 PYX的分子结构及晶胞结构图(灰色:C,蓝色:N,红色: O,白色:H)Fig. 1 Molecular structure and and cell structure of PYX. Gray, carbon; blue, nitrogen; red, oxygen; white, hydrogen. color online

图2 Estane 5703、EPDM和F2311分子模型图Fig.2 Molecular models of Estane 5703, EPDM and F2311

然后建立高聚物与PYX晶面复合体系模型(PYX/Estane 5703、PYX/EPDM和PYX/F2311). 本文以Estane 5703/PYX为例进行说明, 如图3所示. 首先根据粘结剂所占比例(质量比约为5%)来确定所需的聚合物链数, 构建Estane 5703高聚物分子模型. 然后将PYX(011)面扩展为(2×3×2)的超晶胞. 最后将高聚物分子模型与PYX(011)超晶胞构建为PBXs模型. 从±x, ±y, ±z 6个方向压缩晶胞, 使得复合体系达到理论密度. 后续建模方式都采用此方法.

图3 Estane 5703/PYX(011)复合体系的构建过程Fig.3 The building process of Estane 5703/PYX(011) composite system

2.2 MD模拟细节

在MD模拟之前, 首先对模型的初始结构执行几何优化. 对于MD模拟, 首先采用恒温恒压(NPT)系综, Velocity-scaling控温方式和Berendsen控压方式[18], 选取压力为1 GPa, 进行100 ps的MD运算, 确保其模拟的密度曲线收敛并稳定; 然后利用上述所得的构型, 采取相同的NPT系综和控温控压方式, 更改压力为0.0001 GPa, 进行100 ps的MD运算, 并观察密度曲线收敛状态, 以获得稳定的模拟结构; 最后以上述模拟结构为初始构型, 采用NVT系综, 步长为1 fs, 温度设置为298 K, 选择Andersen控温方法[19], 采用Atom Based和Ewald 方法[20]分别计算范德华作用和静电作用, 进行300 ps动力学模拟研究, 最后提取轨迹文件的后100个结果进行分析.

3 结果与分析

3.1 结合能分析

PYX(011)晶面和3种聚合物分子之间的结合能(Ebind)根据以下公式计算：

Ebind=-Einter=-[Etotal-(Ecell+Epolymer)]

(1)

Ebind是聚合物与PYX(011)晶面的结合能;Einter是聚合物与PYX(011)晶面的相互作用能;Etotal是由平衡结构求得的PBX体系的单点能;Ecell为PBX平衡体系去掉高聚物后计算主炸药的单点能;Epolymer为PBX平衡体系去掉主炸药后计算所得高聚物的单点能. 最后通过perl脚本提取轨迹文件中最后100个结构的平均能量得到最终结果.

PBX 体系的结合能(Ebind)可以用来评价PYX与粘结剂体系的相容性, 其值越大, 粘结剂与晶胞之间的相互作用越强, 体系相容性越好. 图4表示常温下粘结剂含量为5%的不同粘结剂与PYX超晶胞之间的结合能, 从图中看到, 不同粘结剂与PYX超胞之间的结合能为正值, 说明体系能够稳定存在. 不同粘结剂与PYX(011)晶面之间的结合能值大小顺序为PYX/Estane 5703 > PYX/EPDM > PYX/F2311, 结果表明Estane5703与PYX(011)晶面之间相互作用最强.

表1 PYX优化前后的晶胞参数及误差

图4 PYX/Estane 5703、PYX/EPDM和PYX/F2311复合体系的结合能Fig.4 The binding energies of PYX/Estane 5703, PYX/EPDM and PYX/F2311 composite systems

3.2 内聚能密度(CED)分析

内聚能密度(CED)是范德华力(vdW)和静电力(electrostatic force)之和, 本质上是一种非键力. CED 反映的是PBX体系中分子间相互作用的强弱. 一般来说, 基团的极性越大, 分子间作用力越强, CED值越高. 通过 MD 模拟计算得到了不同 PYX 基 PBXs 模型的 CED 值以及CED值的分布. CED计算公式为：

CED=(Hv-RT)/Vm

(2)

其中Hv是摩尔蒸发热,RT是汽化过程中所做的膨胀功,Vm是摩尔体积.

基于MD模拟的轨迹文件,可以测得PYX/Estane 5703, PYX/EPDM和PYX/F2311体系的CED等各项数据. 表 2列出了体系内聚能密度的详细计算数值. 由表 2可见, 内聚能密度等于范德华力与静电力之和, 其中静电力对内聚能密度的贡献要大于范德华力, 说明体系中的分子间相互作用以静电力为主. 不同粘结剂与PYX之间的静电力大小顺序和内聚能密度大小顺序均为PYX/F2311 > PYX/Estane 5703 > PYX/EPDM.

表2 PYX/高聚物复合体系的内聚能密度及其分量

3.3 径向分布函数分析

径向分布函数是通过计算局部的空间分布来分析材料结构的一种分析方法. 径向分布函数(RDF)为某一指定原子附近出现另外一种原子的几率, 其表达式为:

g(r)=4πr2ρdr

(3)

式中, r为原子间距, ρ为原子数量密度.g(r)为参考原子指定的距离r处其他原子的概率. 它能反映系统中分子和原子的排列和相互作用. 从径向分布函数的峰值位置可以判断相互作用的类型, 而相互作用力的大小可以从峰的高度进行推断. 一般来说, 分子间短程相互作用包括氢键作用(< 3.1 Å)和范德华作用(3.1～5.0 Å), 分子间长程相互作用通常指静电作用(> 5.0 Å). 比较炸药分子与高聚物粘结剂在相同r处的g(r), 可以得到不同粘结剂与炸药分子之间的相互作用的强弱. 为了研究不同体系高聚物粘结剂与炸药分子之间的作用强度, 分别计算了高聚物粘结剂中的H原子与PYX中的O原子之间的径向分布函数.

图 5分别计算了PYX分子中的O原子与3种粘结剂分子中的H原子之间的径向分布函数. 由图 5可知, 3种体系中, 2.0～3.1 Å之间出现了明显的峰值, 说明它们之间存在着氢键相互作用力; 3.1～5.0 Å之间也有峰值出现, 说明它们之间存在着范德华力相互作用; 作用距离大于5.0 Å范围内也有明显的峰值出现, 说明它们之间存在静电相互作用.

图5 PYX中O原子与聚合物分子(Estane 5703, EPDM, F2311)中H原子之间的径向分布函数Fig.5 Radial distribution functions between O atoms in PYX and H atoms in polymer molecules (Estane 5703, EPDM, F2311)

将3个体系3个相互作用区间的峰点和峰值列于表3 , 由表3可知, 根据3个区间的最强峰强判断, 3个体系中氢键强弱顺序、范德华相互作用强弱顺序和静电相互作用强弱顺序均为PYX/F2311 > PYX/Estane5703 > PYX/EPDM, 和内聚能密度分析中的范德华力和静电力大小顺序比较的结果一致, 且3个体系中作用距离大于5.0 Å范围内的峰强都是最大, 说明O原子与H原子之间主要相互作用是静电相互作用, 与内聚能密度分析中的静电力为主要作用的结果一致.

表3 径向分布函数中三种相互作用区间的峰点和峰值

3.3 力学性能分析

在含能材料领域, 力学性能对材料的制造和应用起着非常重要的作用[4, 21]. 材料在外界刺激下产生内应力和应变, 均满足胡克定律, 弹性系数Cij反映了材料在不同位置的弹性效应. 弹性系数矩阵由广义胡克定律方程中相应分量的一阶偏导数得到, 其数目随着系统对称性的增加而减少. 对于各向同性材料, 它们只有两个独立的弹性系数C11和C12.因此, 我们将定义C12=λ和C11-C12=2μ, 并被称其为lame系数. 用最小二乘法拟合拉伸应力-应变曲线, 得到拉伸模量(或杨氏模量)的平均值, 并得到其它有效的各向同性力学性能, 如拉伸模量(E), 剪切模量(G)和体积模量(K).

5、K、G是评价材料刚度和断裂强度的关键指标. 它们的值越大, 则说明材料的刚度越高.图6为常温下纯PYX和添加3种粘结剂(粘结剂含量为5%时)后PBX体系的弹性模量图，表4 为具体数据. 由表4可见, 不同粘结剂对PYX晶体的刚性的影响不同, 其拉伸模量(E)由大到小排序为PYX > PYX/F2311 > PYX/EPDM > PYX/Estane 5703; 体积模量(K)是断裂强度的衡量标准, 其由大到小依次为PYX > PYX/F2311 > PYX/Estane 5703 > PYX/EPDM; 剪切模量(G)由大到小依次为PYX/F2311 > PYX > PYX/EPDM > PYX/Estane 5703. 以上分析表明, 加入不同种类的粘结剂后, PBX体系拉伸强度和断裂强度都得到了改善; 而除了F2311外, 加入EPDM和Estane 5703两种粘结剂后, 复合体系的抗剪切应变能力得到提高. 总之, 加入粘结剂后, PYX的力学性能得到了改善.