劉旻昀，黃彥平，唐 佳，臧金光，趙學斌，黃善仿，何茂剛，張 博

(1.中國核動力研究設(shè)計院 中核核反應堆熱工水力技術(shù)重點實驗室，四川 成都 610213；2.清華大學 工程物理系，北京 100084；3. 西安交通大學 熱流科學與工程教育部重點實驗室，陜西 西安 710049；4.大連理工大學 能源與動力學院，遼寧 大連 116024)

超臨界流體的熱物理性質(zhì)在其臨界點和擬臨界線附近會發(fā)生劇烈變化，由于這種獨特的物性畸變特性，萃取和化學合成領(lǐng)域可通過細微的溫度壓力調(diào)節(jié)使物質(zhì)的溶解度發(fā)生幾個數(shù)量級的突變，以達到使分離物析出或控制反應的目的[1]。在超臨界二氧化碳布雷頓循環(huán)中，通常將壓縮機進口溫度設(shè)置在擬臨界溫度附近，利用密度的畸變特性降低壓縮功、提高系統(tǒng)效率[2-3]。然而，在超臨界水動力循環(huán)或超臨界碳氫燃料熱管理系統(tǒng)中，超臨界流體的物性畸變特性還可能誘發(fā)傳熱惡化使得壁面過熱失效[4]。因此，為保障和優(yōu)化超臨界流體在能源動力、傳質(zhì)萃取、化學反應、環(huán)境治理等領(lǐng)域的應用，必須對超臨界流體物性畸變特性進行深入研究。目前，關(guān)于物性畸變機理的認識尚未統(tǒng)一。一種觀點[5-7]認為，擬臨界現(xiàn)象與亞臨界氣液相變物性變化特點相近。在擬臨界現(xiàn)象發(fā)生時，超臨界流體內(nèi)部存在類氣結(jié)構(gòu)與類液結(jié)構(gòu)之間的轉(zhuǎn)變，轉(zhuǎn)變過程仍可用超臨界氣液共存和擬沸騰來描述。另一種觀點[8-9]將擬臨界現(xiàn)象視作臨界現(xiàn)象的延伸，強調(diào)體系內(nèi)部漲落的影響。第三種觀點[10]則不考慮擬臨界現(xiàn)象本身的特殊性，在研究相關(guān)問題時看作物性劇烈變化的單相流體。但上述3種觀點均沒有嚴格的理論支撐，同時也無法給出超臨界流體發(fā)生物性畸變的根本原因。

從描述短程相互作用的分子作用力，到描述長程關(guān)聯(lián)的有序結(jié)構(gòu)，再到宏觀上表現(xiàn)出的均一相態(tài)，超臨界流體物性畸變特性存在于微觀、介觀和宏觀的不同尺度上。在微觀尺度上，光散射和中子散射實驗[11]為研究體系的微觀結(jié)構(gòu)變化特點提供了有力手段，而分子動力學[7,12]等模擬手段則作為一種更廉價而精細的補充。在介觀尺度上，對于臨界現(xiàn)象，標度理論[13]和重整化群理論[14-15]從標度不變性出發(fā)，用粗?；乃枷虢⒘送暾睦碚撁枋?；而對于擬臨界區(qū)物性畸變，尚未有理論提出。在宏觀尺度上，關(guān)于超臨界流體傳熱異化的研究[16]已有很多，但在擬臨界區(qū)物性計算方面，如今應用最為廣泛的NIST數(shù)據(jù)庫[17]也有較大偏差。針對上述問題，本文從微觀、介觀和宏觀尺度對超臨界流體的物性畸變特性進行全面分析。采用分子動力學模擬復現(xiàn)和研究超臨界二氧化碳在擬臨界區(qū)的物性畸變特點；提出一種考慮外場的氣液相變模型，補充擬臨界現(xiàn)象認識上的理論空白，同時結(jié)合標度理論，推導部分物性在趨于臨界點時的發(fā)散規(guī)律；開展高精度的二氧化碳密度、定壓比熱、黏度物性測量實驗，指出現(xiàn)有計算模型的缺陷。

