1 引言

奥克托今（HMX）是目前军工领域使用最广泛、综合性能最好的二代炸药之一，通常采用醋酐法及其改进工艺进行合成制备^[1-2]。在合成制备的过程中，难免会有残留的黑索今（RDX）杂质掺杂在HMX晶体中。

杂质的存在会影响晶体生长速率，使晶体的结晶形貌发生变化或导致晶体产生内部缺陷。这不仅会影响晶体的热性能，还会影响晶体的分子弛豫，同时缺陷也会使能量局域化而形成“热点”，进而影响感度和力学性能^[3-4]。徐姣等^[5]研究了RDX和太安（PETN）掺杂碳黑和碳纳米管后，其激光起爆感度的变化，结果表明掺杂后的RDX和PETN激光感度明显升高，且在掺杂质量分数为1%时感度最高。徐金江等^[6]研究了金属离子杂质对HMX结晶特性和热性能的影响。结果表明，掺杂金属离子会改变HMX的结晶面貌，产生微观缺陷，改变晶体品质，进而影响HMX的热性能和安全性，其中Zn²⁺的影响程度最大。李泽生等^[7]将六硝基六氮杂异伍兹烷（CL-20）中的C原子替换成N原子设计了6种CL-20的衍生物，并利用密度泛函理论综合研究了其分子结构、生成热、爆轰性能、热稳定性和感度。结果表明，不同数量N原子的替代掺杂对CL-20的性能影响有优有劣，可以通过该方法改善CL-20的爆轰性能。Bhattacharia S K等^[8-9]研究了PETN晶体被卟啉及其同系物掺杂后的晶体形貌以及升华和热解特性。目前，HMX的同系物RDX杂质及其掺杂率对HMX性能影响的研究未见报道。

为此，本研究采用替代掺杂的形式建立了掺杂有不同数量RDX分子的β-HMX超晶胞，采用分子动力学方法和理论计算得到了不同模型的键连双原子作用能、内聚能密度、溶度参数、爆轰参数和力学参数，分析比较了不同模型的感度、与氟橡胶（F₂₃₁₁）的相容性、爆轰性能和力学性能。

2 计算方法 2.1 模型搭建 2.1.1 HMX晶胞模型搭建

β-HMX（为叙述方便，以下称HMX）属于单斜晶系，空间群为P2₁/C，晶格参数为a=6.54 Å，b=11.05 Å，c=8.70 Å，α=γ=90.00°，β=124.30°，每个单晶胞中包含两个HMX分子^[10]。HMX的单晶胞模型如图 1a所示。在建立HMX超晶胞之前，首先采用Smart算法对HMX单晶胞进行了能量和几何优化，优化结构如图 1b所示，而后将晶格参数分别沿X、Y、Z三个轴向延伸4、3、4倍，建成了HMX超晶胞模型，共96个HMX分子，如图 2所示。

2.1.2 HMX晶胞掺杂模型搭建

在HMX超晶胞模型（为叙述方便，认为此模型掺杂率为0）的基础上，通过随机移除4、8、12和16个HMX分子，并在相应位置添加相同数量的RDX分子建立掺杂率（质量分数）分别为4.17%、8.33%、12.50%和16.67%的四种HMX缺陷模型。特以掺杂率为16.67%的HMX缺陷模型（图 3）为例进行展示。其中，标黄的分子为RDX。

2.1.3 F₂₃₁₁晶胞模型搭建

HMX单独成型能力差，其在使用过程中可通过添加氟橡胶等高分子材料制成高聚物粘结炸药，因此有必要探究RDX的掺杂对HMX与氟橡胶相容性的影响。选择的氟橡胶（F₂₃₁₁）是由偏二氟乙烯和三氟氯乙烯按摩尔比1:1组成的。利用Material Studio（MS）^[11]软件中的Visualizer模块搭建了聚合度为4的F₂₃₁₁分子链，端基用H和F饱和，如图 4a所示，而后利用Amorphous Cell模块构建了含用6个F₂₃₁₁分子链的无定形晶胞，密度设置为1.85 g·cm^-3，如图 4b所示。

