超临界水拟临界点的热膨胀系数计算研究

马栋梁 周 涛 齐 实 陈 杰 夏榜样 肖泽军1（华北电力大学 核科学与工程学院 北京 102206）

2（华北电力大学 核热工安全与标准化研究所 北京 102206）

3（非能动核能安全技术北京市重点实验室 北京 102206）

4（中国核动力研究设计院 核反应堆系统设计技术重点实验室 成都 610041）

5（中国核动力研究设计院 核反应堆热工水力技术重点实验室 成都 610041）

超临界水拟临界点的热膨胀系数计算研究

马栋梁1,2,3周 涛1,2,3齐 实1,2,3陈 杰1,2,3夏榜样4,5肖泽军4,51（华北电力大学 核科学与工程学院 北京 102206）

2（华北电力大学 核热工安全与标准化研究所 北京 102206）

3（非能动核能安全技术北京市重点实验室 北京 102206）

4（中国核动力研究设计院 核反应堆系统设计技术重点实验室 成都 610041）

5（中国核动力研究设计院 核反应堆热工水力技术重点实验室 成都 610041）

在临界点附近时，超临界水的热膨胀系数很大，随着压力的升高，拟临界点的热膨胀系数迅速下降，拟临界点的热膨胀系数是计算临界相变转换数的基础，对于分析超临界压力下的流动不稳定性非常重要。因此，计算超临界水在拟临界点的热膨胀系数，对于了解和掌握超临界水堆中能量的转换或热量传递非常重要。运用MATLAB曲线拟合工具箱，对超临界水的拟临界点的膨胀系数进行了拟合回归分析。拟合得出了超临界水的拟临界点的热膨胀系数的计算公式。该公式具有结构简单易于计算的特点，最大绝对误差为0.20 K-1，最大的相对误差为0.19%，计算精度满足工业研究与分析的要求。

超临界水，拟临界点，膨胀系数，拟合公式

在超临界状态下，液态水和汽态水是没有明显分界点的。在超临界区，用拟临界阶段状态点（拟临界点）作为区分拟液态和拟汽态的分界点。在拟临界点时，对应的温度值为拟临界点温度值，对应此时的定压比热容达到最大值。在超临界水从液态向汽态逐渐转变的过程中，其热力学的各种物性参数会有着各种较为奇特的变化趋势。所以在这些物性参数中，热膨胀系数作为一个导出变量，是一个重要的参数。该参数在超临界核反应堆的流动不稳定性及换热分析计算中具有重要的应用。Ambrosini等[1]在分析超临界压力下的流动不稳定的时候，提出了临界相变转换数(Trans-pseudocritical Number, NTPC)的概念。在临界相变转换数的计算公式中，拟临界点热膨胀系数是其中需要计算的一个重要变量。拟临界点的膨胀系数，即为在一定超临界压力下，对应在拟临界温度点时的热膨胀系数。值得一提的是，“超临界热膨胀系数”跟“拟临界点的热膨胀系数”两者概念是不同的。“超临界热膨胀系数”是在一定压力下，随着温度变化时由连续变化的一组数值组成。而“拟临界点的热膨胀系数”是“超临界热膨胀系数”中的峰值点[2]。计算超临界水在拟临界点的热膨胀系数，对于理解和掌握超临界水堆中能量的转换或热量的传递非常重要。在能源动力工程中的其他各个领域分析研究中，超临界水的拟临界点的热膨胀系数的计算亦有着重要而广泛的应用。但是在大量的能源与动力类的文献[2-7]中，超临界水的拟临界点的热膨胀系数的具体计算数值却鲜有报道。依据国际水和水蒸汽性质协会(International Association for the Properties of Water and Steam, IAPWS)提供的1997年工业用计算标准(IAPWS-IF97)[8]，利用液态水和汽态水的热膨胀系数差值计算的方法，可以直接计算得到在不同超临界压力下的液态水和汽态水的热膨胀系数。确定拟临界点温度之后，从而可以相应得到拟临界点的热膨胀系数。但是直接计算时，计算工作量往往过大，且计算过程较为繁琐。超临界水的拟临界点的热膨胀系数计算模型的提出，能够简单迅速地计算得到在各种不同压力下的拟临界点的热膨胀系数值，且精度范围满足超临界水的拟临界点的热膨胀系数分析的要求。该公式的提出，为进一步编程分析计算超临界水的相变转换数等参数具有重要意义，为便捷快速地量化分析超临界水的流动换热特性奠定了基础，有利于更好地对超临界水进行流动不稳定性分析。

