1 引言

在工程应用中, 硝酸酯增塑的聚氨酯弹性体存在硝酸酯增塑剂的迁移现象, 直接影响着含能材料的使用性能, 如硝酸酯增塑聚醚推进剂(NEPE)的能量性能、力学性能、燃烧性能、贮存性能及推进剂与衬层的界面粘接性能^[1]。一般情况下, 采用浸渍法、气相色谱、高效液相色谱仪和等离子发射光谱等手段分析推进剂和衬层中硝酸酯增塑剂的含量, 可以了解硝酸酯增塑剂在推进剂和衬层中的迁移情况^[2-4], 但实验测定方法需要较长的周期, 且具有一定的危险性。

在进行高分子材料的应用研究时, 运用分子模拟技术不但可以观察到实验难以测定的现象, 还能缩短研究周期, 降低研究成本^[5]。针对固体推进剂所用粘合剂中增塑剂的迁移和扩散现象, 微观尺度上的分子动力学模拟计算获得了较好的结果, 如端羟基聚丁二烯中非极性增塑剂的迁移和扩散^[6-7]、聚氨酯中硝化甘油的扩散性能^[8]等。然而, 分子动力学模拟计算方法仅适用于含有数万个原子的体系, 可观察的时间范围也仅为皮秒至纳秒, 在研究实际问题时难以反映其全貌。最近十几年蓬勃发展起来的耗散粒子动力学DPD(dissipative particle dynamics)属于介观模拟方法, 其模拟体系最多可以包含10⁶个原子, 在微观尺度的快速动力学和宏观尺度的慢速热力学之间架起了联系的桥梁, 能研究微秒范围内聚合物-溶剂体系的相行为^[9-10]。但在含能材料领域, 有关DPD研究的文献报道较少, 如NEPE推进剂粘合剂与增塑剂间相分离^[11]、TATB基高聚物粘结炸药(PBX)的介观结构及其演变过程^[12]等, 而有关固体推进剂或混合炸药所用粘合剂与增塑剂迁移和扩散的DPD研究尚未见到公开的文献报道。

鉴于此, 本研究运用分子模拟软件Materials Studio, 针对聚氨酯弹性体中丙三醇三硝酸酯(NG)、1, 2, 4-丁三醇三硝酸酯(BTTN)的迁移问题进行DPD模拟, 并将介观模拟计算结果与文献报道的微观模拟计算结果相对照, 检验其准确性, 并在更大的时间、空间尺度上探索预聚物分子量、温度和硝酸酯含量对NG和BTTN在聚氨酯弹性体中扩散行为的影响, 为新型高能固体推进剂和新型混合炸药的研制提供参考。

2 分子模型的建立及分子动力学计算

分子模拟计算面临的首要问题便是如何搭建合适的分子模型。在微观尺度上进行分子动力学计算时, 搭建聚合物的三维交联网络已经有了较为成熟的方法^[13], 而介观尺度下的交联网络模型搭建方法尚未见到公开文献报道。本研究参考DPD相关文献^[14], 进行一些探索。

本研究所用预聚物端羟基聚氧化乙烯(PEG)、固化剂缩二脲三异氰酸酯(N-100)及聚氨酯的分子结构式如Scheme 1所示, NG和BTTN的分子结构式如Scheme 2所示。

按照文献[14]的方法, 将聚氨酯分子的结构式划分为三种片段, 如Scheme 3所示。

在Materials Studio的Visualizer模块中分别搭建Scheme 3所示的片段A、B(PEG分子量取5500)、C和两种硝酸酯增塑剂NG、BTTN的分子模型, 进行几何优化。具体条件为: COMPASS力场, Materials Studio内置的Smart minimizer优化方法, 范德华作用选择Atom Based, 静电作用选择Ewald, 位能计算采用球形截断法, cutoff取12.5 Å, spline width取1 Å, buffer width取0.5 Å。

随后, 在Amorphous Cell模块中搭建其三维周期边界模型。以片段B为例, 其三维周期箱的尺寸为20.28 Å×20.28 Å×20.28 Å, 箱内共有878个原子(图 1)。在Forcite Plus模块中对优化并弛豫的无定形体系进行分子动力学模拟计算, 以恒温恒压系综(NPT)、Anderson恒压器、Berendson恒温器、温度298 K、压力10^-4 GPa、时间步长1 fs, 对优化并弛豫的无定形体系进行200 ps的密度优化, 使之接近体系的真实密度。然后, 选取恒温恒容系综(NVT), 按照上述条件, 在真空下进行总模拟时间为200 ps的分子动力学计算, 每0.1 ps记录一次全轨迹文件。

