潮湿环境下金红石型纳米TiO2掺杂聚乙烯的分子动力学模拟

张 涛， 周浩翰， 陈 敏， 孙金莉

（1. 三峡大学 湖北省输电线路工程技术研究中心，湖北 宜昌 443002；2. 国网湖北省电力有限公司电力科学研究院，湖北 武汉 430077）

0 引 言

聚乙烯（PE）具有电绝缘性高、耐环境应力开裂性能好、耐热性能优良等优点，目前作为电力系统高压电缆的主要材料广泛应用于高压电力传输和分配线路中。随着社会和经济高质量发展，工业用电和居民用电需求急剧增加，电力电缆的输送效率和安全保障变得极为重要，然而聚乙烯在运行过程中会遭受热老化、电老化、机械老化以及水树枝老化等因素的协同作用逐渐引发劣化，从而导致聚乙烯各方面性能下降。长时间的热老化以及水分子的扩散运动还会破环聚乙烯的内部结构，进而导致聚乙烯材料不耐热、容易产生变形[1-2]。

纳米技术的快速发展为解决这类问题提供了新的思路。T J LEWIS[3]在1994 年最先提出了纳米电介质的理论，随后纳米电介质被誉为“未来的绝缘材料”[4]。纳米材料具有独特的量子效应、尺寸效应和表面效应，能够全方面提升复合材料的性能[5-7]。其中金红石型纳米TiO2以其独特的化学性能、体积效应、量子尺寸效应和宏观量子隧道效应等新特性在一众金属氧化物纳米颗粒中脱颖而出，作为现今世界上晶体结构排列完好的填料之一，许多学者对TiO2/PE 复合材料的性能进行了研究。刘岩等[8]选用经过NDZ-101处理的TiO2和低密度聚乙烯（LDPE）共混，研究了TiO2不同填充量对LDPE 材料热性能、力学性能以及微观形貌的影响，结果表明，填充TiO2后复合材料的力学性能提高，同时提高了LDPE 的热稳定性，促进了其结晶，但降低了其结晶程度。宫斌等[9]研究了纳米TiO2掺杂对低密度聚乙烯空间电荷行为的影响，结果表明，掺杂纳米TiO2改变了聚乙烯空间电荷的分布和聚合物的微观形貌，使负电荷的注入增强，还产生了大量的深陷阱，有利于阻止电树枝的引发和生长。但是因为纳米粒子尺寸极小和纳米界面的重要性[10]，许多实验无法从宏观层面去观察纳米粒子的作用机理，而分子模拟技术在微观层面观察纳米粒子有着独特的优势与潜力[11]。

凡朋等[12]采用分子模拟技术研究了SiO2/LDPE以及水分子复合体系的微观特性，结果表明，复合体系的玻璃化转变温度随着水分的加入而降低，纳米SiO2的加入增强了复合模型的热稳定性能。张艳芳等[13]采用分子动力学模拟研究了PE/PUF 复合材料中水分的扩散，结果表明，掺杂的聚氨酯（PUF）浓度增大，会导致复合材料的自由体积分数减小，尺寸密度增大，且加快水分的扩散。然而对于纳米TiO2掺杂聚乙烯后，温度和水分子的协同作用对复合体系的微观影响鲜有研究。

本文基于分子动力学（MD）模拟方法，将复合体系的玻璃化转变温度、力学性能、分子链运动以及水分子扩散系数作为核心，研究不同水分子质量分数和不同温度下聚乙烯金红石型纳米TiO2复合体系的微观性能，希望研究结果可为聚乙烯的热力学性能研究和水树枝现象的抑制等方面提供参考。

1 模型构建与模拟方法

