NG/TMETN 混合增塑PEG 体系的力学、感度及相容性能

陈思彤，董可海，夏 成，王 鑫，裴立冠，唐岩辉

（1. 海军航空大学 岸防兵学院，山东 烟台 264001；2. 海军航空大学 航空基础学院，山东 烟台 264001）

1 引言

增塑剂是固体推进剂的重要组成部分，主要作用是改善推进剂的力学性能，还可以提供一定的能量、影响导弹的弹道性能。高能固体推进剂NEPE 因比冲大、能量高而闻名，其常用增塑剂为硝化甘油（NG），而NG 极为敏感，特别是对撞击、冲击等刺激［1］，这导致了NEPE 推进剂的危险品等级为1.1 级、有整体爆轰的危险［2］；随着大量高集成、高价值武器装备的发展和使用，“安全性”问题变得日益突出，所以降低NEPE 推进剂感度的研究十分必要和重要。

已有研究表明，含能钝感增塑剂三羟甲基乙烷三硝酸酯（TMETN）可以起到与NG 类似的增塑作用、且其感度远低于NG［3-4］；最早法国［5］将TMETN 应用于改善双基推进剂的感度；后来美国［6］曾用低能硝酸酯改性NEPE 推进剂、有效地降低了感度；也有研究［7］用它改性缩水甘油叠氮聚醚推进剂（GAP），结果发现该推进剂具有良好的力学性能；相比之下，我国对含能钝感增塑剂TMETN 的研究还较少，特别是将它代替或部分代替NG、降低NEPE 推进剂感度的研究则更少。随着计算机技术的快速发展，分子模拟法成为材料设计、性能评估的快速有效手段，它不需要大量的时间和原材料，还可以为实验提供指导和依据、避免实验的盲目性［8-9］，基于以上背景，本工作采用分子模拟法，研究了TMETN/NG 混合增塑聚乙二醇（PEG）体系的力学性能、感度及相容性，为含能钝感推进剂的发展提供依据。

2 建模过程与模拟方法

2.1 MD 模拟细节与平衡判断

2.1.1 模拟细节

用Material Studio 中Visualizer 模块建立聚合度为30 的PEG 分子［10-11］及增塑剂NG、TMETN 分子，利用Forcite 模块Smart 方法对分子进行几何优化，优化过程选择Compass 力场［11］、Ultra-fine 标准。

为了使各组分均匀随机分散、充分相互作用，利用Amorphous Cell 模 块，在 尺 寸 为200Å ×200Å ×200Å［11-13］的周期箱中，根据增塑比为2.5［14］，建立模型PEG/NG、PEG/TMETN、PEG/NG/TMETN-1、2、3（NG：TMETN 质量比分别为3：1、1：1、1：3），五种体系内分子的数量及原子总数见表1。构建的模型如图1 所示

表1 五种体系的组成Table 1 The compositions of five systems

图1 PEG 和增塑剂的分子构型Fig. 1 Molecular configurations of PEG and plasticizers

将初始周期箱进行5000 步几何优化、及20 ps 正则系综（NVT）的动力学平衡，然后逐步压缩至理论密度［10］（按体积进行平均），每次压缩后都进行20 ps NVT 系综的动力学平衡；将最终达到理论密度的周期箱进行几何优化至收敛，然后进行等温等压系综NPT退火处理，退火条件为［15］：初始温度298 K，中间温度698 K，每次退火结束后自动进行几何优化，在结果中选择能量最低的一帧作为最优构型［16］；然后，对最优构型进行1 nsNPT 系综的密度优化，温度298 K，控温器为nose，压力为标准大气压，控压器为Berendsen，再采用NVT 系综进行1 ns 的平衡计算和1 ns 的动力学计算（用于数据采集）［17］，时间步长均为1 fs，为了保证引发键键长的数据量足够多，数据采集时令每0.01 ps（10 步）输出一帧［18］，引发键键长和径向分布函数均采集1 ns 的，因MD 计算力学性能时间过长，所以力学性能采集最后200 ps 的（已平衡）。

