分子动力学法研究掺杂缺陷对HMX/NQ 共晶炸药性能的影响

杭贵云，余文力，王 涛，王金涛，苗 爽

(火箭军工程大学，陕西 西安 710025)

引 言

HMX是一种综合性能优异的硝胺类高能炸药，但因其感度较高，导致其应用受到限制。硝基胍(NQ)是一种钝感炸药，安全性较好，但威力与能量密度较低。丁雄[1]研究了HMX/NQ共晶炸药的性能，结果表明，共晶炸药的感度低于HMX，且共晶具有较高的能量密度。因此，HMX/NQ共晶炸药有望成为一种新型的高能炸药。

HMX通常采用改进乙酸酐法制备，在产物中会存在杂质成分RDX。因此，在HMX/NQ共晶炸药制备过程中，RDX也会进入其中，导致晶体中存在掺杂缺陷，可能会影响炸药的性能，如稳定性、感度、能量特性和力学性能等，从而进一步影响武器弹药的安全性与作战性能。因此，研究晶体缺陷对炸药性能的影响具有重要意义。LaBarbera等[2]研究了缺陷对含能材料“热点”形成的影响；Luscher等[3]研究了晶体缺陷对α-RDX感度的影响；Xue等[4]研究了位错缺陷对RDX冲击感度的影响；彭亚晶等[5]研究了空位缺陷对RDX性能的影响。以往研究的晶体缺陷类型主要侧重于空位与位错，而研究掺杂缺陷的报道相对较少。

本研究分别建立了“完美”型与含有不同掺杂率缺陷的HMX/NQ共晶炸药模型。在Materials Studio 7.0(简称MS)软件中[6]，采用分子动力学方法预测了模型的稳定性、感度、能量特性和力学性能，综合评估了晶体缺陷对炸药性能的影响，以期为炸药的性能评估提供理论参考。

1 计算模型与计算方法

1.1 单个分子的建立

HMX/NQ共晶炸药由HMX与NQ组成，摩尔比为1∶1。在MS软件中，分别建立HMX与NQ的单个分子模型，如图1所示。

图1 HMX与NQ的分子模型Fig.1 Molecule models of HMX and NQ

1.2 HMX/NQ共晶炸药原始模型的建立

本研究中使用的HMX/NQ共晶炸药结构为文献[1]给出的理论预测结构。HMX/NQ共晶炸药属于单斜晶系，空间群为P21，晶格参数为a=2.264nm，b=0.420nm，c=0.798nm，α=90.00°，β=103.00°，γ=90.00°，单个晶胞中包含2个HMX与2个NQ分子[1]。HMX/NQ共晶炸药的单个晶胞模型如图2(a)所示，而后将单个晶胞模型扩展为36(3×4×3)的超晶胞模型，其中包含72个HMX分子与72个NQ分子，一共144个分子，如图2(b)所示。为了便于与含有晶体缺陷的模型进行比较，将“完美”型晶体模型标记为Model-1。

1.3 HMX/NQ共晶炸药缺陷模型的建立

用2个RDX分子替换“完美”型晶体中的2个HMX分子(图3(a)中标记为黄色)，得到掺杂率为1.39%的缺陷晶体模型，如图3(b)所示，将所得缺陷模型标记为Model-2。

图3 “完美”型晶体模型(Model-1)与掺杂率为1.39%的 缺陷晶体模型(Model-2)Fig.3 Perfect crystal model(Model-1) and defective crystal model with adulteration ratio of 1.39%(Model-2)

在Model-2的基础上，继续用2个RDX分子取代模型中的2个HMX分子，得到掺杂率为2.78%的缺陷晶体模型，标记为Model-3。

采用同样的方法，在Model-3的基础上，继续用4个RDX分子替换4个HMX分子，得到掺杂率为5.56%的缺陷晶体模型，标记为Model-4。

类似地，分别用12个、16个、24个RDX分子替换初始模型中同等数量的HMX分子，得到掺杂率分别为8.34%、11.11%和16.67%的缺陷晶体模型，分别标记为Model-5、Model-6和Model-7。

1.4 计算条件设置