聚乙烯的分子结构是两相结构，由结晶区和无定形区组成。结晶区的微晶体排列规整又紧密，相互之间还存在较大的结合力，在宏观性能上表现出较高的机械强度，相较之下无定形区中的微晶体不规则聚集排列，容易被破坏产生劣化，因此构建模型时优先考虑无定形区。通常而言，聚乙烯因具有疏水性可以很好地隔绝水分子渗透到分子内部，但在潮湿环境下，受电缆运行中电场、温度场和水分子扩散能力的影响，水分子会经过无定形区进入分子内部。此外，聚乙烯材料在生产过程中难免会混入杂质，且老化作用也影响着聚乙烯分子内部结构，这些因素会导致电缆故障损坏。因此本文在构建模型时，为了对比研究金红石型纳米TiO2和水分子对聚乙烯微观性能的影响，分别建立不同类型的复合体系。

1.1 模型的构建

根据模型需求，综合考虑实际状况和模拟时间，本文采用Materials Studio 20.1 软件构建了纯聚乙烯（PE）模型、TiO2/PE 模型、不同水分子质量分数的PE 模型和TiO2/PE 模型。文献[14]研究表明国内500 kV 电力电缆绝缘的基料为低密度聚乙烯，但LDPE 在合成中由于工艺的差异会产生不同的分子链结构、分子量分布等，进而影响性能。为了使LDPE 适用于500 kV 电缆绝缘材料的生产应用，需要减少聚乙烯的链端作用，使分子量多集中于低分子量部分，并且分布更为集中，因此先使用软件的Build 模块构建一条聚合度为150 的聚乙烯单链，再使用软件的Amorphous Cell 模块对4 条相同的单链进行无定形区的复合体系模型构建。文献[15]研究了交联聚乙烯电缆中TiO2/LDPE纳米复合材料的空间电荷特性，其选取的纳米TiO2粒径为25 nm，掺杂的质量分数分别为0.5%、1.0%、3.0%、5.0%。考虑到整体模型的尺寸和纳米TiO2在模型中所占的质量分数和自由体积分数，本研究将纳米TiO2的半径设置为5 Å（质量分数为10%）。在实际中，纳米TiO2表面存在不饱和的化学键，会与空气中的水分子和氧气发生反应来达到稳定状态。为了模拟实际状况，需要对纳米TiO2进行表面处理，即对不饱和的Ti原子添加氢氧根和对不饱和的O原子添加H原子。文献[16]指出聚乙烯材料内的水分子质量分数可达1%以上，因此考虑水分子在模型中所占的质量分数和自由体积分数，同时为了对比分析不同水分子质量分数下复合模型的微观特性变化，将模型中的水分子质量分数设置为0.5%和1%，一共构建6 个复合体系模型如图1 所示，依次为PE、PE+0.5%H2O、PE+1%H2O、TiO2/PE、TiO2/PE+0.5%H2O、TiO2/PE+1%H2O。文献[17]中高压直流电缆用低密度聚乙烯基料的密度均为0.92 g/cm3，而本文搭建的PE 模型密度为0.922 4 g/cm3，与实际基料密度的相对误差为0.26%，说明搭建的模型具有可行性。

1.2 模拟方法

在进行分子动力学模拟计算之前，要确保所构建的复合体系模型处于相对稳定的初始状态，使模型内部的空穴分布接近于实际状况[18]。首先，使用软件里Forcite模块中的Geometry Optimization功能面板对复合体系模型进行结构优化，然后在正则系综（NVT）下对复合体系模型进行5次从100～500 K的循环退火处理，过程中温度间隔为50 K，每个循环周期的模拟时间为10 ps。经过上述处理后，复合体系模型中产生的局部极小值的势阱会被消减，模型结构更加稳定。之后采用循环退火处理后能量值最小的体系模型来模拟复合材料的转变过程，即最小能量值体系模型先进行500 K、200 ps 的NVT系综的动力学模拟，再进行时间为200 ps 的等温等压系综（NPT）的动力学模拟获取此温度下的数据参数，随后以每次递减50 K 降低模型的热力学温度至100 K，在每一次递减后的温度下使用相同体系模型进行相同时间的NVT 系综和NPT 系综的动力学模拟。在此过程中，积分步长均设置为1 fs，压强设置为0.000 1 GPa，体系内各原子的数据参数每隔200 fs 采集一次。考虑到体系模型在不同温度下重新平衡的问题，选取模拟结果的后100 ps 数据参数用于分析。整个模拟方法流程图如图2所示。