2.1.2 平衡判定

首先采用NPT 系综进行了密度优化、然后采用NVT 系综进行了平衡和生产计算，分别进行两种系综平衡的判断。各模拟体系的原子数量类似，以其中一个体系（PEG/NG/TMETN-3）为例，说明体系已达到了平衡：首先，NPT 系综平衡的依据为“体系的密度变化较小”，图2a 为1 ns NPT 系综计算过程中的密度变化情况，红色线为Forcite 模块分析工具得到的密度变化趋势，红色线显示第900 ps 的密度值为1.3912、第950 ps 为1.3912、第1000 ps 为1.3915，可见密度在模拟后期已经达到了动态平衡；NVT 系综平衡的依据为“体系的温度波动较小（大概在5%～10%）、能量波动较小”［11］，如图2b 和图2c 可见，在1 ns NVT 系综生产计算中体系已经是平衡状态，可以证明采集的数据是有效的。

2.2 Meso Dyn 模拟细节

介观动力学模拟是将多个原子或者分子简化成一个珠子，忽略珠子内部原子的相互作用，仅考虑珠子间的排斥力，故而实现了更大尺度的模拟；在进行模拟之前，首先要确定每种分子的粗粒化拓扑结构及各珠子间的相互作用参数。

采用式（1）将高分子链PEG 转化为介观珠子［19］：

式中，Nmeso为粗粒化珠子数量；n 为聚合度；C∞为高分子链的特征比（Synthia 模块计算为4.98）。

一个珠子的质量为219.12（单体的相对分子质量乘以极限特征比），所以PEG 分子（n=30）要采用6 个珠子代替，每个增塑剂用1 个珠子代替。两种珠子的相互排斥参数ν-1εij由公式（2）得到：

式中，ν-1εij为相互排斥参数，kJ·mol-1；R 为气体常数，8.314 J/（mol·K）-1；T为温度，K；χij为Flory-Huggins参数：

图2 平衡参量随时间的变化曲线Fig. 2 Variables vs. simulation time

式中，Δδ 为两物质的溶度参数差值，（J·cm-3）1/2，可由分子动力学模拟MD 模拟计算得到；Vr为参比体积，取两种珠子摩尔体积的平均值，cm3。

采用MD 模拟计算物质间的溶度参数差值：将单分子优化模型构建成纯物质周期箱，PEG 周期箱中包含4条PEG 分子链，NG、TMETN 周期箱中分别包含50 个、39 个分子，对周期箱进行几何优化和NPT 系综退火处理，因每个周期箱仅包含一种分子，体系较容易达到平衡，所以进行800psNPT 系综的密度优化、400 ps NVT系综的动力学计算（后200 ps 用于采集溶度参数）。如图3a 所示，NPT 系综后期已经达到平衡、仅取最后一帧进 行NVT 计 算，如 图3b 和 图3c 所 示，NVT 系 综 后200 ps已经达到平衡，统计的溶度参数结果见表2。

图3 平衡参量随时间的变化曲线Fig. 3 Variables vs. simulation time

PEG 的溶度参数是与分子量相关的，PEG 分子量越大、其溶度参数越小［20］，所以本模拟计算的溶度参数值略大于文献中高分子量PEG 的，这是合理的；此外，NG 和TMETN 模拟得到的密度均略大于真实密度，这与计算所采用的NPT 系综有关，NPT 系综是根据设定的力场自动寻找平衡状态（密度可变）、且最终达到的是动态平衡；此外，NG 溶度参数的计算结果与文献值吻合较好，且各物质的溶度参数是在同一力场条件、参数设置下得到的，故可以进行比较和共同使用［8］。

