聚合物黏结炸药（polymer bonded explosive，PBX）是一种高度填充的颗粒复合材料，主要由含能单质炸药晶粒和聚合物黏结剂以及其他增塑剂、钝感剂等混合而成^[1-3]。PBX具有低感度、高爆速、良好的物理安定性等特点，广泛应用于火箭推进剂、弹头和采矿等军事和民用领域。相对于单质炸药晶粒，黏结剂的弹性模量相对较低且易于变形，在低速撞击、拉伸和压缩等外部刺激下，复杂的热-力-化学耦合反应可能沿着炸药晶体与黏结剂之间的界面产生，从而影响PBX的稳定性^[4-5]。因此，研究炸药晶体与黏结剂之间的界面行为对于PBX的设计和应用具有重要意义。

目前，人们对PBX的力学行为、结构变形和断裂机理开展了大量研究。例如：Yeager等^[6]通过中子反射仪观察到在HMX晶体与黏结剂Estane之间存在一个明显的界面混合层（厚度约6 nm），该混合层可能影响黏结剂与炸药之间的界面附着力；Rae等^[7]通过巴西实验研究了PBX的拉伸断裂行为，发现微裂纹通常在大颗粒的边界形成，并垂直于加载方向扩展，导致PBX断裂，在裂纹处黏结剂被拉成纤维细丝，这被认为会影响炸药的拉伸强度；Xiao等^[8]通过动态拉伸实验研究了RDX基PBX的拉伸性能，发现其拉伸行为表现出明显的高应变率依赖性，炸药颗粒与黏结剂之间的界面脱粘是主要的失效模式。然而，仅通过实验很难完全揭示PBX在拉伸加载下的力学响应，如应力局部化、微观结构演变等。

数值模拟方法是研究PBX力学行为的有效手段，近年来，人们针对含能晶体、聚合物黏结剂及损伤界面提出了不同的本构模型。Wang等^[9]考虑了由于黏结剂破裂和脱粘引起的速率相关损伤，为黏结剂与PBX界面建立了损伤黏弹性模型，并模拟了PBX9501的动态力学响应，发现PBX9501的拉伸强度与应变速率有关，且断裂方式为黏结剂脱粘。Dai等^[10]使用内聚力模型（cohesive zone modeling，CZM）研究了PBX9501在单轴准静态拉伸下的失效演化，发现在整个拉伸过程中，拉伸应变起源于较软的黏结剂，拉伸应变的持续累积导致黏结剂脱粘。分子动力学（molecular dynamics，MD）方法作为一种广泛使用的数值模拟技术，可以揭示界面局部结构演化和微观力学行为^[11-14]。Xiao等^[15-17]应用COMPASS力场计算了一系列HMX基PBX炸药的弹性力学性能，提出了多种黏结剂（Estane5703、PEG、PVDF、PCTFE等）对PBX结构和力学性能的影响。Wang等^[18]通过MD模拟研究了应变率对RDX基PBX的动态力学响应和失效机制的影响，发现PBX的初始损伤主要发生在黏结剂区域。Lv等^[19]使用MD模拟研究了温度和应变率对TATB-F2314界面力学性能及微观结构的影响，发现黏结剂F2314会扩散进入TATB层，导致TATB-F2314界面处出现“内嵌状”界面混合层，并且随着温度的升高，界面混合层厚度逐渐增大。目前，通过MD方法，在原子尺度上对含能晶体-黏结剂界面的微观变形机理和力学性能的研究较少，主要原因是难以建立高颗粒填充且能真实反映PBX细观特性的结构模型。

本研究通过MATLAB编程建立复杂的PBX炸药模型，即含氟聚合物F2311包覆在多个以HMX为基的炸药颗粒外形成HMX-F2311界面模型，采用MD方法模拟单轴拉伸下HMX-F2311炸药的力学性能，通过分析HMX-F2311的拉伸应力-应变曲线、界面微观结构演变和能量分布，探讨应变率对PBX力学性能的影响。

1.   力场和模型的建立

1.1   力　场

HMX-F2311的力场由HMX势、F2311势以及HMX与F2311之间的界面相互作用势组成，因此HMX-F2311系统的总能量可写为

