基于分子動力學(xué)模擬的近臨界區(qū)CO2折射率異常波動分析

孫 輝，章立新，楊其國，劉婧楠，高 明

(1.上海理工大學(xué) 能源與動力工程學(xué)院,上海 200093;2.上海市動力工程多相流動與傳熱重點實驗室，上海 200093)

超臨界二氧化碳(S-CO2)動力循環(huán)系統(tǒng)因其優(yōu)良特性在發(fā)電產(chǎn)業(yè)應(yīng)用前景十分廣闊[1]。超臨界流體表現(xiàn)出介于氣體與液體之間的中間特性，被廣泛應(yīng)用在工業(yè)技術(shù)領(lǐng)域中，成為重要的研究課題。對超臨界流體的研究主要集中在宏觀熱力學(xué)性質(zhì)上，也可以通過研究一些動態(tài)參數(shù)(如吸收系數(shù)αv)或靜態(tài)參數(shù)(如折射率n和介電常數(shù)ε)來分析流體的微觀特性[2]。折射率與密度、濃度間的關(guān)系是利用光學(xué)方法測量物性的重要基礎(chǔ)。亞臨界區(qū)氣體介質(zhì)的密度與折射率間的關(guān)系可直接由Gladstone and Dale經(jīng)驗公式表示，而對于超臨界流體，因其物性的非線性變化，已有大量學(xué)者證明Lorentz-Lorenz公式在描述其折射率與密度間的關(guān)系時精度最高[3-4]。

由于CO2的臨界狀態(tài)(臨界溫度為304.13 K，臨界壓力為7.38 MPa)易于實現(xiàn)，很多學(xué)者對其折射率與物性間的關(guān)系進行了研究[2, 5-7]，其中光波長范圍為447.1～667.8 nm，壓力最高達240 MPa，溫度最高達373.15 K。大部分實驗是通過波長為632.8 nm的激光結(jié)合干涉法根據(jù)光斑條紋來獲取CO2折射率。根據(jù)實驗結(jié)果擬合出的洛倫茲常數(shù)，結(jié)合Lorenz-Lorentz公式可以計算出流體的密度；也可以通過折射率的變化間接反映超臨界流體的相組成、溶質(zhì)濃度和界面張力等性質(zhì)。然而，上述實驗中的工況都遠遠偏離臨界狀態(tài)。夏國正[8]和尹德益[9]采用光屏兩點投射法對近臨界區(qū)CO2折射率進行研究，發(fā)現(xiàn)當溫度為304.25 K時，壓力升至6.77 MPa后光斑飄忽不定，無法獲得折射率的準確值，壓力升至9 MPa后光斑恢復(fù)平穩(wěn)；當溫度為309.65 K時，壓力升至7 MPa后因光線散射無法獲得折射率，而壓力升至10 MPa后光斑平穩(wěn)。因此，CO2在近臨界區(qū)的折射率數(shù)據(jù)仍存在較大空白。孫輝等[10]利用神經(jīng)網(wǎng)絡(luò)技術(shù)預(yù)測CO2的折射率，結(jié)果表明該方法在亞臨界區(qū)和超臨界區(qū)預(yù)測折射率的精度高，反演出的密度與美國國家標準與技術(shù)研究院(NIST)數(shù)據(jù)庫計算值相比偏差小，但是在近臨界區(qū)，由于訓(xùn)練集折射率的空白導(dǎo)致出現(xiàn)異常的非平滑曲線。從宏觀來看，CO2在近臨界區(qū)因壓力的細微變化導(dǎo)致光斑舞動無法聚焦，這是因為在近臨界區(qū)超臨界流體的密度波動大，但此時并不能直觀看到其微觀結(jié)構(gòu)究竟發(fā)生了怎樣的異?，F(xiàn)象。分子動力學(xué)模擬技術(shù)很早就被用來研究團簇體和大分子等微觀結(jié)構(gòu)，而目前使用最廣泛的方法是分子動力(MD)模擬方法。在極端條件下，MD模擬中選擇合適的力場模型可以避免復(fù)雜的實驗，高精度地還原真實的分子微觀狀態(tài)，從而研究物質(zhì)的輸運特性和結(jié)構(gòu)特征等。Bolmatov等[11]發(fā)現(xiàn)利用MD模擬方法結(jié)合液體聲子理論可以解釋超臨界流體新的熱力學(xué)邊界線。Sohrevardi等[12]利用MD模擬方法分別研究了CO2、丙酮和水在一定工況下的輸運特性，結(jié)果表明該方法能有效地揭示混合物的輸運特性。