1 超临界水热膨胀系数计算模型

1.1 计算范围

IAPSW-IF97公式将液态水和汽态水的整个有效区域划分为5个子区域，如图1所示。

图1 IAPWS-IF97的分区Fig.1 IAPWS-IF97 subregion.

在图1中的分区计算公式中，第三区为超临界区，为需要计算热膨胀系数的区域。其中划分第三区和第二区的边界计算公式为：

超临界区的计算公式适用范围如下：623.15K≤T≤T(P)Eq.(1)，Ps(T)≤P≤100MPa。

1.2 定义公式

根据工程热力学的定义，液态水和汽态水的热膨胀系数是在等压情况下比体积随着温度的变化率。计算公式如下：

式中：pβ为热膨胀系数，K-1；υ为比体积，m3·kg-1；T为热力学温度，K。

在超临界区域范围内，第三区给出的基本方程形式是Helmholtz比自由能方程(,T)fρ（对应的无量纲形式的Helmholtz比自由能方程为(,)φ δ τ）。在超临界区域中，(,)φ δ τ以折算密度δ和折算温度τ为独立的变量。压力的关系表达式如下：

式中：ρ为密度，kg·m-3；φδ为φ( δ ,τ)对折算密度的偏导数；基准密度ρc=322 kg·m-3。

1.3 计算公式

因为在超临界区域中给出的基本方程不是直接以P和T为独立变量的，所以，要使用隐含导数求导的方法，将P看成是相对于T独立的变量，然后将压力的关系方程两边同时对T求偏导数，即：

解上述方程可得：

将式(6)带入热膨胀系数的定义式(3)中，得到热膨胀系数的计算公式如下：

式中：δτφ为(,)φ δ τ对折算密度和折算温度的二阶偏导数；δδφ为(,)φ δ τ对折算密度的二阶偏导数。

2 超临界水热膨胀系数参数分析

2.1 不同压力的热膨胀系数

当水在超临界状态下，液态水和汽态水是没有明显分界点的。当水处于温度较低时，可以认为此时水的状态是和亚临界以下的液态水是一致的。当水处于温度较高时，可以认为此时水的状态是和亚临界以下的汽态水是一致的。但是在从液态水向汽态水过渡的过程中，此时处于一种液态水和汽态水的混合状态。在一定的压力条件下，在某一临界温度点附近时，此时液态水和汽态水的热膨胀系数达到最大值，之后迅速减小。超临界水在不同压力条件下的热膨胀系数变化趋势如图2所示。

图2 不同压力下超临界水的热膨胀系数变化趋势Fig.2 Thermal expansion coefficient of supercritical water under different pressure.

由图2可见，在压力一定的条件下，超临界水在从液态水向汽态水状态转变的过程中，热膨胀系数首先迅速升高，达到峰值之后又迅速下降。超临界水的热膨胀系数所对应的温度在拟临界温度点附近，此时的热膨胀系数为拟临界点的热膨胀系数。随着压力的升高，拟临界点的热膨胀系数的峰值逐渐降低，同时对应的拟临界温度值却在逐渐升高。

通过对式(3)进行差值计算，选择其中热膨胀系数的峰值，从而得出不同压力下的拟临界点的温度值和拟临界点的热膨胀系数的数值，如表1所示。

