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

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

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

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

1 分子動(dòng)力學(xué)模擬

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

1.1 計(jì)算模型

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

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

1.2 計(jì)算結(jié)果

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

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

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

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

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

2 考慮外場(chǎng)的氣液相變模型

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

2.1 考慮外場(chǎng)的朗道模型

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

Φ(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ī)項(xiàng)。

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

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

(2)

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

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

由式(1)的熱力學(xué)勢(shì)表達(dá)式可推導(dǎo)出流體的定壓比熱cp為：

cp=T(?s/?T)P

(3)

(4)

由式(4)可知，定壓比熱cp由常規(guī)項(xiàng)cp0和擬臨界增強(qiáng)項(xiàng)兩部分構(gòu)成。在擬臨界區(qū)，密度隨溫度的劇烈變化導(dǎo)致擬臨界增強(qiáng)項(xiàng)的增加。圖4以23 MPa下的超臨界水對(duì)式(4)進(jìn)行了驗(yàn)證?？煽闯觯诮咏R界密度時(shí)，定壓比熱的畸變與?ρ/?T的變化高度相關(guān)。

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

2.2 考慮外場(chǎng)的標(biāo)度理論

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

(5)

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

(6)

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

(7)

3 超臨界二氧化碳高精度物性測(cè)量實(shí)驗(yàn)

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

表1 物性測(cè)量實(shí)驗(yàn)的不確定度

a——定壓比熱實(shí)驗(yàn)測(cè)量系統(tǒng)；b——黏度和密度實(shí)驗(yàn)測(cè)量系統(tǒng)

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

表2 一般區(qū)域內(nèi)實(shí)驗(yàn)測(cè)量結(jié)果與計(jì)算結(jié)果的相對(duì)偏差

表3 擬臨界區(qū)部分實(shí)驗(yàn)測(cè)量結(jié)果與計(jì)算結(jié)果的偏差

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

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

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

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

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

(8)

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

5 結(jié)論

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