超臨界二氧化碳類液-類氣區(qū)邊界線數(shù)值分析*

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

1) (上海理工大學能源與動力工程學院，上海 200093)

2) (上海市動力工程多相流動與傳熱重點實驗室，上海 200093)

超臨界二氧化碳(S-CO2)因在萃取、沉淀、熱力循環(huán)及化學反應等方面有著十分廣闊的應用前景，逐漸成為學術界的重要研究課題.由于在近臨界區(qū)，可以觀察到隨溫度或壓力變化出現(xiàn)大量的物性異變現(xiàn)象，使得各國學者對流體臨界點附近區(qū)域的研究產(chǎn)生了濃厚興趣.隨著分子動力學模擬技術的快速發(fā)展，該技術可輔助傳統(tǒng)實驗方法用于研究近臨界流體的相關物性.為確定S-CO2 在近臨界區(qū)Widom 線范圍及類液-類氣區(qū)的分子結構特征，本文通過分子動力學模擬技術結合聚類分析，研究了溫度和壓力范圍分別在300—350 K和5.5—18.5 MPa 下，CO2 密度時間序列變異系數(shù)及偏度同Widom 線和類液-類氣區(qū)間的關系.結果表明:S-CO2在近臨界區(qū)Widom 線的確定可通過連接密度時間序列曲線變異系數(shù)極大值點來確定，Widom 線沿著臨界點開始延伸直到350 K 時停止;S-CO2 類液區(qū)和類氣區(qū)的分子分布結構可以用數(shù)密度分布的偏度來區(qū)分，偏度在類氣態(tài)時為正值，在類液態(tài)時為負值，而在Widom 線上達到最大值.

1 引言

自19 世紀Charles Cagniard de Tour 發(fā)現(xiàn)超臨界流體以來，超臨界流體已被廣泛研究并應用于不同行業(yè)中.由于超臨界二氧化碳(supercritical carbon dioxide，S-CO2)在萃取、沉淀、熱力循環(huán)及化學反應等方面中有著十分廣闊的應用前景，逐漸成為學術界重要研究課題[1，2].通常，流體超臨界區(qū)被定義為溫度和壓力均高于其臨界點的區(qū)域，流體物性隨著等溫線和等壓線連續(xù)變化，但不會伴隨著液氣相變.然而，近年來學術界對流體臨界點附近區(qū)域的研究產(chǎn)生了濃厚興趣，因為在近臨界區(qū)，可以觀察到隨溫度或壓力變化出現(xiàn)大量的物性異變現(xiàn)象(即所謂的臨界行為)[3].

現(xiàn)有研究指出，臨界點附近二氧化碳的物性變化將顯著降低循環(huán)中的壓縮功，同時對渦輪機設計的影響不可忽略[4].Clarke 等[5]在研究中就強調(diào)物性，如比熱容和密度等的微小變化會顯著影響熱交換器的設計和性能.眾所周知，許多熱力學函數(shù)，包括比熱容、等溫壓縮和熱膨脹等系數(shù)，在臨界點附近出現(xiàn)最大值.有學者指出不同參數(shù)的最大值在(P，T)平面上彼此相距不遠，并引入了單條超臨界區(qū)最大值線，稱其為Widom 線[6]，該術語最初僅用于表示相關長度的異常.在壓力-溫度空間，不同的物性對應著不同的Widom 線，所以這些Widom 線最終形成一個楔形的Widom 區(qū)域，指向臨界點，標志著一些物性(如壓縮性、比熱容、密度、熱膨脹性、聲速)表現(xiàn)出異常行為的區(qū)域.關于超臨界流體Widom 區(qū)，已經(jīng)開展了大量的研究.Nishikawa 等[7-9]就對超臨界流體進行了系列小角度X 光散射實驗，并指出相關長度和密度波動沿著氣液共存曲線延伸到超臨界區(qū)域形成了一個“脊”線，該“脊”線被定義為密度波動在等溫線上達到最大波動尺度的軌跡;Simeoni 等[10]通過非彈性X 光散射和分子動力學模擬，揭示了氬氣在超臨界區(qū)可被分為類液和類氣兩個區(qū)域;Brazhkin等[11]對氬和氖的熱膨脹系數(shù)、可壓縮性和密度漲落進行了分子動力學模擬，指出單一的Widom 線終止于偏離臨界溫度10%的超臨界區(qū);Sedunov[12]利用分子動力學模擬技術提出一種能量平均密度團簇數(shù)的研究方法，分析了S-CO2中類氣體和類液體之間的軟結構轉(zhuǎn)變機理;Bolmatov 等[13]通過衍射測量，結合分子模擬技術，解釋了超臨界流體中的結構交叉現(xiàn)象，并發(fā)現(xiàn)了溫度-壓力圖上的新熱力學邊界，即Frenkel 線，該線與Widom 線的區(qū)別是，前者與超臨界流體中基本微觀結構和動力學特征的交叉相關，而后者由熱力學響應函數(shù)的最大值的軌跡定義;Mareev 等[14，15]就近臨界流體非線性折射率行為進行了研究，表示團簇的形成是導致非線性折射率產(chǎn)生的主要原因;在我們之前的研究中，也曾對S-CO2在Widom 線對應的溫度壓力范圍內(nèi)的折射率進行了預測，并對其光學異變進行了微觀解釋[16，17].綜上所述，Widom 線將流體超臨界區(qū)分成類液區(qū)和類氣區(qū)，隨著分子模擬技術的快速發(fā)展，可輔助傳統(tǒng)實驗方法用于研究近臨界流體的物性異變現(xiàn)象.