通过表1计算得出超临界水的拟临界点的热膨胀系数，可以直观地看出在各种压力下的拟临界点的热膨胀系数数值的大小，可以供相关研究人员对超临界压力下的水的拟临界点的热膨胀系数进行查询参考使用。

表1 不同压力对应的拟临界点的热膨胀系数Table1 Pseudo-critical thermal expansion coefficient under different pressure.

2.2 拟临界点的热膨胀系数变化趋势

在不同压力下的超临界水拟临界温度值和拟临界点的热膨胀系数的变化趋势情况如图3所示。

图3 拟临界温度和拟临界点的热膨胀系数变化趋势Fig.3 Trend of pseudo-critical temperature and pseudo-critical thermal expansion coefficient.

通过图3可以看出，随着压力升高，超临界水的拟临界温度值逐渐升高，而超临界水的拟临界点的热膨胀系数是迅速下降的，下降速度呈现为指数级别的速度下降。当下降到一定程度之后，拟临界点的热膨胀系数变化趋势趋于平缓，说明达到一定压力之后，随着温度的变化，超临界水的热膨胀变化是比较小的。

3 拟合公式的确定及验证

3.1 公式的确定

由于在表1中计算所得的各压力下的拟临界点的热膨胀系数是通过逐点温度差值计算热膨胀系数而得到的，计算工作量较大，且工作较为繁琐。所以根据表1中计算得到的拟临界点的热膨胀系数的原始数据，利用MATLAB拟合数据工具箱对超临界水的拟临界点的热膨胀系数值进行了拟合分析。通过拟合分析得到了如下两种公式结构的拟合公式，计算结果的回归情况较为理想。拟合公式分别如下：

式中：c1=0.327(0.317,0.337)；c2=-7.112(-7.342, -6.882)；d1=-44(-44.01,-43.98)；d2=484(483.7, 484.3)。其中括弧内的系数参数范围为回归系数为95%置信区间的界限参数。

3.2 公式的验证

拟合公式一和拟合公式二计算拟临界膨胀系数的计算结果及相对误差和绝对误差如表2、3所示。

为了将表2、3中的拟合公式计算结果进行直观对比，表2、3中直接计算结果和拟合公式一和拟合公式二得出的结果对比情况如图4所示。

通过图4中拟合公式中的计算结果和直接计算得出的拟临界点的热膨胀系数的对比可以看出，拟合公式一在压力较低时与实际计算得出的数值误差较小，但是压力超过23MPa之后拟合公式一得出的计算结果普遍比实际值偏高。拟合公式二的误差平方和(Sum of the Squared Errors, SSE)跟公式一的误差平方和相比较小，均方根误差(Root Mean Square Error, RMSE)也比公式一的均方根误差要小很多。拟合公式二的回归系数为1，表明拟合结果的实际趋势和实际计算结果的趋势是完全线性相关的。但是在压力低于23MPa的压力范围内，通过拟合公式二计算得出的结果跟拟合公式一相比偏高较多。

表2 拟合公式一和拟临界点的热膨胀系数的结果对比Table2 Contrast of pseudo-critical thermal expansion coefficient with the result of fitting Formula 1.

表3 拟合公式二和拟临界点的热膨胀系数的结果对比Table3 Contrast of pseudo-critical thermal expansion coefficient with the result of fitting Formula 2.

图4 直接计算结果和拟合公式一(a)、拟合公式二(b)计算结果的对比Fig.4 Contrast of directly calculation result with the result of fitting Formula 1 (a) and fitting Formula 2 (b).

3.3 误差分析

拟合公式一和拟合公式二的计算结果的绝对误差和相对误差分布情况如图5所示。通过图5(b)在公式一和公式二的相对误差分布情况也可以看出，在压力较低时公式一的相对误差较小，在压力较高时公式二的相对误差较小。所以综合考虑后确定，当超临界压力小于等于22.5MPa时，建议选用公式一；当超临界压力大于等于23.0MPa时，建议选用公式二。当超临界压力介于22.5MPa和23.0MPa之间时，则采用拟合公式一和拟合公式二组合的方式计算。具体组合过渡形式见§3.4的最终公式。组合计算形式的采用，保证了拟临界点的热膨胀系数的计算在不同压力区间变化时的计算结果具有连续光滑性。

