基于BWR狀態方程的超臨界二氧化碳主要熱力參數推算

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

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

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

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

1 BWR狀態方程及驗證

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

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

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

表1 BWR狀態方程經驗參數Table 1 Empirical parameter of BWR equation of state

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

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

2 SCO2主要熱力參數推算

2.1 基準點與積分路徑

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

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

2.2 內能

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

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

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

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

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

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

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

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

2.3 焓

焓是由開口系工質流動的特殊性而產生的一個新的狀態參數，是熱力系的熱力學能與流動功之和。第二焓方程[14]見式（10）。

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

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

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

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

2.4 熵

熵的物理意義是體系混亂程度的度量，是構成體系的大量微觀離子集體表現出來的性質。第一熵方程見式（16）。

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

2.5 實際比熱容

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

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

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

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

3 結果與討論

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

由于CO2在臨界點附近物性的劇烈變化，提取溫度300～320 K，壓力7～8 MPa范圍內數據，為近臨界區結果比較分析。

3.1 內能

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

3.2 焓

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

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

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

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

3.3 熵

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

3.4 實際比熱容

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

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

圖4 關鍵熱力區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.

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

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

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

4 結論

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

0.52%。

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

符 號 說 明

A0，B0，C0，a，a0，a1，a2，a3，b，c，α，γ經驗參數

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

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

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

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

h焓，kJ/kg

p壓力，Pa

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

s熵，kJ/(kg·K)

T溫度，K

u內能，kJ/kg

v比體積，cm3/kmol

Z壓縮因子

ε誤差

下標

cal 計算值

ref 參考值

0 理想狀態

1 狀態1

2 狀態2