陈思彤，董可海，王 鑫，裴立冠，孔令泽，夏 成

（1. 海军工程大学兵器工程学院，湖北 武汉 430033；2. 海军航空大学岸防兵学院，山东 烟台 264001）

1 引言

硝酸酯增塑聚醚（NEPE）固体推进剂的危险品等级为1.1 级，当受到冲击波和撞击刺激时可能发生爆轰［1］，其中增塑剂硝化甘油（NG）的敏感性较高、并且存在一定程度的迁移［2］；增塑剂的大量迁移会引起“推进剂⁃衬层”界面的脱粘、影响导弹的弹道性能、增加发动机点火后爆炸的可能性等［3］。

目前，关于组分迁移的研究方法主要有实验法［4］和分子模拟法［5］。实验法主要是采用联合试件、利用色谱分析或者X 射线光电子能谱（XPS）分析等手段［6-7］，虽然能够得到最接近于真实条件的迁移量，但存在操作复杂、XPS 对于硝酸酯的检测不准确等不足［8］；近年来分子模拟法快速发展［9］，成为预测扩散性能、研究扩散微观机理的有效手段，具有耗时短、花费低等优点［10-11］。

对于新型含能钝感增塑剂三羟甲基乙烷三硝酸酯（TMETN）［12-13］，法国、日本对其应用于双基推进剂［14］、聚叠氮缩水甘油醚（GAP）推进剂中的性能进行了研究［15］，而将TMETN 应用于NEPE 推进剂的研究还比较少。为此，本研究采用分子模拟法比较三种增塑剂TMETN、NG、1，2，4⁃丁三醇三硝酸酯（BTTN）在NEPE推进剂黏合剂预聚体PEG/N⁃100 中的扩散性能及机理，并讨论温度和增塑比对TMETN 扩散性能的影响，为钝感高能推进剂的发展提供理论依据。

2 模型构建与模拟方法

2.1 分子动力学建模及模拟细节

真实推进剂的成分和结构都较为复杂，鉴于分子模拟（MD）时间和空间尺度的限制，考虑到固化过程是N⁃100 中异氰酸酯基与PEG 中羟基发生反应，所以构建了高分子预聚体网络的一部分（PEG 分子量为6002［16］），如图1 所示。

图1 PEG/N⁃100 聚合物链片段Fig.1 Fragments of PEG/N⁃100 polymer chains

图2 为优化前和优化后的高分子链构型，可见优化过程使高分子链从直链结构变得缠绕卷曲，同时能量也从8188.83 kJ·mol-1降低到-9202.22 kJ·mol-1，因为真实物质的能量总是最低的，所以优化作用明显使得分子构型更接近于实际［17］。

图2 优化前和优化后的预聚物构型Fig.2 The configuration of prepolymer before and after opti⁃mization

PEG 与增塑剂的质量比（增塑比）约为2.5，为减小链端效应［18］，将1 条PEG/N⁃100 高分子链与相应数量的硝酸酯构成周期箱，NG、BTTN、TMETN 的分子数量分别为：65、61、58 个，初始密度按照三种物质的体积平均值进行设定；对周期箱进行smart 方法的几何优化，再进行10 次NPT 系综的退火处理，退火的初始温度为298 K、中间温度为598 K、最终温度为298 K［17］，每次退火后再进行几何优化，并选取能量较低的构型进行1 ns NPT 系综的密度优化，时间步长均为1 fs，然后进行625 ps NVT 系综动力学计算、其中后400 ps用于采集均方位移随时间变化曲线，模拟控温器和控压器均采用Bederson，有文献证实Bederson 控温器对均方位移产生的影响较小［19］。

2.2 介观建模及模拟细节

在进行介观动力学模拟之前，首先要确定各分子所包含的高斯链段数量、各高斯链段间的相互作用参数。高分子预聚物的高斯链段数量有如下公式［20］：

式中，Nmeso为高斯链段数量；n为预聚物的聚合度，此处为136（PEG 分子量为6002）；C∞为全原子高分子链的特征比。

