1 引言

太安(PETN, 季戊四醇四硝酸酯, $\text{C}_{5}\text{H}_{8}\text{O}_{12}\text{N}_{4}$)是典型的硝酸酯类炸药, 具有较高化学稳定性和优良的爆轰性能, 广泛应用于导爆索、雷管、传爆药柱和某些小口径炮弹的装药, 其分子结构和晶体原胞结构如图 1所示。关于PETN晶体结构和性能的经典分子动力学(MD)模拟研究虽已较多^[1-14], 但至今未见对其感度进行理论判别和对其力学性能进行系统预测的MD文献。

感度即安全性，是高能物质最重要性能之一。先前的微观理论判据主要是通过量子化学计算求得的, 用于高能分子较多^[15]。对于成千上万原子的超晶胞或复合材料, 近年来已有关于HMX/AP^{[16, 17]}、HMX基PBX^{[18, 19]}、RDX完美和缺陷晶体及RDX基PBX^{[20, 21]}感度判别的工作发表。我们所建议的引发键最大键长($L_{\text{max}}$)和引发键键连双原子作用能($E_{\text{A}—\text{B}}$)等判据能否对单体炸药适用, 有待考察验证。

力学性能也是高(含)能材料最重要的性能之一。过去的十年, 本课题组已开辟了高能复合材料(PBXs、发射药和推进剂)弹性力学性能的预测研究^[22-30], 主要都是基于静态力学分析。而收集MD平衡后的原子运动全轨迹加以分析的波动法^[31], 虽然耗机时较多, 但对力学性能的预测相对要准确可靠, 故本研究试以波动法求得PETN晶体的力学性能。

在物质结构-性能微观理论研究中, 经典MD模拟较量子力学(QM)计算有一个优点, 就是能求得庞大复杂体系结构、性能随温度递变的统计平均规律。因为研究易爆燃物质的感度和力学性能及其温度效应极为重要, 故而本实验以PETN晶体为研究对象, 选择构建其(4×3×4)超晶胞, 并选择其较为钝感的(100)面^{[3, 13, 32]}加以切割分面, 将所得这两种模型, 在5种不同温度(195, 245, 295, 345，395 K)和NPT系综下, 进行MD模拟研究, 重点探讨结构、能量和力学性能及其温度效应, 寻求和核准微观感度理论判据。

2 模型搭建和模拟细节 2.1 力场选择和模型搭建

本研究之所以选择COMPASS力场^[33], 是因为该力场能较好地适用于凝聚态, 包括硝酸酯及其相关物^[34-36]。以PETN的X射线衍射晶体数据^[37-38]为依据, 构建PETN(4×3×4)超晶胞并置于具周期性边界条件的周期箱中; 每个周期箱中含96个PETN分子, 共2784个原子, 即为PETN晶体模型。沿PETN(4×3×4)超晶胞(100)晶面方向进行“切割”, 并使Z轴平行晶轴$c$向量, 同时垂直于(100)晶面, 晶面(100)面积为28.14×26.84 Å$^{2}$, $c$方向上取$OC$=58.16 Å, 其上真空层高度设置为0。于是构成含144个PETN分子即4176个原子的PETN(100)周期箱。PETN晶体和PETN(100)的初始模型见图 2。

2.2 模拟方法和细节

将搭建好的PETN和PETN(100)两种晶体模型分别在COMPASS^[33]力场、NPT系综下进行MD周期性模拟研究。模拟过程中温度和压力的控制采用Anderson ^[39]和Parrinello^[40]方法, 范德华(vdW)和静电作用(Coulomb)分别用Atom-Based^[41]和Ewald^[42]方法, 截断半径取9.5×10$^{-10}$ m, 并进行截断尾部校正。温度依次选取195, 245, 295, 345, 395 K。初始分子运动速度按Maxwell-Boltzman分布确定; 在周期性边界条件和时间平均等效于系综平均基本假设之上求解牛顿运动方程; 积分采用Verlet方法; 时间步长取1 fs。总模拟步数取2 ns, 其中前1 ns用于热力学平衡, 后1 ns用于统计分析。通过温度-时间和能量-时间曲线以及平衡后晶胞参数比较, 确证PETN和PETN(100)晶体在各温度下均已达到平衡。作为示例, 图 3给出295 K PETN和PETN(100)经MD模拟所得平衡结构。

3 结果与讨论 3.1 $\mathbf{PETN}$晶体的晶胞参数