在得到各物质的溶度参数后，利用式（2）、式（3）计算珠子间的排斥参数，相同种类珠子间的排斥参数为0，不同种类的对应结果列于表3。再根据五个体系中PEG、NG、TMETN 的质量比，在Meso Dyn 模块中设置相关参数，三维周期箱尺寸为32.0 nm×32.0 nm×32.0 nm，温度为298 K，步长为20 ns，总时间为1000 μs。

表2 MD 模拟值与文献值的比较Table 2 Comparison between MD simulation values and literature values

表3 粗粒化珠子的相互排斥参数Table 3 Exclusion parameters between beads

3 结果与讨论

3.1 力学性能

力学模量是用来表征材料抵抗弹性变形能力的指标，塑性和断裂性质也与模量相关联，模量越大，材料的刚度、硬度、屈服强度就越大［10-11］；柯西压（C12-C44）和K/G（体积模量与剪切模量的比值）均用来表征材料的韧性（延展性），材料的延展性越好，其加工性能越好，区别在于前者考虑的是“塑性变形程度”、而后者考虑的是“断裂面形貌”［23-25］；对混合体系NVT 计算得到的后200 ps 轨迹进行力学性能计算，采用恒应变法，每次应变分6 步，最大应变幅为0.003，计算前先对结构进行几何优化，最终得到体系的弹性系数和力学模量见表4。

在图4 中，为了方便表示，各体系按照TMETN 含量增加进行编号，将PEG/NG 体系命名为第1 组，PEG/NG/TMETN-1 为第2 组，以此类推。

由 表4 和 图4 可 见，相 比 于PEG/NG 体 系，PEG/TMETN 体系的杨氏模量和剪切模量均较小、但体积模量较大、柯西压和K/G 也较大（且均为正值），这说明PEG/TMETN 体系更“柔软”、延展性更好，而PEG/NG体系的刚度、硬度、拉伸强度都更大；有研究表明：“与粘合剂相容性较好的增塑剂可以使体系的模量下降、柔性增强、力学性能得到改善”［26］，通过下文计算发现TMETN 与PEG 的相容性优于NG 与PEG 的，所以力学性能的计算结果符合已有的研究规律。

表4 五种体系的弹性系数和力学性能Table 4 The elastic constants and mechanical properties of five systems

图4 体系的杨氏模量、体积模量、剪切模量随TMETN 含量增加的变化Fig. 4 Dependence of young，bulk and shear modulus of the system on the content of TMETN

相比于单一增塑体系，混合增塑体系的杨氏模量、剪切模量均较大，但体积模量却较小，这说明混合增塑体系的刚度和硬度、拉伸强度均较大，这是因为混合增塑剂减小了粘合剂的各向异性［27］；同时，混合增塑体系的柯西压和K/G 均较小、但都是正值，这说明其延展性虽然较单一增塑体系的差、但仍然具有良好的加工性能；由表4 和图4 可见，在三种混合增塑体系中，力学模量的大小顺序均为：PEG/NG/TMETN-2＞PEG/NG/TMETN-1＞PEG/NG/TMETN-3，可见PEG/NG/TMETN-3体系的力学性能最优，并且值得注意的是，虽然混合增塑体系的模量普遍高于单一增塑体系的，但是PEG/NG/TMETN-3（第4 组）的剪切模量、体积模量却低于PEG/NG（第1 组），杨氏模量略高于PEG/NG 的，这说明当TMETN 含量较高时，体系的力学性能得到了改善。

综上所述，在五种增塑体系中，PEG/TMETN 体系的力学性能和加工性能最好，但是众所周知TMETN的能量较NG 低，为了保证NEPE 推进剂的能量水平，进一步考虑TMETN 部分替代NG 的情况，PEG/NG/TMETN-3 体系具有较好的力学性能。

有研究表明，力学性能依赖于体系的组成、物质结构和相互作用等［28］。所以，氢键在一定程度上影响着体系的力学性能［12］，分子间相互作用通常分为氢键（径向分布函数峰r＝2.6～3.1 Å）和范德华作用（r＝3.1～5.0 Å），通过分析预聚物与增塑剂间的径向分布函数，发现PEG 中的氢原子与增塑剂中的端部氧原子存在氢键作用。

