熊先闖，劉正先，趙 明，李孝檢

（天津大學 機械工程學院，天津 300350）

在壓縮機的氣動設計和流場模擬中，氣體的熱力參數(shù)計算是關鍵步驟之一，具體包括壓力、溫度、比體積（v）、焓、內(nèi)能和熵等參數(shù)的計算。當壓縮工質(zhì)為實際氣體時，如碳捕獲與存儲、超臨界二氧化碳（SCO2）布雷頓動力循環(huán)中的壓縮系統(tǒng)，由于存在實際氣體效應，工質(zhì)熱力性質(zhì)顯著區(qū)別于理想氣體[1]。如何快速準確地計算實際氣體熱力參數(shù)是開展以實際氣體為工質(zhì)的壓縮系統(tǒng)熱力設計及數(shù)值模擬的重要前提。SCO2是一種典型的實際氣體。它具有較大的壓縮系數(shù)、比熱容和密度以及較小的黏性。在SCO2布雷頓循環(huán)中，流體工質(zhì)較大的壓縮系數(shù)使得壓縮機消耗的功率較小[2-3]；又因工質(zhì)高密度的性質(zhì)，熱力循環(huán)的幾何結(jié)構(gòu)變得更為緊湊，僅為蒸氣循環(huán)的十分之一。此外，CO2的臨界溫度較低，相對于水等流體工質(zhì)更容易達到臨界狀態(tài)。因此，以SCO2為工質(zhì)的熱力循環(huán)將有潛力取代傳統(tǒng)水蒸氣朗肯循環(huán)。實際氣體熱力參數(shù)計算的主要任務是依據(jù)基本熱力參數(shù)（壓力、溫度、比體積或密度）推算出其他主要熱力參數(shù)（內(nèi)能、焓、熵等）及其衍生參數(shù)。描述基本熱力參數(shù)之間函數(shù)關系的方程叫做狀態(tài)方程。常見的狀態(tài)方程有范德瓦爾方程、RK方程和PR方程等。不同狀態(tài)方程適用的工質(zhì)類型、溫度、壓力范圍各不相同[4]。如RK方程對非極性或輕微極性氣體pVT性質(zhì)的計算結(jié)果符合較好，密度最大適用值可略大于臨界密度，但如果含氫鍵，則會導致很大的誤差。Span-Wagner（SW）模型是一種專門針對SCO2建立的以實驗數(shù)據(jù)標定的半經(jīng)驗解析模型[5]。對于以SCO2為工質(zhì)的熱力系統(tǒng)，如果直接采用SW解析模型，那么計算過程將變得十分繁瑣，眾多的經(jīng)驗參數(shù)及高階非線性特性使得它幾乎無法直接推導出其他衍生熱力參數(shù)的解析表達式。一種更普遍的計算方法是采用高階樣條曲線分段插值實驗數(shù)據(jù)庫以得到工質(zhì)的物性參數(shù)表格[6-7]。NIST（美國國家標準與技術(shù)研究院）的REFPROP數(shù)據(jù)庫，即為基于SW解析模型分段插值的物性表[8-10]。然而，調(diào)用REFPROP數(shù)據(jù)庫生成寬范圍、高精度的物性表，需要消耗較多的時間與內(nèi)存。更重要的是分段插值方法難以得到基本熱力參數(shù)精度較高的一階或二階偏導數(shù)信息，因而無法精確獲得相關的衍生熱力參數(shù)。Benedict-Webb-Rubin（BWR）狀態(tài)方程[11]是精度最高的通用狀態(tài)方程之一，它考慮了分子的聚集，能夠較為準確地描述SCO2的基本熱力參數(shù)。BWR方程具有物理意義明確、解析表達式簡潔等優(yōu)勢?；谠摖顟B(tài)方程及熱力學普遍關系式（熱力學第一定律、麥克斯韋關系式等）可推導出基本熱力參數(shù)及其衍生參數(shù)快速準確的解析計算模型。解析計算模型對于實際氣體熱物理特性的研究具有重要理論意義，同時對于以SCO2為工質(zhì)的熱力過程設計與計算具有重要的工程應用價值。

本工作從熱力基本關系式出發(fā)，結(jié)合BWR狀態(tài)方程，采用了以標準狀態(tài)為基準點，理想狀態(tài)為中間過程的積分路徑，推導了SCO2主要的熱力參數(shù)快速計算解析表達式；并通過與REFPROP數(shù)據(jù)對比驗證了解析模型在不同熱力區(qū)的可靠性。

1 BWR狀態(tài)方程及驗證

狀態(tài)方程的選擇很重要，直接影響導出的物性模型精度，BWR方程能準確描述超臨界流體的pVT行為，具體形式見式（1）。

