基于分子动力学模拟的过冷态水气液界面特性研究*

苗威 吕孝飞 赵会军

（1.常州大学 江苏省油气储运技术重点实验室，江苏 常州 213164；2.常州大学 石油工程学院，江苏 常州 213164）

0 引言

“十四五”规划纲要提出，力争2030年前实现碳达峰，2060年前实现碳中和，并强调加强天然气的开发利用。天然气中通常含有砂砾、水以及CO2和H2S等酸性气体，会对输运管道以及设备造成腐蚀破坏，因此在输送前通常会进行脱水脱硫等处理。冷凝分离是脱水工艺中一个重要的手段，超音速旋流分离器是一种新型的脱水装置［1］，水蒸气在超音速分离器的Laval喷管内膨胀降温，但越过饱和线时水蒸气并不立即凝结，而是继续膨胀，使温度继续降低，在过冷的温度下发生凝结。水蒸气从气态变为液态需要经过凝结成核与液滴生长2个阶段［2］，对这2个阶段进行理论分析时，液滴形成过程中的表面张力是其中非常重要的参数。

目前，过冷液体表面张力的研究方法有理论研究、实验研究以及分子模拟等。理论研究主要有Gibbs界面热力学理论、Laplace理论以及Van der Waals密度梯度理论。过冷区域的水分子极易形成冰晶，所以过冷区域的实验研究很少达到极低的温度。对于过冷区域的热力学性质尤其是拐点的存在仍存在争议。本文采用经典分子动力学方法对水的气液界面特性进行模拟，研究不同温度平衡状态下的密度分布，以及表面张力等界面性质。观察表面张力的温度变化趋势，在常温以及高温区域与IAPWS以及实验数据进行对比，并对深度过冷水表面张力的第2个拐点（SIP）进行讨论。

1 模拟方法

1.1 势能模型

水的分子势能模型的相关研究已经进行了几十年，SPC/E模型在表现热力性质方面与实验值符合较好，因此本文采用SPC/E模型进行模拟。2个水分子之间的相互作用势表示为短程Lennard-Jones势和长程静电Coulomb作用势之和：

表1 SPC/E模型势函数的分子参数

1.2 模拟细节

本次分子动力学（MD）模拟采用Lammps软件。模拟过程在30Å×90Å×30Å的盒子中进行，在中间30Å×30Å×30Å的盒子中创建1 000个分子，模拟单元如图1所示。

分子模拟采用NVT系统，三维周期性边界，系统包括1 000个分子，短程作用势LJ势的截断距离为10Å，长程库伦作用无截断。采用Ewald求和方法计算长程作用力能够有效地计算离子间的相互作用，且在周期性边界条件的模拟过程中可以体现长程作用势的效果。先对系统进行最小化处理，使系统趋于平衡，再采用Nose-Hoover控温方法将温度稳定控制在恒温，模拟步长设为1fs，模拟100万步，最终使体系达到该温度下的平衡状态。

表面张力通过计算系统压力张量获得

采用分子动力学方法计算系统内部的压力张量：

2 结果分析和讨论

2.1 密度分布

为了观察不同温度下过冷态水的气液界面以及表面层的形成，本文计算了不同温度下沿Y轴方向的密度分布。将系统从Y方向切为60个1.5Å厚度的薄片，每1 000步统计每个薄片的质量密度取均值并选用最后1 000步的结果作为密度分布。结果表明：密度在Y方向上有连续快速下降的区域。图2统计了240 K、270 K和300 K温度下Y轴方向密度快速下降区域附近的质量密度。

从图2可以看出，质量密度存在快速下降的区域，这意味着出现了一个气相与液相的界面。液相区域的密度在某个值附近持续振荡，这种振荡幅度在气液界面附近达到最大，在液相区域较小。

随着温度升高，液相分子与气相分子运动更加活跃，因此液相密度会随着温度的升高而减小。将每个温度下处在液相的薄片密度的均值作为密度振荡所围绕的值，统计值如图3。从图3可以看出，液相密度在280 K附近达到最大值，在270~280 K的温度区间内，密度随着温度升高而增加，而在其他的区域保持密度随着温度升高而减小的趋势。

2.2 过冷态水表面张力

为了验证本文模型的准确性，分别对过冷区域和常温区域的气液界面的表面张力进行了模拟，得到220~370 K温度范围内的表面张力如图4所示。

首先，将本文计算获得的常温区域（270~370 K）的表面张力数值和国际水和水蒸气协会IAPWS提出的表面张力计算公式进行比较，验证本文模型的准确性，IAPWS公式形式如下

从图4中可以看出，随着温度的不断减小，表面张力逐渐增大。在常温区，即温度在270~370 K区间时，本文模拟值和IAPWS公式以及VARGAFTIK N B［4］的实验结果比较一致；在温度为250 K左右的低温区域，本文模拟结果和VINS V［5］的实验结果也比较符合，验证了本文模型的准确性。

在220~250 K的过冷区域，本文模拟值稍大于实验值，而且表面张力的温度变化趋势发生改变，对温度的依赖性随着温度的降低而增强，与PT Hacker的观点一致，出现了第2个拐点（SIP），认为第2个拐点出现在250 K附近。

在220~250 K的过冷温度区域内对表面张力进行拟合，公式如下：

3 结论

本文采用分子动力学方法对220~370 K温度范围内的气液界面特性进行了模拟，研究了气液界面处的表面性质以及系统的密度分布，结论如下：

1）液相的密度围绕某个值在不断振荡，这样的振荡在液相幅度较小，在气液界面附近幅度较大。液相的密度在280 K附近达到最大值。

2）表面张力模拟结果在常温区域与实验值以及IAPWS公式计算值符合良好，证明了本文模型的准确性。

3）在220~250 K的过冷区域，表面张力温度曲线的斜率在250 K以下区域出现变化，出现了第2个拐点（SIP）。