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

1模擬方法

1.1模擬體系的建立

采用直角坐標系，模擬盒子如圖1所示，液相位于模擬盒子的中央，汽相分別處于液相的上下兩側，整個(gè)模擬體系中有兩個(gè)汽液界面。模擬盒子在x、y方向的長(cháng)度為L(cháng)x=Ly=L，在z方向的長(cháng)度為L(cháng)z=3L。

圖1模擬盒子的示意圖

對于水的分子動(dòng)力學(xué)模擬研究，采用的勢能模型有很多，如SPC、SPC/E、TIP3P、TIP4P、TIP5P等。本文采用SPC剛體勢能模型，假設只有不同水分子的O原子之間存在短程L－J勢能，不同水分子的H原子之間以及H原子和O原子之間存在長(cháng)程靜電勢能。水分子的總勢能由短程L－J勢能和長(cháng)程靜電勢能兩部分組成，如式(1)所示。SPC模型的勢能參數如表1所示，其中qH和qO分別為水分子中H原子和O原子所帶電荷，rOH為H原子與O原子之間的鍵長(cháng)，θ為兩個(gè)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;為長(cháng)程靜電勢能，kJ/mol;為短程L－J勢能，kJ/mol;N為模擬分子個(gè)數;n為每個(gè)水分子內受靜電作用的作用點(diǎn)數量;i、j為模擬系統內2個(gè)不同的水分子;a、b為分子受靜電作用的作用點(diǎn);為i分子中a作用點(diǎn)所帶電量，C;為j分子中b作用點(diǎn)所帶電量，C;為i分子中a作用點(diǎn)與j分子中b作用點(diǎn)之間的距離，m;εR為真空中介電常數，εR=8.854×10－12F/m;i分子和j分子兩個(gè)O原子之間的距離，m;σO為O原子之間L－J勢能的尺度參數，m;εO為O原子之間L－J勢能的能量參數，kJ/mol。

對于長(cháng)程靜電勢能，采用作用場(chǎng)法。為避免L－J勢能和靜電勢能在邊界處發(fā)生截斷而不連續，導致Hamiltonian函數不守恒問(wèn)題。采用移位法來(lái)處理兩種勢能，如方程(2)和(3)所示。

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

1.2模擬細節

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

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

模擬計算程序是由本課題組采用Fortran語(yǔ)言編寫(xiě)的，其模擬流程如圖2所示。模擬運算中所涉及到的方程如式(5)～(13)所示］。

圖2模擬流程簡(jiǎn)圖

式中:U(k)為第k個(gè)切片的勢能，Uij(k)為i、j分子在第k個(gè)切片內的勢能，nk為第k個(gè)切片的分子數，Vs1為切片的體積，ρ(k)為第k個(gè)切片的數密度，rij為i分子和j分子之間的距離，xij、yij、zij為rij分別在x、y、z方向上的分量，、、分別為i分子中的a原子和j分子中的b原子之間的距離在x、y、z方向上的分量，U()為勢函數U()對的導數，PN(k)、PT(k)分別為第k個(gè)切片的法向應力和切向應力，γ(k)為第k個(gè)切片的局部界面張力，Δz為切片厚度，γ為汽液界面張力，〈〉為系統統計平均，ρV、ρL分別為汽相主體、液相主體密度，NL、NV分別為液相、汽相切片數，UV、UL分別為汽相主體、液相主體勢能(L－J勢能、靜電勢能、總勢能)，z(k)為第k個(gè)切片的位置，z0為Gibbs汽液界面的位置，d為汽液界面厚度。在統計切片內法向應力和切向應力時(shí)，若相互作用的原子a，b均在同一切片內，則計算全部作用;若相互作用原子只有一個(gè)原子在某一切片內，則計算一半作用。