表 1给出PETN晶体295 K下NPT-MD模拟所得晶胞参数和密度。为方便比较, 表中还给出了实验值^[38]和Sorescu等的理论计算结果^{[1, 10, 43]}。由表 1可见, 本工作模拟所得晶胞参数与前人的实验和计算值吻合较好, 晶胞参数相差很小或相等; 与实验密度相对误差也只有0.6%。由此进一步表明COMPASS力场对PETN晶体有较好的适用性, 也表明本研究的MD模拟确已达到平衡, 所得平衡结构可信。

3.2 感度与引发键最大键长的关系

MD模拟能提供各种体系在广义平衡结构下的键长统计分布, 包括在不同温度下的键长分布, 这是MD很有意义的优势所在。PETN的引发键一般认为是分子中的O—$\text{NO}_{2}$键^{[3-4, 9]}。图 4给出295 K PETN和PETN(100)中O—$\text{NO}_{2}$引发键的键长分布。表 2示出PETN和PETN(100)中O—$\text{NO}_{2}$键的最可几键长($L_{\text{prob}}$)、平均键长($L_{\text{ave}}$)和最大键长($L_{\text{max}}$)。图 5示出PETN和PETN(100)晶体中O—$\text{NO}_{2}$键最大键长与温度的关系。

由图 4可见, PETN和PETN(100)的引发键键长分布均呈近似对称的高斯型分布, 最可几键长和平均键长较为接近。PETN晶体引发键键长在1.350~1.480Å范围的占97%, PETN(100)在1.345~1.485Å范围的占98%。由表 2可见, 它们在295 K下的$L_{\text{ave}}$均为1.414Å, 且与PETN的实验键长(1.404Å)^[37]较为相近, 相对误差小于1%;由表 2和图 5可见, 随着温度的升高, 两种模型中的$L_{\text{max}}$单调增大, 与感度随温度升高而增大的实验事实相一致。虽然$L_{\text{max}}$在键长分布中出现的几率只有10$^{-8}$, 但具有$L_{\text{max}}$时分子能量很高, 非常活泼。根据只有“活化”分子才能发生化学反应的碰撞理论, 具有$L_{\text{max}}$的活泼分子中的引发键最容易断裂, 进而引发分解和起爆。这与先前用引发键的$L_{\text{max}}$关联不同温度下的热和撞击感度的研究^{[16, 17, 19-21]}一致。由表 2和图 5还可见, PETN晶体中$L_{\text{max}}$在245 K(1.619Å)和295 K(1.620Å)以及在345 K(1.661Å)和395K(1.662Å)之间的递增幅度很小, 随温度递变规律不甚明显。而在PETN(100)模型体系中, $L_{\text{max}}$随温度的单调递增很显著, 符合通常规律。这就表明, 切割分面进行MD研究, 能提供较好的确定性和规律性。

3.3 感度与引发键连双原子作用能的关系

定义硝酸酯类体系中引发键(O—$\text{NO}_{2}$)连双原子作用能($E_{\text{O}-\text{N}}$), 为在COMPASS力场作用下, 经MD模拟所得平衡体系的总能量$E_{\text{T}}$, 与固定所有引发键中O和N原子求得的体系总能量$E_{\text{T}}$′的差值, 再除以体系中所含O—$\text{NO}_{2}$键的数目$n$, 即$E_{\text{O}-\text{N}}$=($E_{\text{T}}$-$E_{\text{T}}$′)/$n$。

表 3和表 4给出PETN和PETN(100)晶体在不同温度下的$E_{\text{O}-\text{N}}$及其分量, 图 6示出$E_{\text{O}-\text{N}}$随温度的变化。由表 3、表 4和图 6可见, 随着温度的升高, PETN和PETN(100)晶体中$E_{\text{O}-\text{N}}$单调递减, 这与感度随温度升高而升高的实验事实相符。由此表明, MD模拟所得引发键连双原子作用能的统计平均值$E_{\text{O}-\text{N}}$, 亦可作为高能晶体热和撞击感度的理论判据。由表 3和表 4还可见, $E_{\text{O}-\text{N}}$的大小主要由原体系与固定引发键连O与N原子后体系之间的静电力的差值所决定, 其中vdW力和内势能的贡献很小。从图 6可显见, 在各不同温度下, PET(100)晶体模型的$E_{\text{O}-\text{N}}$均比PETN晶体的大, 这与早先认为PETN(100)晶面钝感的实验^[32]和近期量子MD以及反应性力场下的MD结果^{[3, 13]}一致。由此表明, 切割分面尤其是选取较钝感的(100)面进行MD研究能反映PETN晶体的客观实际状况。