分别对HMX/NQ共晶炸药的“完美”型晶体模型和含有掺杂缺陷的晶体模型进行能量最小化处理，对其结构进行优化，而后进行分子动力学计算，其中温度设置为295K，压强设置为0.0001GPa，选择NPT系综与COMPASS力场[7-8]。选择COMPASS力场是因为该力场采用从头算法，适用于凝聚相和不同类型物质相互作用的研究，能在较大范围内对处于孤立体系和凝聚态体系的多种物质的性能进行准确预测。温度采用Andersen控温方法[9]，压强采用Parrinello控压方法[10]，范德华力(vdW)的计算采用Atom-based方法，静电作用的计算采用Ewald方法，时间步长为1fs，总模拟步数为2×105步，其中前105步用于热力学平衡计算，后105步用于统计分析。模拟过程中，每103fs保存一次轨迹，共得100帧轨迹文件。

2 结果与讨论

2.1 平衡判别与平衡体系

在提取数值模拟计算结果时，需要让混合体系达到热力学平衡状态，此时必须同时满足温度平衡和能量平衡。通常认为当温度与能量波动范围为5%～10%时，体系已经达到平衡状态。以Model-4为例，模拟过程中混合体系的温度变化曲线如图4(a)所示，能量变化曲线如图4(b)所示。

从图4可以看出，模拟计算初期，温度和能量均有所上升，且波动幅度较大。随着时间的推移，温度和能量的波动幅度逐渐减小，最终温度波动幅度约为±15K，能量波动幅度约为±5%，偏差相对较小，表明混合体系已达到热力学平衡状态。对于其他的晶体模型，分子动力学计算时，均以温度平衡和能量平衡来判别混合体系是否达到平衡状态。

图4 混合体系的温度和能量随时间的变化曲线Fig.4 Change curves of temperature and energy of mixed system with time

2.2 结合能

对于HMX/NQ共晶及其掺杂缺陷晶体，结合能(Eb)的计算公式如下：

Eb=-Einter=-[Etotal-(ENQ+Eother)]

(1)

式中：Eb为结合能；Einter为分子之间的相互作用力；Etotal为混合体系达到平衡状态时对应的总能量；ENQ为去掉体系中的其他组分后，NQ分子对应的总能量；Eother为去掉NQ分子后，体系中的HMX与RDX分子对应的总能量。上述单位均为kJ/mol。

计算得到Model-1～Model-7的结合能分别为335.7、331.4、320.6、318.3、315.5、310.8、298.6kJ/mol。从以上数据可以看出，“完美”型晶体(Model-1)的结合能最大，表明HMX与NQ分子间的相互作用力最强，炸药的稳定性最好。缺陷晶体的结合能均有不同程度减小，且随掺杂率的增加，结合能逐渐减小。当掺杂率为1.39%(Model-2)时，结合能为331.4kJ/mol；当掺杂率为16.67%(Model-7)时，结合能为298.6kJ/mol，结合能减小幅度为1.28%～11.05%。结合能减小，说明分子之间相互作用力的强度减弱，炸药的稳定性降低。结合能减小的原因可能是由于晶体缺陷的影响，炸药的结构遭到破坏，分子的排列方式发生变化，导致分子间的相互作用力减弱，结合能减小。

2.3 感度

依据“热点”理论[11]与“引发键”思想[12]，同时参考以往的研究工作[13-16]，本研究选用引发键键长、引发键键连双原子作用能和内聚能密度作为判据，从而预测炸药的感度大小。

2.3.1 引发键键长

在HMX/NQ共晶及其缺陷晶体中，HMX所占的比例最大，对外界的刺激最敏感。当受到外界刺激时，HMX更容易发生化学反应，使共晶炸药发生分解或爆炸。HMX的引发键为分子中N-NO2基团上的N—N键[12,17-18]，因此选择HMX分子中的N—NO2键作为引发键来预测体系的感度。

以Model-2为例，图5给出了平衡状态时，体系中引发键(N—NO2键)的键长分布情况。表1给出了平衡状态下，不同模型中引发键的最可几键长(Lprob)、平均键长(Lave)和最大键长(Lmax)。

图5 引发键的键长分布(Model-2)Fig.5 Trigger bond length distribution of Model-2

ModelLprob/nmLave/nmLmax/nm10.13830.13820.153220.13820.13830.153930.13840.13840.154540.13830.13830.155950.13840.13830.157860.13830.13840.159470.13850.13840.1613

