水汽液界面特性的平衡分子動(dòng)力學(xué)模擬

王寶和，王蘊(yùn)博，李 群，夏良志，王 剛

(1.大連理工大學(xué)化工學(xué)院，遼寧大連 116024;2.大連理工大學(xué)化機(jī)學(xué)院，遼寧大連 116024)

水是許多化學(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)注［1－3］。本文擬采用SPC模型，對(duì)水汽液界面特性進(jìn)行平衡分子動(dòng)力學(xué)模擬研究，探討溫度、截?cái)喟霃?、模擬分子數(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等［4－5］。本文采用SPC剛體勢(shì)能模型，假設(shè)只有不同水分子的O原子之間存在短程L－J勢(shì)能，不同水分子的H原子之間以及H原子和O原子之間存在長(zhǎng)程靜電勢(shì)能。水分子的總勢(shì)能由短程L－J勢(shì)能和長(zhǎng)程靜電勢(shì)能兩部分組成，如式(1)所示［2，3］。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)所示［2－3，6］。

式中:rc為截?cái)喟霃?，m;U為校正后的勢(shì)能，kJ/mol;Uc為截?cái)喟霃教幍膭?shì)能，kJ/mol;εS為環(huán)境介電常數(shù)，通常取 εS= ∞［2，3，6］，因此，式(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語(yǔ)言編寫的，其模擬流程如圖2所示。模擬運(yùn)算中所涉及到的方程如式(5) ～ (13)所示［2－3，6］。

圖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ì)算一半作用。

2 模擬結(jié)果與討論

2.1 溫度對(duì)密度分布的影響

在模擬分子數(shù)N=256和截?cái)喟霃絩c=0.9498 nm 的條件下，當(dāng)溫度 T=400、450、500、550 和610 K時(shí)，模擬得到的密度分布如圖3所示。統(tǒng)計(jì)得到的汽相主體密度ρV、液相主體密度ρL及汽液界面厚度d如表2所示。由圖3及表2可見，汽相主體密度和汽液界面厚度隨溫度的提高而增加，而液相主體密度隨溫度的提高而減小。

圖3 溫度對(duì)水密度分布的影響

表2 不同溫度下水的汽相主體密度、液相主體密度及界面厚度的模擬結(jié)果

液相主體密度與汽相主體密度之差(ρL－ρV)與溫度T的關(guān)系如圖4所示。可見，液、汽相主體密度之差隨溫度的升高而降低;從理論上講，在臨界點(diǎn)處，其差值應(yīng)該趨近于零，這與圖3所示的規(guī)律一致。液、汽相主體密度之差與溫度的關(guān)系可以擬合成式(14)的形式。

式中水臨界溫度Tc=647.3 K，利用表2數(shù)據(jù)對(duì)式(14)進(jìn)行擬合，得到參數(shù) ρ0=1545.8 kg/m3，指數(shù)因子x=0.5516。

圖4 液、汽相主體密度之差與溫度之間的關(guān)系

2.2 溫度對(duì)界面張力的影響

在模擬分子數(shù)N=256和截?cái)喟霃絩c=0.9498 nm 的條件下，當(dāng)溫度 T=400、450、500、550 和610 K時(shí)，水汽液界面張力的模擬結(jié)果見表3。

表3 不同溫度下汽液界面張力模擬結(jié)果

圖5為局部界面張力的分布曲線(500 K)。由圖5可見，汽相主體的局部界面張力基本為零;從汽相主體向液相主體的過渡過程中，界面張力值逐漸增加，在汽液界面區(qū)達(dá)到峰值;在液相主體又在零值附近波動(dòng)。水汽液界面張力模擬值隨溫度變化規(guī)律如圖6所示。

圖5 局部界面張力的模擬分布曲線(500 K)

圖6 界面張力與溫度的關(guān)系曲線

由圖6可以看出，隨著溫度的提高，界面張力降低，模擬值與實(shí)驗(yàn)值之間的誤差逐漸減小。界面張力與溫度的關(guān)系可以擬合得到方程(15)。

將表3的數(shù)據(jù)對(duì)式(15)進(jìn)行擬合，得到的參數(shù)γ0=254.3 mN·m－1，指數(shù)因子 k=1.305。

2.3 溫度對(duì)勢(shì)能分布的影響

在模擬分子數(shù)N=256和截?cái)喟霃絩c=0.9498 nm 的條件下，當(dāng)溫度 T=400、450、500、550 和610 K時(shí)，汽相主體總勢(shì)能UV、液相主體總勢(shì)能UL及總勢(shì)能勢(shì)阱深度ΔU的模擬結(jié)果如表4所示。圖7為水分子的勢(shì)能分布曲線(500 K)，圖8為液相主體區(qū)域的勢(shì)能隨溫度的變化趨勢(shì)。

