毛志紅 包福兵 余 霞

(中國計量學(xué)院計量測試工程學(xué)院 杭州 310018)

1 引言

氣液界面熱力學(xué)行為一直是相變傳熱傳質(zhì)研究的重點。由于氣液界面層很薄，并且界面不固定，使得采用理論分析和實驗方法研究界面特性十分困難。近些年來，越來越多的學(xué)者采用數(shù)值模擬，特別是采用分子動力學(xué)方法來研究氣液相變界面特性［1］。

Trokhymchuk 和 Alejandre［2］采用分子動力學(xué)方法計算得到了Lennard－Jones流體氬氣液相平衡時的密度、法向與切向壓力張量分量以及界面張力等的分布曲線，總結(jié)了氣液相平衡時氣液相和界面的密度、壓力張量分量等的變化規(guī)律。Frezzotti［3］等人采用分子動力學(xué)方法系統(tǒng)建立了氣液界面模型，指出非平衡態(tài)時，在氣體與氣液界面接觸處形成努森層，描述了純質(zhì)和多組分流體蒸發(fā)凝結(jié)時努森層的作用，并分析了邊界條件對氣液界面性質(zhì)的影響。Lee［4］采用分子動力學(xué)方法，從基于維里定理、Irving－Kirkwood模型和Harasima模型的計算方法出發(fā)，模擬得到了氬流體在氣液平界面處的壓力張量，并對不同方法得到的結(jié)果進行了對比，發(fā)現(xiàn)采用這3種方法模擬得到的切向壓力張量分量和表面張力基本相同，但通過Irving－Kirkwood模型得到的法向壓力張量分量在低溫時幾乎不變，與另外兩者有很大不同。最近，Chilukoti［5］等人對丁烷、己烷、癸烷等在界面附近的自擴散和界面結(jié)構(gòu)進行了分子動力學(xué)模擬，發(fā)現(xiàn)在溫度相同時，這3種烷在界面處的排序和分子鏈形類似。

中國國內(nèi)學(xué)者對氣液界面微觀性質(zhì)的研究開始得較遲。過增元院士研究小組［6］首次采用分子動力學(xué)模擬方法研究了相界面的形成和變化問題，提出氣液界面過渡區(qū)的寬度與系統(tǒng)的溫度有關(guān)，并研究了簡單氬流體的氣液界面的微觀特性，分析了界面區(qū)的熱力學(xué)性質(zhì)，得到氣液界面的密度、壓力張量分量及局部溫度的分布等。陶文銓院士研究小組［7］研究了氣液界面表面張力的大小與系統(tǒng)尺度的關(guān)系，分析了界面溫度和入射分子數(shù)對界面凝結(jié)過程的影響。曾丹苓教授研究小組［8］系統(tǒng)研究了納米流體的輸運特性，并探究了納米液滴在蒸發(fā)過程中界面現(xiàn)象的物理機理。

但是到目前為止，中國國內(nèi)外關(guān)于氣液界面性質(zhì)的研究大多針對單原子分子，而實際工程中的氣體，往往是多原子分子，例如空分精餾過程中氧氣和氮氣，對這類分子的研究較少。即使有，也多是采用單原子模型來研究，這種模擬方法與實際分子結(jié)構(gòu)存在一定誤差。氧氣作為典型的雙原子分子，具有代表性。隨著鋼鐵、化工、國防、醫(yī)療的發(fā)展，國內(nèi)各行業(yè)對氧氣的需求量越來越大。氧氣主要采用深冷法制取，為提高工業(yè)生產(chǎn)效率，改善生活，對氧的低溫?zé)崃W(xué)性質(zhì)的研究提出了迫切需求。本文采用分子動力學(xué)方法，從Lennard－Jones勢能函數(shù)模型出發(fā)，考慮氧分子中原子間鍵的作用，研究雙原子氧的氣液界面性質(zhì)。

2 模擬方法

采用分子動力學(xué)方法來模擬雙原子氧的氣液界面性質(zhì)，氧分子由2個氧原子通過鍵連接構(gòu)成。氧原子間的勢能采用平移后的Lennard－Jones 12－6勢能函數(shù)模型，其表達式為:

式中:ε為勢阱深度;σ為原子作用直徑;rij為i和j兩原子間距;rc為原子作用截斷半徑，在截斷半徑處，兩原子間勢能為零。

氧分子在研究中一般被看作是一種具有彈性的分子，氧分子中原子間鍵的作用常常采用一種諧振勢函數(shù)表示:

式中:K為彈性常數(shù);xe為平衡鍵長，氧分子內(nèi)鍵長［9］xe=1.208 ?。