式中：${E}_{\mathrm{H}\mathrm{M}\mathrm{X}}$和${E}_{\mathrm{F}2311}$分别为HMX炸药晶体和黏结剂F2311的能量，${E}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}$为炸药与黏结剂之间的界面相互作用能。HMX所用力场是由Smith提出^[20-22]的完整炸药力场，包含分子内和分子间相互作用。黏结剂F2311所用力场为适合聚合物添加剂的COMPASS力场^[23-24]，可准确描述聚合物的链状分子结构。界面相互作用能${E}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}$则用本课题组前期工作中得到的力场^[25]描述，该力场采用Morse势描述界面原子间的吸引力，采用指数衰减函数描述排斥力

式中：${\phi }_{\mathrm{M}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{e}}$、${\phi }_{\mathrm{r}\mathrm{e}\mathrm{p}}$分别为界面原子间的吸引力和排斥力，$r$ 为原子间距离，${D}_{0}$、${r}_{0}$、$\alpha$为Morse力场参数，$A$、$\,\rho$为指数力场参数。由此，将3个不同的势组合成一个描述复杂PBX系统的完整力场，并由LAMMPS程序调用进行MD模拟。

1.2   初始模型

首先，使用MATLAB程序建立F2311包覆HMX的PBX界面初始模型，如图1(a)所示。模型中x、y方向各有4个HMX晶体颗粒，z方向有2个HMX晶体颗粒，黏结剂F2311均匀包裹在颗粒的边缘。模型尺寸约为200 Å×200 Å×100 Å，原子数约为280000。

图  1  HMX-F2311界面的初始模型(a)以及优化后HMX-F2311界面的稳定结构(b)

Figure  1.  Initial model of HMX-F2311 interface (a) and the stable structure of HMX-F2311 interface after optimization (b)

下载: 全尺寸图片 幻灯片

HMX-F2311界面模型的优化过程分为4个步骤：能量最小化、优化原子坐标、退火消除残余应力以及获得平衡构象。首先，对HMX-F2311界面模型进行几何弛豫，采用力收敛公差为10⁻⁸ kcal/(mol·Å)的共轭梯度法进行能量最小化。其次，在NPT系综下进行MD模拟，以优化原子坐标，压力为标准大气压，温度为300 K，步长为0.1 fs，模拟时间为40 ps。接着，为降低系统的能量，消除残余应力，对模型进行退火处理。将体系温度迅速上升至350 K，弛豫50 ps后，降至300 K，继续弛豫50 ps。退火过程在NPT系综下进行，模型尺寸将发生明显变化。最后，采用NVT系综对体系进行结构优化，模拟时间为20 ps。优化后的HMX-F2311界面结构见图1(b)。NVT优化过程中的压力和能量随时间的变化曲线如图2所示。系统的压力和能量最终趋于稳定，上下波动小于5%，表明此时HMX-F2311结构达到平衡。将模拟得到的平衡结构作为初始构型，进行不同应变率下单轴拉伸的MD模拟。

图  2  优化过程中压力和能量随时间的变化

Figure  2.  Pressure and energy fluctuations vs. time during optimization

下载: 全尺寸图片 幻灯片

1.3   拉伸加载

本研究共计算了4组不同应变率下HMX-F2311的单轴拉伸加载过程，应变率分别为0.5×10¹⁰、1.0×10¹⁰、2.5×10¹⁰ 和5.0×10¹⁰ s⁻¹。在NVT系综中分别运行了8.0×10⁵、5.0×10⁵、3.0×10⁵和1.5×10⁵步MD模拟，以执行拉伸过程。整个系统以恒定的应变率沿x轴均匀拉伸，拉伸过程是连续的。因此，整个系统在每一步都处于非平衡状态。在加载过程中，每个步骤都保存了拉伸方向上的压力${p}_{x}$和模拟盒尺寸${L}_{x}$，设${L}_{x0}$为拉伸变形前沿拉伸方向的模拟盒尺寸（$\varepsilon =0$），则拉伸应变$\varepsilon$和应力$\sigma$分别为

2.   动态拉伸下应变率对PBX界面的影响

