水是许多化学反应过程廉价的反应溶剂，也是化工生产过程常用的工质。汽液界面行为是研究水相变传热问题的基础。目前，工程上许多有关水蒸发、水蒸气冷凝、加热干燥等相变传热数据仍主要依赖于实验。随着分子模拟技术的发展，采用分子动力学模拟方法，从分子水平揭示水汽液界面特性的研究，引起了国内外许多学者的极大关注。本文拟采用SPC模型，对水汽液界面特性进行平衡分子动力学模拟研究，探讨温度、截断半径、模拟分子数对水汽液界面特性的影响规律。

1模拟方法

1.1模拟体系的建立

采用直角坐标系，模拟盒子如图1所示，液相位于模拟盒子的中央，汽相分别处于液相的上下两侧，整个模拟体系中有两个汽液界面。模拟盒子在x、y方向的长度为Lx=Ly=L，在z方向的长度为Lz=3L。

图1模拟盒子的示意图

对于水的分子动力学模拟研究，采用的势能模型有很多，如SPC、SPC/E、TIP3P、TIP4P、TIP5P等。本文采用SPC刚体势能模型，假设只有不同水分子的O原子之间存在短程L－J势能，不同水分子的H原子之间以及H原子和O原子之间存在长程静电势能。水分子的总势能由短程L－J势能和长程静电势能两部分组成，如式(1)所示。SPC模型的势能参数如表1所示，其中qH和qO分别为水分子中H原子和O原子所带电荷，rOH为H原子与O原子之间的键长，θ为两个O—H键之间的角度(即键角)，σO为O原子之间L－J势能的尺度参数，εO为O原子之间L－J势能的能量参数，e为基本电荷(1e=1.6×10－19C)，kB为Boltzmann常数(kB=1.3806×10－23J/K)。

表1 SPC模型的参数值

式中:US为总势能，kJ/mol;为长程静电势能，kJ/mol;为短程L－J势能，kJ/mol;N为模拟分子个数;n为每个水分子内受静电作用的作用点数量;i、j为模拟系统内2个不同的水分子;a、b为分子受静电作用的作用点;为i分子中a作用点所带电量，C;为j分子中b作用点所带电量，C;为i分子中a作用点与j分子中b作用点之间的距离，m;εR为真空中介电常数，εR=8.854×10－12F/m;i分子和j分子两个O原子之间的距离，m;σO为O原子之间L－J势能的尺度参数，m;εO为O原子之间L－J势能的能量参数，kJ/mol。

对于长程静电势能，采用作用场法。为避免L－J势能和静电势能在边界处发生截断而不连续，导致Hamiltonian函数不守恒问题。采用移位法来处理两种势能，如方程(2)和(3)所示。

式中:rc为截断半径，m;U为校正后的势能，kJ/mol;Uc为截断半径处的势能，kJ/mol;εS为环境介电常数，通常取εS=∞，因此，式(3)可以简化为方程(4)。

1.2模拟细节

初始时刻，水分子初始位置为各分子的质心以面心立方晶格(FCC)均匀排列在边长为L的液相模拟盒中，液相区上下两侧的汽相区为真空。水分子质心(即O原子所在位置)为分子坐标的原点，H和O原子均在xy平面上，其中一个H原子位于x轴的正方向上，另一个H原子位于xy平面的第二象限区，O和H的位置矢量分别为rO(0，0，0)，rH(0.3159σO，0，0)，rH(－0.1053σO，0.2978σO，0)。水分子初始平动速度由随机数发生器随机给定，初始转动速度为0。

在模拟过程中，对物理量进行无量纲化处理;x、y、z三个方向均采用周期性边界条件;保证系统的体积V、温度T和模拟分子数N保持不变，采用Woodcock变标度恒温法实现系统恒温;不断对体系质心进行矫正，使之处于坐标原点;将模拟盒子沿z方向划分为300个等厚度的薄片;模拟时间步长为0.8fs，总模拟步数为60万步，其中前20万步用于使系统达到平衡，后40万步用于统计界面特性参数。

模拟计算程序是由本课题组采用Fortran语言编写的，其模拟流程如图2所示。模拟运算中所涉及到的方程如式(5)～(13)所示］。

图2模拟流程简图

式中:U(k)为第k个切片的势能，Uij(k)为i、j分子在第k个切片内的势能，nk为第k个切片的分子数，Vs1为切片的体积，ρ(k)为第k个切片的数密度，rij为i分子和j分子之间的距离，xij、yij、zij为rij分别在x、y、z方向上的分量，、、分别为i分子中的a原子和j分子中的b原子之间的距离在x、y、z方向上的分量，U()为势函数U()对的导数，PN(k)、PT(k)分别为第k个切片的法向应力和切向应力，γ(k)为第k个切片的局部界面张力，Δz为切片厚度，γ为汽液界面张力，〈〉为系统统计平均，ρV、ρL分别为汽相主体、液相主体密度，NL、NV分别为液相、汽相切片数，UV、UL分别为汽相主体、液相主体势能(L－J势能、静电势能、总势能)，z(k)为第k个切片的位置，z0为Gibbs汽液界面的位置，d为汽液界面厚度。在统计切片内法向应力和切向应力时，若相互作用的原子a，b均在同一切片内，则计算全部作用;若相互作用原子只有一个原子在某一切片内，则计算一半作用。