图5 拟合公式一与拟合公式二的绝对误差(a)和相对误差(b)对比Fig.5 Absolute error contrast (a) and relative error contrast (b) between fitting Formula 1 and fitting Formula 2.

3.4 公式最终形式

经过§3.3的误差分析，最终确定的公式形式如式(13)所示：

用IAPWS-IF97标准计算拟临界点的热膨胀系数时，需要首先计算已知压力下对应的拟临界点温度，进而通过压力和温度，计算出该压力和温度下的热膨胀系数，即拟临界点的热膨胀系数。而采用新拟合的式(13)计算拟临界点的热膨胀系数时，由于拟临界点的热膨胀系数是压力的函数，所以仅知道压力，即可通过拟合公式求出该超临界压力下对应拟临界点的热膨胀系数。对比该两种不同的计算方法，运用本文公式可以加快反应堆程序的计算速率，进而减少计算时间。

4 结语

根据液态水和汽态水的热膨胀系数计算定义，直接计算得出了超临界水在各个压力下的拟临界点的热膨胀系数。然后根据直接计算得出的各个压力下的拟临界点的热膨胀系数，运用MATLAB曲线拟合工具箱，对超临界水的拟临界点的热膨胀系数进行了拟合回归分析，并得出了两组计算超临界水的拟临界点的热膨胀系数的计算公式。

1) 超临界水拟临界点的热膨胀系数，在临界点附近时达到最大值，随着压力的上升，拟临界点的热膨胀系数迅速下降，当压力超过25.0MPa之后，拟临界点的热膨胀系数的变化趋于平缓。

2) 当超临界压力低于22.5MPa时，建议选用拟合公式一计算拟临界膨胀系数；当超临界压力大于等于23.0MPa时，建议选用拟合公式二计算拟临界膨胀系数。当压力介于22.5MPa和23.0MPa之间时，采用上述两公式组合的方式进行计算。最终推荐计算公式形式如式(13)所示。最终计算公式最大的绝对误差为0.201K-1，最大的相对误差为0.19%，满足计算分析拟临界点的热膨胀系数的要求，且具有结构简单易于编程计算的特点。

1 Ambrosini W, Sharabi M. Dimensionless parameters in stability analysis of heated channels with fluids at supercritical pressures[J]. Nuclear Engineering and Design, 2008, 238: 1917-1929. DOI: 10.1016/j. nucengdes.2007.09.008

2 Pioro I L, Khartabil H F, Duffey R B. Heat transfer to supercritical fluids flowing in channels-empirical correlations (survey)[J]. Nuclear Engineering and Design, 2004, 230: 69-91. DOI: 10.1016/j.nucengdes.2003.10.010

3 李佳, 史秀敏. 基于IAPWS-IF97水和蒸汽体膨胀系数和压缩系数[J]. 燃气与热力, 2012, 4(4): A01-A04. DOI: 10.3969/j.issn.1000-4416.2012.04.001 LI Jia, SHI Xiumin. Volumetric expansion coefficient and compressibility coefficient of water and steam based on IAPWS-IF97[J]. Gas and Heat, 2012, 4(4): A01-A04. DOI: 10.3969/j.issn.1000-4416.2012.04.001

4 傅晟威, 周翀, 杨燕华, 等. 跨临界瞬态数值模拟中的高精度拟临界温度计算[J]. 上海交通大学学报, 2011, 3(3): 413-417 FU Shengwei, ZHOU Chong, YANG Yanhua, et al. Calculation of high-precision pseudo-critical temperatures for trans-critical transient analysis[J]. Journal of Shanghai Jiaotong University, 2011, 3(3): 413-417