1 分子動力學模擬

本文采用Material Studio軟件中的分子動力學模擬功能，研究二氧化碳流體體系跨越擬臨界點時的微觀結(jié)構(gòu)特性。

1.1 計算模型

力場模型為COMPASS力場[18]，其在單分子或凝聚態(tài)物質(zhì)結(jié)構(gòu)、熱力學性質(zhì)等方面的預測能力已得到廣泛證實。模擬體系為立方體盒子，在x、y、z方向上均采用周期性邊界條件以消除邊界效應。根據(jù)分子數(shù)量、截斷半徑無關(guān)性、時間步長無關(guān)性的檢驗結(jié)果，同時盡可能減小計算資源的消耗，設(shè)置二氧化碳分子數(shù)為256，截斷半徑為1.55×10-9m，時間步長設(shè)為7×10-16s。

在模擬過程的前3×10-10s，保持體系的粒子數(shù)、體系體積及溫度不變(即正則系綜)，使得體系達到遲豫平衡。此時，模擬結(jié)果的統(tǒng)計量收斂到一個準穩(wěn)態(tài)，而與分子的初始排布方式無關(guān)。然后，保持體系的粒子數(shù)、體系壓力和溫度不變，對體系的密度、雙體分布函數(shù)、配位數(shù)、靜態(tài)結(jié)構(gòu)因子等物理量進行統(tǒng)計。模擬過程中的溫度控制方法為Berendsen控溫法(即恒溫熱浴耦合法)，壓力控制方法為隨機碰撞法。對于分子之間的非鍵相互作用，采用Group-Based加和方法進行求解。另外，還可根據(jù)能量波動方法求得定壓比熱，通過統(tǒng)計體系內(nèi)粒子數(shù)目的波動求得等溫壓縮系數(shù)等[19]。

1.2 計算結(jié)果

圖1為5 MPa下二氧化碳密度的分子動力學模擬結(jié)果與NIST數(shù)據(jù)庫值之間的偏差?？煽闯觯肿觿恿W模擬能準確預測流體的密度，同時很好地捕捉到該壓力下的氣液相變點，最大相對偏差僅為5.46%。進一步對7.38 MPa和16 MPa下的二氧化碳密度進行了模擬和對比，結(jié)果表明：對于遠離臨界點的狀態(tài)，分子動力學模擬方法在預測流體密度時精度較高，模擬值與數(shù)據(jù)庫值的相對偏差小于5%；在臨界點附近計算偏差顯著上升，但仍有相同的變化趨勢。因此，本文所搭建的分子動力學模型可作為定性研究物性畸變現(xiàn)象的手段。

圖1 5 MPa下二氧化碳密度模擬結(jié)果與NIST數(shù)據(jù)庫值的對比

以雙體分布函數(shù)g(r)為工具，研究體系的微觀特性。作為分子排列有序性的一種度量，g(r)表示以一個任意指定的粒子為中心，在距離r處找到其他粒子的概率。圖2為7.38 MPa下T=280～350 K的二氧化碳體系的C-C原子雙體分布函數(shù)的變化情況，作為對比，圖2中虛線分別為280 K、0.1 MPa下二氧化碳和氦氣的g(r)曲線。從圖2可看出，隨溫度的升高，雙體分布函數(shù)的第1峰與第2峰之間的峰谷逐漸升高，說明體系的有序程度不斷減弱、分子排列變得稀疏。當溫度高于310 K后，g(r)曲線中僅存在第1峰，而第1峰之外，與280 K、0.1 MPa下的氣態(tài)二氧化碳情形十分相似。這說明體系中僅存在嚴格的近程有序，體系的微觀狀態(tài)更接近于氣態(tài)。但與氦氣情形下的近程無序且長程無序不同，二氧化碳在氣態(tài)條件下仍是近程有序的，這實際上反映了二氧化碳分子不為0的四極矩之間的相互作用。當溫度小于304.13 K或大于305 K時，體系的g(r)曲線基本重合；而在304.13～305 K的溫度區(qū)間內(nèi)，曲線特征發(fā)生顯著的轉(zhuǎn)變。