图2 模拟方法流程图Fig.2 Flow chart of simulation method

整个模拟过程采用COMPASSⅢ力场，此力场能够准确模拟各种分子和聚合物的气相性质（结构、构象、振动等）和凝聚相性质（状态方程、内聚能等），已被广泛应用于有机分子和无机分子构成的复合体系[19]。如文献[20]采用Nose-Hoover 法对温度进行控制，文献[21]采用Berendsen 压浴法对压强进行控制，文献[22]采用Velocity Verlet 跳蛙算法对牛顿运动方程进行求解。对于非键之间的相互作用，可采用Atom Based 方法[23]对范德华作用力进行计算，采用Ewald法[24]对静电作用力进行计算。

2 模拟结果及分析

2.1 玻璃化转变温度

玻璃化转变是高分子聚合物固有的性质，是其形态转变的宏观体现，它直接影响到材料的热力学性能，因此玻璃化转变温度（Tg）是高分子物理研究的热点内容。在分子动力学模拟中确定体系模型Tg值的方法是比体积-温度关系法[25]，其中比体积为密度的倒数。在相同的模拟方法和设置条件下，模型的分子链会随着模拟温度的变化产生运动变化，从而导致比体积-温度关系曲线呈现不同的趋势。研究通过比体积-温度关系法求得模型的玻璃化转变温度。待模拟计算完成后，取每一温度下模型的平均密度，然后以该密度的倒数即比体积作为纵坐标，以对应的温度点作为横坐标，便可获得比体积-温度关系曲线图。再对此曲线进行线性拟合，得到曲线的拐点即为该体系模型的Tg值，如图3 和图4所示。从图3 和图4 可以看出，水质量分数分别为0、0.5%、1.0%时，PE 复合体系的Tg值分别为252、238、218 K，而TiO2/PE 复合体系的Tg值分别为285、277、264 K，掺杂纳米TiO2提升了PE 体系的Tg值。文献[26]通过Materials Studio软件模拟得到PE的Tg值为255 K，而本研究得到的Tg值与之相差3 K，这是由于本软件构建的复合模型具有随机性，模型在构建时有不同的初始分布形态，同时参数设置的不同与模拟过程是否彻底也会影响模拟结果，所以两者间存在误差是可以接受的。

图3 PE复合体系比体积-温度曲线Fig.3 The specific volume-temperature curves of PE composite system

图4 TiO2/PE复合体系比体积-温度曲线Fig.4 The specific volume-temperature curves of PE+TiO2 composite system

温度升高，所有复合体系的比体积会随之增大，在Tg值附近出现转折，即模型正在从玻璃态转向高弹态，分子链开始运动。水分子质量分数提高，PE 复合体系的Tg值随之降低，而掺杂TiO2后TiO2/PE复合体系的Tg值相较于PE复合体系的Tg值有了明显的提升。这是因为水分子在复合体系中占有体积较小，受外界因素影响较大但不与聚乙烯链产生氢键效应，容易产生位移使分子间的间隔变大，从而为聚乙烯分子链提供活动的空间，进而降低了复合体系的玻璃化转变温度；掺杂TiO2后，纳米TiO2的表面活性效应增加了复合体系的热稳定性，从而提高了复合体系的玻璃化转变温度。玻璃化转变温度的变化反映了材料的热稳定性，由此可知，水分子质量分数的提高会降低复合体系的热稳定性，而掺杂金红石型纳米TiO2会有效提升复合体系的热稳定性。

2.2 力学性能

在分子动力学模拟计算中，可以采用静态常应变法计算体系模型的力学性能，即对平衡后的体系施加微小的应力，使其产生形变，从而获得应力-应变的关系。通常应力-应变关系满足胡克定律，如式（1）所示。

式（1）中：σi为应力；Cij为弹性刚度常数；εj为应变。