5 李精精，周涛，段军, 等. 基于遗传神经网络方法的流动不稳定起始点研究[J]. 核动力工程, 2014, 4(2): 63-66 LI Jingjing, ZHOU Tao, DUAN Jun, et al. Study on onset of flow instability by genetic neural network[J]. Nuclear Power Engineering, 2014, 4(2): 63-66

6 陈玮玮, 方贤德, 商辉, 等. 超临界压力下竖直管内水的传热关系式研究[J]. 工程热物理学报, 2016, 1(1): 104-110 CHEN Weiwei, FANG Xiande, SHANG Hui, et al. Investigation of heat transfer correlations of water under supercritical pressure in vertical tubes[J]. Journal of Engineering Thermo Physics, 2016, 1(1): 104-110

7 Liao Y X, Lucas D, Krepper E, et al. Flashing evaporation under different pressure levels[J]. Nuclear Engineering and Design, 2013, 265: 801-813. DOI: 10.1016/j.nucengdes.2013.09.027

8 Wagner W, Kruse A. 水和蒸汽的性质[M]. 项红卫, 译.北京: 科学出版社, 2003 Wagner W, Kruse A. Properties of water and steam[M]. XIANG Hongwei, Tr. Beijing: Science Press, 2003

Calculation research on pseudo-critical thermal expansion coefficient of supercritical water

MA Dongliang1,2,3ZHOU Tao1,2,3QI Shi1,2,3CHEN Jie1,2,3XIA Bangyang4,5XIAO Zejun4,51(School of Nuclear Science and Engineering, North China Electric Power University, Beijing 102206, China)
2(Institute of Nuclear Thermal-hydraulic Safety and Standardization, North China Electric Power University, Beijing 102206, China)
3(Beijing Key Laboratory of Passive Safety Technology for Nuclear Energy, Beijing 102206, China)
4(Key Laboratory of Reactor System Design Technology, Nuclear Power Institute of China, Chengdu 610041, China)
5(Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China, Chengdu 610041, China)

Background: The supercritical water thermal expansion coefficient is very large when the reactor is approaching the critical point, but drops rapidly with the increase of pressure. The pseudo-critical thermal expansion coefficient is the basis for calculation of the critical phase transition transformation number which is of great importance in the analysis of flow instability of supercritical pressure water, therefore its calculation is very important to understand and master the conversion of energy or heat transfer in supercritical water reactor. Purpose: This study aims to derive a fitting formula for direct calculation of the pseudo-critical expansion coefficient. It is useful for thermal-hydraulic program analysis. Methods: First of all, the pseudo-critical temperature need to be calculated, then a regression analysis for the pseudo-critical thermal expansion coefficient of supercritical water is fitted by using

Supercritical water, Pseudo-critical, Thermal expansion coefficient, Fitting formula

MA Dongliang, male, born in 1982, graduated from North China Electric Power University with a master’s degree in 2008, doctoral student, focusing on nuclear thermal hydraulic and safety

ZHOU Tao, E-mail: zhoutao@ncepu.edu.cn

TL99

10.11889/j.0253-3219.2016.hjs.39.120601

中核核反应堆热工水力技术重点实验室基金(No.20130901)、核反应堆系统设计技术重点实验室基金(No.2015BJ0151)资助

马栋梁，男，1982年出生，2008年于华北电力大学获硕士学位，博士研究生，研究方向为核热工水力与安全

周涛，E-mail: zhoutao@ncepu.edu.cn

Supported by Nuclear Reactor Thermal Hydraulic Technology Key Laboratory Fund (No.20130901), Key Laboratory of Nuclear Reactor System Design Technology (No.2015BJ0151)

2016-10-12，

2016-11-10

MATLAB curve fitting toolbox. The fitting calculation formula about the pseudo-critical thermal expansion coefficient is derived from those fitting data in the form of combination formula. Results: The calculation result is continuous smoothness with maximum absolute error no more than 0.20 K-1, and the biggest relative error is 0.19%. Conclusion: The formula is simple in structure and easy to calculate with high efficiency, and the calculation precision meets the requirements of industry research and analysis.