王寶和 ，雷廣平,2，王 維

(1.大連理工大學 化工學院，遼寧 大連 116024 ； 2.中北大學 機械與動力工程學院，山西 太原 030051 ; 3.大連理工大學 化工機械與安全學院，遼寧 大連 116024)

0 前言

眾所周知，高壓氣體的熱力學性質(zhì)明顯背離了理想氣體，因此，高壓氣體的狀態(tài)通常可以用逸度或逸度系數(shù)等熱力學參數(shù)來描述。對于高壓氣-液平衡等過程，采用普通的實驗方法難以實現(xiàn)，利用分子模擬的方法進行逸度或逸度系數(shù)的計算已經(jīng)成為可能。MIZAN等[1]利用分子動力學(MD)模擬，研究了超臨界水中過氧化氫自由基的逸度系數(shù)。WIERZCHOWSKI等[2]采用恒溫恒壓蒙特卡羅(MC)分子模擬法，研究了三種水模型的水蒸氣逸度系數(shù)。本文擬采用分子動力學模擬技術，探討了氣-液平衡狀態(tài)下飽和氬氣體壓力和逸度等熱力學性質(zhì)的計算。

1 模擬方法

1.1 模擬體系的建立

采用直角坐標系，模擬盒子如圖1所示，液相位于模擬盒的中部，氣相分別處于模擬盒的上下兩側(cè)，整個模擬體系中有兩個氣-液界面。對于氬流體，分子i和j間的相互作用勢能函數(shù)U(rij)如式(1)所示[3]。

(1)

式中：rij為分子i和j之間的距離，ε為L-J能量參數(shù)，σ為L-J尺度參數(shù)。對于氬流體，σ=3.396 7×10-10m，ε=1.615 3×10-21J[3]。

圖1 模擬盒示意圖

1.2 模擬細節(jié)

模擬體系由4 000個氬分子組成。初始時刻，飽和氬液體按照面心立方(FCC)方式排布于模擬盒子中部，上下兩側(cè)為真空區(qū)域，作為氣相區(qū)。粒子初始速度隨機給定。為使體系不發(fā)生宏觀運動，模擬過程中不斷調(diào)整質(zhì)心位置及體系總動量，使質(zhì)心位置處于坐標原點，體系總動量為零。截斷半徑rc=4σ(σ為氬分子的尺度參數(shù))。時間步長Δt=4.331×10-15s，z方向劃分300薄層用于統(tǒng)計密度分布。模擬中，x、y、z方向均采用周期性邊界條件。系統(tǒng)采用NVT系綜，采用Woodcock變標度法實現(xiàn)恒溫控制，采用Velocity-Verlet法來求解牛頓運動方程。模擬過程共運算13萬動力學步，前3萬步用來實現(xiàn)氣-液平衡。待系統(tǒng)平衡后，開始統(tǒng)計體系的熱力學性質(zhì)。本文所用模擬程序均為本課題組采用Fortran語言自行編寫。算法流程如圖2所示。

2 結(jié)果與討論

2.1 飽和氬氣體的摩爾體積

采用MD模擬方法，計算得到氬流體在不同溫度和氣-液平衡條件下，液相主體密度(ρL)、氣相主體密度(ρG)及飽和氬氣體的摩爾體積(Vm)，如表1所示；密度分布曲線如圖3所示。

(注：stepe為平衡步數(shù)；steps為總運算步數(shù))

T/KρL/kg·L-1ρG/kg·L-1Vm/cm3·mol-1104.31.2690.0241 696.09113.31.2020.0401 011.56123.21.1180.068611.04131.51.0240.099413.38

圖3 不同溫度下的密度分布曲線

由表1和圖3可見，隨著溫度的升高，液相主體密度和飽和氬氣體的摩爾體積逐漸降低，而汽相主體密度逐漸增大。

2.2 飽和氬氣體壓力的計算

2.2.1RK方程法

本文采用Ridlich-Kwang狀態(tài)方程(2)(以下簡稱為RK方程)[4]，來計算飽和氬氣體的壓力(pRK)。

(2)

式中：R為通用氣體常數(shù)，T為溫度，Vm為飽和氬氣體的摩爾體積，a、b為參數(shù)，可根據(jù)流體的臨界溫度(Tc)和臨界壓力(pc)，利用方程(3)計算得到[4]。

(3a)

b=0.086 64R2Tc/pc

(3b)

根據(jù)表1中飽和氬氣體摩爾體積的MD模擬值，就可以利用方程(3)和(2)，計算出不同溫度下的飽和氬氣體壓力(pRK)，如表2所示。

表2 不同溫度下的pRK 值

2.2.2維里方程法

利用維里方程，采用分子動力學模擬的方法，計算得到不同溫度下飽和氬氣體壓力(pMD，詳細計算步驟參見文獻)，如表3所示[4]。

表3 不同溫度下的pMD值

2.2.3兩種方法計算結(jié)果的比較

在氣-液平衡條件下，飽和氬氣體壓力的實驗值(pexp)，如表4所示[6]。

表4 不同溫度下的pexp值

圖4 不同方法計算得到的壓力值比較

圖4給出了不同方法得到的壓力值的比較?？梢?，三種方法得到的飽和氣體壓力均隨著溫度的升高而增大。在整個模擬溫度范圍內(nèi)，維里方程法及RK方程法得到的壓力值均略低于實驗值。在低溫條件下，pMD與pRK比較接近；在高溫條件下，pRK略小于pMD。整體而言，維里方程法計算得到的壓力值更接近于實驗值。

2.3 飽和氬氣體逸度的計算