式中，p為壓強，105Pa；T為溫度，K；v為比體積，cm3/kmol；R為通用氣體常數(shù)，數(shù)值為8.314（cm3·MPa）/（kmol·K）-1；A0，B0，C0，a，b，c，α，γ均為經(jīng)驗參數(shù)。

在CO2溫度范圍300～600 K、壓力范圍7～30 MPa下，對應的經(jīng)驗參數(shù)值由文獻[12]得到，見表1。

表1 BWR狀態(tài)方程經(jīng)驗參數(shù)Table 1 Empirical parameter of BWR equation of state

為了驗證擬合方程和參數(shù)質(zhì)量的好壞，需要進行回代檢驗。在溫度300～600 K、壓力7～30 MPa內(nèi)取樣，通過REFPROP數(shù)據(jù)得到v，并進行相應的單位轉(zhuǎn)換。將溫度和v代入式（2），通過MATLAB編程計算，最后得到壓力計算值（pcal）與REFPROP數(shù)據(jù)（pref）進行比較分析。pcal與pref在整個熱力學區(qū)的相對誤差在2.6%以內(nèi)，平均相對誤差僅為0.51%。參數(shù)質(zhì)量及擬合方程對pVT行為描述較好，利用賦予這些參數(shù)值的BWR方程，進一步推算熱力性質(zhì)上是可靠的。本工作所有參考數(shù)據(jù)都是通過REFPROP軟件查詢得到的數(shù)據(jù)。壓縮因子（Z）主要被定義為實際體積與理想體積之比，實際氣體的狀態(tài)方程引入了Z，簡化形式見式（2）。結(jié)合式（1）與式（2），得到Z具體表達式（3）。

為了驗證式（3）的精確性，將Z計算值與參考值進行比較。由BWR方程推算的Z最大誤差不超過2.6%，全局平均誤差0.52%，誤差分布也比較平緩，無明顯偏差點。

2 SCO2主要熱力參數(shù)推算

2.1 基準點與積分路徑

實際氣體的熱力性質(zhì)為狀態(tài)量，與過程無關。但是推導方法應用了熱力普遍關系式，需要通過一定路徑進行積分以及合適的基準點差值計算，才能保證結(jié)果的可靠性[13]。

以標準狀態(tài)1（T1=298 K，p1=105Pa）為基準點，基準點熱力學參數(shù)已知。積分路徑常用的是從標準狀態(tài)出發(fā)，先定壓升溫至待定狀態(tài)溫度T2，再定溫升壓至待定壓力p2，但是由于積分量中涉及到比熱項，是溫度的函數(shù)，在不同壓力下的表達式也不盡相同。為了避免被積函數(shù)的不定性，將比熱函數(shù)設定為理想狀態(tài)。先計算狀態(tài)1等溫膨脹到理想狀態(tài)（p→0，v→∞），再對理想狀態(tài)下的升溫過程（T1→T2）進行計算，最后等溫升壓到待求狀態(tài)2（T2，p2）。三個過程量疊加即為狀態(tài)2與狀態(tài)1的差值，狀態(tài)1數(shù)據(jù)可由REFPROP直接查得。

2.2 內(nèi)能

內(nèi)能代表了儲存于熱力體系內(nèi)部的能量，第一熱力學能方程[14]見式（4）。

式中，u為內(nèi)能，kJ/kg；cv為恒容摩爾比熱容，kJ/（kg·K）。

對式（4）從狀態(tài)1到狀態(tài)2的路徑進行積分，并計算得到差值，見式（5）。

則式（5）中第一、三項中被積函數(shù)見式（7）。

第二項中被積函數(shù)為理想恒容摩爾比熱容（cv∞）。理想恒壓摩爾比熱容（cp0）計算式見式（8）。

式中，參數(shù)通過實驗數(shù)據(jù)擬合得到a0= 1.782×102，a1= 8.677×10-1，a2= -8.465×10-4，a3= 3.667×10-7。

由邁爾公式（9）得到cv∞。

被積函數(shù)已全部轉(zhuǎn)換成T和v的具體表達式，代入式（5）積分求得內(nèi)能差值，而u1為基準點內(nèi)能（已知），即可得到待定狀態(tài)2內(nèi)能的快速計算式，經(jīng)過換算最終單位為kJ/kg。

2.3 焓

焓是由開口系工質(zhì)流動的特殊性而產(chǎn)生的一個新的狀態(tài)參數(shù)，是熱力系的熱力學能與流動功之和。第二焓方程[14]見式（10）。