基于NVT系綜，計算區(qū)域如圖1所示，左右兩側(cè)為氣相區(qū)，中間區(qū)域為液相區(qū)，系統(tǒng)尺寸為55?×165?×55?，3個方向都采用周期性邊界條件。為了提高計算效率，采用元胞鏈表法來計算粒子間作用力，采用Verlet蛙跳法積分粒子的位置和速度。計算中，粒子數(shù)N=3 000，時間步長取為Δt=1 fs。一開始所有分子均勻分布在液相區(qū)域，根據(jù)溫度設(shè)定初始速度。根據(jù)文獻結(jié)果及對計算數(shù)據(jù)的分析，系統(tǒng)經(jīng)過20萬步時達到平衡，再統(tǒng)計之后的180萬步的結(jié)果進行分析。

為了統(tǒng)計界面性質(zhì)，沿y方向把計算區(qū)域平均分成Ns=165層。根據(jù)每層中原子數(shù)和速度，統(tǒng)計得到每層體相的密度、飽和壓強及表面張力，獲得平衡態(tài)氣相飽和壓強、氣相及液相密度、界面區(qū)表面張力等性質(zhì)。每層厚度為Δy的分子數(shù)密度統(tǒng)計公式:

圖1 氣液相變系統(tǒng)示意圖Fig.1 Schematic diagram of gas-liquid phase change system

式中:Vs=Lx×Δy×Lz。

根據(jù)Kirkwood－Buff對表面張力的定義，結(jié)合本文建立的氣液界面模擬體系，表面張力γ可以通過式(4)來計算:

其中:PN和PT分別為界面處的法向壓力張量分量和切向壓力張量分量。

式中:［］表示求總體均值，ρ(y)、T(y)、V(y)分別表示局部數(shù)密度、溫度和體積;dφ(r)/d r為勢函數(shù)對r的導(dǎo)數(shù);kB為波爾茲曼常數(shù)，1.38×10－23J/K。

3 氧原子參數(shù)確定

采用分子動力學(xué)方法模擬氣液界面性質(zhì)過程中，勢能參數(shù)對計算結(jié)果影響很大。以往有關(guān)氣液界面性質(zhì)研究中，通往把氧分子看成一個整體，忽略了分子中原子的振動和取向?qū)φw的影響，而采用雙原子模型研究氧氣液界面性質(zhì)的文獻并不多見。Arora［9］采用雙原子模型研究了氧分子在碳納米管中的運動中，使用了如下的參數(shù):ε/kB=52 K、σ =2.99 ?、xe=1.208 ?、K=11.77 N/cm，該參數(shù)能夠較好地描述氧分子與碳納米管的相互作用。但是，當(dāng)采用上述雙原子模型和參數(shù)模擬存在相變現(xiàn)象的氧氣的氣液界面時，所得結(jié)果與實驗數(shù)據(jù)［10］相比誤差較大。

分子間的引力和斥力作用對氣液相變影響巨大。引力越大，界面中的分子越難擺脫液相區(qū)分子的束縛進入氣相區(qū)，而氣相區(qū)分子越容易被界面中的分子捕捉而進入液相;斥力越大，界面中的分子越容易受到液相區(qū)分子的斥力作用而成為氣相分子，而氣相區(qū)分子越難進入界面區(qū)進而進入液相。勢阱深度ε代表了分子間作用力的強弱，對氣液界面性質(zhì)有較大影響，因此，為簡化計算，模擬中，根據(jù)實驗測得的氣液相密度，調(diào)整勢阱深度ε，而保持其他參數(shù)不變。根據(jù)典型溫度T=80 K，截斷半徑R=4.5σ，對勢阱深度 ε/kB=52、57.97、59.36、59.94、60.17、60.29、61.6 K共7種情況下模擬結(jié)果與實驗數(shù)據(jù)(ρg=1.468 4 kg/m3，ρl=1 190.5 kg/m3)［10］進行了比較，具體見表1。

表1 80 K時，氧在不同ε下氣液相密度與實驗值的比較Table 1 Gas-liquid density of oxygen compared with experimental value under differentεat 80 K

圖2 氧氣氣液相密度隨溫度的變化Fig.2 Gas-liquid density of oxygen under different temperatures

圖3 氧氣的飽和壓強隨溫度的變化Fig.3 Saturation pressure of oxygen under different temperatures

從表1中可以看出，隨著ε/kB的不斷增大，ρg不斷減小，ρl不斷增大。當(dāng) ε/kB=60.17 K時，氣液相密度誤差在2.5%以內(nèi)，與實驗值十分接近，因此，在后續(xù)計算中都采用這個值。圖2和圖3分別給出了系統(tǒng)達到平衡時氧的氣液相密度隨溫度的變化以及氧的飽和壓強隨溫度的變化曲線，并與實驗值進行了比較。從圖中可以看出，隨著溫度升高，液相區(qū)密度逐漸減小，氣相區(qū)密度逐漸增大。計算中，當(dāng)溫度達到154.58 K時，系統(tǒng)密度處于均勻狀態(tài)，定義此時溫度為氧的臨界溫度，模擬值與實驗值符合很好，驗證了雙原子氧模擬程序的正確性以及確定的勢能參數(shù)的可靠性。

