大工況范圍內(nèi)超臨界區(qū)二氧化碳物性快速計算方法研究

文_彭兆睿 鄭秋云 張信榮

1 北京大學工學院 2 湖南文理學院 3 北京市城市熱管理工程技術(shù)研究中心

傳統(tǒng)氟利昂類及烷烴類制冷劑因其較高的臭氧破壞潛能（ODP）和溫室效應潛能（GWP）逐漸在制冷供熱領(lǐng)域被限制使用，而天然工質(zhì)二氧化碳（CO2）由于具有環(huán)保、無毒、不可燃的特性逐步受到關(guān)注。CO2臨界溫度低，容易在常規(guī)工況范圍內(nèi)跨越臨界點達到超臨界狀態(tài)，利用超臨界CO2無氣液相變、低粘度、高密度及高可壓縮性的特點，可將其應用在熱泵制熱制冷循環(huán)、熱發(fā)電循環(huán)等場景中并獲得較好性能，在“雙碳”目標下顯示出巨大潛力。

在實際蒸汽壓縮循環(huán)或發(fā)電循環(huán)運行過程中，實時獲取物性參數(shù)來評估并優(yōu)化循環(huán)性能是降低能耗的重要步驟，其中關(guān)鍵在于根據(jù)測點溫壓及時、快速、高精度計算出物性參數(shù)。利用狀態(tài)方程方法進行計算時，精度較高但計算量大，隱式方程則需要迭代進一步增大了計算成本，耗費計算時間。擬合關(guān)聯(lián)式法以既有物性庫為數(shù)據(jù)源，不基于理論原理而是采用統(tǒng)計學方法擬合顯式關(guān)聯(lián)式，具有簡化計算、降低成本的優(yōu)勢，在傳統(tǒng)制冷劑（如R12、R22等）應用中已獲得有效證明。對于超臨界CO2，因為其物性在臨界點及擬臨界線附近劇烈非線性變化的特點，增加了擬合的難度和擬合關(guān)聯(lián)式的復雜性，既有研究存在擬合關(guān)聯(lián)式過于復雜、非線性系數(shù)過多、分區(qū)復雜等特點，增加計算成本；同時，不同物性關(guān)聯(lián)式較少，針對各個物性需要進一步發(fā)展。對此，本文提出了一種新的擬合思路，以焓值為例，可獲得大工況范圍內(nèi)適用的非線性系數(shù)少、計算成本低、精度高的顯式關(guān)聯(lián)式，可為相關(guān)工業(yè)應用提供參考。

1 物性范圍

本文針對超臨界區(qū)域內(nèi)的焓值建立物性庫，數(shù)據(jù)涵蓋壓力7.38～22.18MPa，溫度220～550K的大工況范圍，為了準確描述靠近臨界壓力及擬臨界溫度附近的非線性趨勢，數(shù)據(jù)變化梯度較大處所取數(shù)據(jù)點更為密集，對應數(shù)據(jù)點如圖1所示。本文數(shù)據(jù)來源為美國國家標準技術(shù)研究所（NIST）開發(fā)的REFPROP。由圖1可以看出，遠離擬臨界線（比熱最大時對應的溫度構(gòu)成的曲線）時的焓值隨溫度變化呈現(xiàn)出線性趨勢，擬臨界線左側(cè)遠端和右側(cè)遠端具有不同的斜率，跨越擬臨界線附近區(qū)域為兩種斜率的過渡區(qū)域；越靠近臨界壓力，兩遠端斜率相差越大，對應過渡區(qū)域斜率越大，非線性程度越大，導致超臨界區(qū)焓值關(guān)聯(lián)式擬合難度增大。

圖1 本文擬合所取的焓值數(shù)據(jù)點

2 擬合方法

本文本質(zhì)上屬于非線性數(shù)據(jù)擬合問題，由于擬合目標在于獲得在大工況范圍內(nèi)適用的計算成本低的擬合關(guān)聯(lián)式，要求不僅是關(guān)聯(lián)式精度高，而且要關(guān)聯(lián)式分區(qū)少及系數(shù)少，因此采用適當?shù)臄M合關(guān)聯(lián)式形式至關(guān)重要，經(jīng)過對Logistic方程等經(jīng)典S形曲線方程的測試，尚無滿足要求的簡單方程形式。對此，本文非線性擬合方法的基本思想在于，針對強烈非線性區(qū)段數(shù)據(jù)，降低數(shù)據(jù)病態(tài)以及降低數(shù)據(jù)非線性程度，進一步優(yōu)化擬合關(guān)聯(lián)式非線性系數(shù)獲得適宜形式。

首先將超臨界壓力下二氧化碳焓值數(shù)據(jù)及對應的壓力、溫度數(shù)據(jù)進行數(shù)域壓縮，通過數(shù)據(jù)重整，將各參數(shù)數(shù)據(jù)集中在量級為1的正區(qū)間，避免因數(shù)據(jù)量級差異及量級較高造成擬合求解時的系數(shù)矩陣病態(tài)，影響擬合效果。重整結(jié)果如圖2所示，其中Tc和hc分別表示臨界溫度和對應的臨界焓值，為簡明表示，后續(xù)h*=h/hc,T*=T/Tc，P*=P/Pc。