3 结果与讨论 3.1 微观模拟

对图 1所示的三维周期箱进行分子动力学计算, 所得势能-时间曲线和温度-时间曲线如图 2所示。由图 2可以看出, 在100~200 ps的区间内, 其势能和温度波动较小, 说明体系达到了平衡。分析相对应的轨迹文件, 可以采集各个组分的溶度参数(δ)和密度(ρ)(见表 1)。此外, 由分子量(M_n)和密度(ρ), 可以算出摩尔体积(V)、各组分与NG摩尔体积的比值(R)。这些数据也一并列入表 1中。

3.2 粗粒化的介观模型

为了进行介观模拟(DPD)计算, 必须对上述微观分子模型进行粗粒化处理。按照参考文献[14]的方法, 将1个NG分子视为1个DPD珠子, 根据其余各部分的摩尔体积比值(即表 1中的R), 可以将Scheme 3中的片段A、B、C分别视为1、35、1个DPD珠子, 1个BTTN分子也视为1个DPD珠子。由此, 即可搭建单条聚氨酯分子链的粗粒化模型(见图 3a, 其中的蓝色珠子为片段A, 棕色珠子为片段B, 灰色珠子为片段C)。通过粗粒化处理, 一个NG、BTTN分子分别由20、23个原子缩减为一个DPD珠子, 而一条聚氨酯分子链则从图 4中的878个原子缩减为图 3a的108个DPD珠子。由表 1可知, NG的摩尔体积为140.48 cm³, 一个NG分子的体积即为233.36 ų。按照文献[15]报道的方法, 计算出一个DPD珠子的半径为:

$ {R_{\rm{c}}} = \sqrt[3]{{3 \times 233.36}} = 8.88Å $

(1)

随后, 根据实际体系中聚氨酯和硝酸酯的质量分数(NG和BTTN均为6%, 聚氨酯为88%), 可以搭建聚氨酯/NG/BTTN混合体系的三维周期边界条件模型(见图 3b, 其中的红色珠子为NG, 绿色珠子为BTTN), 三维周期箱的尺寸为100×100×100R_c³, 即888 Å×888 Å×888 Å, 箱内共有15000个DPD珠子。从图 3b可以看出, 多条聚氨酯分子链与NG、BTTN混合在一起, 构成了硝酸酯増塑聚氨酯体系的粗粒化介观模型。

3.3 介观模拟的相关参数

获得表 1中各个组分的溶度参数、摩尔体积后, 由式(2)可以得到组分之间的Flory-Huggins相互作用参数χ_ij^[16], 结果见表 2。

$ {\chi _{ij}} = \frac{{{{({\delta _i}-{\delta _j})}^2}{V_{\rm{r}}}}}{{{\rm{R}}T}} $

(2)

