水(shui)是(shi)許多(duo)化(hua)(hua)學(xue)反應過(guo)程廉價(jia)的(de)(de)(de)(de)(de)反應溶劑(ji)，也是(shi)化(hua)(hua)工生產過(guo)程常用(yong)(yong)的(de)(de)(de)(de)(de)工質。汽(qi)液界面行為是(shi)研究水(shui)相變傳(chuan)熱(re)問題的(de)(de)(de)(de)(de)基礎。目(mu)前，工程上許多(duo)有關(guan)水(shui)蒸發、水(shui)蒸氣冷凝、加熱(re)干燥等相變傳(chuan)熱(re)數(shu)據仍主要(yao)依賴于實驗(yan)。隨著分(fen)子模(mo)(mo)(mo)擬(ni)技(ji)術的(de)(de)(de)(de)(de)發展，采用(yong)(yong)分(fen)子動(dong)力學(xue)模(mo)(mo)(mo)擬(ni)方法(fa)，從(cong)分(fen)子水(shui)平揭示水(shui)汽(qi)液界面特(te)性(xing)的(de)(de)(de)(de)(de)研究，引起了國內外許多(duo)學(xue)者的(de)(de)(de)(de)(de)極(ji)大(da)關(guan)注。本文(wen)擬(ni)采用(yong)(yong)SPC模(mo)(mo)(mo)型，對水(shui)汽(qi)液界面特(te)性(xing)進(jin)行平衡分(fen)子動(dong)力學(xue)模(mo)(mo)(mo)擬(ni)研究，探討溫度(du)、截斷半徑、模(mo)(mo)(mo)擬(ni)分(fen)子數(shu)對水(shui)汽(qi)液界面特(te)性(xing)的(de)(de)(de)(de)(de)影響(xiang)規(gui)律。

1模擬方法

1.1模擬體系的建(jian)立

采(cai)用直(zhi)角(jiao)坐標系(xi)，模擬盒子如圖(tu)1所示，液(ye)相位于模擬盒子的中央，汽相分別處于液(ye)相的上下兩(liang)側(ce)，整個模擬體系(xi)中有(you)兩(liang)個汽液(ye)界面。模擬盒子在x、y方(fang)向的長度為(wei)Lx=Ly=L，在z方(fang)向的長度為(wei)Lz=3L。

圖1模擬盒子(zi)的(de)示意圖

對于(yu)水(shui)的(de)(de)(de)(de)分子(zi)(zi)(zi)(zi)(zi)(zi)(zi)動力學模(mo)擬研究，采用(yong)的(de)(de)(de)(de)勢能(neng)(neng)模(mo)型(xing)有很多(duo)，如SPC、SPC/E、TIP3P、TIP4P、TIP5P等。本文采用(yong)SPC剛體(ti)勢能(neng)(neng)模(mo)型(xing)，假設只有不同(tong)(tong)水(shui)分子(zi)(zi)(zi)(zi)(zi)(zi)(zi)的(de)(de)(de)(de)O原(yuan)(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)之(zhi)間存(cun)在(zai)短程(cheng)L－J勢能(neng)(neng)，不同(tong)(tong)水(shui)分子(zi)(zi)(zi)(zi)(zi)(zi)(zi)的(de)(de)(de)(de)H原(yuan)(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)之(zhi)間以及H原(yuan)(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)和(he)O原(yuan)(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)之(zhi)間存(cun)在(zai)長程(cheng)靜電勢能(neng)(neng)。水(shui)分子(zi)(zi)(zi)(zi)(zi)(zi)(zi)的(de)(de)(de)(de)總勢能(neng)(neng)由短程(cheng)L－J勢能(neng)(neng)和(he)長程(cheng)靜電勢能(neng)(neng)兩部分組成(cheng)，如式(1)所(suo)示。SPC模(mo)型(xing)的(de)(de)(de)(de)勢能(neng)(neng)參(can)(can)數(shu)如表1所(suo)示，其中(zhong)qH和(he)qO分別為水(shui)分子(zi)(zi)(zi)(zi)(zi)(zi)(zi)中(zhong)H原(yuan)(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)和(he)O原(yuan)(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)所(suo)帶電荷，rOH為H原(yuan)(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)與O原(yuan)(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)之(zhi)間的(de)(de)(de)(de)鍵長，θ為兩個(ge)O—H鍵之(zhi)間的(de)(de)(de)(de)角(jiao)度(即鍵角(jiao))，σO為O原(yuan)(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)之(zhi)間L－J勢能(neng)(neng)的(de)(de)(de)(de)尺度參(can)(can)數(shu)，εO為O原(yuan)(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)之(zhi)間L－J勢能(neng)(neng)的(de)(de)(de)(de)能(neng)(neng)量參(can)(can)數(shu)，e為基(ji)本電荷(1e=1.6×10－19C)，kB為Boltzmann常(chang)數(shu)(kB=1.3806×10－23J/K)。

