定子橡胶在一定程度上影响着潜油螺杆泵采油系统中的使用寿命，而定子橡胶在其加工和使用过程中因受到温度和化学介质等作用，会逐渐出现粘连、变硬、龟裂等由于物理性能失衡导致的一系列现象.因此，在潜油螺杆泵行业中使用的橡胶制品应具有良好的耐高温、耐低温、耐油与耐臭氧老化性能.氢化丁腈橡胶(HNBR)是一种分子主链高度饱和的橡胶材料，主要通过丁腈橡胶(NBR)氢化制得[1]，即将NBR链段上的丁二烯单元进行选择性氢化，将C—C不饱和双键加氢生成单键，既保持了原有NBR的耐油、耐介质和耐磨性，又兼具优良的耐热、耐候、耐氧、耐化学介质腐蚀性、高拉伸强度和长期压缩永久变形较小的综合性能，特别是其热老化性能得以极大提升.HNBR分子链中的丙烯腈单元提供了材料优异的耐油性和高拉伸强度规格稳定性，可用于生产重复性优异的聚合物.然而，高玻璃化转变温度(Tg)和不良的机械性能阻碍了其应用.为了改善HNBR的Tg和机械性能，常用的方法是添加适量的合适增塑剂.普通增塑剂邻苯二甲酸二丁酯(DBP)[2]可以显着降低HNBR的分子间相互作用，从而可以降低Tg并改善机械性能.因此，有必要研究DBP对HNBR性能的影响.物理上稳定的分散体需要考虑由溶度参数(δ)判断的增塑剂与粘合剂之间的互溶性.聚合物的性能可以通过关键指标Tg进行评估，该指标的值决定工艺温度范围.根据Fox-Flory的自由体积理论，大多使用自由体积(VF)和密度(ρ)方法来计算Tg，即温度对应于密度随温度变化曲线的调整点.此外，机械性能对于聚合物的制备、加工和应用也很重要.
许多研究人员通过MD方法对聚合物的混溶性、玻璃化转变温度和力学性能进行了模拟，所得结果与实验结果具有很好的一致性，证实了模拟的可靠性[3].分子模拟可用来计算聚合物的物理和机械性能[4].
1 分子建模过程
构建HNBR单链[1]，由丙烯腈、顺式1，4-丁二烯、反式1，4-丁二烯、1，2-丁二烯和丁烷单体组成的重复单元结构并将其首尾连接，各基体组分如表1所示.
表1 单链各基体结构及配比
Tab.1 Structure and ratio of single chain in each matrix
名称比例结构丙烯腈0.32顺式1,4丁二烯0.18反式1,4-丁二烯0.391,2-丁二烯0.09丁烷0.02
单链重复数暂定50，构建后的HNBR单链构成及DBP分子模型分别如图1～2所示.对HNBR单链和DBP分子分别在Forcite模块中利用Smart Minimizer功能对能量进行最小化，当能量达到最低时，用Amorphous Cell模块分别建立了纯HNBR，以及DBP体积分数为0.1，0.15，0.2，0.25的HNBR/DBP二元体系的混合物模型.绘制出85HNBR/15DBP的AC模型，如图3所示.
图1 HNBR单链
Fig.1 HNBR single chain
图2 DBP分子模型
Fig.2 DBP molecular model
图3 85HNBR/15DBP的AC模型
Fig.3 AC model for 85HNBR/15DBP
通过分子力学软件Material Studio7.0的Forcite模块对混合物模型进行优化.对于DBP较小的分子结构，在273 K下进行优化可以得到最优的几何构型.然而，对于HNBR复杂的高分子结构，在273 K下优化得不到最优的几何构型，因此，将模拟温度设置为298 K.选用NVT系综、COMPASS力场[5]对高分子体系进行MD模拟，时间为500 ps，在0.1 MPa下进行退火处理，将温度从100 K升温至远高于HNBR玻璃化温度的500 K再降温至100 K，经过5次循环，体系在高温下经历多次松弛，得到能量最低的最优结构[6].
2 溶度参数的确定
溶度参数(δ)的定义为内聚能密度(CED)的平方根，已广泛用于研究材料的物理混溶性，它是指每单位体积蒸发一摩尔液体所需的能量，其数学表达式[7]为
(1)
=(δ1-δ2)2φ1φ2
(2)
式中：EV为汽化过程中内部能量的变化；HV为蒸发的焓；Vm为摩尔数在蒸发温度下液体的体积；R和T分别为气体常数和绝对温度；φ1、φ2分别为组分1、2的体积分数.由式(1)、(2)可得，δ1、δ2的值越接近，则ΔHV越小；ΔHV越小，体系相容性越好.因此，溶度参数差对各组分间的相容性预测有一定参考价值.
基本溶解度理论表明，分子溶质的溶解度与CED有关.可以通过MD模拟计算出的δ来预测混合材料的相溶性.对HNBR和DBP的平衡AC结构在NPT系综下进行500 ps的MD模拟，然后在NVT系综下进行再一次500 ps的MD模拟来计算溶度参数.在模拟过程中，均采用周期性边界条件.温度控制选择Andersen恒温器，压力控制则选择Berendsen恒压器.选用Verlet速度时间积分方法(步长为1 fs)，用于积分牛顿运动方程.范德华相互作用选取Atom Based，截断半径为19.5 Å.静电相互作用选择Ewald方法，精度为0.000 1 kcal/mol[8].模拟结果为不同聚合度下HNBR与DBP的溶度参数，如图4所示.表2为HNBR与DBP溶度参数范围.
图4 HNBR与DBP溶度参数
Fig.4 Solubility parameters of HNBR and DBP
表2 HNBR与DBP溶度参数范围
Tab.2 Ranges of solubility parameters for HNBR and DBP (J·cm-3)1/2
材料δexpδMDHNBR20.27～20.5320.4DBP19.73～19.8119.7
对于高分子体系，分子间如果不存在强极性基团，若溶度参数之差的绝对值小于1.3～2.1(J·cm-3)1/2，则混合体系相容[9].显然，本文中二者的溶度参数绝对值远小于1.3(J·cm-3)1/2，故二者相容.由图4可知，HNBR的溶度参数在聚合度达到30之后保持在20.4(J·cm-3)1/2，即HNBR单体链长设置为50合理.
3 体系相容性预测
径向分布函数g(r)是反映材料微观结构的特征物理量，它是在距离原子一定距离处找到另一个原子的概率与在完全均匀分布的情况下找到一个原子的概率之间的比[10].g(r)反映出原子排布的有序性问题.如果近程出现一个峰值，峰高且陡，表明有序性较高，属于晶体体系，原子紧密相连.如果出现小于0.3 nm的峰，则表示分子链有序性较低，属于非晶结构即无定型结构.分子间径向分布函数表示原子自由散落之间的相互作用方式，氢键作用范围为0.26～0.31 nm，范德华力作用范围为0.31～0.5 nm[11].
HNBR、DBP物质的分子间和分子内的径向分布函数如图5～6所示.由图5可知，HNBR分子在0.3 nm直到平稳过程中没有出现强峰，说明结构为无定型结构，分子间主要作用方式为范德华力.由图6可知，HNBR与DBP二者分子内作用力的径向分布函数曲线分布较为接近，依据相似相容原理，可以证明二者相容性较优.
图5 HNBR与DBP分子间作用力
Fig.5 Inter-molecular forces between HNBR and DBP
图6 HNBR与DBP共混物的分子内作用力
Fig.6 Intra-molecular forces of HNBR and DBP blend
4 玻璃化转变温度
玻璃化转变温度(Tg)是聚合物的重要物理性能，决定了其应用温度范围.获取Tg的方法是在NPT系综下设置平衡混合系统，然后使用脚本完成计算.本质上，玻璃化转变是聚合物材料中的常见现象.当聚合物经历玻璃化转变时，其各种物理性质都会发生剧烈变化.目前，还没有关于玻璃化转变的完整理论来解释其实验事实.在玻璃化转变中，摩尔体积根据温度变化而变化.Fox和Flory提出了自由体积理论来解释玻璃化转变，在Tg以下，随着温度的升高，聚合物的自由体积或空体积的变化很小，但却在玻璃化状态时自由体积出现突变[12]；共聚物的密度(ρ)随温度的拐点也是求得Tg的一种方法，本文根据不同物质含量的二元混合体系的摩尔体积随温度变化的图像得出其玻璃化转变温度.在NPT系综下将模拟温度设定在一定温度范围内(200～600 K).用最小二乘法分别调整玻璃化转变温度以上和以下的温度范围的数据点，得到两条直线，两条线的交点为玻璃化转变温度Tg.
在NPT下进行了500 ps的MD模拟，采用平衡体积分数为0、0.10、0.15、0.20和0.25的HNBR/DBP二元体系的AC结构.随后进行高温退火处理，使用NPT系综在大气压下以20 K的间隔将温度从600 K降至200 K.选择Andersen恒温器和Berendsen恒压器分别控制温度和压力.通过Ewald方法计算静电相互作用，并且通过基于原子的方法计算范德华相互作用.重复该步骤，直到温度远低于HNBR/DBP的Tg.从平均偏摩尔体积随温度的变化曲线中观察到的斜率转折点被认为是HNBR/DBP的Tg，计算结果如图7所示.
图7 HNBR与DBP共混物的玻璃化转变温度
Fig.7 Glass transition temperature of HNBR/DBP blend
在Tg以下，密度随温度升高而降低，因为分子链的能量不足以克服在一定温度范围内的旋转能垒.一旦温度超过一定范围，由于系统处于高能量状态，平均偏摩尔体积的提升速度会加快.此外，由于DBP降低了HNBR之间的相互作用，因此，Tg随着DBP体积分数的增加而降低.
5 力学性能
力学性能主要是由于聚合物系统在不同的弹性性能测试模式下表现出的性能.定义如下特征模量：杨氏模量(E)、体积模量(K)和剪切模量(G).材料的应力σ和应变张量ε之间的最一般关系遵循广义的胡克定律，其表达式为
σi=Cijεj
(3)
式中，C为弹性系数的对称(6×6)矩阵.矩阵的对称性将系数减小到21，以描述各向异性材料的应力和应变之间的关系，并且当材料中存在弹性对称表面时，数量将继续减少.基于罗伊斯均值方法，可以通过下式计算有效的体积模量和剪切模量，即
KR=[S11+S11+S11+2(S12+S23+S31)]-1
(4)
GR=5[4(S11+S22+S33)-4(S12+S23+S31)+
(S44+S55+S66)]
(5)
式中：下标R表示罗伊斯平均值；柔度矩阵S为弹性系数矩阵C的逆矩阵，即S=C-1，可以通过MD模拟来计算.将DBP体积分数为0、0.10、0.15、0.20、0.25的平衡HNBR/DBP的AC模型在NPT系综进行进一步的动力学计算，持续时间为500 ps，然后通过NVT系综的MD模拟执行500 ps的平衡结构计算HNBR/DBP的K和G.获得K和G后，计算拉伸模量(E)和泊松比(v)，即
E=2G(1+v)=3K(1-2v)
(6)
机械性能是聚合物最重要的性能之一，主要由聚合物增塑剂决定.拉伸模量(E)与拉伸弹性有关；体积模量(K)测量聚合物的不可压缩性；剪切模量(G)是指材料的刚度，即E、K的值表明HNBR/DBP在受到外力作用时易变形的程度.泊松比(v)反映了材料的横向变形.此外，K/G和柯西压力是描述材料延展性的另外两个主要参数.它们与材料承受应力时抵抗裂纹形成的能力有关.K/G表示塑性范围(拉伸伸长率)，因此K/G的高值与延展性相关，而K/G的低值与脆性相关.K/G值越高，延展性越好.对于柯西压力，负值表示脆性，而正值表示延展性.为了获得上述参数，在298 K下进行了HNBR的力学性能模量计算，结果如表3所示.为了研究DBP对HNBR力学性能的影响，在298 K下对DBP体积分数为0、0.10、0.15、0.20、和0.25的HNBR/DBP体系进行了力学性能模量分析，如表3所示.
表3 不同条件下HNBR的有效各向同性力学性能
Tab.3 Effective isotropic mechanical properties of HNBR under different conditions
体积分数%C11C12C13C22C23C33C44C55C6603.601.191.473.971.162.040.981.331.76102.731.031.192.871.082.870.911.291.58153.311.171.994.112.043.831.021.651.60204.051.311.144.622.252.590.901.331.62253.601.881.044.692.333.121.141.301.21
表4列出了298 K下HNBR/DBP系统的弹性系数和机械性能.从C12～C44的正值可以得出在298 K下的HNBR/DBP共混物是易延展材料的结论.此外，HNBR/DBP体系的泊松比均约为0.285，表明HNBR/DBP体系具有塑性体的性质(0.2<v<0.4).图8为HNBR机械模量.当DBP体积分数为0.10时，E、K、G值即相关机械模量值达到最低；当DBP体积分数在0.10～0.20时，虽然其E、K、G呈上升趋势，但其拉伸模量、体积模量及剪切模量均低于0%HNBR；但DBP体积分数大于0.20时，其体积模量K已超过0%HNBR.这种变化表明，HNBR/DBP系统的弹性模量E在DBP体积分数为0.1时增加，同时其体积模量K与剪切模量G的降低使其可塑性增强，刚性和脆性下降，机械性能得到改善.随着DBP体积分数的增加，K/G值增加.因此，添加一定量DBP可以提高HNBR的延展性.
表4 不同条件下HNBR的弹性系数
Tab.4 Elastic coefficients of HNBR under different conditions
体积分数%EGPaKGPaGGPaK/GC12～C44v02.931.991.171.700.210.25102.241.590.891.790.120.26152.631.881.041.810.150.28202.742.211.062.080.420.29252.562.450.972.530.740.31
图8 HNBR机械模量
Fig.8 Mechanical modulus of HNBR
6 结 论
本文通过分子动力学(MD)模拟预测HNBR/DBP二元体系的混溶性、玻璃化转变温度和机械性能，得出如下结论：通过MD模拟计算了HNBR和DBP的径向分布函数，以预测其相容性，表明HNBR和DBP是可混溶的；根据平均偏摩尔体积-温度图得出0、0.10、0.15、0.20和0.25的HNBR/DBP玻璃化转变温度变化曲线，结果表明，Tg随DBP体积分数的增加而降低；分析了弹性常数，例如体积模量、剪切模量、拉伸模量、柯西压力、K/G及泊松比v，结果表明，适当添加DBP可以提高HNBR的延展性并降低其脆性，可以促进HNBR的应用.
[1] 吴丝竹，胡涛，周浩，等.不同氢化度氢化丁腈橡胶的结构表征与分子模拟[J].高等学校化学学报，2013，34(4)：1027-1032.
(WU Si-zhu，HU Tao，ZHOU Hao，et al.Characte-rization and molecular simulation of HNBR with diffe-rent hydrogenation degree[J].Chemical Journal of Chinese Universities，2013，34(4)：1027-1032.)
[2] Li J，Jin S H，Lan G C，et al.Molecular dynamics si-mulations on miscibility，glass transition temperature and mechanical properties of PMMA/DBP binary system[J].Journal of Molecular Graphics and Modelling，2018，84：182-188.
[3] 潘艳秋，刘易杭，徐瑞松，等.PEK-C膜热解交联过程反应机理[J].沈阳工业大学学报，2020，42(1)：40-46.
(PAN Yan-qiu，LIU Yi-hang，XU Rui-song，et al.Reaction mechanism of pyrolysis crosslinking process for PEK-C based membrane[J].Journal of Shenyang University of Technology，2020，42(1)：40-46.)
[4] 尚海鑫，褚金奎，高佳丽.SU-8光刻胶玻璃化转变温度及力学性能的分子动力学模拟[J].高分子材料科学与工程，2014，30(4)：95-98.
(SHANG Hai-xin，CHU Jin-kui，GAO Jia-li.Molecular dynamics simulation on the mechanical properties and glass transition temperature of SU-8 photoresis[J].Polymer Materials Science & Engineering，2014，30(4)：95-98.)
[5] 陈正隆，徐为人，汤立达.分子模拟的理论与实践[M].北京：化学工业出版社，2007.
(CHEN Zheng-long，XU Wei-ren，TANG Li-da.Theory and practice of molecular simulation[M].Beijing：Chemical Industry Press，2007.)
[6] 马泽坤.三元复合驱采出水稳定性机理研究[D].北京：中国矿业大学，2018.
(MA Ze-kun.Study on the stability mechanism of produced water from ASP flooding[D].Beijing：China University of Mining and Technology，2018.)
[7] Erlebach A，Ott T，Otzen C，et al.Thermodynamic compatibility of actives encapsulated into PEG-PLA nanoparticles：in Silico predictions and experimental verification[J].Journal of Computational Chemistry，2016，37(24)：2220-2227.
[8] 张光磊，陈一泊，贾伟韬，等.SiO-2气凝胶纳米多孔结构的自组装过程模拟[J].人工晶体学报，2017，46(12)：2515-2520.
(ZHANG Guang-lei，CHEN Yi-bo，JIA Wei-tao，et al.Self-assembly simulation of nano-porous structure of silica aerogel[J].Journal of Synthetic Crystals，2017，46(12)：2515-2520.)
[9] 居学海，叶财超，徐司雨.含能材料的量子化学计算与分子动力学模拟综述[J].火炸药学报，2012，35(2)：1-8.
(JU Xue-hai，YE Cai-chao，XU Si-yu.Overview on quantum chemical computing and molecular dynamic simulations of energetic materials[J].Chinese Journal of Explosives & Propellants，2012，35(2)：1-8.)
[10] 杨俊清.高能有机叠氮增塑剂的设计与计算研究[D].南京：南京理工大学，2017.
(YANG Jun-qing.Design and computational study of energetic organic azido plasticizers[D].Nanjing：Nanjing University of Science and Technology，2017.)
[11] 孙会刚.水分对油纸绝缘热老化及寿命的影响与热老化程度表征研究[D].重庆：重庆大学，2011.
(SUN Hui-gang.Study on influence of moisture on thermal aging & life and thermal aging condition eva-luation of oil-paper insulation[D].Chongqing：Chongqing University，2011.)
[12] 付一政，刘亚青，张丽燕，等.不同构型聚丙烯的玻璃化转变温度的分子模拟[J].分子科学学报，2009，25(10)：53-56.
(FU Yi-zheng，LIU Ya-qing，ZHANG Li-yan，et al.Molecular simulation of the glass transition of diffe-rent configuration of polypropylene[J].Journal of Molecular Science，2009，25(10)：53-56.)