由Synthia 模块得到PEG 的C∞是4.98，故PEG 高斯链段的质量为219.12（特征比乘以单体的相对分子质量），所以1 条PEG 分子链等于27 个高斯链段、1 个增塑剂分子等于1 个高斯链段。

在MesoDyn 中，两种高斯链段的相互作用ν⁃1εij表示为：

式中，ν-1εij单位为kJ·mol-1；R为气体常数，8.314 J·mol-1·K-1；T为温度，K；χij为Flory⁃Huggins 参数：式中，Δδ是两物质的溶度参数差值，（J·cm-3）1/2；Vr取两种高斯链段摩尔体积的平均值，cm3。

根据文献［21］中MD 计算的各物质溶度参数，计算得到两种高斯链段间相互作用参数ν-1εij、并列于表1。在Meso Dyn 模块中根据增塑比设置珠子的混合比例，三维周期箱尺寸为32.0 nm×32.0 nm×32.0 nm，温度为298 K，步长为20 ns，总时间为1000 μs。

表1 增塑剂、PEG 的溶度参数及高斯链段间的相互作用参数Table 1 Solubility parameters of plasticizers and PEG，as well as the interaction parameters among the Gauss chains

3 结果与讨论

3.1 三种增塑剂的扩散性能

首先，应阐述“扩散”与“迁移”存在的联系与区别，扩散是发生在推进剂内部，而迁移是指推进剂内部的组分扩散至界面处、并进一步扩散到其它层次中（比如衬层），扩散是迁移的前提、可以在一定程度上预测迁移性能，但迁移较扩散更为复杂、还涉及到迁移组分与衬层材料的相互作用等，本节仅比较TMETN、NG、BTTN 三种增塑剂在预聚体内部的扩散性能、不考虑衬层种类对增塑剂扩散的影响。

采用MD 计算扩散系数的常用方法是Einstein 关系式［22］：通过MD 计算出均方位移⁃时间斜率m，即可得到扩散系数D。

采集后400 ps 增塑剂的均方位移，统计长度采用200 ps，原点间隔采用1 帧（0.01 ps），得到MSD随时间的变化曲线，如图3 所示。

图3 三种增塑剂的均方位移随时间变化曲线Fig.3 Mean⁃square displacement（MSD）⁃t curves of plasti⁃cizers

由表2 可见，NG、BTTN 在PEG/N⁃100 中扩散系数的数量级为10-11、TMETN 扩散系数的数量级为10-12，将本次模拟的结果与类似的文献模拟结果进行比较：李红霞［23］采用MD 计算了DOS在HTPB 中的扩散系数为2×10-11，王晓［14］采用MD 计算NG 在PET/MDI、TDI、IPDI 中的扩散系数分别为6.5×10-13、1.25×10-12、1.11×10-11；本文研究的预聚体、固化剂与文献中的稍所不同，但模拟结果的数量级在同一范围内，说明本模拟与文献结果的吻合性较好。由表2 可见，三种增塑剂在预聚体PEG/N⁃100 中扩散系数的大小顺序为：NG＞BTTN＞TMETN，结果表明NG 最容易在黏合剂基体中发生扩散，BTTN 次之，TMETN 最难发生扩散。

表2 三种增塑剂均方位移⁃时间曲线的斜率及其对应的扩散系数Table 2 Slope of MSD⁃t curves and the diffusion coefficient of three plasticizers

3.2 三种增塑剂的扩散机理

据文献［14，24-25］，影响增塑剂扩散的微观因素主要有：增塑剂与预聚体的分子间相互作用、体系的自由体积分数和增塑剂分子的尺寸，下面从这三个角度分析增塑剂在PEG/N⁃100 预聚体中的微观扩散机理。

3.2.1 增塑剂与预聚体的分子间相互作用

径向分布函数g（r）可以反映体系的相互作用类型（峰值位置）及大小（曲线高低）［26］，在混合物的径向分布函数曲线中，异类分子间的相互作用曲线高于同类分子间的相互作用曲线，则体系相容性较好［27，29］。