綜上所述，雖然關(guān)于S-CO2光學(xué)性質(zhì)的研究成果較多，但是其研究工況都遠離臨界區(qū)，對于近臨界區(qū)折射率的異常波動未曾從微觀層面給出直觀的解釋。因此，筆者利用MD模擬方法，通過徑向分布函數(shù)(RDF)揭示近臨界狀態(tài)下CO2的微觀結(jié)構(gòu)特性，通過可視化軟件展現(xiàn)近臨界流體的聚類現(xiàn)象，并對密度時間序列曲線進行分析，可以直觀地反映不同溫度下近臨界壓力附近密度的劇烈波動。

1 MD模擬方法

CO2力場模型較多，已有學(xué)者對其模擬結(jié)果的絕對偏差和相對偏差進行了對比。結(jié)果表明，不同的力場模型對物理參數(shù)的預(yù)測精度不同。因此，在研究CO2物理參數(shù)時應(yīng)選擇合適的力場模型，以獲得最佳預(yù)測效果[13]。比較各模型發(fā)現(xiàn)，EPM2、Zhang和SAFT-γ 3種模型使用最多且模擬表現(xiàn)優(yōu)異。在本文研究中，密度對折射率發(fā)生異?，F(xiàn)象的影響較大，考慮到不同模擬所參考的狀態(tài)方程或?qū)嶒瀰?shù)不一，因此有必要對上述3種模型進行比較。

1.1 模型與勢函數(shù)

本文模型使用的勢函數(shù)有Lennard-Jones勢和廣義Mie勢，后者的表達式為：

(1)

式中：E和rij分別為原子i與原子j之間的勢能和距離；εij為勢能阱深度，表明原子間相互作用的強弱；σij為原子間相互作用勢為0的有限距離，與原子直徑有關(guān)；m和b分別為斥力指數(shù)和吸引力指數(shù)；C為常數(shù)。

C是給定指數(shù)的函數(shù)，定義為：

(2)

當m=12，b=6時，C=4，式(1)表示標準的12-6 Lennard-Jones勢，是描述雙原子短程相互作用最常見的勢之一。表1中給出了所需的CO2模型參數(shù)。異對原子的相互作用參數(shù)可根據(jù)洛倫茲-貝特洛組合規(guī)則(Lorentz-Berthelot combining rules)計算[14]。

(3)

式中：σii和σjj為相同原子間相互作用勢為0的有限距離；εii和εjj為相同原子的勢能阱深度。

對于庫侖長程相互作用，采用基于原子的電荷模型，其表達式為：

(4)

式中：qi和qj分別為2個原子的電荷；ε為真空介電常數(shù)。

對于Zhang和EPM2模型，由于其屬于剛體模型，因此只考慮原子間能量。

表1 力場模型參數(shù)

1.2 物理模型

通過自編程建立了包含1 000個CO2分子，邊長為10 nm的立方體盒子(見圖1)，并將系統(tǒng)設(shè)定為等溫等壓系綜(NPT)，所有的模擬都是在LAMMPS開源軟件中進行的。系統(tǒng)在x、y、z3個方向上均采用周期性邊界條件。模擬過程中，對單位點模型采用 Nose-Hoover非哈密頓運動方程，而對剛體模型則采用Kamberaj描述的算法。設(shè)定時間步長為1 fs，將運動方程與Velocity-Verlet算法相結(jié)合，在各時間步長下更新原子的位置和速度。在每個溫度和壓力下，共運行500萬次，前300萬次使系統(tǒng)達到平衡狀態(tài)，后200萬次的結(jié)果每100 fs輸出一次，作為最終各模型的比較數(shù)據(jù)。潛在的截斷值設(shè)為4σ(σ為尺寸參數(shù))。密度輸出后取平均值，并與REFPROP軟件給出的密度進行比較，從而對各模型進行評價分析。選取密度平均相對誤差最小的模型進行后續(xù)分析。

圖1 物理模型

由于筆者主要研究CO2在近臨界區(qū)折射率異?，F(xiàn)象的影響因素，于是在溫度T為300 K、305 K和310 K，壓力p為5.5～15.5 MPa狀態(tài)下對3種物理模型進行了模擬，結(jié)果如圖2所示。

(a) 各模型預(yù)測密度與NIST預(yù)測密度的對比

(b) 不同溫度下的平均相對誤差

