孫輝 劉婧楠? 章立新 楊其國(guó) 高明

1) (上海理工大學(xué)能源與動(dòng)力工程學(xué)院，上海 200093)

2) (上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室，上海 200093)

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

1 引言

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

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

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

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

圖1 為物理模型示意圖.模擬體系是各邊長(zhǎng)為10 nm 的立方體盒子，內(nèi)有1000 個(gè)CO2分子.系統(tǒng)沿著x，y，z三個(gè)方向均采用周期性邊界條件.模擬過程中對(duì)整個(gè)系統(tǒng)施加NPT 系綜，溫度和壓力的控制采用Nose-Hoover 恒溫器[20]，并采用Kamberaj 等[21]描述的算法動(dòng)態(tài)更新坐標(biāo)和速度.設(shè)定時(shí)間步長(zhǎng)為1 fs，將運(yùn)動(dòng)方程和velocity-Verlet算法相結(jié)合，在各時(shí)間步長(zhǎng)下更新原子的位置和速度.在每個(gè)溫度和壓力下，共運(yùn)行1000 萬時(shí)間步，前500 萬時(shí)間步使系統(tǒng)達(dá)到平衡狀態(tài)，后500 萬時(shí)間步的結(jié)果每100 fs 輸出一次.對(duì)于每一個(gè)分子，分子動(dòng)力學(xué)模擬只計(jì)算以該分子為中心，以4σ為半徑的球內(nèi)所有粒子與它的相互作用，從而使計(jì)算更加精確，本文潛在的截?cái)嘀祬⒖糀imoli 等在文獻(xiàn)[22]中的數(shù)值，設(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(對(duì)于溫度340 K 和350 K 兩種工況，壓力范圍分別為1.01Pc—2.5Pc，1.15Pc—2.64Pc).本文所作分子動(dòng)力學(xué)模擬采用標(biāo)準(zhǔn)的12-6 Lennard-Jones 勢(shì)能模型，表達(dá)式為：

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

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

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

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

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

3 結(jié)果與討論

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

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

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

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

3.1 S-CO2 的Widom 線確定

在表示一組數(shù)據(jù)離散程度的各項(xiàng)指標(biāo)中，實(shí)驗(yàn)標(biāo)準(zhǔn)差(Sx)是最重要且最常用的，表達(dá)式如下:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

利用Ovito 軟件可以直觀地展示分子在不同壓力下的聚類情況，設(shè)定分子間截?cái)嗑嚯x為3.5 ?(一個(gè)CO2分子的大小約為3.2 ?)，在此范圍內(nèi)的所有分子稱為一個(gè)簇，用相同顏色表示.圖8 給出了溫度為305 K，壓力分別為6.5，7.5，15.5 MPa 下分子聚類示意圖，并在穩(wěn)定后對(duì)20 個(gè)時(shí)間步平均得到了團(tuán)簇個(gè)數(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.

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

4 結(jié)論

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

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

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