逸度的計算方法很多，既可以通過熱力學數(shù)據(jù)計算獲得，也可以通過氣體狀態(tài)方程計算得到。本文采用RK狀態(tài)方程結(jié)合分子動力學模擬數(shù)據(jù)及實驗數(shù)據(jù)，來計算飽和氣體逸度(以下簡稱為RK方程法)。此外，本文還提出一種直接通過分子動力學模擬，來計算飽和氣體逸度的新方法(以下簡稱為虛擬理想氣體法)。

2.3.1RK方程法

在恒溫條件下，純組分逸度與pVT之間的關系，如方程(4)所示[6]。

(4)

式中：f為逸度，p為壓力。將RK方程(2)代入方程(4)，積分整理可以得方程(5)。

(5)

式中：Z為壓縮因子，可利用方程(6)求得；Bp及A/B為參數(shù)，可由方程(7)計算得到。

pVm=ZRT

(6)

Bp=pb/RT

(7a)

A/B=a/bRT1.5

(7b)

由上述方程計算得到不同溫度下的參數(shù)A/B、Bp及壓縮因子Z值，如表5所示。可見，在氣-液平衡條件下，飽和氬氣體的壓縮因子隨溫度的升高而降低，這說明隨著溫度的升高，飽和氣體與理想氣體的偏差會越來越大。

表5 不同溫度下參數(shù)A/B、Bp及壓縮因子Z值

根據(jù)表1的飽和氬氣體的摩爾體積及表2的壓力值，可以利用方程(5)和(6)求得不同溫度下的逸度(fRK)，如表6所示。為便于比較，表6還列出了相應溫度下的壓力值(pRK)。

表6 不同溫度下的pRK及fRK值

2.3.2虛擬理想氣體法

本文提出一種直接利用分子動力學模擬，來計算飽和氣體逸度的新方法，即：將飽和氣體逸度看作是相同體積內(nèi)含有的與飽和氣體具有相同總能量的虛擬理想氣體的壓力。

眾所周知，理想氣體僅具有動能，即其總能量等于動能；而真實氣體的總能量包括動能與勢能。為將飽和氣體轉(zhuǎn)化為相同體積下的虛擬理想氣體，本文將氣體總能量看作理想氣體分子運動的虛擬動能。則相同汽體體積中應含有的理想氣體分子數(shù)(NV)可由式(8)來計算。

(8)

式中：Et為飽和氣體總能量，Ek與Ep分別為飽和氣體所具有的動能與勢能，kB為波爾茲曼常數(shù)。

氣體中理想氣體的壓力可根據(jù)理想氣體狀態(tài)方程計算，且該壓力即視為飽和氣體的逸度(fMD)，如式(9)所示。

fMDV=NVkBT

(9)

式中：V為氣體體積。

根據(jù)上述方程計算得到的不同溫度下的飽和氣體的動能(Ek)、勢能(Ep)、氣體體積(V)、氣體分子數(shù)(N)、虛擬理想氣體分子數(shù)(Nv)，如表7所示?？梢姡S著溫度的升高，氣體中的分子數(shù)增大，這是由于溫度升高分子獲得更多動能，從而更容易掙脫液體對它的束縛而進入氣相區(qū)。

表7 不同溫度下的氣體動能、勢能、氣體體積、氣體分子數(shù)及虛擬理想氣體分子數(shù)

采用虛擬理想氣體法計算得到的飽和氣體逸度如表8所示。為便于比較，表8還列出了相應溫度下的壓力模擬值(pMD)。

根據(jù)方程(6)及(9)，可以得到逸度系數(shù)、壓縮因子及虛擬理想氣體分子數(shù)之間的關系，如方程(10)所示。

表8 不同溫度下的pMD與fMD值

(10)

式中：φ為逸度系數(shù)。

根據(jù)表5及7可以計算得到Nv/(ZN)值，由表8可以計算出φ值，進而可以得到兩者之間的相對誤差，如表9所示?？梢姡S著溫度的升高，Nv/(ZN)與逸度系數(shù)(φ)之間的相對誤差逐漸增大。

表9 不同溫度下Nv/(ZN)、φ及兩者相對誤差

2.3.3不同方法計算結(jié)果的比較

將飽和氣體壓力實驗值及摩爾體積實驗值，利用方程(5)和(6)便可得到實驗壓力所對應的逸度(fexp)。不同溫度下的fexp如表10所示。為便于比較，表10中還列出了相應溫度下的壓力模擬值pexp。

表10 不同溫度下pexp與fexp值

圖5 不同方法所得逸度隨飽和氣體壓力的變化曲線

圖5給出了不同方法計算得到的逸度值與相應飽和氣體壓力的關系曲線?？梢?，隨著飽和氣體壓力的升高，氣體逸度與相應飽和壓力的偏離程度增大，這說明隨著溫度的升高，飽和氣體與理想氣體偏離程度加大；在整個模擬范圍內(nèi)三種方法所得的逸度基本一致。虛擬理想氣體法所得逸度fMD與fexp相比，其值稍微偏高，但相對誤差在15%以內(nèi)，而fRK與fexp則符合得更好，相對誤差<5%。

3 結(jié)論

采用分子動力學模擬技術，對飽和氬氣體的熱力學性質(zhì)進行了模擬計算，得到以下結(jié)論：①隨著溫度的升高，液相主體密度和飽和氬氣體的摩爾體積逐漸降低，汽相主體密度逐漸增大。②維里方程法和RK方程法計算得到的飽和氬氣體壓力與實驗值基本一致。③虛擬理想氣體法及RK方程法計算得到的飽和氬氣體逸度值及實驗值三者基本一致。