PEG/增塑剂共混物相容性的分子动力学模拟和介观模拟

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

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

1 引言

增塑剂用于改善推进剂力学性能，并提供一定的能量，如果增塑剂和粘合剂的相容性不好，将会导致增塑剂迁移，引起推进剂力学性能变差、内部能量流失、界面脱粘、衬层绝热层不耐烧蚀［1-2］、甚至还会影响推进剂的安定性和感度［3］，所以，增塑剂与其他组分的相容性是评价推进剂性能的首要指标之一。随着计算机软硬件的发展，很多研究者采用分子模拟法探究混合体系的相容性［4-5］，此法较实验法简单易行［6］，还可以从微观角度研究相容性的本质，但多数研究集中于HTPB 推进剂［7-9］。

高能固体推进剂NEPE常用的增塑剂为硝化甘油（NG），因其感度较高导致推进剂的危险品等级为1.1级［10］，并且存在一定程度的增塑剂迁移现象［11］。有研究表明：含能钝感增塑剂TMETN可与NG起到类似的增塑作用，且TMETN感度更低［12］；法国已应用TMETN降低了双基推进剂的感度［13］；美国也已将低能硝酸酯增塑剂［14］用于战术助推的钝感NEPE推进剂；日本 Yoshio Oyumi［15］发现 TMETN增塑的GAP/AN推进剂具有良好力学性能，TMETN有助于提高燃速。而国内关于降低NEPE推进剂感度的研究还较少，特别是将TMETN应用于战术NEPE推进剂的研究则更少。

本文立足于以上背景，从相容性能的角度研究了TMETN应用于NEPE推进剂的可行性。采用分子动力学模拟得到溶度参数和径向分布函数，对聚乙二醇（PEG）与增塑剂NG、BTTN、TMETN的相容性进行了预测和比较，并通过结合能分析相容性的本质；然后把微观结果转化为介观模拟的输入参数，对比了NG/PEG、BTTN/PEG、TMETN/PEG三种体系的介观形态随时间演变过程，为低易损NEPE推进剂的研制提供依据。

2 模型构建与模拟方法

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

利用Material Studio中Visualizer模块建立分子量为6002的PEG分子及NG、BTTN、TMETN分子（分子结构式见图1），利用Forcite模块对分子进行几何优化，几何优化是为了寻找局部最优构型，优化过程均选择Compass力场［16］；然后进行若干次退火处理、直到分子的能量稳定在最小值，退火可以使分子跨过局部能垒、寻找全局最优构型，退火条件为初始温度为298 K，中间温度为 598 K［17］，每次退火后自动进行几何优化，这是在全局最优构型中进一步寻找，从中选取能量最低的一帧构型作为输出的最优单链构型，为下一步建立无定形结构做好准备。

图1 PEG、NG、BTTN、TMETN的分子结构式Fig.1 Chemical structures of PEG、NG、BTTN、TMETN

2.1.1 纯物质体系模拟细节

为减小尺寸效应，确保无定形分子模型中有1000个以上原子［7］，PEG、NG、BTTN、TMETN 无定形分子模型中包含的分子数量及物质的密度见表1。

表1 PEG、NG、BTTN、TMETN 4种纯物质无定形分子模型参数Table 1 The parameters of pure substance（PEG、NG、BTTN、TMETN）in amorphous cells

每种物质均采用Amorphous cell模块构建10个构型，从中筛选能量最低者，进行几何优化和若干次退火处理，为保证纯物质密度保持为设定值，采用微正则系综（NVE）［18］进行退火，退火后进行几何优化，从若干祯中筛选能量最低者进行分子动力学模拟［17］。

动力学模拟用Andersen控温法［19］，范德华和静电力的计算分别用 Atom-based［20］和 Ewald 方法，非键截断半径为 0.95 nm，时间步长 1 fs［7］，其它为默认，先进行300p正则系综（NVT）的动力学松弛计算，再进行400 ps的NVE系综动力学计算，后50 ps已经平衡、用于采集溶度参数及径向分布函数数据；此外，对PEG进行了1 ns的模拟，1 ns得到的溶度参数20.809与400 ps的结果20.941比较接近，证实了400 ps模拟时间确实已足够长。