表1 SPC模型(xing)的參(can)數值

式中:US為(wei)(wei)總勢(shi)能(neng)(neng)，kJ/mol;為(wei)(wei)長程靜電(dian)勢(shi)能(neng)(neng)，kJ/mol;為(wei)(wei)短程L－J勢(shi)能(neng)(neng)，kJ/mol;N為(wei)(wei)模擬分(fen)(fen)(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)個(ge)數;n為(wei)(wei)每(mei)個(ge)水分(fen)(fen)(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)內受靜電(dian)作(zuo)(zuo)(zuo)(zuo)用(yong)的作(zuo)(zuo)(zuo)(zuo)用(yong)點(dian)數量;i、j為(wei)(wei)模擬系統(tong)內2個(ge)不同的水分(fen)(fen)(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi);a、b為(wei)(wei)分(fen)(fen)(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)受靜電(dian)作(zuo)(zuo)(zuo)(zuo)用(yong)的作(zuo)(zuo)(zuo)(zuo)用(yong)點(dian);為(wei)(wei)i分(fen)(fen)(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)中a作(zuo)(zuo)(zuo)(zuo)用(yong)點(dian)所帶電(dian)量，C;為(wei)(wei)j分(fen)(fen)(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)中b作(zuo)(zuo)(zuo)(zuo)用(yong)點(dian)所帶電(dian)量，C;為(wei)(wei)i分(fen)(fen)(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)中a作(zuo)(zuo)(zuo)(zuo)用(yong)點(dian)與j分(fen)(fen)(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)中b作(zuo)(zuo)(zuo)(zuo)用(yong)點(dian)之(zhi)間(jian)(jian)(jian)的距離，m;εR為(wei)(wei)真空中介(jie)電(dian)常(chang)數，εR=8.854×10－12F/m;i分(fen)(fen)(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)和j分(fen)(fen)(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)兩個(ge)O原子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)之(zhi)間(jian)(jian)(jian)的距離，m;σO為(wei)(wei)O原子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)之(zhi)間(jian)(jian)(jian)L－J勢(shi)能(neng)(neng)的尺(chi)度(du)參(can)數，m;εO為(wei)(wei)O原子(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)(zi)之(zhi)間(jian)(jian)(jian)L－J勢(shi)能(neng)(neng)的能(neng)(neng)量參(can)數，kJ/mol。

 對于長程靜電勢能，采用(yong)作用(yong)場(chang)法(fa)。為避免(mian)L－J勢能和靜電勢能在邊界處發生截斷而不連續，導致Hamiltonian函數(shu)不守(shou)恒(heng)問題。采用(yong)移位法(fa)來處理兩種勢能，如方程(2)和(3)所(suo)示。

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

1.2模擬細節(jie)

