水是許多化學反應過程廉價的反應溶劑，也是化工生產過程常用的工質。汽液界面行為是研究水相變傳熱問題的基礎。目前，工程上許多有關水蒸發(fā)、水蒸氣冷凝、加熱干燥等相變傳熱數據仍主要依賴于實驗。隨著分子模擬技術的發(fā)展，采用分子動力學模擬方法，從分子水平揭示水汽液界面特性的研究，引起了國內外許多學者的極大關注。本文擬采用SPC模型，對水汽液界面特性進行平衡分子動力學模擬研究，探討溫度、截斷半徑、模擬分子數對水汽液界面特性的影響規(guī)律。

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為模擬系統(tǒng)內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勢能和靜電勢能在邊界處發(fā)生截斷而不連續(xù)，導致Hamiltonian函數不守恒問題。采用移位法來處理兩種勢能，如方程(2)和(3)所示。

式中:rc為截斷半徑，m;U為校正后的勢能，kJ/mol;Uc為截斷半徑處的勢能，kJ/mol;εS為環(huán)境介電常數，通常取εS=∞，因此，式(3)可以簡化為方程(4)。

1.2模擬細節(jié)

初始時刻，水分子初始位置為各分子的質心以面心立方晶格(FCC)均勻排列在邊長為L的液相模擬盒中，液相區(qū)上下兩側的汽相區(qū)為真空。水分子質心(即O原子所在位置)為分子坐標的原點，H和O原子均在xy平面上，其中一個H原子位于x軸的正方向上，另一個H原子位于xy平面的第二象限區(qū)，O和H的位置矢量分別為rO(0，0，0)，rH(0.3159σO，0，0)，rH(－0.1053σO，0.2978σO，0)。水分子初始平動速度由隨機數發(fā)生器隨機給定，初始轉動速度為0。

在模擬過程中，對物理量進行無量綱化處理;x、y、z三個方向均采用周期性邊界條件;保證系統(tǒng)的體積V、溫度T和模擬分子數N保持不變，采用Woodcock變標度恒溫法實現系統(tǒng)恒溫;不斷對體系質心進行矯正，使之處于坐標原點;將模擬盒子沿z方向劃分為300個等厚度的薄片;模擬時間步長為0.8fs，總模擬步數為60萬步，其中前20萬步用于使系統(tǒng)達到平衡，后40萬步用于統(tǒng)計界面特性參數。

模擬計算程序是由本課題組采用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為切片厚度，γ為汽液界面張力，〈〉為系統(tǒng)統(tǒng)計平均，ρV、ρL分別為汽相主體、液相主體密度，NL、NV分別為液相、汽相切片數，UV、UL分別為汽相主體、液相主體勢能(L－J勢能、靜電勢能、總勢能)，z(k)為第k個切片的位置，z0為Gibbs汽液界面的位置，d為汽液界面厚度。在統(tǒng)計切片內法向應力和切向應力時，若相互作用的原子a，b均在同一切片內，則計算全部作用;若相互作用原子只有一個原子在某一切片內，則計算一半作用。