2.1.2 混合体系模拟细节

在保证周期箱内原子数量大于1000的基础上，为了便于比较增塑剂分子与PEG分子间相互作用，将三种增塑剂分子的数量均取50个，周期箱的初始密度按体积比例进行加和［7］，NG/PEG、BTTN/PEG、TMETN/PEG三种混合体系的具体参数设置如表2，无定形分子模型（图2）的构造和优化与2.1.1节纯物质的类似。

表2 混合物无定形分子模型参数Table 2 The parameters of blends in amorphous cells

无定形分子模型要选取重复出现且能量较低的构型，退火后进行400 ps NVT系综的动力学松弛，再采用Berendsen控压方法［7］在等温等压系综中（NPT）进行400 ps的动力学计算，后50 ps用于分析性能。为了证明模拟的时间足够长，又对NG/PEG体系进行了1 ns NPT系综计算，不同模拟时间的径向分布函数结果如图3，因径向分布函数只是分析曲线的相对位置高低、进而判断分子间作用的相对大小，所以400 ps的结果与1 ns的结果很相近、可以用于性能的分析。

2.2 介观建模过程及模拟细节

介观模拟需要确定：每种物质的粗粒化高斯链段数量，及不同粗粒化链段间的相互作用参数。

高分子预聚物的高斯链段数量采用式（1）［8］：

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

采用Synthia模块计算预聚物PEG的C∞为4.98，粗粒化珠子的质量为219.12（用特征比乘以单体分子质量），所以需采用27个粗粒化珠子表示分子量为6002的PEG分子，根据各粗粒化珠子的质量相同，采用一个珠子代替一个增塑剂分子。

此外，不同珠子间的相互作用参数为：

图2 NG/PEG、BTTN/PEG、TMETN/PEG混合物的无定形分子模型Fig.2 Amorphous cells of the NG/PEG、BTTN/PEG、TMETN/PEG

式中，ν-1εij为 MesoDyn 相互作用参数 kJ·mol-1；χij为Flory-Huggins参数（无量纲）；R为气体常数8.314 J·（mol·K）-1，T取 298 K。

在确定分子粗粒化珠子数量及相互作用参数后，在Meso Dyn模块中进行模拟：根据增塑剂与粘合剂预聚物的质量比（增塑比）为3∶1，设置珠子的混合比例，三维周期箱尺寸为32.0 nm×32.0 nm×32.0 nm，温度为298 K，步长为20 ns，总时间为1000 μs，其它为默认。

3 结果与讨论

3.1 MD模拟结果分析

3.1.1 相容性的溶度参数判断

图3 不同模拟时长的NG/PEG混合物中的分子间径向分布函数Fig.3 Intermolecular radial distribution function for NG/PEG with different simulation time

根据高分子溶液理论［21］，恒温和恒压下的溶解过程自发进行的前提是：Gibbs自由能的变化小于零。两组分的溶度参数差值越小，越容易满足此条件、相容性越好，这就是物质的“相似相溶”原理，也可以通过几种物质的极性相近、分子间作用力类型和大小相近、各自的径向分布函数曲线相近，来评价相容性的优劣。

采用分子动力学模拟计算得到PEG和增塑剂的溶度参数结果见表3，δ为本次结算结果，为了便于比较，同时将文献［22-23］的分子动力学计算结果δ′和δ"、文献［24］采用 Hoftyzer-Krevelen（1976）方法的估算值δ′"均列于表3。由表3可见，本次模拟（δ）与文献模拟值（δ′）、（δ"）相近，但与估算值（δ′"）有一定差距，分析原因为所采用的计算方法不同，而本模拟是在相同力场和参数设置下进行的，所以可以考察溶度参数的相对值［4］；并且不同方法得到物质溶度参数的大小顺序是相同的，即通过计算Δδ判断共混物相容性优劣的结果是一致的，由表4可见，PEG与三种增塑剂的相容性优劣顺序为：TMETN/PEG＞BTTN/PEG＞NG/PEG。