2.2 分子动力学模拟

Compass ^[12]力场是一种从头算力场，其中多数力场参数是采用量子力学计算得到的，对单个分子或凝聚相材料性能的模拟和预测具有较高的准确性，尤其是其对HMX超晶胞和RDX超晶胞^[13-14]均进行过成功的模拟计算。因此，本研究选择Compass力场进行分子动力学（MD）模拟。

在HMX超晶胞及其掺杂模型进行MD模拟之前，首先采用Smart算法对初始模型进行结构优化，优化收敛精度设置为4.186×10^-3 kJ·mol^-1·Å^-1。当优化结果显示的最大导数低于0.05时认为优化模型实现了能量极小化，内应力已被平衡。而后，又对优化后的模型进行了模拟退火，从中选取了能量最小的退火模型进行MD模拟。MD模拟选择在恒压恒焓（NPT）系综下进行，模拟温度设置为298 K，精度设置为“fine”。模拟过程中分别采用Anderson^[15]和Berendsen^[16]方法对温度和压力进行控制，范德华（vdW）和静电作用（Coulomb）的加和方法分别采用Atom-based^[17]和Ewald^[18]，截断半径取1.55 nm，并对截断尾部进行矫正。原子运动的初始速度由Boltzmann分布确定，并采用Verlet方法求解牛顿运动方程的积分。对所建模型进行了20 ps的MD模拟，前10 ps对模型体系进行平衡，后10 ps用于统计分析确定能量、力学参数和其它参数。每0.1 ps取样一次，共获得100帧轨迹。对原子运动轨迹的统计分析必须要基于模拟体系达到平衡才有意义。模拟体系平衡的标志是温度和能量达到平衡。通过温度及能量随时间变化的曲线确认HMX超晶胞及其掺杂模型已达到平衡状态。

3 结果讨论 3.1 RDX掺杂对HMX炸药感度的影响

肖继军等^[19-20]利用引发键键连双原子作用能和内聚能密度作为判据评估HMX和RDX的感度，所得结果与实验相符。因此，本研究也选用引发键键连双原子作用能和内聚能密度作为判据，评估RDX的掺杂对HMX炸药感度的影响。

3.1.1 键连双原子作用能

HMX和RDX的引发键都是N—NO₂基团中的N—N键，但是由于RDX在体系中占的质量分数较小，因此以HMX的键连双原子作用能作为整个体系的键连双原子作用能。

引发键键连双原子作用能（E_N—N）指的是引发键的键能，其定义为：

$ {E_{{\rm{N - N}}}} = \frac{{{E_1} - {E_2}}}{n} $

(1)

式中，E₁表示混合体系在平衡状态下系统的总能量，J·mol^-1；E₂表示在平衡状态下，固定HMX中所有的N原子后模型体系的总能量，J·mol^-1；n为体系中HMX分子中包含的N—N键的数量。

按照掺杂率由大到小的顺序，HMX超晶胞及其掺杂模型的键连双原子作用能（E_N—N）分别为175.88，166.35，157.80，148.57，139.52 kJ·mol^-1。由此可以看出，随着RDX掺杂率的增加，HMX超晶胞及其掺杂模型的E_N—N呈下降趋势，下降幅度为9.53~36.36 kJ·mol^-1，且与掺杂率基本呈线性。E_N—N减小，表明引发键束缚N原子的能力减弱，在受到外界激励的作用下，引发键更易发生解离，从而使炸药分解甚至爆炸。因此，随着掺杂率的增加，HMX超晶胞及其掺杂模型的感度升高，安全性能恶化。

3.1.2 内聚能密度