分子動力學模擬常用于模擬氣體、液體和固體分子的性質(zhì)，以高精度還原真實的分子微觀狀態(tài).本文以分子動力學為基礎，用計算機模擬了粒子間相互作用的時間演化過程，所有的分子模擬均是在LAMMPS 開源軟件下進行計算的.值得注意的是，在Widom 線所處的溫度壓力區(qū)域內(nèi)，正常液體和氣體的一些已建立的流量和熱交換計算方法不能使用，因為它們不能處理與這些異變行為有關的突然變化[18].對涉及超臨界狀態(tài)的流動(包括泄漏)和換熱的問題，為了避免在Widom 線所處的溫度壓力區(qū)域內(nèi)使用常規(guī)計算方法，了解這一異變區(qū)域的位置是至關重要的，可以提醒工藝設計人員應注意到異變區(qū)域出現(xiàn)的壓力-溫度范圍[19].因此本文利用分子動力學模擬方法，通過對密度時序曲線、分子聚類等特征進行分析，揭示了S-CO2Widom線對應的溫度壓力范圍、類液-類氣區(qū)結構特征及其判定方法，能夠輔助工藝設計人員或科研人員尋找異變區(qū)域大致位置范圍及對類液-類氣結構進行判定.

2 物理模型及模擬細節(jié)

圖1 為物理模型示意圖.模擬體系是各邊長為10 nm 的立方體盒子，內(nèi)有1000 個CO2分子.系統(tǒng)沿著x，y，z三個方向均采用周期性邊界條件.模擬過程中對整個系統(tǒng)施加NPT 系綜，溫度和壓力的控制采用Nose-Hoover 恒溫器[20]，并采用Kamberaj 等[21]描述的算法動態(tài)更新坐標和速度.設定時間步長為1 fs，將運動方程和velocity-Verlet算法相結合，在各時間步長下更新原子的位置和速度.在每個溫度和壓力下，共運行1000 萬時間步，前500 萬時間步使系統(tǒng)達到平衡狀態(tài)，后500 萬時間步的結果每100 fs 輸出一次.對于每一個分子，分子動力學模擬只計算以該分子為中心，以4σ為半徑的球內(nèi)所有粒子與它的相互作用，從而使計算更加精確，本文潛在的截斷值參考Aimoli 等在文獻[22]中的數(shù)值，設為4σ(σ為尺寸參數(shù)).

圖1 物理模型Fig.1.Physical model.

S-CO2的臨界溫度(Tc)為304.13 K，臨界壓力(Pc)為7.377 Mpa[23].如圖2 所示，本文模擬工況溫度范圍在0.98Tc—1.15Tc，壓力范圍在0.75Pc—2.1Pc(對于溫度340 K 和350 K 兩種工況，壓力范圍分別為1.01Pc—2.5Pc，1.15Pc—2.64Pc).本文所作分子動力學模擬采用標準的12-6 Lennard-Jones 勢能模型，表達式為：

圖2 模擬參數(shù)范圍Fig.2.Range of simulation parameters.

式中，E和rij分別表示原子i和原子j之間的勢能和距離;εij為勢能阱深度，表明粒子間相互作用的強弱;σij表示粒子間相互作用勢為0 的有限距離，與原子直徑有關.本文力場模型選用EPM2 剛性三位點模型，該模型已被學術界廣泛用于S-CO2分子動力學模擬中，并有學者經(jīng)過綜合比較，認為該模型相較于其他模型有著更高的適用性和準確性[22，24，25].表1 中列出了該模型的參數(shù)，其中r0指的是CO2中的C—O 化學鍵的長度.對于異對原子的相互作用參數(shù)，可用洛倫茲-貝特洛組合(Lorentz-Berthelot combining rules)規(guī)則計算[26]，具體為：

表1 EPM2 力場[27]模型參數(shù)Table 1.Force field model parameter.