四种物质的内聚能密度和溶度参数的分量结果见表5。分析表5结果发现，在分子间相互作用中，每种物质内聚能密度的范德华分量（CEDv）均大于静电分量（CEDe），即范德华作用均大于静电作用、占据主要地位，特别是对于聚乙二醇而言；因为非极性分子之间只存在范德华作用，所以范德华力代表着非极性作用、而静电力代表着极性作用［22］，由表5可见，PEG的极性较弱，TMETN与PEG的极性更相近、相容性更好。

表3 溶度参数模拟值与文献值Table 3 Simulated values and literature values of solubility parameters

表4 3种增塑剂与PEG的溶度参数的差值Table 4 The difference of solubility parameters between plasticizers and PEG （J·cm-3）1/2

表5 内聚能密度及其范德华分量、静电分量（CED=CEDv+CEDe）Table 5 Cohesive energy density and the van der Waals、electrostatic components（CED=CEDv+CEDe）J·m-3

此外，运用公式（3）可以将表4所示的溶度参数差值转换为Flory-Huggins参数（χ）：

式中，Δδ为两种物质溶度参数的差值，（J·cm-3）1/2；Vr为参比体积，取两种珠子摩尔体积的平均值，cm3；R为气体常数 8.314 J·（mol·K）-1；T为 298 K。三种增塑剂与PEG的Flory-Huggins参数结果列于表6，再利用公式（2）将其转换为介观动力学模拟的输入参数（ν-1εij），结果列于表 7。

3.1.2 相容性的径向分布函数判断

图4为四种纯物质的分子内、分子间径向分布函数，分子内径向分布函数可体现物质的有序度信息［25］：图 4a中NG-NG、BTTN-BTTN两条曲线进行了向上平移10个单位处理，目的在于方便观察，可见PEG与三种增塑剂均属于无定性结构（仅在0.3 nm以前出现峰）［26］。

表6 增塑剂/PEG的Flory-Huggins参数Table 6 Flory-Huggins parameters of plasticizer/PEG

表7 混合体系粗粒化珠子的相互作用参数Table 7 Exclusion parameters between different beads

图4 纯物质分子内和分子间径向分布函数Fig.4 Intramolecular and intermolecular radial distribution functions for pure substance

分子间径向分布函数可以反映分子间相互作用的类型（峰值位置）及大小（曲线高低）［27］，氢键作用范围为0.26～0.31 nm，范德华作用范围为0.31～0.50 nm。由图4b可见，四种纯物质内的分子间作用均为范德华作用，此外，根据“纯物质的分子间径向分布函数越接近，其混合后的相容性越好［25］”，也可得出体系相容性优劣顺序为：TMETN/PEG＞BTTN/PEG＞NG/PEG；通过分析三种增塑剂分子结构发现，分子内均含有3个硝基，相比于NG分子，BTTN分子中仅多了一个亚甲基，根据“结构决定性质”，可以说明亚甲基削弱了BTTN分子间的相互作用、降低了溶度参数值，而TMETN分子是一个中心碳原子周围连接四部分官能团，其类四面体结构导致分子间作用力最小、溶度参数与PEG最接近，相容性最好。

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

进一步分析混合物分子间的径向分布函数，由图5可见，增塑剂与PEG的相互作用方式均为VdW作用；有文献［28-30］记载，在混合物中，两种分子间的径向分布函数越是高于单一物质自身的分子间径向分布函数，混合物的相容性越好。在NG/PEG体系中，NG与NG分子间作用明显大于NG与PEG的分子间作用，所以体系容易发生相分离、相容性最差；在BTTN/PEG体系中，当r＜0.4 nm时，BTTN与BTTN的分子间作用略大于BTTN与PEG的分子间作用，当r＞0.4 nm时，前者明显大于后者，所以BTTN/PEG体系相容性较差；在TMETN/PEG体系中，TMETN与PEG分子间作用绝大多数大于TMETN与TMETN分子间作用，说明TMETN与PEG最难发生相分离、相容性最好，这与溶度参数的判断结果一致。

3.1.3 相容性的本质分析

相容的本质就是组分间的相互作用，结合能可以表示分子间相互作用的类型和大小［31］，其表达式为：