圖2 7.38 MPa下T=280～350 K時C-C原子間雙體分布函數(shù)對比

對于一定壓力下的超臨界流體，溫度較低時體系的有序程度較高，表現(xiàn)為近程、中程有序而長程無序，與液體類似；溫度較高時體系表現(xiàn)為近程有序而中長程無序，與氣體相似；在中間的某一狹窄溫度區(qū)間內(nèi)，體系的結(jié)構(gòu)發(fā)生劇烈的變化。因此，可進一步將超臨界流體區(qū)劃分為類氣區(qū)、類液區(qū)和擬臨界區(qū)，這與散射實驗[11]的測量結(jié)果也是吻合的。

2 考慮外場的氣液相變模型

分子動力學模擬從微觀視角上給出了擬臨界現(xiàn)象的直觀認知，但仍無法解釋超臨界流體發(fā)生擬臨界現(xiàn)象的原因。為解釋這一問題，本文基于朗道二級相變理論[20-21]和標度理論，提出一種考慮外場的氣液相變模型。

2.1 考慮外場的朗道模型

根據(jù)朗道理論在其他體系的推廣形式[22-24]，如描述液晶體系相變的Landau-de Gennes模型[25]，本文提出了一種新模型。在新模型中，序參量η為(ρ-ρc)/ρc，其中ρ為流體的平均溫度，下標c表示臨界點狀態(tài)。在這種定義下，η在整個流體域均有明確的定義，且在臨界等容線(ρ=ρc)上等于0。此外，朗道理論中熱力學勢Φ的表達式也按實際流體的特性進行了改進。首先定義一組無量綱變量(p,t,η)，其中p=(P-Pc)/Pc為無量綱壓力，而t=(T-Tc)/Tc為無量綱溫度。新的熱力學勢表達式為：

Φ(p,t,η)=Φ0(p,t)+C(-(p-bt)η+atη2+Bη4)

(1)

式中：a、b、B、C均為大于0的系數(shù)，b、C為常數(shù)，a、B僅與壓力相關(guān)；Φ0(p,t)為與流體溫度和壓力相關(guān)的常規(guī)項。

相比于經(jīng)典朗道理論，新模型基于重新定義的序參量η，同時增加了η的一次項。此一次項表示了外場h=p-bt對熱力學勢的作用，也對應于溫度和壓力對于流體狀態(tài)的影響。根據(jù)朗道理論，一個穩(wěn)定的熱力學狀態(tài)(p,t)對應于(?Φ/?η)p,t=0，由式(1)則有：

h=p-bt=2atη+4Bη3

(2)

因此，對于臨界等容線η=0，流體的壓力與溫度應滿足線性關(guān)系(圖3)。無論是對于水、二氧化碳、氮氣還是氧氣，不同流體的臨界等容線均在很大的溫度壓力范圍內(nèi)保持線性，進而證明了引入的線性外場項的正確性。

圖3 4種不同流體的臨界等容線

由式(1)的熱力學勢表達式可推導出流體的定壓比熱cp為：

cp=T(?s/?T)P

(3)

(4)

由式(4)可知，定壓比熱cp由常規(guī)項cp0和擬臨界增強項兩部分構(gòu)成。在擬臨界區(qū)，密度隨溫度的劇烈變化導致擬臨界增強項的增加。圖4以23 MPa下的超臨界水對式(4)進行了驗證?？煽闯?，在接近臨界密度時，定壓比熱的畸變與?ρ/?T的變化高度相關(guān)。

圖4 臨界等容線附近超臨界水的定壓比熱畸變規(guī)律

2.2 考慮外場的標度理論

除了可在定量上解釋擬臨界區(qū)的物性畸變特性，本文所提出的模型亦適用于亞臨界氣液相變。但對于臨界點鄰域，模型的預測值與實驗測量值存在顯著的偏差。以臨界等溫線和臨界等壓線為例，則：

(5)