對式（10）沿指定路徑積分，從狀態(tài)1到狀態(tài)2焓變見式（11）。

由于d（pv）=pdv+vdp，剩余項一般形式為式（13）。

將式（12）和式（13）合并即可得到式（10）中第一、三項的一般形式，見式（14）。

式（14）按照一般形式拓展即可得到可計算的焓差見式（15）。

2.4 熵

熵的物理意義是體系混亂程度的度量，是構(gòu)成體系的大量微觀離子集體表現(xiàn)出來的性質(zhì)。第一熵方程見式（16）。

沿指定路徑積分，得到狀態(tài)2與狀態(tài)1之間的熵增，見式（17）。

2.5 實際比熱容

氣體在壓力和容積不變的條件下被加熱時的比熱容，分別為恒壓摩爾比熱容（cp）和cv。實際cv定義式為式（18）。

其中快速計算式已由式（5）得到，為T和v雙變量解析式，經(jīng)直接求偏導得式（19）。

通過熱容容差關系（不同下標關系、循環(huán)關系和麥克斯韋關系）整理可得式（20）。

CO2的實際比熱容已具化成T和v的快速計算式，輸入變量即可與REFPROP數(shù)據(jù)對比精確性，換算后最終單位為kJ/（kg·K）。

3 結(jié)果與討論

對上述推導得到的熱力性質(zhì)快速計算式，進行編程得到結(jié)果。為了分析計算結(jié)果在關鍵熱力學區(qū)的可靠性，即溫度T為300～600 K，壓力p為7～30 MPa，比體積v大概在50～700 cm3/kmol之間。此區(qū)域橫跨了臨界點，包括了部分超臨界區(qū)與部分亞臨界區(qū)。取樣繪圖，由于此區(qū)域v跨度較大，而且v不隨壓力線性變化，采用常規(guī)方法采樣無法突顯CO2所處熱力狀態(tài)。故采用設置不同壓力間接得到v的方法，得到對應的v，再將T和v代入快速計算式，得到對應的熱力數(shù)據(jù)。繪制熱力性質(zhì)在不同溫度下隨v變化的關系曲線，并將計算值與REFPROP數(shù)據(jù)進行誤差比較。

由于CO2在臨界點附近物性的劇烈變化，提取溫度300～320 K，壓力7～8 MPa范圍內(nèi)數(shù)據(jù)，為近臨界區(qū)結(jié)果比較分析。

3.1 內(nèi)能

圖1 為關鍵熱力區(qū)及近臨界區(qū)的內(nèi)能計算值與參考值隨v變化及誤差分布。由圖1可知，在關鍵熱力學區(qū)，內(nèi)能計算結(jié)果與REFPROP數(shù)據(jù)最大誤差不超過2.7%，平均誤差僅為0.65%；近臨界點附近，最大誤差不超過1.9%，平均誤差0.61%。在相同v條件下，誤差在320 K附近降至最低，在305 K附近會略微升高。在300～320 K區(qū)間，誤差會隨著v的升高而減小。在溫度大于305 K的超臨界區(qū)，誤差基本穩(wěn)定，溫度及v的影響可以忽略，接近零。當v＞100 cm3/kmol時，誤差降至低水平，且受溫度影響變小。說明了內(nèi)能計算結(jié)果的可靠性，對于采用CO2為工質(zhì)的數(shù)值模擬與工程實踐中，內(nèi)能可采用本工作數(shù)據(jù)。而且在SCO2的應用中，結(jié)果具有更好地適用性。

3.2 焓

圖2 為關鍵熱力區(qū)及近臨界區(qū)的焓計算值與參考值隨v變化的誤差分布。

圖1 關鍵熱力區(qū)（a）和近臨界區(qū)（b）內(nèi)能相對誤差Fig.1 Relative error of internal energy(εu) in key thermal area(a) and in near critical area(b).v：specific volume.

圖2 關鍵熱力區(qū)（a）和近臨界區(qū)（b）焓相對誤差Fig.2 Relative error of enthalpy(εh) in key thermal area(a) and in near critical area(b).

由圖2可知，由于焓與內(nèi)能的推導表達式類似，所以整體的數(shù)據(jù)分布差異不大，但也有一定區(qū)別。在關鍵熱力學區(qū)，最大誤差不超過2.5%，平均誤差0.6%；近臨界區(qū)最大誤差不超過2.3%，平均誤差0.58%。計算可知，焓的誤差在可接受范圍內(nèi)，對于不同v區(qū)的精度分布與內(nèi)能接近，即適用范圍也相似。

3.3 熵