为了研究应变率对HMX-F2311界面拉伸变形的影响，在应变率为0.5×10¹⁰、1.0×10¹⁰、2.5×10¹⁰ 和5.0×10¹⁰ s⁻¹下进行了一系列动态拉伸模拟，应力-应变曲线如图3(a)所示。不同应变率下的拉伸应力-应变曲线表现出相似的特征，均可分为3个阶段：弹性阶段（Ⅰ）、塑性阶段（Ⅱ）和破坏阶段（Ⅲ）。在拉伸加载初期，应力随应变的增加而线性增加，即弹性阶段（$\varepsilon =0.025$）；在应力-应变曲线中没有明显的屈服点，应力持续非线性增加至PBX的拉伸强度（$\varepsilon =0.080$），即塑性阶段；达到拉伸强度后，随着应变的增大，应力持续减小，PBX的力学性能逐步劣化失效，为破坏阶段。在模拟中，采用应力-应变曲线中的极值应力作为PBX的拉伸强度，并将弹性阶段的应力-应变曲线斜率作为PBX的弹性模量。如图3(b)和图3(c)所示，拉伸强度和弹性模量均随拉伸应变率的增加而增大，随对数应变率的增加而线性增大。

图  3  HMX-F2311在不同应变率下的应力-应变曲线 (a)以及拉伸强度 (b) 和弹性模量 (c) 与对数应变率的关系

Figure  3.  (a) Stress-strain curves of HMX-F2311 under different tension loading; (b) tension strength vs. logarithmic strain rate; (c) elastic modulus vs. logarithmic strain rate

下载: 全尺寸图片 幻灯片

为了比较不同应变率（0.5×10¹⁰、1.0×10¹⁰、2.5×10¹⁰ 和5.0×10¹⁰ s⁻¹）下PBX的微观结构演变，通过OVITO可视化程序计算了拉伸过程中每个原子的局部应变分布，如图4所示。在低应变率（0.5×10¹⁰、1.0×10¹⁰ s⁻¹）拉伸加载下：当 $\varepsilon =0.05$ 时，靠近黏结剂区域的应变明显高于炸药晶体的内部应变，表明在拉伸加载下原子应变主要集中在较软的黏结剂F2311上；当应力达到拉伸强度后（$\varepsilon=0.10$），炸药晶体与黏结剂F2311界面处出现了部分孔洞和微裂纹；当 $\varepsilon =0.15$ 时，靠近加载侧的界面逐渐开裂，出现黏结剂脱粘现象；随着拉伸的持续，黏结剂的脱粘降低了PBX的承载能力，导致材料的力学性能逐渐劣化；当$\varepsilon =0.25$时，PBX的拉应力降到最低，在靠近加载侧的界面处沿垂直于加载方向产生了一条主裂纹。PBX的失效模式为黏结剂脱粘。

下载: 全尺寸图片 幻灯片

图  4  不同拉伸加载应变率下HMX-F2311界面在x方向的应变分布

Figure  4.  Distribution of strain field in the x direction of HMX-F2311 interface under different tension strain rates

下载: 全尺寸图片 幻灯片

如图5所示，相比于低应变率拉伸加载，在高应变率（2.5×10¹⁰ 和5.0×10¹⁰ s⁻¹）拉伸加载下，PBX微观结构演变整体上差异不大。加载初期，PBX的拉伸原子形变主要集中在软黏结剂附近。随着拉伸的持续，相对较弱的HMX-F2311界面发生脱粘。不同的是，在低应变率拉伸下形成了一条大致垂直于加载方向的主裂纹，且断裂位置靠近加载面一侧。随着拉伸应变率的增加，在远离加载侧的界面处也出现了大量微裂纹，破坏路径分布在整个模型上。当 $\varepsilon =0.20$ 时，低应变率拉伸加载下HMX-F2311靠近加载侧的界面处已经产生了一条裂纹，而高应变率加载下HMX-F2311界面处只存在孔洞，裂纹并未开始扩展，表明随着应变率的增大，HMX-F2311的抗拉性能增强。HMX-F2311在单轴动态拉伸下的断裂失效是由于黏结剂脱粘，与Rae等^[7]通过显微镜观察到的断裂方式相似。

图  5  不同拉伸应变率下HMX-F2311界面在x方向的应变分布