初(chu)始時刻，水分(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)初(chu)始位(wei)置(zhi)為(wei)(wei)(wei)(wei)各分(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)的(de)(de)質(zhi)心(xin)(xin)以面心(xin)(xin)立方(fang)晶格(FCC)均(jun)勻排列在(zai)邊(bian)長為(wei)(wei)(wei)(wei)L的(de)(de)液(ye)相模擬(ni)盒(he)中，液(ye)相區上下兩側的(de)(de)汽相區為(wei)(wei)(wei)(wei)真空。水分(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)質(zhi)心(xin)(xin)(即(ji)O原(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)所在(zai)位(wei)置(zhi))為(wei)(wei)(wei)(wei)分(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)坐(zuo)標的(de)(de)原(yuan)點，H和O原(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)均(jun)在(zai)xy平(ping)面上，其中一個(ge)H原(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)位(wei)于x軸(zhou)的(de)(de)正方(fang)向上，另一個(ge)H原(yuan)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)位(wei)于xy平(ping)面的(de)(de)第二(er)象限區，O和H的(de)(de)位(wei)置(zhi)矢(shi)量分(fen)別為(wei)(wei)(wei)(wei)rO(0，0，0)，rH(0.3159σO，0，0)，rH(－0.1053σO，0.2978σO，0)。水分(fen)子(zi)(zi)(zi)(zi)(zi)(zi)(zi)初(chu)始平(ping)動(dong)速度(du)由隨機數發(fa)生器(qi)隨機給定，初(chu)始轉(zhuan)動(dong)速度(du)為(wei)(wei)(wei)(wei)0。

在模擬過程中(zhong)，對物(wu)理量進行(xing)無量綱化處理;x、y、z三(san)個(ge)方向均采(cai)用(yong)周(zhou)期性邊界(jie)條件;保證系統的體積V、溫(wen)度(du)T和(he)模擬分(fen)(fen)子數N保持不(bu)(bu)變，采(cai)用(yong)Woodcock變標度(du)恒(heng)溫(wen)法實現(xian)系統恒(heng)溫(wen);不(bu)(bu)斷對體系質心進行(xing)矯正，使(shi)之處于(yu)坐標原點;將模擬盒子沿z方向劃分(fen)(fen)為(wei)300個(ge)等厚(hou)度(du)的薄片;模擬時間(jian)步(bu)長為(wei)0.8fs，總模擬步(bu)數為(wei)60萬(wan)(wan)步(bu)，其中(zhong)前20萬(wan)(wan)步(bu)用(yong)于(yu)使(shi)系統達(da)到平(ping)衡，后40萬(wan)(wan)步(bu)用(yong)于(yu)統計界(jie)面(mian)特性參數。

模(mo)擬計算(suan)程(cheng)序(xu)是由本課題組(zu)采用Fortran語言(yan)編(bian)寫的，其(qi)模(mo)擬流程(cheng)如圖2所(suo)示。模(mo)擬運(yun)算(suan)中所(suo)涉及到的方程(cheng)如式(5)～(13)所(suo)示］。

圖(tu)2模擬流程簡圖(tu)