接著進行自適應回歸方程優(yōu)化，通過將擬合函數(shù)拆分為基準函數(shù)和調(diào)整函數(shù)的方法降低數(shù)據(jù)源非線性的影響。根據(jù)數(shù)據(jù)源趨勢，P=22.18MPa下焓值-溫度參數(shù)趨勢接近線性變化如圖2所示，因此以該壓力下擬合得到的擬合函數(shù)f(T*)作為基準函數(shù)，同時定義Pb=22.18MPa為基準壓力。為獲得不同壓力下焓值，這里采用調(diào)整函數(shù)計算不同壓力下焓值與基準壓力下焓值的調(diào)整量，該函數(shù)受壓力和溫度同時主導，為二元擬合函數(shù)g(T*,P*)，調(diào)整量匯成的數(shù)據(jù)點分布如圖3所示。根據(jù)數(shù)據(jù)點趨勢，初步采用類似Planck定律的關(guān)聯(lián)式形式見式（1）、式（2），設(shè)置b1-b8共8個非線性系數(shù)，進一步分別對各系數(shù)進行敏感性分析和擬合效果評價以簡化關(guān)聯(lián)式形式，直至關(guān)聯(lián)式誤差滿足需求。

圖2 數(shù)據(jù)重整后的焓值-溫度-壓力參數(shù)分布及基準函數(shù)、調(diào)整函數(shù)示意圖

另外，根據(jù)圖3的數(shù)據(jù)趨勢還可以發(fā)現(xiàn)，在T*＜0.95區(qū)間內(nèi)，焓值受壓力影響可以忽略不計，因此此段的擬合函數(shù)h*即為基準函數(shù)f(T*)；在T*＞0.95區(qū)間內(nèi)，擬合函數(shù)為基準函數(shù)與調(diào)整函數(shù)疊加之和即h*=f(T*)+g(T*,P*)。

圖3 用于擬合調(diào)整函數(shù)的數(shù)據(jù)點

3 擬合結(jié)果

根據(jù)上述擬合思路，首先基準函數(shù)f(T*)采用三階多項式形式見式（2）。

擬合后系數(shù)分別為

對于調(diào)整函數(shù)g(T*,P*)，對類似Planck定律的關(guān)聯(lián)式形式的各系數(shù)進行敏感性分析和擬合效果評價，簡化系數(shù)后得到如式(3)擬合關(guān)聯(lián)式形式。

擬合后系數(shù)分別為b1=3.8345e10,b2=2839574,b3=-7.5962e10,b4=0.1343。RMSE=0.0571，最大誤差降至約12%（見圖4）。結(jié)果顯示，最大誤差在（擬）臨界溫度附近。該擬合關(guān)聯(lián)式在超臨界區(qū)域大工況范圍內(nèi)將參數(shù)系數(shù)降低至8個，在除0.95＜T*＜1.1的（擬）臨界溫度區(qū)間外誤差均滿足5%的要求，顯示出較高精度。

圖4 擬合關(guān)聯(lián)式結(jié)果（圖中實線，左軸）與標準數(shù)據(jù)（圖中點，左軸）對比及兩者之間相對誤差（圖中虛線，右軸），1＜ P＊＜3，T＊＞0.95

進一步對比通過不同途徑計算焓值所需的時間，圖5顯示了分別調(diào)用NIST，調(diào)用本文擬合關(guān)聯(lián)式，調(diào)用既有文獻關(guān)聯(lián)式的時間及對應誤差。結(jié)果表明，采用本文關(guān)聯(lián)式可以有效降低計算量，計算成本為直接調(diào)用NIST的0.2%，為既有關(guān)聯(lián)式的20%，且誤差MRE僅為既有關(guān)聯(lián)式的2倍，加之本文關(guān)聯(lián)式適用范圍更廣，體現(xiàn)出在應用上更具有優(yōu)勢。

圖5 通過不同途徑計算焓值的時間對比，從左到右依次為：直接調(diào)用NIST，調(diào)用本文擬合關(guān)聯(lián)式，調(diào)用既有文獻關(guān)聯(lián)式，調(diào)用三階多項式（非擬合結(jié)果，僅做對比）

4 結(jié)論

針對超臨界CO2物性非線性變化的特點，本文探究了一種顯式擬合方法，采用基準函數(shù)結(jié)合調(diào)整函數(shù)降低非線性程度的方式，獲得了可適用于大工況范圍內(nèi)的擬合關(guān)聯(lián)式，該關(guān)聯(lián)式非線性系數(shù)少、精度高、計算成本小，本方法得到的關(guān)聯(lián)式平均相對誤差8%，計算速度為傳統(tǒng)關(guān)聯(lián)式的20%，在超臨界CO2熱泵制熱制冷循環(huán)及發(fā)電循環(huán)運行場景中具有一定的應用和推廣潛力。