4 結(jié)果及分析

根據(jù)確定的氧原子勢能參數(shù)，采用分子動力學(xué)方法模擬獲得不同溫度下氧的壓力張量分量分布曲線，圖4給出了60、90、120和140 K時的氣液界面的法向壓力張量分量PN和切向壓力張量分量PT在界面附近的分布。從圖中可以看出，法向壓力張量分量和切向壓力張量分量在氣相區(qū)域和液相區(qū)域的值基本相等，液相區(qū)域壓力張量分量的波動比氣相區(qū)域的波動較大，而在氣液界面處存在較大的勢壘和勢阱，與此相比較，氣液區(qū)域波動可忽略不計。法向壓力張量分量PN先下降再上升，最后在平衡位置波動。切向壓力張量分量PT先下降再在平衡位置波動。

圖4 不同溫度下氧氣壓力張量分量分布曲線Fig.4 Distribution curve of oxygen stress tensor component under different temperature

根據(jù)式(4)對表面張力的定義，可以得到界面附 近的表面張力分布及界面區(qū)的表面張力。圖5給出了不同溫度下界面附近氧氣表面張力的分布。由圖可以看出，在氣相和液相區(qū)域，不存在界面，此時表面張力的值為0。而在氣液界面區(qū)，表面張力值較大，表面張力的值隨著溫度的升高而減小。圖6給出了不同溫度下氧氣的表面張力值，并與實驗值［11］進行了比較。從圖中可以看出，模擬值與實驗值符合的很好。表面張力隨溫度升高而線性減小。當(dāng)溫度達到臨界溫度154.58 K時，表面張力值減小到0，不存在界面區(qū)。

圖5 表面張力在界面附近的分布Fig.5 Distribution of surface tension near the interface

圖6 表面張力隨溫度的變化Fig.6 Variation of surface tensions with temperatures

4 結(jié)論

采用Lennard－Jones勢能模型，對雙原子氧的氣液界面特性進行了分子動力學(xué)模擬研究。通過計算結(jié)果與實驗值的對比，發(fā)現(xiàn)當(dāng)氧原子的ε/kB=60.17 K時，計算得到的氧氣的氣液相密度以及飽和壓強隨飽和溫度的變化曲線都與實驗結(jié)果符合得很好。通過計算發(fā)現(xiàn)，氣相和液相區(qū)域的切向和法相壓力張量分量相等，這些區(qū)域中的表面張力為零。界面區(qū)壓力張量的法向與切向分量存在較大的勢阱和勢壘，界面區(qū)的表面張力隨溫度線性降低，當(dāng)達到臨界溫度時，表面張力減小到0。通過模擬值與實驗值的對比，驗證了雙原子模型及氧原子勢能參數(shù)的正確性，為進一步研究氧原子在各種氧化物中的性質(zhì)提供了幫助。

1 陳進，徐征，王惠齡.Ar－Kr低溫界面?zhèn)鳠徇^程分子動力學(xué)模擬［J］.低溫工程，2005(3):45－46.

2 Trokhymchuk A，Alejandre J.Computer simulations of liquid/vapor interface in Lennard－Jones fluids:Some questions and answers［J］.Journal of Computational Physics，1999，111(18):8510－8523.

3 Frezzotti A.Boundary conditions at the vapor－liquid interface ［J］.Chem.Phys.，2011，23(3):030609－1－030609－9.

4 Lee SH.Pressure Analyses at the Planar Surface of Liquid－Vapor Argon by a Test－Area Molecular Dynamics Simulation［J］.Bulletin of the Korean Chemical Society，2012，33(9):3039－3042.

5 Chilukoti H K，Kikugawa G，Ohara T.A molecular dynamics study on transport properties and structure at the liquid－vapor interfaces of alkanes［J］.International Journal of Heat and Mass Transfer，2013，59:144－154.

6 王遵敬，陳 民，過增元.汽液界面動力學(xué)行為與熱力學(xué)性質(zhì)的分子動力學(xué)研究［J］.工程熱物理學(xué)報，2001，22(1):70－73.

7 熊建銀，陶文銓，何雅玲.氣液界面特性和凝結(jié)過程的分子動力學(xué)研究［C］.中國工程熱物理學(xué)會第十一屆年會，北京:2005.

8 蔡治勇，劉娟芳，劉 朝，等.納米液滴蒸發(fā)過程的分子動力學(xué)模擬［J］.工程熱物理學(xué)報，2008，29(10):1630－1632.

9 Arora G，Sandler SI.Air separation by single wall carbon nanotubes:Mass transport and kinetic selectivity［J］.Journal of Chemical Physics，2006，124:084702.

10 陳國邦，黃永華，包 銳.低溫流體熱物性［M］.北京:國防工業(yè)出版社，2006.

11 Lemmon E W，Penoncello S G.The Surface Tension of Air and Air Component Mixtures［J］.Advances in Cryogenic Engineering，1994，39:1927－1934.