由图4 可见：TMETN/PEG/N⁃100 体系的相容性较好，而在另外两种体系中，增塑剂与增塑剂分子间的相互作用大于增塑剂与PEG/N⁃100 分子间的相互作用（NG 所在的体系更为明显），所以，BTTN 与PEG/N⁃100 的相容性较差、NG 与PEG/N⁃100 的相容性最差；体系相容性的顺序与扩散系数的顺序在逻辑上保持一致，即体系的相容性越好，增塑剂越难发生扩散。

图4 NG/PEG、BTTN/PEG、TMETN/PEG 混合物的分子间径向分布函数Fig.4 Intermolecular radial distribution functions for NG/PEG，BTTN/PEG and TMETN/PEG

相容性就是复杂的分子间相互作用，为了进一步比较分析各体系内分子间相互作用的差异，将混合物中分子间径向分布函数按照增塑剂⁃PEG、PEG⁃PEG、增塑剂⁃增塑剂的作用类型进行处理，结果如图5 所示。

在三个体系中，PEG 的分子数量相同，NG、BTTN、TMETN 的分子数量依次递减，由图5a 可见，虽然体系中TMETN 的分子数量最少，但该体系中PEG分子附近出现TMETN 的概率却最大，故PEG 与TMETN 的结合能力最强，同理可发现，BTTN 与PEG的结合能力居中、NG 与PEG 的结合能力最弱；由图5b 可见，NG/PEG 体系中PEG 的自聚集能力最强，BTTN/PEG 体系中次之、TMETN/PEG 体系中最弱；由图5c 可见，因三个体系中增塑剂分子数量的大小顺序与该径向分布函数曲线的高低一致，故无法判定“NG附近出现NG 概率最大”是由NG 的自聚集能力强引起的、还是由NG 分子数量多引起的。综上，PEG 与TMETN 的结合能力最强、该体系中PEG 的自聚集能力最弱，这也是TMETN 不易发生扩散的微观原因。

图5 三个体系中增塑剂⁃PEG、PEG⁃PEG、增塑剂⁃增塑剂的分子间径向分布函数Fig.5 Intermolecular radial distribution functions of plasticizer⁃PEG，PEG⁃PEG and plasticizer⁃plasticizer in three systems

由于分子间径向分布函数曲线是对各原子间相互作用进行平均而得到的，所以可能抹平原子间的氢键作用，氢键作用的范围为2.6～3.1 nm，3.1 nm 以后为范德华作用的范围，针对于此，对可能存在氢键的增塑剂与PEG/N⁃100 原子对进行径向分布函数分析，最终发现PEG 中的H 原子与增塑剂中的硝基O 原子存在氢键作用。

三个体系中氢键作用的结果见图6 和表3。由图6和表3可见，三个体系均存在氢键作用，且TMETN_O 与PEG/N100_H 的氢键作用最强，峰值比BTTN 体系中的高0.096、比NG 体系中的高0.1456，说明原子间氢键作用越强、增塑剂越难发生扩散。

图6 增塑剂与PEG 间O—H 原子的径向分布函数Fig.6 Radial distribution functions between the oxygen atom of plasticizer and hydrogen atom of PEG

表3 氢键作用的位置及峰值Table 3 Location and peak of hydrogen bond interaction

增塑剂的扩散会导致体系的混乱程度增加、引起相分离，这将对装药的性能产生较大威胁，而相分离的时间和特征区尺寸都超出了MD 的范围、需借助介观动力学模拟（Meso Dyn）来研究。根据体系的增塑比为2.5，参数设置如表1，最终得到三个体系的密度分布如图7 所示，显示颜色的地方代表着密度高于该物质平均密度的位置和水平。

从图7可见，TMETN/PEG体系中的两相分布较均匀，而BTTN/PEG 和NG/PEG 体系发生了相聚集、且后者更为严重，这也说明TMETN在PEG相中发生扩散的能力最弱，而BTTN和NG发生扩散近而聚集的能力更强；此外还发现，PEG易向周期箱四周聚集、增塑剂则是向中间聚集。

图7 三个体系中的浓度分布Fig.7 Concentration distribution in three systems

3.2.2 体系的自由体积分数

体系内的自由体积为高分子链的运动和小分子的扩散提供空间，MD 中自由体积表示为：

其中，Vf是自由体积，Vsp是周期箱体积，Vw是范德华表面占有体积，单位均为Å3。

