雙原子氧氣液界面特性的分子動(dòng)力學(xué)模擬

毛志紅 包福兵 余 霞

(中國(guó)計(jì)量學(xué)院計(jì)量測(cè)試工程學(xué)院 杭州 310018)

1 引言

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

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

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

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

2 模擬方法

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

式中:ε為勢(shì)阱深度;σ為原子作用直徑;rij為i和j兩原子間距;rc為原子作用截?cái)喟霃?，在截?cái)喟霃教?，兩原子間勢(shì)能為零。

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

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

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

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

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

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

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

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

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

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

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

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

表1 80 K時(shí)，氧在不同ε下氣液相密度與實(shí)驗(yàn)值的比較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 氧氣的飽和壓強(qiáng)隨溫度的變化Fig.3 Saturation pressure of oxygen under different temperatures

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

4 結(jié)果及分析

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

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

根據(jù)式(4)對(duì)表面張力的定義，可以得到界面附 近的表面張力分布及界面區(qū)的表面張力。圖5給出了不同溫度下界面附近氧氣表面張力的分布。由圖可以看出，在氣相和液相區(qū)域，不存在界面，此時(shí)表面張力的值為0。而在氣液界面區(qū)，表面張力值較大，表面張力的值隨著溫度的升高而減小。圖6給出了不同溫度下氧氣的表面張力值，并與實(shí)驗(yàn)值［11］進(jìn)行了比較。從圖中可以看出，模擬值與實(shí)驗(yàn)值符合的很好。表面張力隨溫度升高而線性減小。當(dāng)溫度達(dá)到臨界溫度154.58 K時(shí)，表面張力值減小到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勢(shì)能模型，對(duì)雙原子氧的氣液界面特性進(jìn)行了分子動(dòng)力學(xué)模擬研究。通過計(jì)算結(jié)果與實(shí)驗(yàn)值的對(duì)比，發(fā)現(xiàn)當(dāng)氧原子的ε/kB=60.17 K時(shí)，計(jì)算得到的氧氣的氣液相密度以及飽和壓強(qiáng)隨飽和溫度的變化曲線都與實(shí)驗(yàn)結(jié)果符合得很好。通過計(jì)算發(fā)現(xiàn)，氣相和液相區(qū)域的切向和法相壓力張量分量相等，這些區(qū)域中的表面張力為零。界面區(qū)壓力張量的法向與切向分量存在較大的勢(shì)阱和勢(shì)壘，界面區(qū)的表面張力隨溫度線性降低，當(dāng)達(dá)到臨界溫度時(shí)，表面張力減小到0。通過模擬值與實(shí)驗(yàn)值的對(duì)比，驗(yàn)證了雙原子模型及氧原子勢(shì)能參數(shù)的正確性，為進(jìn)一步研究氧原子在各種氧化物中的性質(zhì)提供了幫助。

1 陳進(jìn)，徐征，王惠齡.Ar－Kr低溫界面?zhèn)鳠徇^程分子動(dòng)力學(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 王遵敬，陳 民，過增元.汽液界面動(dòng)力學(xué)行為與熱力學(xué)性質(zhì)的分子動(dòng)力學(xué)研究［J］.工程熱物理學(xué)報(bào)，2001，22(1):70－73.

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

8 蔡治勇，劉娟芳，劉 朝，等.納米液滴蒸發(fā)過程的分子動(dòng)力學(xué)模擬［J］.工程熱物理學(xué)報(bào)，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 陳國(guó)邦，黃永華，包 銳.低溫流體熱物性［M］.北京:國(guó)防工業(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.