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

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

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

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

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

1.3 計(jì)算模型

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

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

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

(5)

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

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

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

1.4 模型驗(yàn)證

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

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

表  3  熔融氫氧化鈉密度與黏度計(jì)算值與實(shí)驗(yàn)值的對(duì)比

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

2 結(jié)果與討論

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

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

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

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

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

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