图5 五种体系中的原子间径向分布函数Fig.5 Radial distribution functions between atoms in the five systems

氢键强弱不仅与体系内硝基的浓度有关，还与增塑剂（硝基）与PEG（氢原子）的浓度配比有关［12］。由图5 可见，五种体系内的氢键作用强弱为：PEG/TMETN＞PEG/NG＞PEG/NG/TMETN-2＞PEG/NG/TMETN-1＞PEG/NG/TMETN-3，单一增塑体系中的氢键作用较混合增塑体系中的更强，而混合增塑比例（NG∶TMETN）为1∶1 时的氢键最强、比例为1∶3 时的氢键最弱。

对于三种混合增塑剂体系，氢键作用的强弱与力学模量的大小顺序一致，这说明对于组成相似的体系，氢键可以较好地反映力学性能；对于两种单一增塑体系，PEG/TMETN 的氢键作用比PEG/NG 的更强，但其力学模量却更小，这说明对于组成不同的体系，力学性能受物质本身的性质影响较大。

3.2 感度性能

3.2.1 不同体系的感度比较

引发键是体系中最脆弱的化学键，当受到外界刺激时，活化分子中的引发键最容易发生断裂、引起爆炸，所以引发键可以用来关联感度，硝酸酯增塑剂的热解和起爆引发键为O—NO2键［28］，图6 为PEG/NG 体系中引发键键长的分布情况。由图6 可见，引发键键长呈现类高斯分布［29］，最可几键长（最大概率键长）出现在平均键长附近，且最大键长和最小键长出现的概率均非常低，下面统计五种增塑体系的最大引发键键长、最可几键长、平均键长及其概率。

图6 PEG/NG 体系中O—NO2的键长分布Fig.6 Distribution of O—NO2 bond length in PEG/NG blends

五种体系内的引发键键长及其概率见表5。由表5可见，体系中最大引发键键长出现的概率为10-8，与文献［30-31］中的数量级相同，虽然其概率很低，但它代表着引发键的最脆弱状态，其值越大、引发键越容易断裂、体系的感度越大［13］。五种增塑体系中最大引发键键长的大小 顺 序 为：PEG/NG＞PEG/NG/TMETN-1＞PEG/NG/TMETN-2＞PEG/NG/TMETN-3＞PEG/TMETN；可见，随着TMETN 含量的增加，体系的感度不断降低，在混合增塑体系中，PEG/NG/TMETN-3体系具有最低的感度。

由表5 可见，各体系的最可几键长差别不大、均略低于平均键长，并且PEG/NG 体系的最可几键长小于PEG/TMETN 体系的，可见不能用最可几键长表征感度；此外，平均引发键键长的大小顺序为：PEG/NG＞PEG/NG/TMETN - 1＞PEG/NG/TMETN - 2=PEG/NG/TMETN-3＞PEG/TMETN，可见随着TMETN 含量增加，体系内引发键键长呈现整体减小的趋势，不过对于感度不同的两种体系2 和3（3 的感度更低），它们的平均键长却相等，所以平均值并不能较好地反应“感度”这一极端情况。

综上所述，在五种增塑体系中，PEG/TMETN 体系的感度最低，在混合增塑体系中，PEG/NG/TMETN-3体系的感度最低。

表5 五种体系内的引发键键长及其概率Table 5 Length and probability of O—NO2 trigger bond in five systems

3.2.2 混合增塑体系的热感度判据

为了更准确地探究规律，对体系PEG/NG/TMETN-1和PEG/NG/TMETN-3 均进行了温度为298，313，328，343 K 的动力学计算，统计了不同温度下的引发键键长数据，结果见表6 和表7。由表6 和表7 可见，随着温度升高，两种体系的最大引发键键长均明显增大，说明体系正在变得敏感，这与事实规律相符；最可几键长的变化规律不明显；平均引发键键长不断增大、但变化太小，这说明引发键整体在变脆弱，但平均键长仍无法明显地反映感度的变化。此外，体系的内聚能密度（CED）是衡量单位体积内每摩尔物质从凝聚态变为气态所需能量的指标，这种状态改变所需的能量越大（CED 越大），表明越难发生、体系越不敏感［17］。