Figure  5.  Distribution of strain field in the x direction of HMX-F2311 interface under different tension strain rates

下载: 全尺寸图片 幻灯片

为了从能量角度理解PBX在单轴拉伸过程中的形变机制，分析了不同拉伸应变率（0.5×10¹⁰和5.0×10¹⁰ s⁻¹）下PBX界面的势能随应变的变化情况，如图6所示。在低应变率（0.5×10¹⁰ s⁻¹）拉伸加载下：当PBX拉伸处于弹性阶段和塑性阶段时，势能随着应变的增大逐渐增大；当达到HMX-F2311的拉伸强度时，势能也达到极值（0.68×10⁴ kcal/mol），随后迅速下降。然而，在高应变率（5.0×10¹⁰ s⁻¹）拉伸加载下，当达到HMX-F2311的拉伸强度后，势能仍然持续升高至极值（3.26×10⁴ kcal/mol），随后下降。随着拉伸应变率的增加，HMX-F2311的势能增加了4倍。

图  6  不同应变率下总势能随应变的变化

Figure  6.  Variations of potential energy with stain under different strain rates

下载: 全尺寸图片 幻灯片

模拟体系的势能一般分为成键相互作用和非键相互作用，其中：成键相互作用主要包括键能、角能、二面角能以及非正常二面角扭转能等，而非键相互作用包括范德华力作用和库仑相互作用。拉伸应变率为0.5×10¹⁰和5.0×10¹⁰ s⁻¹时PBX界面各势能分项随应变的变化曲线如图7所示。为了便于分析，势能以及各能量分项均减去未拉伸时的初始值，其中：ΔE_potential、ΔE_bond、ΔE_angle、ΔE_dihed、ΔE_vdwl、ΔE_imp分别为势能、键能、角能、二面角能、范德华力相互作用能、非正常二面角扭转能的变化量，ΔE_vdwl、ΔE_bond、ΔE_angle、ΔE_dihed以及ΔE_imp均由LAMMPS直接输出。从图7可以看出：单轴拉伸过程中影响势能变化的主要因素是范德华力相互作用，尤其在5×10¹⁰ s⁻¹的高应变率下，范德华力相互作用对势能的变化起决定性作用；成键能对势能的贡献较小，其值随应变率的变化不明显；键能和角能在拉伸过程中对势能有一定的影响，尤其是弹性和塑性阶段；二面角能和非正常二面角扭转能在整个拉伸过程中几乎没有变化。

图  7  不同应变率下各个势能分项随应变的变化

Figure  7.  Variations of each energy sub-item of potential energy with strain under different strain rates

下载: 全尺寸图片 幻灯片

PBX的势能随拉伸应变率的增加而快速增大的原因可能与黏结剂F2311的构象有关。在较低的拉伸应变率下，F2311链有更多的时间进行平衡，并在每一步拉伸后进行滑动、旋转和展开以达到平衡，从而减少界面变形的约束。然而，在高应变率下，F2311链改变原子取向的时间较短，导致HMX-F2311界面处的非平衡状态和更大的变形阻力，从而需要较大的应力才能使PBX界面进一步形变。由此可知，低应变率更有利于黏结剂F2311达到构象平衡。

3.   结　论

采用MD方法模拟了HMX-F2311在不同拉伸速率下的力学性能和断裂机理，得到如下结论。

(1) HMX-F2311的力学性能表现出高应变率依赖性。拉伸强度和弹性模量随着应变率的增加而增大，并与对数应变率呈线性关系。

(2) 在单轴拉伸加载下，应变主要集中在黏结剂F2311附近，并逐渐在界面处形成孔洞和裂纹。PBX的断裂是由于黏结剂F2311的脱粘。在低应变率下形成了一条垂直于加载方向的主裂纹，随着拉伸应变率的增加，破坏路径将分布在整个模型上。

(3) 在单轴拉伸加载过程中，HMX-F2311的势能随拉伸应变率的增加而迅速增大。在高应变率下，范德华力相互作用对拉伸过程中势能的演变起决定作用。

研究结果对于理解PBX界面的微观结构演变及失效机理具有重要意义，也有助于通过建立材料微观结构与性能之间的关系设计和改进PBX复合材料。