基于BWR狀態(tài)方程的超臨界二氧化碳主要熱力參數(shù)推算

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

（天津大學(xué) 機(jī)械工程學(xué)院，天津 300350）

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

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

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

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

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

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

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

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

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

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

2.1 基準(zhǔn)點(diǎn)與積分路徑

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

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

2.2 內(nèi)能

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

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

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

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

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

式中，參數(shù)通過實(shí)驗(yàn)數(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的具體表達(dá)式，代入式（5）積分求得內(nèi)能差值，而u1為基準(zhǔn)點(diǎn)內(nèi)能（已知），即可得到待定狀態(tài)2內(nèi)能的快速計(jì)算式，經(jīng)過換算最終單位為kJ/kg。

2.3 焓

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

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

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

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

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

2.4 熵

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

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

2.5 實(shí)際比熱容

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

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

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

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

3 結(jié)果與討論

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

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

3.1 內(nèi)能

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

3.2 焓

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

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

3.3 熵

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

3.4 實(shí)際比熱容

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

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

圖4 關(guān)鍵熱力區(qū)cp計(jì)算值及參考值（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ì)，本工作的推導(dǎo)結(jié)果與REFPROP數(shù)據(jù)的誤差在高于320 K的溫度都比較小，整體的數(shù)據(jù)分布趨勢符合度高。在臨界溫度附近（300～320 K），物性的變化趨勢明顯，雖然該區(qū)域的物性擬合難度較大，但結(jié)果還比較可觀。熱力性質(zhì)快速解析計(jì)算模型在不同區(qū)域的平均誤差見表2。

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

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

4 結(jié)論

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

0.52%。

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

符 號 說 明

A0，B0，C0，a，a0，a1，a2，a3，b，c，α，γ經(jīng)驗(yàn)參數(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壓縮因子

ε誤差

下標(biāo)

cal 計(jì)算值

ref 參考值

0 理想狀態(tài)

1 狀態(tài)1

2 狀態(tài)2