式中, δ_i、δ_j分别为组分i、j的溶度参数, (J·cm^-3)^0.5; V_r为参比体积(取组分i、j摩尔体积的平均值, cm³; R为气体常数, 8.314 J·(mol·K) ^-1, T为温度, 298 K。

由式(3), 可得到介观模拟计算所需DPD珠子之间的相互作用参数a_ij^[15], 结果见表 3。

$ {a_{ij}} = 25 + 3.50{\chi _{ij}} $

(3)

3.4 DPD计算

介观分子模型搭建好后, 根据表 3中的a_ij数据, 可在介观模拟计算模块Mesocite中构建相应的DPD力场, 并以该力场对图 3b所示的聚氨酯/硝酸酯混合体系的介观三维周期箱进行5000步的几何优化, 条件为: Smart minimizer优化方法, 范德华作用选择Bead Based, 静电作用选择Ewald, cutoff取9.24Å, spline width取0, buffer width取0.5 Å。

随后进行DPD动力学计算, 条件为:温度298 K、步长500 fs、总模拟时间50000 ps、每1 ps记录一次全轨迹文件。由势能-时间曲线(图 4a)和温度-时间曲线(图 4b)可以看出, 其势能和温度波动较小, 说明体系达到了平衡。在40000 ~50000 ps的区间内, 对全轨迹文件进行分析, 得到NG、BTTN在聚氨酯/硝酸酯无定形混合物中的MSD-t(均方位移-时间)曲线(图 5), 经Origin进行线性拟合后得到曲线的斜率。

通过爱因斯坦关系式, 可求得扩散系数D^[13]:

$ D = \mathop {{\rm{lim}}}\limits_{t \to \infty } \frac{{ < |r\left( t \right)-r\left( 0 \right)|{^2} > }}{{6t}} = \mathop {{\rm{lim}}}\limits_{t \to \infty } \frac{{MSD}}{{6t}} $

(4)

式中, r(t)代表t时刻粒子的坐标, 而r(0)为其初始坐标, 尖括号代表对所有原子和初始构型(或时间系综)进行平均, MSD即为均方位移。也就是说, 将介观模拟计算获得的MSD-t曲线的斜率除以6, 即可得到扩散系数D。

3.5 硝酸酯种类对介观扩散系数的影响

扩散系数是表征分子扩散、迁移能力的重要参数, 其含义是沿扩散方向, 单位时间内每单位浓度梯度下垂直通过单位面积所扩散某物质的质量或摩尔数^[13], 其单位为m²·s^-1。针对图 5中的MSD-t曲线(温度298 K、片段B的分子量5500、NG和BTTN的重量百分含量均为6%), 经Origin进行线性拟合后得到曲线的斜率, 由式(4)可算出NG和BTTN在聚氨酯中的扩散系数D, 见表 4。另外, 以NG/聚氨酯体系在微观条件下进行分子动力学计算所获得的扩散系数^[8]也列入表 4中。

表 4中的介观扩散系数与文献所报道相似体系微观尺度下分子动力学模拟计算所获得的扩散系数^[8]均处于10^-12 m²·s^-1这一数量级上, 证明其结果是可信的。

从表 4还可以看出, NG的扩散系数大于BTTN的扩散系数。由图 3a来看, 聚氨酯分子中片段A、C仅分别为2个、1个DPD珠子, 而片段B却为35个DPD珠子, 故片段B在聚氨酯分子与硝酸酯的相互作用中起着绝对的主导作用。从图 2可以看出, BTTN的分子中比NG多了一个亚甲基, 其极性弱于NG, 故BTTN的溶度参数值也小于NG而更加接近于片段B的溶度参数值(见表 1)。根据高分子物理的溶液理论^[16]可知, 若溶剂和聚合物的溶度参数分别表示为δ₁和δ₂, 则根据Hildebrand和Scatchard的论点, 混合热∆H_m可由下式得出

$ \Delta {H_m} = V{\left( {{\delta _1}-{\delta _2}} \right)^2}{\varphi _1}{\varphi _2} $

(5)

式中, V为溶液的总体积, $φ$₁、$φ$₂分别为溶剂和聚合物的体积分数。聚合物与溶剂混溶性主要取决于混合热, ∆H_m越小则二者越容易混溶。根据(5)式, 若δ₁和δ₂相差越小, ∆H_m就越小, 则溶剂与聚合物的混溶性越好。也就是说, 可以根据溶度参数差值∆δ来判断溶剂与聚合物的混溶性。实验结果也证实, 若溶度参数相近则聚合物与溶剂易于相溶^[17]。由于BTTN与片段B之间的溶度参数差值较小, 二者的混溶性比NG-片段B好, 故BTTN在聚氨酯中的迁移系数比NG小。而70 ℃加速老化条件下, 硝酸酯在NEPE推进剂中的含量测试结果表明, BTTN含量随贮存时间的变化幅度小于NG^[1], 说明BTTN在NEPE推进剂中的迁移量确实比NG小。

3.6 预聚物分子量对硝酸酯介观扩散系数的影响

预聚物的分子量对聚氨酯弹性体的性能有很大影响。为考察这一参数对NG和BTTN在聚氨酯中扩散系数的影响, 在温度298 K、NG和BTTN质量分数均为6%的条件下, 选择分子量不同的PEG预聚物(即片段B), 按照上述步骤搭建分子模型并进行DPD计算, 所得NG和BTTN扩散系数D的变化情况如图 6所示。

从图 6可以看出, 随着预聚物PEG分子量的提高, NG和BTTN在聚氨酯中的扩散系数均增大, 且预聚物分子量相同时BTTN的扩散系数依然小于NG的扩散系数。

按照前面提及的模拟方法及步骤进行微观尺度下的分子动力学计算后, 可以得到不同分子量预聚物(即片段B)的溶度参数, 如图 7所示。

由图 7可以看出, 随着分子量增大, 预聚物的溶度参数减小, 导致分子量增大后片段B与NG和BTTN的溶度参数差值增大, 由式(5)可知片段B与NG和BTTN的混溶性变差, 文献报道也证实了这一点^[11]。因此, 硝酸酯更加容易发生迁移, 故二者的扩散系数增大。而相同条件下NG扩散系数大于BTTN的原因, 也依然在于其溶度参数不同, 前已述及。

3.7 温度对硝酸酯介观扩散系数的影响

为研究温度对硝酸酯在聚氨酯中扩散系数的影响, 根据某硝酸酯增塑炸药的实际配方, 以片段B的分子量为5500、NG和BTTN的质量分数均为6%的条件, 分别在5个温度下进行了DPD计算, 所得NG和BTTN扩散系数D的变化情况如图 8所示。

由图 8可以看出, 随着温度的升高, NG和BTTN在聚氨酯中的扩散系数也增大, 且同一温度下NG的扩散系数仍然大于BTTN。硝酸酯在聚氨酯中的迁移, 可理解为硝酸酯与高聚物的相互扩散混合。根据热力学理论, 扩散混合的自由能∆F_m可以表示为

$ \Delta {F_{\rm{m}}} = \Delta {H_{\rm{m}}}-T\Delta S $

(6)

式中, ∆H_m为扩散体系的混合热, J·mol^-1; ∆S为混合熵, J·K^-1; T为绝对温度, K。扩散是混乱度增加的过程, 故∆S＞0。温度越高, ∆F_m越小, 扩散作用越容易进行^[16]。

此外, 有研究者认为, 温度升高后运动单元热运动能量随之升高, 聚合物链段间的活动性增加, 分子间距离增加, 聚合物中的自由体积膨胀、溶胀度升高, 使得增塑剂进行扩散的有效空间增大; 且温度升高后, 增塑剂分子的活动性增加, 加上粘接体系的粘度下降, 使得增塑剂的扩散行为更活跃^[7]。所以, NG和BTTN在聚氨酯中的扩散系数随温度升高而增大。

3.8 硝酸酯含量对其介观扩散系数的影响

在聚氨酯弹性体中提高硝酸酯的含量, 有助于提高固体推进剂或新型混合炸药的能量。为研究NG和BTTN的含量对其在聚氨酯弹性体中扩散系数的影响, 以温度298 K、片段B的分子量5500的条件, 分别在5种硝酸酯含量下(NG和BTTN的质量分数相同)进行了DPD模拟计算, 所得NG和BTTN扩散系数D的变化情况如图 9所示。

从图 9可以看出, NG和BTTN在聚氨酯中的扩散系数随其含量的提高而下降, 而含量相同时NG的扩散系数大于BTTN。这是因为在交联聚合物的溶胀过程中, 存在着2种相反趋势的平衡过程:增塑剂力图渗入高聚物内部, 使高聚物体积膨胀, 引起分子网交联点之间分子链的伸展, 降低高聚物的构象熵值; 构象熵值降低又必然引起分子网的弹性收缩力, 力图使分子网收缩。当弹性体中增塑剂含量超过共容增塑限量后, 增塑剂的增加对溶胀度的影响很小, 几乎保持不变; 过量的增塑剂占据大量空间, 不仅束缚分子网的伸展, 而且促使三维分子网络收缩, 使其自由体积减少, 降低了增塑剂的扩散活动^[7]。此外, 随NG和BTTN含量的增加, NG和BTTN分子间的相互作用增强, 也会降低NG和BTTN在聚氨酯弹性体内孔穴间的跳跃扩散^[8]。所以, NG和BTTN的扩散系数出现了下降的趋势。

4 结论

针对本研究所搭建的粗粒化模型, 经过模拟计算后可以得出以下结论。

(1) DPD计算的介观扩散系数与相似体系分子动力学模拟计算的微观扩散系数均处于10^-12 m²·s^-1这一数量级上, 证明其结果是可信的。

(2) 随着预聚物PEG分子量的提高, NG和BTTN在聚氨酯中的扩散系数增大。

(3) 随着温度的升高, NG和BTTN在聚氨酯中的扩散系数也增大。

(4) NG和BTTN在聚氨酯中的扩散系数随其含量的提高而下降。

(5) 同样条件下, NG的扩散系数大于BTTN。

(6) 介观分子模拟计算结果与文献报道值、理论趋势较为符合, 说明图 3的硝酸酯/聚氨酯粗粒化模型可以用于考察NG和BTTN在聚氨酯中扩散的情况, 这也许是因为迁移、扩散仅仅涉及NG和BTTN分子与聚氨酯分子之间的相溶性所致, 且有待于与实验测定数据进一步对比。如欲考察介观条件下该体系的力学性能等, 还需进一步探索。