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

1模擬方法

1.1模擬體系的建立

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

圖1模擬盒子的示意圖

對(duì)于水的分子動(dòng)力學(xué)模擬研究，采用的勢(shì)能模型有很多，如SPC、SPC/E、TIP3P、TIP4P、TIP5P等。本文采用SPC剛體勢(shì)能模型，假設(shè)只有不同水分子的O原子之間存在短程L－J勢(shì)能，不同水分子的H原子之間以及H原子和O原子之間存在長(zhǎng)程靜電勢(shì)能。水分子的總勢(shì)能由短程L－J勢(shì)能和長(zhǎng)程靜電勢(shì)能兩部分組成，如式(1)所示。SPC模型的勢(shì)能參數(shù)如表1所示，其中qH和qO分別為水分子中H原子和O原子所帶電荷，rOH為H原子與O原子之間的鍵長(zhǎng)，θ為兩個(gè)O—H鍵之間的角度(即鍵角)，σO為O原子之間L－J勢(shì)能的尺度參數(shù)，εO為O原子之間L－J勢(shì)能的能量參數(shù)，e為基本電荷(1e=1.6×10－19C)，kB為Boltzmann常數(shù)(kB=1.3806×10－23J/K)。

表1 SPC模型的參數(shù)值

式中:US為總勢(shì)能，kJ/mol;為長(zhǎng)程靜電勢(shì)能，kJ/mol;為短程L－J勢(shì)能，kJ/mol;N為模擬分子個(gè)數(shù);n為每個(gè)水分子內(nèi)受靜電作用的作用點(diǎn)數(shù)量;i、j為模擬系統(tǒng)內(nèi)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為真空中介電常數(shù)，εR=8.854×10－12F/m;i分子和j分子兩個(gè)O原子之間的距離，m;σO為O原子之間L－J勢(shì)能的尺度參數(shù)，m;εO為O原子之間L－J勢(shì)能的能量參數(shù)，kJ/mol。

對(duì)于長(zhǎng)程靜電勢(shì)能，采用作用場(chǎng)法。為避免L－J勢(shì)能和靜電勢(shì)能在邊界處發(fā)生截?cái)喽贿B續(xù)，導(dǎo)致Hamiltonian函數(shù)不守恒問題。采用移位法來處理兩種勢(shì)能，如方程(2)和(3)所示。

式中:rc為截?cái)喟霃剑琺;U為校正后的勢(shì)能，kJ/mol;Uc為截?cái)喟霃教幍膭?shì)能，kJ/mol;εS為環(huán)境介電常數(shù)，通常取εS=∞，因此，式(3)可以簡(jiǎn)化為方程(4)。

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

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

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

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

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

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