式中(zhong):U(k)為(wei)(wei)(wei)(wei)第(di)(di)(di)(di)k個(ge)切片(pian)(pian)(pian)的(de)(de)勢(shi)能(neng)(neng)，Uij(k)為(wei)(wei)(wei)(wei)i、j分(fen)(fen)子(zi)在(zai)第(di)(di)(di)(di)k個(ge)切片(pian)(pian)(pian)內(nei)的(de)(de)勢(shi)能(neng)(neng)，nk為(wei)(wei)(wei)(wei)第(di)(di)(di)(di)k個(ge)切片(pian)(pian)(pian)的(de)(de)分(fen)(fen)子(zi)數，Vs1為(wei)(wei)(wei)(wei)切片(pian)(pian)(pian)的(de)(de)體積，ρ(k)為(wei)(wei)(wei)(wei)第(di)(di)(di)(di)k個(ge)切片(pian)(pian)(pian)的(de)(de)數密(mi)度(du)，rij為(wei)(wei)(wei)(wei)i分(fen)(fen)子(zi)和(he)(he)j分(fen)(fen)子(zi)之間(jian)的(de)(de)距離，xij、yij、zij為(wei)(wei)(wei)(wei)rij分(fen)(fen)別(bie)(bie)在(zai)x、y、z方向(xiang)(xiang)(xiang)上的(de)(de)分(fen)(fen)量，、、分(fen)(fen)別(bie)(bie)為(wei)(wei)(wei)(wei)i分(fen)(fen)子(zi)中(zhong)的(de)(de)a原(yuan)(yuan)子(zi)和(he)(he)j分(fen)(fen)子(zi)中(zhong)的(de)(de)b原(yuan)(yuan)子(zi)之間(jian)的(de)(de)距離在(zai)x、y、z方向(xiang)(xiang)(xiang)上的(de)(de)分(fen)(fen)量，U()為(wei)(wei)(wei)(wei)勢(shi)函數U()對的(de)(de)導數，PN(k)、PT(k)分(fen)(fen)別(bie)(bie)為(wei)(wei)(wei)(wei)第(di)(di)(di)(di)k個(ge)切片(pian)(pian)(pian)的(de)(de)法向(xiang)(xiang)(xiang)應力和(he)(he)切向(xiang)(xiang)(xiang)應力，γ(k)為(wei)(wei)(wei)(wei)第(di)(di)(di)(di)k個(ge)切片(pian)(pian)(pian)的(de)(de)局(ju)部界(jie)面(mian)張(zhang)(zhang)力，Δz為(wei)(wei)(wei)(wei)切片(pian)(pian)(pian)厚(hou)度(du)，γ為(wei)(wei)(wei)(wei)汽(qi)(qi)液(ye)(ye)界(jie)面(mian)張(zhang)(zhang)力，〈〉為(wei)(wei)(wei)(wei)系統(tong)(tong)統(tong)(tong)計平均，ρV、ρL分(fen)(fen)別(bie)(bie)為(wei)(wei)(wei)(wei)汽(qi)(qi)相(xiang)主體、液(ye)(ye)相(xiang)主體密(mi)度(du)，NL、NV分(fen)(fen)別(bie)(bie)為(wei)(wei)(wei)(wei)液(ye)(ye)相(xiang)、汽(qi)(qi)相(xiang)切片(pian)(pian)(pian)數，UV、UL分(fen)(fen)別(bie)(bie)為(wei)(wei)(wei)(wei)汽(qi)(qi)相(xiang)主體、液(ye)(ye)相(xiang)主體勢(shi)能(neng)(neng)(L－J勢(shi)能(neng)(neng)、靜電勢(shi)能(neng)(neng)、總勢(shi)能(neng)(neng))，z(k)為(wei)(wei)(wei)(wei)第(di)(di)(di)(di)k個(ge)切片(pian)(pian)(pian)的(de)(de)位置，z0為(wei)(wei)(wei)(wei)Gibbs汽(qi)(qi)液(ye)(ye)界(jie)面(mian)的(de)(de)位置，d為(wei)(wei)(wei)(wei)汽(qi)(qi)液(ye)(ye)界(jie)面(mian)厚(hou)度(du)。在(zai)統(tong)(tong)計切片(pian)(pian)(pian)內(nei)法向(xiang)(xiang)(xiang)應力和(he)(he)切向(xiang)(xiang)(xiang)應力時，若相(xiang)互(hu)作用(yong)的(de)(de)原(yuan)(yuan)子(zi)a，b均在(zai)同一(yi)切片(pian)(pian)(pian)內(nei)，則計算全部作用(yong);若相(xiang)互(hu)作用(yong)原(yuan)(yuan)子(zi)只(zhi)有一(yi)個(ge)原(yuan)(yuan)子(zi)在(zai)某一(yi)切片(pian)(pian)(pian)內(nei)，則計算一(yi)半作用(yong)。