3.4 感度与内聚能密度的关系

内聚能密度(CED)是单位体积1 mol凝聚态克服分子间作用力变为气态时所需能量。在MD模拟中CED是vdW力与静电力之和, 即分子的非键力。基于MD模拟轨迹所得不同温度下PETN晶体和PETN(100)的内聚能密度及其分量。见表 5。

由表 5可见, 随温度升高, PETN和PETN(100)的CED、vdW力和静电力均单调递减, 表明它们由晶态变为气态时所需能量变小, 这与温度升高感度增大的实验事实相一致。由此表明, 在一定条件下, PETN晶体的CED也可用于热感度相对大小的理论判据; 且在每个温度下PETN和PETN(100)的CED及其分量数值都惊人地一致, 表明CED的统计平均值不受MD模型的影响。

3.5 力学性能比较

弹性模量是评价材料刚性的指标, 是材料抵抗弹性形变能力的度量^[44]。而弹性性能和塑性性能是相关联的^[45], 通常表示阻止塑性形变能力的硬度和拉伸强度与剪切模量呈正比; 断裂强度与本体模量呈正比; 本体模量与剪切模量的比值($K$/$G$)则用来衡量体系的延展性, 即材料发生形变而不产生裂缝的能力, 其值较大则表明材料延展性较好。柯西压($C_{12}$-$C_{44}$)亦可用来度量体系的延展性^[46], 若其是正值, 则表明材料延展性较好, 若其为负值, 则材料显脆性。在用于判别材料延展性相对大小时, $K$/$G$值与柯西压的区别在于前者对材料延展性能的判别是基于产生塑性形变的程度, 而后者是基于材料断裂面的形貌^[19]。基于PETN和PETN(100)不同温度下MD模拟的原子运动轨迹, 通过波动法分析, 表 6给出所得力学性能。图 7示出PETN(100)的力学性能随温度的变化。

根据广义虎克定律, 反映应力应变关系的弹性系数应有36个; 因$C_{ij}$=$C_{ji}$, 故极端各向异性体只有21个独立的弹性系数。表 6中未列入的数值等于或接近于零, 表明PETN不是极端各向异性弹性体; 随着温度升高, $C_{11}$, $C_{22}$和$C_{33}$组、$C_{12}$, $C_{13}$和$C_{23}$组以及$C_{44}$, $C_{55}$和$C_{66}$组内的数值彼此接近, 表明其接近各向同性弹性体。由表 6和图 7可见, 随着温度升高, PETN和PETN(100)晶体的弹性系数和各弹性模量均单调下降, 说明其各向同性增长和刚性下降, 柔韧性增强, 符合力学性能递变的实验事实。但与感度随温度升高的实验事实不符, 这是因为影响感度的因素很多, 随温度变化, 结构和能量成为主导因素, 而弹性模量降为次要因素所致。另一方面, 随温度升高, PETN(100)的$C_{12}$-$C_{44}$值明显减小, $K$/$G$值也呈下降趋势, 表明PETN(100)基于柯西压和$K$/$G$的延展性随温度升高而降低, 与感度随温度升高而升高的实验事实相符合。这表明, 切割分面的PETN晶体具有稳定一致的力学性能。

4 结论

通过PETN和PETN(100)晶体两种模型在5个温度下的NPT-MD模拟, 得出如下主要结论:

(1) 模拟所得PETN晶胞参数和密度与实验值接近, 表明COMPASS力场适用、所得平衡结构可信。

(2) 模拟求得PETN晶体的引发键(O—$\text{NO}_{2}$)键长分布。平均键长接近实验值, 最大键长($L_{\text{max}}$)随温度升高单调递增, 与感度随温度升高而增大的实验事实相一致。表明引发键$L_{\text{max}}$确可作为高能晶体热和撞击感度相对大小的理论判据。

(3) PETN晶体的引发键连双原子作用能($E_{\text{O}-\text{N}}$)和内聚能密度(CED), 随温度升高而单调递减, 表明$E_{\text{O}-\text{N}}$和CED与感度存在负相关。$E_{\text{O}-\text{N}}$和CED在一定条件下均可作为热和撞击感度相对大小的理论判据。

(4) 由波动法求得PETN晶体的力学性能。弹性系数和弹性模量随温度升高而减小, 表明各向同性递增, 刚性降低。

(5) 切割分面研究PETN晶体, 能获得较确定的结构、相互作用能和力学性能及其温度效应。