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：