自由体积是指除去PEG/N⁃100 分子和增塑剂分子共同占有体积后的剩余体积，采用Atom Volumes& Surfaces 工具对三种增塑体系的自由体积进行计算，探针半径为1.4 Å，因为各体系的总体积大小不一样，为了能够比较三种增塑剂的相对运动空间，定义自由体积分数FFVSim.为自由体积与各周期箱体积的比值：

体系的自由体积分数越大，提供给增塑剂分子运动的相对空间就越大，扩散也就越容易发生。由表4可见，体系内自由体积分数的大小顺序为TMETN/PEG/N⁃100＞NG/PEG/N⁃100＞BTTN/PEG/N⁃100，可见TMETN 具有最大的相对运动空间，但在“自由体积分数”和“TMETN 与PEG/N⁃100 的强分子间作用”的综合影响下，TMETN 的扩散系数还是最小的；此外，NG/PEG/N⁃100 体系中的自由体积分数居中、BTTN/PEG/N⁃100 体系的最小，这说明NG 的相对运动空间较BTTN 的 大，又 由3.2.1 节 知“NG 与PEG/N⁃100 的相容性也较差”，这两种因素的综合作用导致了NG 的扩散系数最大，BTTN 次之。

表4 体系的自由体积分数Table 4 Fractional free volume of systems

3.2.3 增塑剂分子的尺寸

增塑剂分子的大小也会影响扩散性能，分子越小，越容易在聚合物孔洞中跳跃；由增塑剂分子结构可知，三种增塑剂分子的大小顺序为：TMETN＞BTTN＞NG，这与扩散系数的顺序在逻辑上是一致的，即NG 分子最小、最容易在聚合物大分子孔洞中发生跳跃，BTTN次之，TMETN 最难。

综上，从“分子间相互作用”、“体系的自由体积分数”、“增塑剂分子的尺寸”三个方面讨论了增塑剂扩散性能的微观机理，在TMETN/PEG/N⁃100 体系中，体系相容性最好、TMETN 分子尺寸最大，这两方面都加大了扩散的难度，但同时该体系的自由体积分数也最大，这方面会促进扩散的发生，最终在三种因素的综合作用下，TMETN 的扩散系数还是最小的；在NG/PEG/N⁃100 体系中，体系的相容性最差、NG 的分子尺寸最小、体系的自由体积分数较BTTN/PEG/N⁃100 的大，这三种因素均导致NG 的扩散系数比BTTN 的大。

3.3 温度对TMETN 扩散性能的影响

3.3.1 温度对体系相容性的影响

根据推进剂储存、使用及固化等条件，比如北方冬天低温可达0 ℃，常温是25 ℃，固化温度大约为50 ℃，所以选取0、25、40、50 ℃，分析TMETN/PEG/N⁃100 体系中TMETN 的扩散系数，采集的MSD 随时间变化曲线如图8，并以图8 为根据、经过拟合得到扩散系数，将扩散系数列于表5。不同温度下扩散系数的变化如图9 所示。

图8 不同温度下TMETN 均方位移随时间的变化曲线Fig.8 Mean⁃square displacement（MSD）⁃t curves of plasti⁃cizers at various temperatures

由表5 和图9 可见，随着温度升高，扩散系数先缓慢增加、后剧烈增加，即温度越高、扩散系数增加的速率越大，这与文献［23］中的规律一致、也与高温加速老化的规律一致；常温298 K 时扩散系数为低温273 K的1.2 倍左右；固化温度323 K 时的扩散系数大约为常温298 K 时扩散系数的2 倍，这解释了为何增塑剂在固化过程中已发生较大的迁移［30］。

表5 不同温度下增塑剂TMETN 的均方位移⁃时间曲线的斜率及对应的扩散系数Table 5 The slopes of MSD⁃t curves and the diffusion coeffi⁃cient of TMETN at various temperatures

图9 不同温度下扩散系数的变化Fig.9 The diffusion coefficient at various temperatures