對于庫侖成對相互作用，采用基于原子的部分電荷模型，并使用(3)式計算庫侖相互作用:

其中C是能量轉(zhuǎn)換常數(shù)，qi和qj分別表示兩原子電荷，而ε是真空介電常數(shù).

3 結果與討論

盡管超臨界流體的物性隨溫度和壓力的變化而連續(xù)變化，但在臨界點附近出現(xiàn)大多數(shù)熱力學特性的異變行為.在臨界點附近，可以觀察到由吉布斯熱力學勢的二階導數(shù)確定的物性參數(shù)的臨界行為，如壓縮系數(shù)、熱膨脹系數(shù)和定壓比熱等會隨著壓力的變化出現(xiàn)極值.以定壓比熱(Cp)為例，通過NIST 數(shù)據(jù)庫[28]計算了CO2在一段溫度范圍內(nèi)的定壓比熱，結果如圖3 所示.S-CO2在近臨界溫度開始，隨著壓力的增加，定壓比熱會出現(xiàn)極值點，且隨著溫度的增加，極值點所對應的壓力也會變大.這種現(xiàn)象會隨著溫度遠離臨界溫度而逐漸消失.實際上，對于其他的物性參數(shù)，也會出現(xiàn)類似的情況.在P-T圖中連接這些極值點，就會出現(xiàn)不同的Widom 線，構成前文所述的Widom 楔形三角區(qū).

圖3 不同溫度下CO2 的定壓比熱Fig.3.Specific heat of CO2 under constant pressure at different temperatures.

通過分子動力學模擬可以發(fā)現(xiàn)，密度的時間序列曲線波動幅度和上述極值點出現(xiàn)的位置非常接近.以溫度310 K 為例，分子模擬在不同壓力下的密度時間序列曲線如圖4 所示，可知Cp值和密度時間序列曲線波動幅度均在8.5 MPa 時表現(xiàn)出最大值.鑒于時間序列曲線波動幅度與物性異變極大值發(fā)生的壓力范圍一致，因此考慮，可通過密度時間序列曲線振幅，分析Widom 線對應的溫度壓力范圍，并界定類液-類氣區(qū)間.

圖4 密度時間序列曲線Fig.4.Density time series curve.

3.1 S-CO2 的Widom 線確定

在表示一組數(shù)據(jù)離散程度的各項指標中，實驗標準差(Sx)是最重要且最常用的，表達式如下:

式中，Xi表示第i個時間下的密度值;n表示時間序列對應的密度樣本總數(shù);表示密度平均數(shù).表2是本文模擬不同溫度和壓力工況下的密度時間序列曲線標準差，可以看出在不同溫度下，標準差隨著壓力的變化會出現(xiàn)最大值，標準差越大，則表示曲線波動程度越大.

然而，當溫度較低時，壓力的線性變化將會帶來密度的非線性跳躍增加，各組數(shù)據(jù)的測量尺度差異太大，直接比較標準差是不適合的.為消除量級的影響，引入變異系數(shù)Cv，表示標準差Sx同平均數(shù)之間的比值:

變異系數(shù)沒有量綱，其值越大，表示變量的離散程度越大，變量的均衡性越差.圖5 表示了溫度和壓力分別在300—350 K 和5.5—19.5 MPa 范圍內(nèi)S-CO2密度時間序列的變異系數(shù).可以看出，隨著溫度的增加，變異系數(shù)極值點對應的壓力逐漸變大，且與NIST 計算的物性參數(shù)極大值對應壓力一致，在溫度為350 K 時，這種現(xiàn)象逐漸消失.

圖5 密度時間序列曲線變異系數(shù)Fig.5.Coefficient of variation of density time series curve.

連接變異系數(shù)的極值點，可以在P-T圖中確定出一條Widom 線，如圖6 所示.

圖6 Widom 線的確定Fig.6.Determination of Widom Line.

在Widom 線的上半部分稱為類液區(qū)，在Widom線的下半部分稱為類氣區(qū).本文同時利用NIST 對相同工況范圍內(nèi)的S-CO2定容比熱、等溫壓縮系數(shù)、體積膨脹系數(shù)進行了計算，將極大值進行連線后發(fā)現(xiàn)，在310 K 之前，NIST 計算的3 種參數(shù)極大值對應的壓力與本文變異系數(shù)預測值一致，310—350 K 范圍內(nèi)接近.從而進一步證明了本文方法的有效性.

Widom 線像是亞臨界區(qū)氣液邊界線的延長線，終止于350 K.但是，僅通過變異系數(shù)無法判斷類液區(qū)和類氣區(qū)的特征，因此考慮引入了偏度系數(shù).

3.2 S-CO2 類液-類氣區(qū)特征