由表8 和表9 可见，随着温度升高，两体系的内聚能密度及其范德华分量、静电分量均呈现下降趋势，这说明体系状态改变所需要的能量变小，即更容易分解和起爆。综上，由于最大引发键键长和体系的内聚能密度与热感度之间的关联性较强，所以可以将其作为“判断体系热感度”的参考依据。

表6 体系PEG/NG/TMETN-1 在不同温度下的引发键键长及其概率Table 6 Length and probability of O—NO2 trigger bond at different temperature（PEG/NG/TMETN-1）

表7 体系PEG/NG/TMETN-3 在不同温度下的引发键键长及其概率Table 7 Length and probability of O—NO2 trigger bond at different temperature（PEG/NG/TMETN-3）

表8 体系PEG/NG/TMETN-1 不同温度下体系的内聚能密度及其分量Table 8 Cohesive energy density and the components of systems at different temperature（PEG/NG/TMETN-1）

表9 体系PEG/NG/TMETN-3 不同温度下体系的内聚能密度及其分量Table 9 Cohesive energy density and the components of systems at different temperature（PEG/NG/TMETN-3）

3.3 相容性能

当增塑剂与粘合剂的相容性较差时，增塑剂容易在推进剂内部聚集、无法有效减缓粘合剂PEG 的结晶现象［32］，增塑剂还会向衬层、绝热层迁移，导致推进剂内部能量流失、衬层变得不耐烧蚀等［33-34］，所以增塑剂与粘合剂的相容性尤为重要。介观动力学模拟不仅可以得到体系的形态、密度分布情况，还可以通过“有序度参数值”比较“同一体系内各物质、不同体系内同一物质的聚集情况”，找到各物质随TMETN 和NG 配比变化的聚集规律、非常直观有效地比较各体系的相容性；经过1000 μs 的介观动力学模拟，最终得到五种体系中各物质的有序度参数随时间的变化曲线，如图7 所示。

Meso Dyn 是先将各粒子混合成为均匀体系后再开始模拟，由图7 可见，在模拟初期，体系内各粒子的有序度参数均从0 迅速增大、然后缓慢增大、最后稳定，值得注意的是，PEG/TMETN 体系的有序度参数迅速达到稳定、且数值很小，这说明该体系的熵变很小、容易达到新的动态平衡状态；粒子的有序度参数越大，其聚集能力越强、程度越大，进一步统计1000 μs 附近粒子的有序度参数，结果见表10。各粒子的有序度参数随TMETN 含量增加的变化如图8 所示。

表10 五种体系内粒子的有序度参数Table 10 Order parameters of the particles for five systems

图7 不同体系的有序度参数随模拟时间的变化Fig.7 Order parameters vs. simulation time for various systems

图8 各粒子的有序度参数随TMETN 含量增加的变化Fig.8 Order parameters of particles vs.the content of TMETN

由 表10 和 图8 可 见，在PEG/NG 体 系 中，PEG 的有序度参数大于NG 的，这说明此体系中PEG 的聚集能力强于NG 的；在混合增塑体系中，PEG 的聚集能力较NG、TMETN 的都更强；而在PEG/TMETN 体系中，TMETN 的聚集能力却略超过PEG 的，不过这两者的聚集能力都很弱。