圖3 為關鍵熱力區(qū)及近臨界區(qū)的熵計算值與參考值隨v變化的誤差分布。由圖3可知，在關鍵熱力學區(qū)，最大誤差達到了15%，平均誤差5.2%；在近臨界區(qū)最大誤差10%，平均誤差5.3%。熵的計算誤差較內(nèi)能及焓普遍偏大，這與BWR方程本身以及推導方法有關。相同溫度條件下，隨著v的升高，誤差逐漸減小，趨于穩(wěn)定，最低能為1%。相反在v相同情況下，溫度越高，誤差也會隨之增高。當v＞100 cm3/kmol時，誤差降至低水平，變化平緩，且受溫度影響不大。所以對于需要CO2熵的情況下，需要考慮工質(zhì)的工作溫度及v，以及項目對精度的要求。

3.4 實際比熱容

對于實際比熱容，推導過程中涉及到BWR方程的二階偏導，精度相比其他熱力性質(zhì)會偏低，所以不對其近臨界區(qū)做分析，但是給出計算數(shù)據(jù)與參考數(shù)據(jù)具體分布圖可直觀看出趨勢。圖4為關鍵熱力區(qū)cp計算值及參考值和cp及cv的計算相對誤差。由圖4可知，在低溫區(qū)（低于320 K）的差異較大，當v＞150 cm3/kmol時，比熱容值變化較平緩，且計算誤差也極小。但在305 K時，誤差變化劇烈，所以當工況溫度在305 K左右時，需要視精度要求使用。cv與cp平均誤差分別為5.39%，3.28%。

圖3 關鍵熱力區(qū)（a）和近臨界區(qū)（b）熵相對誤差Fig.3 Relative error of entropy(εs) in key thermal area(a) and in near critical area(b).

圖4 關鍵熱力區(qū)cp計算值及參考值（a-b）和cp及cv相對誤差（c-d）Fig.4 Calculated and reference value of cp(a-b) and relative error of cp(c) and cv(d) in key thermal area.

對于上述熱力性質(zhì)，本工作的推導結(jié)果與REFPROP數(shù)據(jù)的誤差在高于320 K的溫度都比較小，整體的數(shù)據(jù)分布趨勢符合度高。在臨界溫度附近（300～320 K），物性的變化趨勢明顯，雖然該區(qū)域的物性擬合難度較大，但結(jié)果還比較可觀。熱力性質(zhì)快速解析計算模型在不同區(qū)域的平均誤差見表2。

表2 快速解析計算模型平均誤差Table 2 Average error of fast analytical calculation model

由表2可知，內(nèi)能、焓、熵、cp及cv在關鍵熱力區(qū)的平均誤差分別為0.65%，0.60%，5.2%，5.39%，3.28%。

4 結(jié)論

1）對溫度300～600 K、壓力7～30 MPa范圍內(nèi)的CO2，首先進行方程參數(shù)質(zhì)量驗證。代入REFPROP數(shù)據(jù)，壓力計算值與參考值在關鍵熱力學區(qū)的平均誤差僅為0.5%，參數(shù)質(zhì)量以及擬合方程對pVT行為描述較好。對實際氣體引入Z，Z的計算結(jié)果與REFPROP數(shù)據(jù)符合良好，平均誤差

0.52%。

2）推算了熱力性質(zhì)的快速計算定量表達式，包括內(nèi)能、焓、熵、比熱容。對結(jié)果進行了全局比較以及近臨界區(qū)的加密取樣比較，平均誤差都在可接受范圍以內(nèi)。內(nèi)能、焓、熵、cp及cv在關鍵熱力區(qū)的平均誤差分別為0.65%，0.60%，5.2%，5.39%，3.28%。對于v＞100 cm3/kmol狀態(tài)，尤其是溫度高于320 K的SCO2，本工作計算結(jié)果均與REFPROP數(shù)據(jù)接近，可適用于高精度計算要求。

符 號 說 明

A0，B0，C0，a，a0，a1，a2，a3，b，c，α，γ經(jīng)驗參數(shù)

cp恒壓摩爾比熱容，kJ/(kg·K)

cp0理想恒壓摩爾比熱容，kJ/(kg·K)

cv恒容摩爾比熱容，kJ/(kg·K)

cv∞理想恒容摩爾比熱容，kJ/(kg·K)

h焓，kJ/kg

p壓力，Pa

R通用氣體常數(shù)，8.314 (cm3·MPa)/(kmol·K)-1

s熵，kJ/(kg·K)

T溫度，K

u內(nèi)能，kJ/kg

v比體積，cm3/kmol

Z壓縮因子

ε誤差

下標

cal 計算值

ref 參考值

0 理想狀態(tài)

1 狀態(tài)1

2 狀態(tài)2