大工况范围内超临界区二氧化碳物性快速计算方法研究

文_彭兆睿 郑秋云 张信荣

1 北京大学工学院 2 湖南文理学院 3 北京市城市热管理工程技术研究中心

传统氟利昂类及烷烃类制冷剂因其较高的臭氧破坏潜能（ODP）和温室效应潜能（GWP）逐渐在制冷供热领域被限制使用，而天然工质二氧化碳（CO2）由于具有环保、无毒、不可燃的特性逐步受到关注。CO2临界温度低，容易在常规工况范围内跨越临界点达到超临界状态，利用超临界CO2无气液相变、低粘度、高密度及高可压缩性的特点，可将其应用在热泵制热制冷循环、热发电循环等场景中并获得较好性能，在“双碳”目标下显示出巨大潜力。

在实际蒸汽压缩循环或发电循环运行过程中，实时获取物性参数来评估并优化循环性能是降低能耗的重要步骤，其中关键在于根据测点温压及时、快速、高精度计算出物性参数。利用状态方程方法进行计算时，精度较高但计算量大，隐式方程则需要迭代进一步增大了计算成本，耗费计算时间。拟合关联式法以既有物性库为数据源，不基于理论原理而是采用统计学方法拟合显式关联式，具有简化计算、降低成本的优势，在传统制冷剂（如R12、R22等）应用中已获得有效证明。对于超临界CO2，因为其物性在临界点及拟临界线附近剧烈非线性变化的特点，增加了拟合的难度和拟合关联式的复杂性，既有研究存在拟合关联式过于复杂、非线性系数过多、分区复杂等特点，增加计算成本；同时，不同物性关联式较少，针对各个物性需要进一步发展。对此，本文提出了一种新的拟合思路，以焓值为例，可获得大工况范围内适用的非线性系数少、计算成本低、精度高的显式关联式，可为相关工业应用提供参考。

1 物性范围

本文针对超临界区域内的焓值建立物性库，数据涵盖压力7.38～22.18MPa，温度220～550K的大工况范围，为了准确描述靠近临界压力及拟临界温度附近的非线性趋势，数据变化梯度较大处所取数据点更为密集，对应数据点如图1所示。本文数据来源为美国国家标准技术研究所（NIST）开发的REFPROP。由图1可以看出，远离拟临界线（比热最大时对应的温度构成的曲线）时的焓值随温度变化呈现出线性趋势，拟临界线左侧远端和右侧远端具有不同的斜率，跨越拟临界线附近区域为两种斜率的过渡区域；越靠近临界压力，两远端斜率相差越大，对应过渡区域斜率越大，非线性程度越大，导致超临界区焓值关联式拟合难度增大。

图1 本文拟合所取的焓值数据点

2 拟合方法

本文本质上属于非线性数据拟合问题，由于拟合目标在于获得在大工况范围内适用的计算成本低的拟合关联式，要求不仅是关联式精度高，而且要关联式分区少及系数少，因此采用适当的拟合关联式形式至关重要，经过对Logistic方程等经典S形曲线方程的测试，尚无满足要求的简单方程形式。对此，本文非线性拟合方法的基本思想在于，针对强烈非线性区段数据，降低数据病态以及降低数据非线性程度，进一步优化拟合关联式非线性系数获得适宜形式。

首先将超临界压力下二氧化碳焓值数据及对应的压力、温度数据进行数域压缩，通过数据重整，将各参数数据集中在量级为1的正区间，避免因数据量级差异及量级较高造成拟合求解时的系数矩阵病态，影响拟合效果。重整结果如图2所示，其中Tc和hc分别表示临界温度和对应的临界焓值，为简明表示，后续h*=h/hc,T*=T/Tc，P*=P/Pc。

接着进行自适应回归方程优化，通过将拟合函数拆分为基准函数和调整函数的方法降低数据源非线性的影响。根据数据源趋势，P=22.18MPa下焓值-温度参数趋势接近线性变化如图2所示，因此以该压力下拟合得到的拟合函数f(T*)作为基准函数，同时定义Pb=22.18MPa为基准压力。为获得不同压力下焓值，这里采用调整函数计算不同压力下焓值与基准压力下焓值的调整量，该函数受压力和温度同时主导，为二元拟合函数g(T*,P*)，调整量汇成的数据点分布如图3所示。根据数据点趋势，初步采用类似Planck定律的关联式形式见式（1）、式（2），设置b1-b8共8个非线性系数，进一步分别对各系数进行敏感性分析和拟合效果评价以简化关联式形式，直至关联式误差满足需求。

图2 数据重整后的焓值-温度-压力参数分布及基准函数、调整函数示意图

另外，根据图3的数据趋势还可以发现，在T*＜0.95区间内，焓值受压力影响可以忽略不计，因此此段的拟合函数h*即为基准函数f(T*)；在T*＞0.95区间内，拟合函数为基准函数与调整函数叠加之和即h*=f(T*)+g(T*,P*)。

图3 用于拟合调整函数的数据点

3 拟合结果

根据上述拟合思路，首先基准函数f(T*)采用三阶多项式形式见式（2）。

拟合后系数分别为

对于调整函数g(T*,P*)，对类似Planck定律的关联式形式的各系数进行敏感性分析和拟合效果评价，简化系数后得到如式(3)拟合关联式形式。

拟合后系数分别为b1=3.8345e10,b2=2839574,b3=-7.5962e10,b4=0.1343。RMSE=0.0571，最大误差降至约12%（见图4）。结果显示，最大误差在（拟）临界温度附近。该拟合关联式在超临界区域大工况范围内将参数系数降低至8个，在除0.95＜T*＜1.1的（拟）临界温度区间外误差均满足5%的要求，显示出较高精度。

图4 拟合关联式结果（图中实线，左轴）与标准数据（图中点，左轴）对比及两者之间相对误差（图中虚线，右轴），1＜ P＊＜3，T＊＞0.95

进一步对比通过不同途径计算焓值所需的时间，图5显示了分别调用NIST，调用本文拟合关联式，调用既有文献关联式的时间及对应误差。结果表明，采用本文关联式可以有效降低计算量，计算成本为直接调用NIST的0.2%，为既有关联式的20%，且误差MRE仅为既有关联式的2倍，加之本文关联式适用范围更广，体现出在应用上更具有优势。

图5 通过不同途径计算焓值的时间对比，从左到右依次为：直接调用NIST，调用本文拟合关联式，调用既有文献关联式，调用三阶多项式（非拟合结果，仅做对比）

4 结论

针对超临界CO2物性非线性变化的特点，本文探究了一种显式拟合方法，采用基准函数结合调整函数降低非线性程度的方式，获得了可适用于大工况范围内的拟合关联式，该关联式非线性系数少、精度高、计算成本小，本方法得到的关联式平均相对误差8%，计算速度为传统关联式的20%，在超临界CO2热泵制热制冷循环及发电循环运行场景中具有一定的应用和推广潜力。