内聚能密度（CED）是指单位体积内1 mol凝聚态物质汽化为气态物质所需要的能量，是分子间相互作用的综合反映，是范德华力和静电力（Electrostatic）之和，本质上是一种非键力，反映的是体系中分子间相互作用的强弱，与体系的感度之间存在一定的关联。炸药的CED值越小，表明炸药受热发生相变所需的能量越小，间接表明同等条件下炸药对热越敏感。其计算公式为：

$ {\rm{CED}} = \frac{{\Delta {H_{\rm{V}}} - RT}}{{{V_{\rm{m}}}}} $

(2)

式中，ΔH_V为摩尔蒸发热，J·mol^-1；R为气体常数，8.314 J·mol^-1·K^-1，T为温度，K。V_m为摩尔体积，cm³·mol^-1。计算得到的HMX超晶胞及其掺杂模型的CED值及其分量范德华力和静电力的值见表 1。

由表 1可知，4种掺杂缺陷模型与HMX超晶胞模型的CED随着掺杂率的增加而减小，这表明随着RDX掺杂率的增加，HMX的热感度上升，安全性降低。

3.2 RDX掺杂对HMX与F₂₃₁₁相容性的影响

两种物质形成相容体系的热力学条件为：

$ \Delta {G_{\rm{M}}} = \Delta {H_{\rm{M}}} - T\Delta {S_{\rm{M}}} < 0 $

(3)

式中，ΔG_M为混合自由能，J·mol^-1；ΔH_M为混合热，J·mol^-1；ΔS_M为混合熵，J·K^-1·mol^-1。由于ΔH_M的值总是大于零，因此两物质是否相容，就需要看ΔS_M的贡献是否足够克服ΔH_M。然而，对于凝聚态物质和高分子材料的混合，熵的增加十分有限，所以，ΔG_M值的正负取决于ΔH_M的大小。

由于CED能够体现物质间的相互作用力，Hildebrand^[21]等引入溶度参数的概念，其定义为内聚能密度的平方根：

$ \delta = \sqrt {\frac{{\Delta E}}{{{V_{\rm{m}}}}}} = \sqrt {\frac{{\Delta {H_{\rm{V}}} - RT}}{{{V_{\rm{m}}}}}} $

(4)

式中，ΔE为体系的内聚能，J·mol^-1；V_m为摩尔体积，cm³·mol^-1；ΔH_V为摩尔蒸发热J·mol^-1，R为气体常数，8.314 J·mol^-1·K^-1，T为温度，K。高分子材料相互作用过程中的焓变与其溶度参数具有下列关系^[22]：

$ \frac{{\Delta {H_{\rm{M}}}}}{{{V_{\rm{m}}}}} = \left( {{\delta _1} - {\delta _2}} \right){\sigma _1}{\sigma _2} $

(5)

式中，σ₁、σ₂为各组分的体积分数。由（5）式可以看出，溶度参数差值（Δδ）可以用来评估组分间相容性，其越趋于0，组分之间的相容性越好。表 2列出了MD计算得到的HMX超晶胞及其掺杂模型和F₂₃₁₁无定形晶胞的溶度参数。

由表 2可以看出，随着掺杂RDX的增多，HMX超晶胞及其掺杂模型同F₂₃₁₁之间的溶度参数差值逐渐变小，即ΔH_M的值更趋于0。这表明随着RDX掺杂率的增加，HMX超晶胞及其掺杂模型同F₂₃₁₁之间的相容性逐渐增强。

3.3 RDX掺杂对HMX炸药爆轰性能的影响

通过修正氮当量法^[23]计算了炸药的爆速（D）和爆压（p），并根据质量加权法^[24]计算了爆热（Q），以此对其能量特性进行评估。

对于C_aH_bO_cN_d炸药，氧平衡系数（OB）计算公式^[25]如下：

$ {\rm{OB}} = \frac{{\left[ {c - \left( {2a + b/2} \right)} \right]}}{{{M_{\rm{r}}}}} \times 16 \times 100\% $