式中，EPEG/plasticizer为混合体系的总能量，Eplasticizer和EPEG分别为增塑剂和PEG的平均单点能。

混合体系内分子间结合能的结果见表8。由表8可见，在增塑剂/PEG体系中，两种分子间的相互作用仅来源于非键作用，在NG/PEG体系中范德华力约占61.64%，静电力约占38.36%，在BTTN/PEG体系中范德华力约占69.11%，静电力约占30.89%，在TMETN/PEG体系中范德华力约占72.67%，静电力约占27.33%，这说明混合物分子相容的本质为：“分子间主要通过范德华作用相结合，并且，增塑剂与PEG的相容性越好，分子间通过非极性作用（范德华作用）相结合的占比越大”。

3.2 Meso Dyn模拟结果分析

分子动力学模拟反映的是微观角度的信息，而高分子预聚物的相分离特征区尺寸和时间都远大于其模拟范围，下面进一步借助介观模拟，忽略原子尺度的信息，得到三种体系的等密度图在1000 μs内的变化规律，如图6所示，为了便于观察，图中只显示了增塑剂相。

由图6可见，在模拟初期，三种增塑剂都均匀分散在PEG相中，随着时间延长，NG/PEG体系中NG相开始聚集，最终形成较大的颗粒相分散于体系之中；相比之下，BTTN/PEG体系相分离发生时间更晚、程度更小；而TMETN/PEG体系中PEG仅发生轻微同相归并，未出现明显相分离。

表8 混合体系分子间的结合能及其源项Table 8 The binding energy between molecules and the components kJ·mol-1

图6 体系中增塑剂相NG、BTTN、TMETN的等密度图随时间变化过程Fig.6 The isosurface of the density fields for plasticizer NG，BTTN和TMETN with time evolution in blends

体系内硝酸酯的有序度参数随时间变化曲线如图 7所示。由图7可见，在模拟初期NG的有序度参数迅速增大，50 μs后增大速率变慢，说明体系的熵已明显变大、发生了相分离，在300 μs后有序度趋于不变，说明体系处于新的稳定状态；BTTN有序度参数增大的速率和程度均低于NG，在200 μs后处于新的稳定状态，说明BTTN发生相分离的程度低于NG；相比于前两者，TMETN的有序度参数基本稳定在0附近，体系熵变化很小、基本处于稳定状态。从体系自由能密度变化也可以看出：NG/PEG体系的自由能密度迅速下降；BTTN/PEG体系的自由能密度稍有下降；而TMETN/PEG体系的自由能密度基本稳定，可见TMETN/PEG体系的相容性最好，其次为BTTN/PEG，最差为NG/PEG体系。

图7 有序度参数和自由能密度随模拟时间变化图Fig.7 Order parameter and free energy density for plasticizer with time evolution

4 结论

通过对PEG/NG，PEG/BTTN，PEG/TMETN共混体系的分子动力学和介观动力学模拟，主要得到以下结论：

（1）分子动力学模拟PEG和三种增塑剂的溶度参数、纯物质及混合物的分子间径向分布函数，结果均表明体系相容性的优劣顺序为：TMETN/PEG＞BTTN/PEG＞NG/PEG，讨论其微观原因为：BTTN中的亚甲基以及TMETN的分子结构削弱了各自的分子间相互作用，减小了与PEG的溶度参数差值。

（2）对混合体系的分子间径向分布函数和结合能分析得到：PEG与三种增塑剂相容本质为“通过非键作用相结合，其中范德华作用占主要比重”，此外，增塑剂的极性越接近于PEG，体系通过范德华作用相结合的比重越大。

（3）增塑剂/PEG体系介观形态演变、等密度图和有序度参数图表明：TMETN/PEG体系未发生相分离、熵变最小，混溶性最好；BTTN/PEG体系发生较轻微相分离、混溶性较好；而NG/PEG体系发生相分离、混溶性较差。

（4）MD和MesoDyn模拟结果表明：含能钝感增塑剂TMETN与粘合剂预聚物PEG的相容性优于增塑剂NG和BTTN，同时TMETN的感度远低于NG，可以考虑将其代替或部分代替NG、用于改善NEPE推进剂的易损性，研制低易损战术武器装备。