进一步分析不同温度下TMETN 中硝基O 原子与PEG 中H 原子间的氢键作用，为了便于比较各温度下氢键的峰值和位置，将273，298，313，323 K 条件下的径向分布函数曲线均置于图10 中。不同温度下氢键作用的位置及峰值见表6。

图10 不同温度下TMETN 与PEG 间H—O 原子的径向分布函数曲线Fig.10 Curves for radial distribution functions between plas⁃ticizer（O）and PEG（H）at various temperatures

由 图10 和 表6 可 见，随 着 温 度 升 高，TMETN 与PEG 原子间氢键作用峰的位置后移、峰值变低，这说明高温使增塑剂与黏合剂间的相互作用变弱，这导致了TMETN 更容易发生扩散。

表6 不同温度下氢键作用的位置及峰值Table 6 Location and peak of hydrogen bond interaction at various temperatures

3.3.2 温度对体系自由体积分数的影响

由表7 可见，随着温度升高，TMETN/PEG/N⁃100体系的自由体积分数变大，说明高温使聚合物的活动性增强、分子间距离增加，这给增塑剂提供了更大的活动空间，在宏观上表现为聚合物的溶胀度提高［23］；并且温度的升高也会增加小分子的动能，即小分子的跳跃能力变强，这均使增塑剂更易发生扩散。

表7 不同温度下体系的自由体积分数Table 7 Fractional free volume for systems at various tem⁃peratures

3.4 增塑比对TMETN 扩散性能的影响

根据NEPE 推进剂中增塑剂的含量，研究了在常温298 K 条件下，增塑比为2.5、2.8、3 时TMETN 的扩散 性 能［31］，得 到MSD⁃t曲 线 拟 合 斜 率 及 扩 散 系数（表8）。

表8 增塑剂TMETN 的均方位移⁃时间曲线斜率及对应的扩散系数Table 8 Slopes of MSD⁃t curves and the diffusion coeffi⁃cients of TMETN

由表8 可见，随着增塑比的增加，TMETN 扩散系数在减小，为了观察增塑比对体系形态分布的影响，采用MesoDyn 进行了1000 μs 的介观动力学模拟，得到TMETN 等密度图及有序度参数。由图11 可见，随着增塑比的增加，组分发生同相归并的程度变小。由图12 可见，TMETN 的有序度参数减小、体系的熵增减小，所以体系相容性变好，这是TMETN 的扩散系数减小的一个原因；此外有研究表明［23］，随着增塑剂含量的增加，大量的增塑剂会导致高分子网络收缩、使得增塑剂在孔隙中发生有效跳跃的难度变大、从而使扩散系数减小，这一点在达到其共容增塑极限后更为明显。

图11 不同增塑比条件下TMETN 的等密度图Fig.11 The isosurface of the density fields for TMETN at vari⁃ous plasticization ratios

图12 不同增塑比条件下TMETN 的有序度参数Fig.12 The order parameters for TMETN at various plastici⁃zation ratios

4 结论

（1）增塑剂在PEG/N⁃100 中的扩散系数大小顺序为：NG＞BTTN＞TMETN，这说明含能钝感增塑剂TMETN 具有较强地抗扩散能力，可以进一步研究将其代替或部分代替NG、应用于NEPE 推进剂中。

（2）从增塑剂与黏合剂的分子间相互作用、体系的自由体积分数和增塑剂分子的尺寸三个方面分析了增塑剂扩散性能的微观机理，在三个增塑体系中，TMETN 与PEG/N⁃100 的分子间相互作用最强、体系的自由体积分数最大、TMETN 分子的尺寸最大，各因素综合作用导致了TMETN 的扩散系数最小。

（3）随着温度升高，TMETN 扩散系数先缓慢增加、后剧烈增加，这与高温加速老化的规律一致；固化温度323 K 时的扩散系数为常温298 K 时的2 倍左右，这在一定程度上解释了“增塑剂在固化过程中已发生较大迁移”的现象；从微观角度来解释扩散系数增大的原因为：高温导致TMETN 与PEG 的原子间氢键作用减弱、体系的自由体积分数变大。

（4）随着增塑比增加，TMETN 扩散系数减小，介观动力学研究表明“体系的有序度参数减小、体系相容性变好”是其中的原因之一。