从图5可以看出，当混合体系达到平衡状态时，引发键(N—NO2键)的键长分布呈近似对称的高斯分布。从表1可以看出，对于不同的模型，最可几键长与平均键长近似相等，并且变化范围很小，表明晶体缺陷对最可几键长与平均键长的影响很小，而最大键长变化很明显，且不同模型之间的差异较大。对于“完美”型晶体模型(Model-1)，最大键长为0.1532nm，对于掺杂率分别为1.39%(Model-2)、16.67%(Model-7)的晶体模型，最大键长分别为0.1539、0.1613nm，增大幅度为0.46%～5.29%。最大键长增大，说明引发键的强度减小，预示含能材料的感度升高，安全性降低。此外，从表1中还可以看出，随着掺杂率的增加，最大键长逐渐增大，表明炸药的感度随着缺陷数量的增加而逐渐升高。

2.3.2 引发键键连双原子作用能

引发键键连双原子作用能(EN-N)的计算公式如下[13-16]：

(2)

式中：E1为混合体系达到热力学平衡状态时对应的总能量(kJ/mol)；E2为固定HMX中所有的N原子后，混合体系达到平衡状态时的总能量(kJ/mol)；n为体系中HMX分子中包含的N—NO2键的数量。

通过计算，得到Model-1～Model-7的键连双原子作用能分别为149.75、148.81、146.52、142.04、139.18、132.5、123.94kJ/mol。从以上数据可以看出，“完美”型晶体(Model-1)的键连双原子作用能最大，而缺陷晶体的键连双原子作用能均有不同程度的减小。其中，当掺杂率为1.39%(Model-2)时，键连双原子作用能为148.81kJ/mol；当掺杂率为16.67%(Model-7)时，键连双原子作用能为123.94kJ/mol；键连双原子作用能减小的幅度为0.63%～17.24%。键连双原子作用能减小，说明引发键的强度减弱，键断裂时需要的能量降低，即含能材料的感度升高、安全性降低。此外，键连双原子作用能的数据也表明，随着掺杂率的增加，即缺陷数量的增加，键连双原子作用能逐渐减小，预示炸药的感度逐渐升高，安全性逐渐减弱。因此，晶体缺陷会对炸药的安全性产生不利影响。在缺陷晶体中，键连双原子作用能减小的原因可能是由于晶体结构与分子排列方式发生改变，NQ的钝感效应降低，HMX的活性增强，从而使引发键键能减小、感度升高。

2.3.3 内聚能密度

内聚能密度(CED)的计算公式如下：

(3)

式中：CED为内聚能密度，kJ/cm3；HV为摩尔蒸发热，kJ/mol；RT为物质气化时所作的膨胀功，kJ/mol；Vm为摩尔体积，cm3/mol。

通过分子动力学模拟，得到不同模型的内聚能密度、范德华力与静电力，结果如表2所示。

从表2可以看出，“完美”型晶体的内聚能密度、范德华力与静电力最大。在缺陷晶体中，内聚能密度的最大值与最小值分别为0.832和0.748kJ/cm3；与“完美”型晶体相比，内聚能密度减小幅度为0.83%～10.85%。内聚能密度减小，表明炸药由凝聚态变为气态时吸收的能量减小，预示炸药的感度升高、安全性降低。此外，随着掺杂率的增加，内聚能密度逐渐减小，表明炸药的感度随着晶体缺陷的增加而逐渐升高。内聚能密度减小的原因可能是晶体结构遭到破坏，分子排列方式发生改变，从而导致分子间的范德华力与静电力作用强度减弱，内聚能密度减小、感度升高。

表2 不同模型的内聚能密度、范德华力与静电力

2.4 爆轰性能

本研究采用修正氮当量法[19]计算炸药的爆轰参数并预测其能量特性。

对于由C-H-O-N 4种元素组成、分子式为CaHbOcNd的炸药，氧平衡系数(OB)的计算公式如下：

(4)

式中：a、b、c分别为炸药分子中包含的C原子、H原子与O原子的数目；Mr为炸药的相对分子质量。

对于混合炸药，氧平衡系数的计算公式如下：

OB=∑wiOBi

(5)

式中：wi为混合炸药中第i种组分所占的质量分数，%；OBi为第i种组分对应的氧平衡系数。

爆轰参数(D、p)的计算公式如下：

(6)

式中：D为炸药的爆速，m/s；p为炸药的爆压，GPa；ρ为炸药的密度，g/cm3；∑Nch为炸药的修正氮当量；pi为1 mol炸药爆炸时生成第i种爆轰产物的摩尔数；Npi为第i种爆轰产物的氮当量系数；BK为炸药分子中第K种化学键出现的次数；NBK为炸药分子中第K种化学键的氮当量系数；Gj为炸药分子中第j种基团出现的次数；NGj为炸药分子中第j种基团的氮当量系数。公式(6)中各变量的计算方法详见参考文献[20]。