随着体系内TMETN 含量的增加，PEG、NG 的有序度参数均减小，并且PEG 的值是先平缓减小、后剧烈减小，而NG 的值则是一直剧烈减小，这说明TMETN 有利于改善PEG、NG 的聚集情况，并且TMETN 的含量增加时，这种改善能力迅速变强；而TMETN 自身的有序度参数却先缓慢增加、完全替代NG 时又有所减小，这是因为TMETN 的含量增大、分子间作用增强、导致同相聚集的可能性变大，当完全替代NG 时体系的相容性最好、聚集能力最弱。此外，随着TMETN 含量的增加，增塑剂和PEG 的有序度参数变得更接近，这说明它们的聚集能力变得相当。

进一步观察各体系介观形态随时间的演变过程（图9），在图9 中，红色点为PEG 相、蓝色点为NG 相、绿色点位TMETN 相，显示颜色的地方代表着密度高于该物质平均密度的位置和水平。

由图9 可见，各体系在模拟初期的分布都是较为均匀的，但随着时间的延长，有的体系开始出现相分离；在同一时刻，PEG/NG 体系的相分离最为严重，如图9b 至图9e 所示。随着TMETN 含量的增加，体系在同一时刻的相分离程度减弱、分布变得越来越均匀；当TMETN 完全替代NG 时（见图9e），该体系未出现明显的相分离，这说明TMETN 的加入确实改善了体系的相容性；如图9b 至图9d 所示，在混合增塑体系中，由相分离的边界可见，相比于NG 相（蓝色点），TMETN相的边界（绿色点）更加“接近于”PEG 相的边界（红色点），这也说明了TMETN 与PEG 的相互作用更强、相容性更好；此外，在混合增塑体系中，两种增塑剂相混合得较为均匀，说明这两者的相容性也很好。

综合力学性能、感度和相容性，PEG/TMETN 体系的各方面性能均最佳，不过TMETN 的能量较NG 低，如果再考虑能量水平的话，混合增塑剂则综合了NG的高能量和TMETN 的优异性能，研究表明部分掺入TMETN 可以改善体系的感度和相容性、TMETN 含量较高时（PEG/NG/TMETN-3）体系的力学性能较好。

4 结论

采用分子模拟法，对体系PEG/NG、PEG/TMETN、PEG/NG/TMETN-1、2、3（TMETN 的含量递增）进行了力学性能、感度、相容性的模拟和分析，主要得到如下结论。

（1）在力学性能和加工性能方面，PEG/TMETN 体系最优，混合增塑体系较单一增塑体系略差，但当TMETN 含量较高（PEG/NG/TMETN-3）时，体系的力学性能得到了改善、并具有良好的加工性能。

（2）对于组成相似的三种混合增塑体系，氢键作用的强弱与力学模量的大小顺序一致，这说明氢键可以较好地反应混合增塑体系的力学性能，而对于组成不同的两种单一增塑体系，物质本身的性质、结构等因素对力学性能的影响更大。

（3）在感度方面，随着TMETN 含量的增加，体系内的最大引发键键长不断减小，即感度不断降低，PEG/TMETN 体系具有最低感度；最大引发键键长和体系的内聚能密度均与体系的热感度关联性较强，可将其作为反映体系热感度的参考依据；平均引发键键长可以有效地表示引发键整体的变化趋势、但无法明显地反应体系的感度。

（4）在相容性方面，通过对粒子有序度参数和介观体系密度分布的分析，发现在PEG/NG 体系及混合增塑体系中，PEG 相的聚集能力均大于增塑剂相的聚集能力、体系均出现相分离，随着TMETN 含量的增加，PEG、NG 各自的聚集情况迅速得到改善、体系相分离程度不断减弱，在PEG/TMETN 体系中，粒子的聚集能力最小、体系未出现明显的相分离。

（5）在设计的五种体系中，PEG/TMETN 体系的力学性能、延展性、感度、相容性均最优；如果考虑到TMETN 的能量较低，“混合增塑剂”综合了“NG 高能量和TMETN 优异性能”、具有较大地发展前景，在三种混合增塑体系中，PEG/NG/TMETN-3 体系具有最优的力学性能、感度及相容性。

图9 五种体系中密度分布随时间的演变过程Fig.9 Evolution of density distribution over time for five systems