因为聚乙烯接近于各向同性材料，所以可以将模拟中的复合体系视为各向同性材料[27]，因此刚度矩阵C可简化为式（2）。

式（2）中，λ和μ是Lamé 常数，其计算公式分别为式（3）、式（4）。

复合体系杨氏模量（E）、剪切模量（G）的计算公式分别为式（5）、式（6）。

由文献[28]可知，聚乙烯在-60℃下仍可保持良好的力学性能，但使用温度在80～110 ℃。因此对稳定的复合体系在150～350 K 温度下进行静态常应变模拟，得到不同温度下各体系的力学模量如图5～6 所示。从图5～6 可以看出，随着温度升高，所有复合体系的杨氏模量和剪切模量出现下降趋势，并且PE 复合体系在200～350 K 温度范围内的下降幅度增大，TiO2/PE 复合体系在250～350 K 的下降幅度增大。由于温度的升高会使分子链更加剧烈运动，导致原子间结合不稳定，从而降低了复合体系的力学性能。再则复合体系处在玻璃化转变温度区间时，正在从玻璃态转向高弹态，分子链运动随着温度升高更加剧烈，容易产生断裂，导致复合体系的力学性能在玻璃化转变温度区间及后续有效温度范围内大幅下降。随着水分子质量分数增加，所有复合体系的杨氏模量和剪切模量随之降低。这是水分子在复合体系中起到了增塑作用的缘故，水分子插入体系的分子之间，削弱了分子间的作用力，从而导致复合体系的力学性能下降。

图5 不同温度下各体系的杨氏模量Fig.5 Young's modulus of each system at different temperatures

图7为300 K 下PE 模型和TiO2/PE 模型的力学模量。从图7 可以看出，300 K 下，PE 模型的杨氏模量和剪切模量分别为2.317 3 GPa 和0.849 2 GPa，TiO2/PE 模型的杨氏模量和剪切模量分别为12.239 3 GPa 和3.724 8 GPa，杨氏模量和剪切模量分别提升了428.17%和338.62%。这是由于掺杂金红石型纳米TiO2后，纳米TiO2的超微尺寸使其表面活性增大，减少了分子链的张力，对复合体系的缺陷进行了有效的填补，并最大限度地抑制了复合体系内部的活性基团，从而大幅度提高了复合体系的力学性能。

图7 300 K温度下不同模型的力学性能Fig.7 Mechanical properties of different models at 300 K

2.3 自由体积占比

根据T G FOX 等[29]的自由体积理论，将材料内部分子之间的空穴称为自由体积（VF），被分子链所占据的空间称为占有体积（VO），自由体积为分子链的运动提供了活动空间，因此复合体系模型的热学性能和力学性能变化可以采用自由体积理论来解释。由于不同聚乙烯体系具有不同的体积参数，不能直接对各体系的自由体积进行比较，所以采用自由体积分数（αFFV）表示自由体积相对大小，其计算公式为式（7）。

本研究采用软件中的Atom Volumes & Surfaces功能面板测试在150～350 K 温度下稳定的复合体系的自由体积分数。表1 是不同温度下TiO2/PE+1%H2O 模型的自由体积分数，表2 是300 K 温度下所有复合体系的自由体积分数。

表1 不同温度下TiO2/PE+1%H2O模型自由体积占比Tab.1 Free volume fraction of TiO2/PE+1%H2O model at different temperature

表2 300 K温度下所有复合体系的自由体积分数Tab.2 Free volume fraction of all composite systems at 300 K

从表1可以看出，随着温度升高，复合体系的自由体积和自由体积分数增大，在250～350 K 温度范围内增大趋势上升。由于温度升高分子链运动增强，运动空间增大，增大趋势上升表示复合体系处在玻璃化转变范围内，分子链运动更加剧烈，使得自由体积增大，从而降低了复合体系的力学性能。从表2可以看出，随着水分子质量分数增加，复合体系的自由体积和自由体积分数增大。这是由于水分子的增塑作用，加之聚乙烯分子链间无氢键效应对分子链运动产生的空间阻碍很小，导致自由体积分数增大，因此水分子的加入会降低复合体系的玻璃化转变温度和力学性能。掺杂金红石型纳米TiO2后，复合体系的自由体积和自由体积分数降低。这种变化缘于纳米TiO2占据了复合体系中的一部分空间，阻碍了分子链的运动，导致复合体系自由体积下降，自由体积分数随之降低，从而改善了复合体系的热稳定性和力学性能。