對應的臨界指數(shù)δ=3與實驗測量的結(jié)果δ=4.8±0.02不相符。針對此問題，標度理論和重整化群理論的結(jié)論為：在臨界點附近，體系內(nèi)部漲落的影響不可忽略，導致忽略了漲落的朗道理論不再適用。此時不能遵循平均場近似，需考慮體系內(nèi)部劇烈漲落對熱力學勢的影響，因此有必要對模型在臨界點鄰域進行修正。根據(jù)標度理論，趨于臨界點時體系的關(guān)聯(lián)長度發(fā)散，此時系統(tǒng)具有標度不變性，熱力學函數(shù)可轉(zhuǎn)換為系統(tǒng)參數(shù)的廣義齊次函數(shù)f。對于外場影響下的序參量，則：

(6)

式中：±分別對應η>0和η<0；β為熱膨脹系數(shù)。

(7)

3 超臨界二氧化碳高精度物性測量實驗

在宏觀尺度上，目前超臨界二氧化碳擬臨界區(qū)熱物理性質(zhì)實驗數(shù)據(jù)稀少，且不同來源實驗數(shù)據(jù)間的一致性也較差。因此，本文聯(lián)合西安交通大學和大連理工大學開展了兩組“背靠背“的物性測量實驗進行驗證。以二氧化碳為實驗工質(zhì)，分別采用U型管振蕩法[26]、流動型量熱法[27-28]、毛細管法[29]，完成25～500 ℃溫度范圍、7～14 MPa壓力范圍內(nèi)的密度、定壓比熱和黏度的測量，共得到數(shù)據(jù)點1 940個，圖5為實驗裝置示意圖。表1列出在置信因子k=2、置信度為95%的條件下兩組實驗中壓力、溫度、定壓比熱及黏度的測量不確定度。

表1 物性測量實驗的不確定度

a——定壓比熱實驗測量系統(tǒng)；b——黏度和密度實驗測量系統(tǒng)

從“背靠背”驗證的角度來看，在相同工況下，兩組密度測量結(jié)果的平均相對偏差僅為0.127%，吻合較好。將實驗結(jié)果與當前應用最為廣泛的Span-Wagner方程[30]和Vesovic黏度計算模型[31]的計算結(jié)果進行對比，對于物性變化較為平緩的一般區(qū)域，實驗結(jié)果與計算結(jié)果的偏差很小(表2)。但對于物性劇烈變化的擬臨界區(qū)，實驗結(jié)果與計算結(jié)果的偏差顯著增大，且表現(xiàn)出物性畸變越劇烈則偏差越大的趨勢(表3)?，F(xiàn)有的擬臨界區(qū)數(shù)據(jù)點仍較少，數(shù)據(jù)點壓力偏高且溫度間隔過大，實驗精度需進一步評估和優(yōu)化。

表2 一般區(qū)域內(nèi)實驗測量結(jié)果與計算結(jié)果的相對偏差

表3 擬臨界區(qū)部分實驗測量結(jié)果與計算結(jié)果的偏差

現(xiàn)有的二氧化碳熱物理性質(zhì)理論模型在擬臨界區(qū)之外計算精度較好，但在近臨界區(qū)仍存在缺陷，計算精度需實驗驗證并進一步提高。對于超臨界二氧化碳布雷頓循環(huán)，為了獲得更高的系統(tǒng)效率，往往需要盡可能將背壓設(shè)置得低(如7.4 MPa)，利用這種畸變特性來減小壓縮功。因此，掌握二氧化碳熱物性在近/超臨界區(qū)域精確的實驗數(shù)據(jù)及其變化規(guī)律是一個關(guān)鍵問題，需要在實驗技術(shù)、理論模型等方面進一步研究。

4 超臨界流體區(qū)的嚴格劃分