根据修正氮当量理论，通过计算得到不同模型的爆轰参数，结果如表3所示，其中炸药的密度可以直接从分子动力学结果中提取得到。

表3 不同模型的密度与爆轰参数

从表3可以看出，在不同的模型中，“完美”型晶体的密度、爆速与爆压最大，分别为1.798g/cm3、8664m/s和34.00GPa，表明炸药的能量密度最高。在缺陷晶体中，由于RDX的能量密度小于HMX，因此缺陷晶体的氧平衡系数、密度与爆轰参数均呈现出逐渐减小的变化趋势且掺杂率越高，密度与爆轰参数越小。与“完美”型晶体相比，缺陷晶体的密度、爆速与爆压减小幅度分别为0.89%～7.06%、0.68%～5.41%、1.85%～14.18%。密度与爆轰参数减小，表明炸药的威力减小，能量密度降低，因此晶体缺陷会对炸药的能量特性产生不利影响。

2.5 力学性能

力学参数可由弹性系数求得，满足广义胡克定律，其表达式为[21-22]：

σi=Cijεj

(7)

式中：σ为应力；ε为应变；Cij(i，j=1，2，…，6)为弹性系数，满足Cij=Cji，因此独立的弹性常数只有21个，对于完全的各向同性体，独立的弹性常数只有2个(C11，C22)。

体积模量(K)与剪切模量(G)的计算公式如下：

KR=[S11+S22+S33+2(S12+S23+S31)]-1

(8)

GR= 15[4(S11+S22+S33)-4(S12+S23+S31)+

3(S44+S55+S66)]-1

(9)

式中：下标R表示Reuss平均。柔量系数矩阵S=[Sij]等于弹性系数矩阵C的逆矩阵，即S=C-1。

力学参数之间存在如下关系：

E=2G(1+γ)=3K(1-2γ)

(10)

基于上述公式，可以求得拉伸模量(E)与泊松比(γ)的表达式：

(11)

(12)

通过计算，得到不同模型的弹性系数与力学性能参数，结果如表4与图6所示。

表4 HMX/NQ共晶“完美”型模型与缺陷模型的力学性能

图6 不同模型的力学参数Fig.6 Mechanical parameters of different models

从表4与图6可看出，“完美”型模型(Model-1)的拉伸模量、体积模量与剪切模量的值最大，分别为13.389、8.225和5.490GPa，而柯西压(C12-C44)的值最小(0.209GPa)，表明炸药的刚性与硬度最强，而延展性较差。由于晶体缺陷的影响，E、K、G的值减小，而柯西压增大。在缺陷晶体中，Model-2(掺杂率为1.39%)的E、K、G最大，分别为13.283、8.139和5.408GPa，柯西压为0.317GPa；而Model-7(掺杂率为16.67%)对应模量最小，分别为9.021、5.652和3.655GPa，柯西压为1.996GPa。与“完美”型模型相比，E、K、G减小的幅度分别为0.106～4.368GPa、0.086～2.573GPa和0.082～1.835GPa，而柯西压增大幅度为0.108～1.787GPa。模量减小，表明炸药的刚性与硬度降低；柯西压增大，表明炸药的延展性增强。在缺陷晶体中，模量减小的原因可能是由于晶体结构遭到破坏，分子排列方式发生了变化，从而使晶体的对称性破坏，抵抗变形的能力减弱，在外界作用下，晶体更容易产生形变。

3 结 论

(1)与“完美”型晶体相比，缺陷晶体的结合能减小，分子之间的相互作用力减弱，炸药的稳定性降低。随着缺陷数量的增加，结合能逐渐减小，炸药的稳定性逐渐减弱。

(2)缺陷晶体的引发键键长增大，而键连双原子作用能与内聚能密度减小，表明炸药的感度升高，安全性降低。

(3)由于RDX的能量密度小于HMX，因此缺陷晶体的密度与爆轰参数减小，威力与能量密度降低。随着晶体缺陷数量的增加，炸药的能量特性逐渐减弱。

(4)在缺陷晶体中，由于晶体结构遭到破坏，炸药的拉伸模量、体积模量与剪切模量减小，柯西压增大，表明炸药的刚性与硬度降低、延展性增强。在外界作用下，炸药更容易产生形变。