偏度系數(shù)k是統(tǒng)計數(shù)據(jù)分布偏斜方向和程度的度量[29]，表示曲線相對于均值不對稱程度的特征數(shù)，表達式為：

其中E表示期望值.如果數(shù)據(jù)對稱，偏度等于0;若偏度大于0，則表示均值右邊的數(shù)據(jù)更為分散，表明數(shù)據(jù)右偏;若偏度小于0，則表示數(shù)據(jù)左偏.表3 列出了本文模擬工況下的偏度值.

對照變異系數(shù)來看，偏度和類液-類氣的區(qū)分有著很大的關聯(lián).如圖5 中，當溫度為305 K 時，變異系數(shù)極大值對應的壓力為7.5 MPa，在圖6中，當壓力高于此數(shù)值時，均表現(xiàn)為類液區(qū)，而在表3 中此區(qū)域?qū)钠染鶠樨撝?因此從總體看，當壓力小于變異系數(shù)極大值對應的壓力時，偏度值為正數(shù)，即右偏，相對于正態(tài)分布密度曲線，其峰值靠左，此時小于均值的數(shù)據(jù)比大于均值的數(shù)據(jù)多，總體密度偏小，此時流體偏氣態(tài).當壓力大于變異系數(shù)極大值對應的壓力時，偏度值為負數(shù)，即左偏，相對于正態(tài)分布密度曲線，其峰值靠右，此時大于均值的數(shù)據(jù)比小于均值的數(shù)據(jù)多，總體密度偏大，流體偏液態(tài).

表3 密度時間序列曲線偏度Table 3.Skewness of density time series curve.

圖7 給出了溫度分別為305，330，350 K 下不同壓力處的偏度，可以看出，當溫度為350 K 時，偏度的正負區(qū)分消失.

圖7 NIST 計算結果同變異系數(shù)比較Fig.7.NIST results were compared with coefficient of variation.

利用Ovito 軟件可以直觀地展示分子在不同壓力下的聚類情況，設定分子間截斷距離為3.5 ?(一個CO2分子的大小約為3.2 ?)，在此范圍內(nèi)的所有分子稱為一個簇，用相同顏色表示.圖8 給出了溫度為305 K，壓力分別為6.5，7.5，15.5 MPa 下分子聚類示意圖，并在穩(wěn)定后對20 個時間步平均得到了團簇個數(shù)(其余溫度下的聚類情況與之類似).圖9 為分子聚類在 P-T 圖中的表示.

圖8 不同壓力下偏度Fig.8.Different pressure of skewness.

圖9 分子聚類在P-T 圖中的表示Fig.9.Representation of molecular clustering in a P-T plot.

結合偏度可知，當壓力為6.5 MPa 時(虛線最低點)，此時偏度為正值，分子共產(chǎn)生641 個團簇，其中最大團簇包含68 個CO2分子，分子間空隙較大，表現(xiàn)為類氣體特征;當壓力為7.5 MPa 時(虛線中間點)，偏度達到最大值，分子團簇數(shù)明顯減少，僅包含493 個團簇，最大團簇包含130 個CO2分子，密度分布的不對稱性增加，但分子間的空隙依舊較大，近似于類氣體狀態(tài)，但實際上正處于類液到類氣的過渡狀態(tài)，即偽沸騰狀態(tài);當壓力為15.5 MPa 時(虛線最高點)，此時偏度為負值，分子共產(chǎn)生8 個團簇，其中最大團簇包含2973 個CO2分子，分子間空隙減小，表現(xiàn)為類液體特征.

4 結論

本文通過NIST 數(shù)據(jù)庫計算了CO2的相關物性，發(fā)現(xiàn)由吉布斯熱力學勢的二階導數(shù)確定的物性參數(shù)，在不同溫度下臨界行為對應的壓力值，同分子動力學模擬中密度時間序列曲線波動幅度相關，提出可以通過時間序列曲線的變異系數(shù)來確定S-CO2的Widom 線范圍;通過分析偏度并結合聚類分析，提出可通過偏度值對S-CO2類液區(qū)和類氣區(qū)進行界定，得到以下結論.

1) S-CO2在近臨界區(qū)Widom 線的確定可通過連接時間序列曲線變異系數(shù)極大值點來確定，它沿著臨界點開始延伸直到350 K 時停止;由變異系數(shù)確定的Widom 線可近似替代不同參數(shù)在此范圍內(nèi)異常值的連線，但僅通過變異系數(shù)無法判別類液區(qū)和類氣區(qū)結構特征.

2) S-CO2類液區(qū)和類氣區(qū)的分子分布結構可以用密度時間序列的偏度來區(qū)分.它在類氣態(tài)時為正，在類液態(tài)時為負，在Widom 線上達到最大值.