1.2 模拟细节

本研究基于LAMMPS 软件开展分子动力学模拟研究。模拟采用的混合势函数包括Buckingham 势(截断半径为10 Å)、L-J 12-6 势(截断半径为 10 Å)以及长程库仑相互作用(截断半径为 15 Å)。分子内键长与键角通过SHAKE 算法进行约束。长程库仑力采用PPPM(particle-particle particle-mesh)方法计算，精度设置为 10-5。为避免液膜整体平移对应力计算产生影响，模拟过程中每25 000 步对体系动量进行重整化以消除整体漂移。初始构型经控压(NPT)弛豫后获得平衡结构；生产计算采用等容过程(NVT)，时间步长均为1.0 fs。

相较于均匀体系，混合物界面性质的计算需排除界面间的相互干涉。文献明确了界面体系的建构准则，即模拟盒子沿液膜法向设为长方体构型，液相体系模拟应包含2 000~5 000 个粒子，液膜厚度不低于30 Å。本文针对氢氧化钠液膜体系的模拟，选取5 000 个氢氧化钠分子为研究对象，初始模拟盒子为 140 Å ×100 Å ×100 Å 的长方体 (140 Å 为液膜法向尺寸)，在控压弛豫过程中保持盒体高度 x 不变，通过 y 、 z 平面的距离调控体系压力，弛豫后模拟盒子尺寸变为 140 Å ×80 Å ×80 Å，对应液膜厚度为40 Å。为探究水分子含量对熔融氢氧化钠表面张力的影响，还分别构建了含250、125 和30 个水分子的液膜体系。水分子扩散系数的计算体系具有各向同性特征，因此采用正方体模拟盒子建模，体系由5 000 个氢氧化钠分子和250 个水分子组成。初始模拟盒子尺寸为 80 Å ×80 Å ×80 Å，控压收缩后，尺寸变为 57 Å ×57 Å ×57 Å。图1 为两个体系在生产计算前的初始构型。所有模拟体系均采用周期性边界条件。

图 1 用于生产计算的初始构型

a--弛豫后的液膜体系(计算表面张力)；b--弛豫后的体相体系(计算扩散系数)

1.3 计算模型

基于Green-Kubo公式[23]计算黏度（η）：

其中：V为体系体积；KB为玻尔兹曼常数；T为热力学温度；t为时间，取t=0为初始时刻；m为原子质量；τ为对应t时刻（式（3））的时间延迟；α、β为方向；riα和viα为第i个分子在α方向上的位移和速度；rijα为粒子i和j在α方向的位置矢量分量；fijβ为粒子i和j在β方向的相互作用力分量；〈·〉表示体系内所有原子的平均。

本研究采用体系应力张量法[24]计算表面张力系数。该方法利用体系内法向应力和切向应力之间的差异性计算表面张力。对于以z方向为界面法向的体系，应力与表面张力系数的关系如下：

(5)

其中：σxx、σyy和σzz为应力张量的分量；Lx、Ly和Lz为模拟盒子的边长。

表面张力系数是表面张力的定量表征，反映的是单位长度界面上分子间的作用力。本文采用应力张量法计算得到的即为表面张力系数，后续对表面张力变化规律的分析均以此为基础。

需要指出的是，在高温条件下，水分子可能从氢氧化钠液膜中逸出至真空区域，若对整个模拟盒子进行积分统计，将导致计算结果失真，因此仅对液膜所在区域进行积分处理。由于水分子在氢氧化钠液膜界面发生了聚集，为合理纳入其对液膜性质的贡献，本文以液膜法向方向上的总分子数密度作为划分液膜区的判据：将数密度高于液膜中心密度10%的区域定义为液膜，以确保积分范围涵盖液膜主体及其附近水分子的影响。

1.4 模型验证

为验证势场精度，对比了不同温度下熔融氢氧化钠密度、黏度的计算值与实验值的误差。根据文献[25]的实验可得知密度ρ（kg/m3）和黏度η（mPa·s）与温度T（K）的关系式为：

根据式（6）、（7）可知，氢氧化钠的密度与温度呈负线性关系，黏度随温度呈指数衰减。熔融氢氧化钠密度和黏度及其与实验值的对比如表3所列。

表  3  熔融氢氧化钠密度与黏度计算值与实验值的对比

由表3可知，熔融氢氧化钠密度计算值与实验值的相对误差小于5%，黏度计算值与实验值的相对误差小于9%，该误差处于可接受区间。可见，Buckingham势能准确描述熔融氢氧化钠密度和黏度性质，对该体系具有良好的适用性。

2 结果与讨论

2.1 温度和水分子对熔融氢氧化钠表面张力的影响

表面张力系数是液体界面分子间相互作用的宏观表现。不同体系液膜的表面张力系数如图2 所示。由图2 可见，对于纯熔融氢氧化体系，600~750 K 温度区间，表面张力的计算值与实验值吻合良好；750 K 以上的高温区，二者出现偏差，计算值高于实验值。但总体来看，该混合力场可以预测熔融氢氧化钠表面张力系数的变化趋势。与纯氢氧化钠体系对比，含水体系的表面张力系数随水分子含量的升高呈梯度式降低。600 K 时，5% 含水体系的表面张力系数仅为54 mN/m，相较于纯氢氧化钠体系(143 mN/m)降幅超过 60%，这表明在低温区间，水分子对液膜张力有不可忽视的影响。所有含水体系的表面张力系数均随温度的升高持续上升，当体系温度接近900 K 时，不同含水体系的表面张力系数趋于一致，水分子对表面张力的影响显著减弱。

图 2 不同体系液膜的表面张力系数

为揭示水分子影响液膜界面张力的微观机制，本研究采用径向分布函数gAB(r)（RDF）评价各粒子间的相互作用强度。RDF是描述系统中i、j两粒子间距离为r处相对分布的函数。其定义为：

其中：ρA和ρB分别为粒子A和粒子B的平均密度；NAB（r）为在距离r范围内，粒子A周围出现粒子B的数目；V为体系体积；nA和nB分别为粒子A和粒子B的总数。

通过RDF曲线特征峰的高度与峰位，可判断粒子间的相互作用强度与空间分布规律：特征峰越高，表明粒子间相互作用越强；峰位则反映体系在不同分子间距下的结构有序性。对RDF曲线积分可计算以A粒子为中心的B粒子的平均配位数MAB：