1.2 模擬細(xì)節(jié)

本研究基于LAMMPS 軟件開展分子動力學(xué)模擬研究。模擬采用的混合勢函數(shù)包括Buckingham 勢(截斷半徑為10 ?)、L-J 12-6 勢(截斷半徑為 10 ?)以及長程庫侖相互作用(截斷半徑為 15 ?)。分子內(nèi)鍵長與鍵角通過SHAKE 算法進行約束。長程庫侖力采用PPPM(particle-particle particle-mesh)方法計算，精度設(shè)置為 10-5。為避免液膜整體平移對應(yīng)力計算產(chǎn)生影響，模擬過程中每25 000 步對體系動量進行重整化以消除整體漂移。初始構(gòu)型經(jīng)控壓(NPT)弛豫后獲得平衡結(jié)構(gòu)；生產(chǎn)計算采用等容過程(NVT)，時間步長均為1.0 fs。

相較于均勻體系，混合物界面性質(zhì)的計算需排除界面間的相互干涉。文獻明確了界面體系的建構(gòu)準(zhǔn)則，即模擬盒子沿液膜法向設(shè)為長方體構(gòu)型，液相體系模擬應(yīng)包含2 000~5 000 個粒子，液膜厚度不低于30 ?。本文針對氫氧化鈉液膜體系的模擬，選取5 000 個氫氧化鈉分子為研究對象，初始模擬盒子為 140 ? ×100 ? ×100 ? 的長方體 (140 ? 為液膜法向尺寸)，在控壓弛豫過程中保持盒體高度 x 不變，通過 y 、 z 平面的距離調(diào)控體系壓力，弛豫后模擬盒子尺寸變?yōu)?140 ? ×80 ? ×80 ?，對應(yīng)液膜厚度為40 ?。為探究水分子含量對熔融氫氧化鈉表面張力的影響，還分別構(gòu)建了含250、125 和30 個水分子的液膜體系。水分子擴散系數(shù)的計算體系具有各向同性特征，因此采用正方體模擬盒子建模，體系由5 000 個氫氧化鈉分子和250 個水分子組成。初始模擬盒子尺寸為 80 ? ×80 ? ×80 ?，控壓收縮后，尺寸變?yōu)?57 ? ×57 ? ×57 ?。圖1 為兩個體系在生產(chǎn)計算前的初始構(gòu)型。所有模擬體系均采用周期性邊界條件。

圖 1 用于生產(chǎn)計算的初始構(gòu)型

a--弛豫后的液膜體系(計算表面張力)；b--弛豫后的體相體系(計算擴散系數(shù))

1.3 計算模型

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

其中：V為體系體積；KB為玻爾茲曼常數(shù)；T為熱力學(xué)溫度；t為時間，取t=0為初始時刻；m為原子質(zhì)量；τ為對應(yīng)t時刻（式（3））的時間延遲；α、β為方向；riα和viα為第i個分子在α方向上的位移和速度；rijα為粒子i和j在α方向的位置矢量分量；fijβ為粒子i和j在β方向的相互作用力分量；〈·〉表示體系內(nèi)所有原子的平均。

本研究采用體系應(yīng)力張量法[24]計算表面張力系數(shù)。該方法利用體系內(nèi)法向應(yīng)力和切向應(yīng)力之間的差異性計算表面張力。對于以z方向為界面法向的體系，應(yīng)力與表面張力系數(shù)的關(guān)系如下：

(5)

其中：σxx、σyy和σzz為應(yīng)力張量的分量；Lx、Ly和Lz為模擬盒子的邊長。

表面張力系數(shù)是表面張力的定量表征，反映的是單位長度界面上分子間的作用力。本文采用應(yīng)力張量法計算得到的即為表面張力系數(shù)，后續(xù)對表面張力變化規(guī)律的分析均以此為基礎(chǔ)。

需要指出的是，在高溫條件下，水分子可能從氫氧化鈉液膜中逸出至真空區(qū)域，若對整個模擬盒子進行積分統(tǒng)計，將導(dǎo)致計算結(jié)果失真，因此僅對液膜所在區(qū)域進行積分處理。由于水分子在氫氧化鈉液膜界面發(fā)生了聚集，為合理納入其對液膜性質(zhì)的貢獻，本文以液膜法向方向上的總分子數(shù)密度作為劃分液膜區(qū)的判據(jù)：將數(shù)密度高于液膜中心密度10%的區(qū)域定義為液膜，以確保積分范圍涵蓋液膜主體及其附近水分子的影響。

1.4 模型驗證

為驗證勢場精度，對比了不同溫度下熔融氫氧化鈉密度、黏度的計算值與實驗值的誤差。根據(jù)文獻[25]的實驗可得知密度ρ（kg/m3）和黏度η（mPa·s）與溫度T（K）的關(guān)系式為：

根據(jù)式（6）、（7）可知，氫氧化鈉的密度與溫度呈負(fù)線性關(guān)系，黏度隨溫度呈指數(shù)衰減。熔融氫氧化鈉密度和黏度及其與實驗值的對比如表3所列。

表  3  熔融氫氧化鈉密度與黏度計算值與實驗值的對比

由表3可知，熔融氫氧化鈉密度計算值與實驗值的相對誤差小于5%，黏度計算值與實驗值的相對誤差小于9%，該誤差處于可接受區(qū)間。可見，Buckingham勢能準(zhǔn)確描述熔融氫氧化鈉密度和黏度性質(zhì)，對該體系具有良好的適用性。

2 結(jié)果與討論

2.1 溫度和水分子對熔融氫氧化鈉表面張力的影響

表面張力系數(shù)是液體界面分子間相互作用的宏觀表現(xiàn)。不同體系液膜的表面張力系數(shù)如圖2 所示。由圖2 可見，對于純?nèi)廴跉溲趸w系，600~750 K 溫度區(qū)間，表面張力的計算值與實驗值吻合良好；750 K 以上的高溫區(qū)，二者出現(xiàn)偏差，計算值高于實驗值。但總體來看，該混合力場可以預(yù)測熔融氫氧化鈉表面張力系數(shù)的變化趨勢。與純氫氧化鈉體系對比，含水體系的表面張力系數(shù)隨水分子含量的升高呈梯度式降低。600 K 時，5% 含水體系的表面張力系數(shù)僅為54 mN/m，相較于純氫氧化鈉體系(143 mN/m)降幅超過 60%，這表明在低溫區(qū)間，水分子對液膜張力有不可忽視的影響。所有含水體系的表面張力系數(shù)均隨溫度的升高持續(xù)上升，當(dāng)體系溫度接近900 K 時，不同含水體系的表面張力系數(shù)趨于一致，水分子對表面張力的影響顯著減弱。

圖 2 不同體系液膜的表面張力系數(shù)

為揭示水分子影響液膜界面張力的微觀機制，本研究采用徑向分布函數(shù)gAB(r)（RDF）評價各粒子間的相互作用強度。RDF是描述系統(tǒng)中i、j兩粒子間距離為r處相對分布的函數(shù)。其定義為：

其中：ρA和ρB分別為粒子A和粒子B的平均密度；NAB（r）為在距離r范圍內(nèi)，粒子A周圍出現(xiàn)粒子B的數(shù)目；V為體系體積；nA和nB分別為粒子A和粒子B的總數(shù)。

通過RDF曲線特征峰的高度與峰位，可判斷粒子間的相互作用強度與空間分布規(guī)律：特征峰越高，表明粒子間相互作用越強；峰位則反映體系在不同分子間距下的結(jié)構(gòu)有序性。對RDF曲線積分可計算以A粒子為中心的B粒子的平均配位數(shù)MAB：