(6)

式中，a，b，c分别为炸药分子中C、H、O三种原子的数目，M_r为炸药的摩尔质量，g·mol^-1。

根据修正氮当量理论，爆速和爆压的计算公式如下：

$ \left\{ \begin{array}{l} D = \left( {690 + 1160\rho } \right)\Sigma {N_{{\rm{ch}}}}\\ p = 1.106{\left( {\rho \Sigma {N_{{\rm{ch}}}}} \right)^2} - 0.84\\ \Sigma {N_{{\rm{ch}}}} = \frac{{100}}{{{M_{\rm{r}}}}}\left( {{p_i}{N_{pi}} + \Sigma {B_K}{N_{BK}} + \Sigma {G_j}{N_{Gj}}} \right) \end{array} \right. $

(7)

式中，D为炸药的爆速，m·s^-1；p为炸药的爆压，GPa；ρ为炸药的密度，g·cm^-3；∑N_ch为炸药的修正氮当量；p_i为1 mol炸药生成第i种爆轰产物的摩尔数；N_pi为第i种爆轰产物的氮当量系数；B_K为炸药分子中第K种化学键出现的次数；N_BK为炸药分子中第K种化学键的氮当量系数；G_j为炸药分子中第j种基团出现的次数；N_Gj为炸药分子中第j种基团的氮当量系数。

关于修正氮当量法的具体计算过程以及p_i，N_pi，B_K，N_BK，G_j，N_Gj等参数的来源可以通过文献[25]获得。

根据质量加权法，混合炸药的计算公式如下：

$ Q = N{\omega _i}{Q_i} $

(8)

式中，ω_i是混合炸药中第i种组分的质量分数；Q_i是混合炸药中第i种组分的爆热，kJ·kg^-1。采用Meyer R等^[26]实验测得的HMX和RDX的爆热（6197 kJ·kg^-1和6322 kJ·kg^-1），对HMX超晶胞模型和掺杂模型的爆热进行预测。

计算得到不同模型的爆轰参数见表 3。

由表 3可见，随着RDX掺杂率的增加，HMX超晶胞及其掺杂模型的OB没有变化，主要是因为HMX和RDX是同系物，具有相同的最简分子式。HMX超晶胞及其掺杂模型的ρ、D和p随着掺杂率的增加而减小，减小幅度分别为1.12%~5.59%、0.84%~4.19%和2.27%~11.14%。对于这三个爆轰参数的变化，本质上是ρ的变化引起的，因为HMX和RDX的修正氮当量系数是一样的。HMX超晶胞及其掺杂模型的Q随着掺杂率的增加呈上升趋势，主要原因在于RDX的爆热大于HMX的爆热。爆热的增加幅度为0.06%~0.26%，变化较为微小，主要原因是HMX和RDX之间爆热的差值比较小，而且RDX掺杂的数量也较少。因此，可以得出RDX掺杂缺陷对HMX炸药爆轰性能的影响主要体现在降低其密度、爆速和爆压，提升其爆热。

3.4 力学性能

根据广义虎克定律^[27]，通过最小二乘法拟合弹性系数得出平均的拉伸应力应变，获得体积模量（K_b）和剪切模量（G_s）。虎克定律、体积模量和剪切模量的计算公式如下所示：

$ {\sigma _i} = {C_{ij}}{\varepsilon _j} $

(9)

$ {K_{{\rm{b}}\;{\rm{R}}}} = {\left[ {{S_{11}} + {S_{22}} + {S_{33}} + 2\left( {{S_{12}} + {S_{23}} + {S_{31}}} \right)} \right]^{ - 1}} $

(10)

$ {G_{{\rm{s}}\;{\rm{R}}}} = 15{\left[ {4\left( {{S_{11}} + {S_{22}} + {S_{33}}} \right) - 4\left( {{S_{12}} + {S_{23}} + {S_{31}}} \right) + 3\left( {{S_{44}} + {S_{55}} + {S_{66}}} \right)} \right]^{ - 1}} $

(11)

式中，σ为应力，Pa；ε为应变；K_b为体积模量，Pa；G_s为剪切模量，Pa；下标R为Reuss平均；C_ij（i，j=1，2，……，6）为弹性系数矩阵；S_ij（i，j=1，2，……，6）为柔量系数矩阵，等于C_ij的逆矩阵，即S=C^-1。

力学参数之间具有相互联系，关系式如下所示：

$ {E_{\rm{t}}} = 2{G_{\rm{s}}}\left( {1 + \mu } \right) = 3{K_{\rm{b}}}\left( {1 - 2\mu } \right) $

(12)

式中，E_t为弹性模量，Pa；μ是泊松比。

根据公式（12）推导获得弹性模量以及泊松比的计算公式如下所示：

$ {E_{\rm{t}}} = \frac{{9{G_{\rm{s}}}{K_{\rm{b}}}}}{{3{K_{\rm{b}}} + {G_{\rm{s}}}}} $

(13)

$ \mu = \frac{{3{K_{\rm{b}}} - 2{G_{\rm{s}}}}}{{2\left( {3{K_{\rm{b}}} + {G_{\rm{s}}}} \right)}} $

(14)

根据公式（9）~（14），求得HMX超晶胞及其掺杂模型的力学性能参数，结果如表 4所示。从表 4可以看出，与HMX超晶胞模型相比，掺杂缺陷模型的弹性模量（E_t）、体积模量（K_b）和剪切模量（G_s）均随掺杂缺陷率的增加而减小，减小幅度分别为1.04，1.63，2.57，3.63 GPa；0.58，1.03，1.51，1.73 GPa；0.42，0.64，1.02，1.45 GPa；柯西压（C₁₂-C₄₄），体积模量与剪切模量的比值（K_b/G_s）均随掺杂缺陷率的增加而增加，增加幅度分别为0.35，0.54，1.54和2.69 GPa以及0.11，0.15，0.30和0.64。这表明RDX杂质的存在使HMX的刚性、硬度和断裂强度性能减弱，柔性和延展性变好，在外界作用下，炸药更容易变形。

4 结论

采用分子动力学方法，探究了HMX超晶胞及其掺杂模型的感度、与F₂₃₁₁的相容性、爆轰性能以及力学性能，探讨了RDX掺杂缺陷对HMX炸药性能的影响情况。结果表明：

（1）掺杂缺陷导致HMX炸药的键连双原子作用能和内聚能密度减小，减小幅度分别为9.53~36.36 kJ·mol^-1，0.028~0.135 kJ·cm^-3，表明HMX炸药感度增大，安全性减弱。随着掺杂率的增加，感度逐渐增大，安全性逐渐减弱。

（2）随着掺杂缺陷浓度的增加，HMX炸药与F₂₃₁₁的溶解度参数的差值逐渐减小，减小幅度分别为0.51，0.89，1.78，2.32 J^1/2·cm^-3/2，表明HMX炸药与F₂₃₁₁的相容性增强。

（3）掺杂缺陷导致HMX炸药的密度、爆速和爆压分别减小1.12%~5.59%、0.84%~4.19%和2.27%~11.14%；爆热略微增加，增幅为0.06%~0.26%。但其对炸药的氧平衡没有影响。随着掺杂率的增加，受影响的爆轰参数的变化程度逐渐增大。

（4）掺杂缺陷使得HMX炸药的弹性模量、体积模量和剪切模量减小，减小幅度分别为1.04~3.63 GPa、0.58~1.73 GPa和0.42~1.45 GPa；柯西压和体积模量与剪切模量的比值增大，增幅分别为0.35~2.69 GPa和0.11~0.64。随着掺杂率的增加变化程度逐渐增加，表明炸药的抗变形能力减弱，柔性和延展性变好，力学性能整体变差。