2.4 分子链运动能力

聚乙烯体系中的原子总是处在不断运动的状态中，原子运动会导致聚乙烯分子链发生整体运动，从而影响聚乙烯复合材料的热学性能和力学性能。均方位移（αMSD）可以有效描述分子链的运动能力，其计算公式为式（8）。

式（8）中：ri(t)和ri(0)代表体系中任意原子i在t时刻和初始时刻的位移矢量；< >代表对组内所有原子平均。

对稳定的复合体系模型进行均方位移分析，得到不同温度下PE 模型和TiO2/PE 模型中PE 分子链的均方位移曲线，如图8所示。

图8 PE分子链的均方位移曲线Fig.8 MSD curves of PE molecular chain

从图8 可以看出，随着温度升高，两个模型中PE 分子链的均方位移增大，表明温度升高会加剧PE 分子链的运动。TiO2/PE 模型中PE 分子链的均方位移小于PE 模型中PE 分子链的均方位移，表明掺杂金红石型纳米TiO2后，纳米TiO2占据了复合体系的部分自由体积，减少了PE 分子链的运动空间，从而增强了复合体系的稳定性。此外，两个模型中PE 分子链的均方位移均出现突变现象，PE 模型中PE 分子链的均方位移在250～300 K 和400～450 K区间发生突变，而TiO2/PE 模型中PE 分子链的均方位移在300～400 K 和400～500 K 区间发生突变。第一次突变是由于复合体系从玻璃态转变为高弹态，分子链运动阻碍变小，运动空间增大，均方位移发生了突变。第二次突变可能是由于复合体系从高弹态转变为粘流态，分子链运动趋近于流体状态，导致均方位移发生突变。这种突变现象验证了分子链运动与玻璃化转变过程的相关性。

2.5 水分子扩散

水分子在复合体系中可以自由扩散，其扩散能力可以通过扩散系数进行描述。在分子动力学模拟计算中，水分子的扩散系数（D）可由Einstein 关系式[30]求解，如式（9）所示。

式（9）中，N代表扩散的原子数。

因为模拟的时间足够长，且模拟过程中已经对N个扩散水分子的αMSD进行了平均操作，所以水分子的扩散系数计算公式可以简化为式（10）。

式（10）中，k为对水分子的均方位移曲线图进行最小二乘线性拟合得到的斜率。

以TiO2/PE+1%H2O 模型为例，作出在不同温度下水分子的均方位移曲线图，如图9 所示。从图9可以看出，随着温度升高，水分子的αMSD值和曲线斜率均上升，100～350 K 温度区间内水分子的扩散运动比较稳定，而400～500 K 温度区间内水分子的扩散运动具有一定的随机性。

图9 TiO2/PE+1%H2O模型水分子的均方位移曲线Fig.9 MSD curves of water molecules in TiO2/PE+1%H2O model

通过均方位移曲线图求出不同温度下各模型水分子的扩散系数，如表3 所示。从表3 可以看出，水分子的运动随着温度升高而加剧，各模型水分子的扩散系数逐渐增大。水分子质量分数的增加，也会使各模型水分子的扩散系数增大，这是由于水分子的增塑作用造成的。而掺杂金红石型纳米TiO2后，各模型水分子的扩散系数大幅降低，主要有两方面的原因：一方面，纳米TiO2占据了复合体系的部分自由体积，阻碍了水分子在其中的运动；另一方面，纳米TiO2与水分子可以形成氢键，减缓了水分子的运动，从而降低了水分子的扩散系数。另外，从表3 还可以看出，PE+H2O 复合体系和TiO2/PE+H2O 复合体系中水分子在温度达到350 K 后扩散运动都变得更加剧烈，这个变化也验证了玻璃化转变温度的有效性。