表4 不同溫度下水分子勢(shì)能的模擬結(jié)果

圖7 水分子的勢(shì)能分布曲線(500 K)

圖8 液相主體區(qū)域的勢(shì)能隨溫度的變化趨勢(shì)

前已述及，水的勢(shì)能分為L(zhǎng)－J勢(shì)能和靜電勢(shì)能。由圖7可以看出，L－J勢(shì)能均為正值，在液相區(qū)形成勢(shì)壘，勢(shì)壘高度ΔULJ為液相主體L－J勢(shì)能與汽相主體L－J勢(shì)能之差;靜電勢(shì)能均為負(fù)值，在液相區(qū)形成勢(shì)阱，勢(shì)阱深度ΔUe為汽相主體靜電勢(shì)能與液相主體靜電勢(shì)能之差;由于靜電勢(shì)能起主導(dǎo)作用，總勢(shì)能也為負(fù)值，同樣在液相區(qū)形成勢(shì)阱，分子之間主要為吸引作用。從圖8可以看出，汽相主體勢(shì)能作用不明顯，勢(shì)壘高度隨溫度升高而降低，液相主體勢(shì)能的勢(shì)阱深度隨體系溫度的升高而減小。

2.4 模擬分子數(shù)對(duì)模擬結(jié)果的影響

在溫度500 K和截?cái)喟霃絩c=0.9498 nm的條件下，當(dāng)模擬分子數(shù) N=108、256、500和864時(shí)，模擬得到的密度分布見圖9。統(tǒng)計(jì)得到的汽相主體密度ρV、液相主體密度ρL及汽液界面厚度d見表5。

圖9 水分子數(shù)對(duì)密度分布的影響

表5 不同水分子數(shù)下界面性質(zhì)的模擬結(jié)果

由表5和圖9可見，隨著模擬分子數(shù)的增加，液相主體密度有所增加，液相主體區(qū)域?qū)挾燃哟?，汽液界面厚度稍有增大，汽相主體密度有所波動(dòng)。

2.5 截?cái)喟霃綄?duì)模擬結(jié)果的影響

在溫度為500 K和模擬分子數(shù)為864的條件下，當(dāng)截?cái)喟霃?rc=0.7915、0.9498、1.2660 nm時(shí)，模擬得到的密度分布如圖10所示。統(tǒng)計(jì)平均得到的汽相主體密度ρV、液相主體密度ρL及汽液界面厚度d如表6所示。從表6和圖10可以看出，隨著截?cái)喟霃降脑黾?，液相主體密度增大，汽相主體密度減小，汽液界面厚度變化不大。

表6 不同截?cái)喟霃较碌乃航缑嫘再|(zhì)模擬結(jié)果

圖10 截?cái)喟霃綄?duì)密度分布的影響

3 結(jié)論

采用SPC模型，對(duì)水汽液界面特性的分子動(dòng)力學(xué)模擬研究結(jié)果表明，隨著溫度的升高，汽相主體密度增加，汽液界面厚度增大，液相主體密度降低，界面張力逐漸減小，液相主體區(qū)域勢(shì)能的勢(shì)阱深度也逐漸降低。隨著模擬分子數(shù)的增加，液相主體密度增加，汽液界面厚度稍有增大。隨著截?cái)喟霃降脑黾?，液相主體密度增加，汽液界面厚度變化不大。

［1］李印實(shí)，何雅玲，孫 杰，等.不同水分子模型凝結(jié)系數(shù)的分子動(dòng)力學(xué)模擬對(duì)比研究［J］.西安交通大學(xué)學(xué)報(bào)，2006，40(11):1272 －1275.

［2］陳正隆，徐為人，湯立達(dá).分子模擬的理論與實(shí)踐［M］.北京:化學(xué)工業(yè)出版社，2007.

［3］王蘊(yùn)博.Mg(OH)2粉體干燥動(dòng)力學(xué)及水分子動(dòng)力學(xué)模擬［D］.大連:大連理工大學(xué)，2012.

［4］Jorgensen W L，Chandrasekhar J，Madura J D，et al.Comparison of simple potential functions for simulating liquid water［J］.J Chem Phys，1983，79(2):926 －935.

［5］Mahoney M W，Jorgensen W L.A five－site model for liquid water and the reproduction of the density anomaly by rigid，nonpolarizable potential function［J］.J Chem Phys，2000，112(20):8910 －8922.

［6］Rapaport D C.The art of molecular dynamics simulation［M］.Cambridge:Cambridge University Press，1995.