由圖2可知，不同工況下3種模型預(yù)測所得的CO2密度的變化趨勢與NIST預(yù)測值一致，總體偏差較小，但存在極個別偏差較大的數(shù)據(jù)點，恰在近臨界區(qū)。從圖2(a)可以看出，EPM2模型預(yù)測值與NIST預(yù)測值最接近。從圖2(b)可以看出，EPM2模型在溫度300～310 K、壓力5.5～15.5 MPa內(nèi)的平均相對誤差最小。Zhang、SAFT-γ和EPM2模型的平均相對誤差分別為8.87%、9.66%和6.34%。因此，選用EPM2模型進行分析。

2 結(jié)果及分析

2.1 徑向分布函數(shù)分析

徑向分布函數(shù)gm(r)為：

(5)

式中：N為系統(tǒng)的分子數(shù)目；dN為距參考分子中心距離r到r+dr間的分子數(shù)目；ρs為系統(tǒng)的密度。

徑向分布函數(shù)表示系統(tǒng)區(qū)域密度與平均密度之比，可以描述系統(tǒng)的結(jié)構(gòu)性質(zhì)。徑向分布函數(shù)通常會出現(xiàn)高低不同的峰谷，這是因為隨著r的變化，區(qū)域密度與平均密度的比值隨之改變，一般來說在距離參考分子較遠處區(qū)域密度與平均密度相同，即徑向分布函數(shù)值在1附近波動。通過計算得到不同壓力下徑向分布函數(shù)隨溫度的變化趨勢，如圖3所示。

(a) p=5.5 MPa

(b) p=7.5 MPa

(c) p=8.5 MPa

(d) p=9.5 MPa

(e) p=10.5 MPa

(f) p=12.5 MPa

由圖3(a)可知，當壓力低于臨界值時，隨著溫度的升高，徑向分布函數(shù)的波峰逐漸減小，這是由于CO2逐漸由液態(tài)過渡到氣態(tài)，系統(tǒng)的區(qū)域密度逐漸減小。由圖3(b)可知，壓力為7.5 MPa時，隨著溫度的升高，徑向分布函數(shù)的峰值先升后降，差別很大。300 K時的峰值急劇下降，表現(xiàn)出類氣體的異常行為，導(dǎo)致這一現(xiàn)象的實際原因是溫度和壓力均接近臨界值，分子間形成部分團簇，而在聚類情況下，系統(tǒng)的區(qū)域密度與平均密度十分接近，因此其峰值只有1.25。而后隨著溫度的繼續(xù)升高，系統(tǒng)逐漸過渡到超臨界狀態(tài)，展現(xiàn)出類液體現(xiàn)象。同理，由圖3(c)可知，當壓力為8.5 MPa，溫度分別為300 K和305 K時，系統(tǒng)表現(xiàn)出聚類特征，直至溫度升至310 K時才過渡到類液體的超臨界區(qū)域。由圖3(d)～圖3(f)可知，當壓力遠離臨界值時，溫度的變化對結(jié)構(gòu)特征的影響不大，超臨界流體為高度的壓縮性流體，密度逐漸增大，整個系統(tǒng)形成巨大團簇。

聚類分析可以通過OVITO軟件直觀地展現(xiàn)出來，以溫度為305 K，壓力分別為7.5 MPa、10.5 MPa為例，畫出系統(tǒng)平衡后某一時刻的粒子位置圖(見圖4)，設(shè)定分子間截斷距離為3.5 ?(1個CO2分子的大小約為2.5 ?)，在此范圍內(nèi)的分子稱為1個簇。

圖4(a)給出了壓力為7.5 MPa時的分子聚集狀態(tài)，其最大團簇含有80個分子，最小團簇含有1個分子，共形成446種團簇。圖4(b)中壓力為10.5 MPa時只含有24種團簇。實際上，這些聚集體、團簇的形成會改變CO2的光學(xué)性質(zhì)(如線性和非線性折射率)[15]。在聚類出現(xiàn)的情況下，超臨界流體的Lorentz-Lorenz公式[16]為：

(6)

實際上，在強光場中，折射率n還與光強有關(guān)，式(6)中n=n0+n2E2,n0為線性折射率,當非線性折射率n2=0時，n=n0。若存在單個分子和團簇共同存在的二元混合物，其非線性折射率n2可表示為：

(7)

式中：RLL為洛倫茲常數(shù)；β為單個分子的極化率；β2為包含多個分子團簇的極化率；mmol為分子的質(zhì)量；s為團簇中的分子數(shù)；N、N1和N2分別為分子總聚集度、單個分子聚集度和團簇聚集度，N=N1+N2s；γ為形狀因子，對于任意的簇，0.5<γ<1。

