重油体系黏度特性的分子动力学研究

1. 研究背景

润滑油基础油被广泛应用于机械加工领域，其化学成分包括高沸点、高分子量烃类和非烃类混合物，分子中碳原子数量一般在20-40之间。这些重质液体对于核燃料循环设施至关重要，可能会影响当地的中性环境。此外，分子结构的变化会影响重油的动态特性，使建模过程具有挑战性。在这项工作中，作者构建了特殊的石蜡油分子模型，使用MedeA-LAMMPS材料设计平台构建了分子集合，并计算了体系的静态和动态性质(密度，粘度，粘度指数以及扩散系数)。

2. 结果与讨论

2.1. 模型构建

如图1 (a)(b)所示，作者利用Amorphous Materials Builder模块，将四种类型的烷烃分子链按其在重油中的真实比例分布构建了重油体系的超胞模型。

图1 (a) 四种烷烃分子链(顺时针)：线性支链烷烃、环烷烃、芳香烷烃和二环烷烃；(b) 重油体系的超胞模型

2.2. 力场选择

作者选用了专门描述有机小分子以及聚合物运动的COMPASS力场来研究重油体的分子动力学过程。通过Forcefield—Read/Choose将力场导入并匹配到结构中。如图2所示，为在COMPASS力场中描述有机分子骨架振动的各种键合相互作用，包括键伸缩振动能、键角振动能、二面角扭转能以及面外振动能。

图2  COMPASS力场中的相互作用能项

2.3. 平衡模拟

作者将构建好的重油体系超胞模型在313K温度下做NPT系综模拟，如图3所示，在400ps后，体系的温度、压力以及能量到基本收敛，到达平衡状态。如表示1所示，体系在温度293K的平衡密度为0.86g/cm3，与实验值的误差在5%以内。

图3 重油体系在313K下的NPT平衡模拟

表1 体系的平衡密度

2.4. 非平衡模拟

作者通过平衡模拟得到了体系的静态性质(密度)，对于体系动态性质(粘度)的计算，作者使用了非平衡模拟方法，在体系中施加了剪切应力用于计算各种流体运输性质。如图4所示，体系被施加恒定的剪切应变速率，在MedeA-LAMMPS-Viscosity模块中利用随时间变化的剪切应力来计算粘度。粘度由公式(1)给出。

    (1)

图4 体系对剪切应力τ的响应，u为剪切速率，y为层的垂直方向。

如图5(a)所示，利用上述方法，作者得到了313K和373K温度下粘度随剪切速率的变化关系。并通过粘度-剪切速率曲线外推得到零剪切粘度(外推线与y轴的交点截距)，表明温度升高会降低重油体系的粘度大小。如图5(b)所示，将313K和373K温度下计算的粘度与实验值对比，并分析了误差，发现误差在10%以内。表明这种非平衡方法计算粘度的鲁棒性。

图5 (a) 313K和373K下重油体系的粘度随剪切速率的变化图；(b) 313K和373K温度下粘度计算值与实验值对比的误差分析

如图6所示，作者最后利用MedeA-LAMMPS-Diffusion模块，通过计算两个氢离子均方位移(MSD)研究了氢气分子在重油体系以及水中的扩散行为，得到了氢气分别在重油体系和水中的扩散系数，见表2。可以看出氢气在水中的扩散系数比在重油中的扩散系数高两个数量级。与粘度计算的结果一致。

图6 氢原子扩散的均方位移(MSD)

表2 氢气在水和重油中的扩散系数

3. 结果与展望

在这项工作中，建立了石蜡油的分子模型并用于分子动力学模拟，并计算了体系的静态性质和动态性质(密度、粘度、扩散系数)，通过与实验值的对比以及误差分析，表明了分子动力学模拟用于重油运输性质计算的鲁棒性。

参考文献：

Manring C A, Hawaii A I. Assessment of thermal neutron scattering in a heavy paraffinic molecular material [J]. Annals of Nuclear Energy, 2019, 128(JUN.): 140-147.

使用MedeA模块：

• MedeA-Environment

• MedeA-LAMMPS

• Amorphous Materials Builder

• MedeA-Forcefield

• MedeA-LAMMPS-Viscosity

• MedeA-LAMMPS-Diffusion