通過微觀、介觀、宏觀3個不同尺度的研究分析，得到超臨界流體發(fā)生物性畸變的機理。分子動力學模擬的結(jié)果表明，超臨界流體在擬臨界區(qū)發(fā)生物性畸變時，同時伴隨類氣-類液結(jié)構(gòu)之間的轉(zhuǎn)變，而在漲落主導的臨界點鄰域內(nèi)，流體分子之間存在長程關(guān)聯(lián)。長程關(guān)聯(lián)帶來的體系自相似性為采用粗?；椒ń硕壤碚撎峁┝艘罁?jù)，但也導致傳統(tǒng)連續(xù)體假設(shè)的失效。因此，研究處于臨界點鄰域的流體必須從更為基礎(chǔ)的假設(shè)出發(fā)，采用分子動力學模擬或格子玻爾茲曼方法[32]。隨著遠離臨界點，漲落迅速衰減以至于可忽略，此時基于平均場思想的朗道理論變得適用。同時由于不存在長程關(guān)聯(lián)，傳統(tǒng)的宏觀描述(如Navier-Stokes方程)同樣成立，但需注意物性畸變對于流動、傳熱等宏觀行為的影響。因此，超臨界流體區(qū)域可進一步劃分為4個子區(qū)域(圖6)：臨界點鄰域、類氣區(qū)、類液區(qū)和擬臨界區(qū)。

圖6 臨界等容線附近超臨界水的定壓比熱畸變規(guī)律

關(guān)于各區(qū)域的劃分，首先可估算臨界點鄰域的范圍。當沿臨界等容線趨近臨界點時，關(guān)聯(lián)長度ξ的發(fā)散遵循[33]：

(8)

當|t|足夠大時，流體分子之間無長程關(guān)聯(lián)，而只有分子間作用力決定的短程關(guān)聯(lián)，因此ξ≈ξ0≈10-10m。只有當t≈10-6時，關(guān)聯(lián)長度才會增長到與流體層厚度相當?shù)摩蘭量級。對于二氧化碳，t≈10-6對應于T-Tc≈3×10-4K。因此，臨界點鄰域的范圍是很小的。對于大部分的超臨界流體應用場景而言，傳統(tǒng)的宏觀描述均適用。對于其他3個區(qū)域，Banuti[34]推廣了亞臨界條件下的克拉伯龍方程，得到關(guān)于超臨界流體的劃分方法；Kurganov等[35]則從流體壓縮性的評估出發(fā)，以相對膨脹功Eq=pβ/ρcp作為劃分依據(jù)?？偟貋碚f，擬臨界區(qū)的劃分更多是經(jīng)驗性的，因此兩種劃分方法都可行。在實際工程中，為了使用方便，也可根據(jù)第2節(jié)中的結(jié)論，以無量綱密度作為劃分擬臨界區(qū)的簡單判據(jù)。

5 結(jié)論

本文通過微觀、介觀和宏觀3個不同尺度的分析論證，發(fā)現(xiàn)了導致超臨界流體物性畸變的不同機制。根據(jù)影響機制的不同，在相圖上，超臨界流體區(qū)可進一步細分為臨界點鄰域、擬臨界區(qū)、類氣區(qū)和類液區(qū)，得到如下結(jié)論：1) 分子動力學模擬的結(jié)果表明，在跨越擬臨界線時，超臨界流體的微觀構(gòu)型會發(fā)生類氣-類液之間的轉(zhuǎn)變；2) 在臨界點鄰域內(nèi)，體系的熱力學函數(shù)、漲落和關(guān)聯(lián)長度趨于臨界點時按冪律發(fā)散，此時考慮外場的標度理論可很好地描述，但傳統(tǒng)的宏觀描述和平均場理論均失效；3) 臨界點鄰域的范圍很小，以無量綱溫度t=(T-Tc)/Tc表示時，近似滿足t<10-6；4) 對于絕大部分工程應用場景，超臨界流體內(nèi)部的漲落影響均可忽略，此時考慮外場的氣液相變模型能很好地描述物性畸變的規(guī)律；5) 無量綱密度|η|<0.25可作為劃分擬臨界區(qū)的簡單判據(jù)，對于二氧化碳，對應密度為350～585 kg/m3；6) 超臨界二氧化碳的物性測量實驗結(jié)果表明，現(xiàn)有的物性模型在描述擬臨界區(qū)物性畸變特性方面仍有不足，但對于擬臨界區(qū)之外的一般區(qū)域已有足夠高的計算精度。