(a) p=7.5 MPa

(b) p=10.5 MPa

圖5給出了不同溫度下的團簇種類個數(shù)與壓力的關(guān)系。在不同溫度下，隨著壓力的變化，團簇種類個數(shù)差異很大。當壓力處于亞臨界狀態(tài)時，團簇種類少，幾乎是單個分子成為1個個體，非線性折射率幾乎為0，但是隨著壓力的升高，尤其在臨界值附近，聚簇和解聚同時發(fā)生，團簇種類和數(shù)量波動劇烈，單個分子和疏密不同的團簇多元共存，團簇密度和系統(tǒng)整體密度劇烈變化，折射率也呈現(xiàn)出強烈的非線性變化。從OVITO軟件分析來看，團簇種類個數(shù)在不同瞬態(tài)波動較大，隨著壓力繼續(xù)升高，聚簇動力增加，團簇種類個數(shù)明顯減少，系統(tǒng)形成1個巨型團簇(圖4(b)中分子聚集程度較圖4(a)密，說明在10.5 MPa壓力下，團簇種類個數(shù)減少)。此時團簇內(nèi)的總分子數(shù)目幾乎不變，在Widom線(即超臨界區(qū)類液體與類氣體的分界線)[17-18]附近達到極大值。當壓力高于Widom線時，隨著壓力的升高，形成巨晶層，折射率的非線性也不明顯了。

圖5 不同溫度下的團簇種類個數(shù)

根據(jù)費馬最小時間原理，在兩點之間所有的可能路徑中，光會選擇耗時最少的路徑[19]。射線在介質(zhì)中每點的速度是c0/n(c0為真空中的光速)，因此射線需要尋找一條優(yōu)先通過折射率小的路徑。從結(jié)構(gòu)特性來看，近臨界區(qū)團簇種類繁雜，變化劇烈，射線每次通過路徑所需的時間不同，所以從光屏來看，會產(chǎn)生無規(guī)則的長條狀。

2.2 密度時間序列曲線分析

密度時間序列曲線可由Lammps直接輸出，通過其波動規(guī)律來描述混沌和類隨機特性[20]。筆者主要通過時序時間來直觀展現(xiàn)近臨界區(qū)CO2密度的劇烈波動，從而解釋近臨界區(qū)CO2異常的光學(xué)現(xiàn)象。圖6為不同壓力下的密度時間序列曲線圖。從圖6可以看出，隨著溫度變化，不同壓力下的密度波動程度也隨之變化。

(a) 300 K

(b) 305 K

(c) 310 K

當溫度低于臨界值時，密度會在臨界壓力處劇烈波動；當溫度超過臨界值時，密度會在近臨界壓力處劇烈波動，這與上述對徑向分布函數(shù)的分析結(jié)論一致。從密度時間序列曲線來看，在遠高于臨界壓力和低于臨界壓力處，密度的波動逐漸呈現(xiàn)出規(guī)則化。這是由于前者生成了穩(wěn)定的團簇體，而后者是穩(wěn)定的單分子，而在近臨界區(qū)，分子聚簇和解聚變化劇烈，團簇種類繁多，所以密度波動很大。

3 結(jié) 論

(1) 通過徑向分布函數(shù)的分析和聚類現(xiàn)象的可視化，從微觀層面揭示了CO2折射率在近臨界區(qū)出現(xiàn)異常的原因，即由于近臨界區(qū)CO2的單分子和團簇體多元混合，團簇密度波動，進而影響光對最短路徑的找尋，從而產(chǎn)生非線性折射率。這種情況將隨著壓力的升高而消失，因為在遠離臨界壓力后，系統(tǒng)內(nèi)形成穩(wěn)定的巨型團簇，不會出現(xiàn)明顯的非線性折射。

(2) 在近臨界區(qū)CO2密度時間序列曲線出現(xiàn)混沌，系統(tǒng)密度發(fā)生劇烈波動，根據(jù)Lorentz-Lorenz公式，密度的波動會直接引起折射率的劇烈波動。

(3) 綜合結(jié)構(gòu)特征和密度時間序列曲線來看，在測量近臨界區(qū)CO2折射率時，可以對瞬態(tài)的折射率進行抓拍。類比密度時間序列曲線，只要瞬態(tài)圖形足夠多，就可進一步分析求解出平均折射率，亦可以根據(jù)Lorentz-Lorenz公式對近臨界區(qū)的折射率進行預(yù)測。