表3 不同温度下各模型水分子的扩散系数Tab.3 Diffusion coefficients of water molecule in each model at different temperatures

2.6 径向分布函数

径向分布函数（radial distribution function,RDF）是描述粒子密度距参考原子距离的函数，其可以描述体系中的氢键效应。图10 中定义黑色小球为参考原子，白色小球为粒子，假设与黑色小球中心距离由r增至（r+dr）的粒子数目为dN，则径向分布函数g(r)表示为式（11）。

图10 径向分布示意图Fig.10 Radial distribution diagram

式（11）中：ρ为体系中的粒子密度；N为体系中的粒子数目。

分子的有序度可由分子内的径向分布函数曲线图表示，分子之间的相互作用情况可由分子之间的径向分布函数曲线图表示。图11 为TiO2/PE+1%H2O 模型中纳米TiO2与水分子之间的径向分布函数曲线图。从图11可以看出，在纳米TiO2与水分子之间的分子间距超过5 Å 之后曲线趋于平缓，因此只取5 Å之前的数据进行分析。曲线图在2～3 Å出现了明显的峰值，这个范围处于氢键作用力的范围，说明在这个范围内纳米TiO2与水分子之间发生反应形成了氢键，进一步验证了纳米TiO2与水分子的氢键效应。另外在4～5 Å 出现了细微的峰值，这个范围处于范德华作用力的范围，说明在这个范围内纳米TiO2与水分子之间也发生了反应，但由于范德华作用力相对于氢键作用力较弱，在分析纳米TiO2抑制水分子扩散运动时可以忽略不计。

图11 纳米TiO2与水分子之间的RDF曲线图Fig.11 RDF curve between nano TiO2 and water molecules

图12为TiO2/PE+1%H2O 模型中PE 与水分子之间的径向分布函数曲线图。由图12 可知，曲线在3～10 Å 之间没有出现明显的峰值，通过与图11 对比，可以说明在3～10 Å 范围内PE 与水分子之间没有产生氢键效应，也没有范德华作用力的影响。

图12 PE与水分子之间的RDF曲线图Fig.12 RDF curve between PE and water molecules

3 结 论

（1）由比体积-温度关系法确定了各体系的玻璃化转变温度，PE、PE+0.5%H2O、PE+1%H2O、TiO2/PE、TiO2/PE+0.5%H2O、TiO2/PE+1%H2O 各体系Tg值分别为252、238、218、285、277、264 K。掺杂纳米TiO2提升了聚乙烯材料的热稳定性。水分子的加入降低了各体系的Tg值，且随着水分子质量分数的增加，下降幅度也越大，导致聚乙烯材料的热稳定性降低。

（2）随着水分子质量分数的增加，各体系的杨氏模量和剪切模量下降。随着温度升高，各体系的力学模量下降，即水分子加入和温度升高均会降低聚乙烯材料的力学性能。而掺杂纳米TiO2将PE 体系的杨氏模量和剪切模量分别提升了428.17%和338.62%，使聚乙烯材料的力学性能有效提高。

（3）随着温度升高，PE 分子链的均方位移显著增大，且在有效范围内发生两次突变，表明温度升高加剧了PE分子链的运动，并与玻璃化转变保持一致性，破坏了聚乙烯材料的稳定性。掺杂纳米TiO2后PE 分子链的均方位移比未掺杂时低，表明纳米TiO2抑制了PE分子链的运动，增强了聚乙烯材料的稳定性。更优异的热力学稳定性使复合材料具有应用于更恶劣环境的可能性。

（4）纳米TiO2的掺杂会占据复合体系的自由体积，加之纳米TiO2与水分子之间具有氢键效应，致使水分子的扩散系数骤降。表明掺杂金红石型纳米TiO2在提高PE 材料耐水树枝能力和使用寿命上起着主导作用。