水和水蒸气物性参数计算程序开发

蓝文鸿， 于新国

(华北电力大学 核科学与工程学院，北京 102206)

符号说明：

ρ——密度，kg/m3

h——比焓，kJ/kg

T——热力学温度，K

w——当地声速，m/s

p——压强，MPa

s——比熵，kJ/(kg·K)

u——热力学能，kJ/kg

v——比体积，m3/kg

cp——比定压热容，kJ/(kg·K)

cv——比定容热容，kJ/(kg·K)

精确快速计算出水和水蒸气的物性参数对核电站的设计和热工水力分析程序开发至关重要。许多高校和研究所均有自己独立开发的水和水蒸气物性参数计算程序[1]，这些程序大多基于C++和Fortran语言编写[2-3]，能对少量数据进行精确计算，精度要求满足目前工业应用需要。但现有的IAPWS-IF97标准对水和水蒸气物性参数的计算主要侧重于优化迭代速度，降低迭代次数[4]，较少涉及到因程序截断误差引起的计算结果不准确的分析。同时现有程序也没有对计算结果进行校核及进一步的检验，导致由于精度要求的提高，其截断误差随着迭代次数的增加而积累，从而引起结果不准确的问题。此外，目前的程序大多采用手动少量数据输入进行运算的模式，不具备对大规模的数据进行批量导入、计算和导出的功能。

过去基于工程语言编写的程序存在计算精度和运算速度不足的问题，而Matlab语言简单易懂，能够以图形化界面直观展示大规模计算结果的数据变化规律，同时其对截断误差的处理优于C++和Fortran。在程序结构上，笔者设计了更合理的计算步骤和程序逻辑，基于导出式分析解与迭代式数值解相互校核检验结果的思路提出反馈校核函数，保证程序在不同精度要求下仍能稳定快速输出计算结果。在数据输入、输出上设计了相关接口，能快速输出多种数据格式，为大规模计算提供可能。水和水蒸气物性参数计算程序基于IAPWS-IF97标准[5]进行编写，区域划分采用IAPWS-IF97标准推荐划分方式，其计算参数覆盖范围为0≤p≤100 MPa、273.15 K≤T≤2 273.15 K。

1 理论推导及分析

根据IAPWS-IF97标准，结合文献[6]给定的单位质量工质的亥姆霍兹自由能f=u-Ts和单位质量工质的吉布斯自由能g=h-Ts，确定水和水蒸气状态，依据水和水蒸气的状态将参数待计算范围划分为5个区域[5](见图1)，分别对应：常规水区、过热蒸汽区、临界区、饱和水线和高温低压区。在临界区，依据临界水的状态划分了26个二级区域[7-8]进行计算。采用吉布斯自由能函数g(p,T)和亥姆霍兹自由能函数f(ρ,T)为基本方程，具体公式如下:

f(ρ,T)=u(ρ,T)-T×s(ρ,T)

(1)

g(p,T)=h(p,T)-T×s(p,T)

(2)

通过g(p,T)和f(ρ,T)的导出式即可计算基本的水和水蒸气物性参数。

图1 区域和对应方程Fig.1 Division of regions and corresponding equations of IAPWS-IF97

对于简单可压缩的纯物质系统，任意一个状态参数都可由另外2个相互独立的状态参数的特征函数表示。通过对特征函数应用全微分得到麦克斯韦关系式：

(3)

根据麦克斯韦关系式和拉普拉斯方程推导得出各物性参数，结果见表1和表2。

2 程序设计与功能

2.1 程序设计

采用偏导导出式计算时，需要逐次迭代解高阶偏微分方程，虽然计算结果精确，但计算时间过长，不利于程序的快速响应。因此，本文程序没有采用偏导导出式的解法，而是采用自拟合的迭代式的数值解法。对IAPWS-IF97标准中的拟合迭代式进行修正，得出本程序各区域的最终迭代式。

表1 g(p,T)推导得出的物性参数Tab.1 Physical properties obtained by partial derivative of g(p,T)

本文程序的主要工作原理是基于迭代式的运算

表2 f(ρ,T)推导得出的物性参数Tab.2 Physical properties obtained by partial derivative of f(ρ,T)

求得水和水蒸气物性参数。主程序的流程图见图2。

图2 程序流程图Fig.2 Program flow chart

程序输入的独立参数为(p,T)，当程序读入输入参数(p,T)时，会先利用区域判别算法对参数所在区域进行判别,再将待求参数输入到对应子区域进行运算。当程序在Matlab软件下运行时，激活反馈校核算法，对迭代式计算结果与偏导导出式计算结果进行比较后输出计算结果。当需要快速响应或不在Matlab软件下运行时，将不会使用反馈校核算法。

反馈校核算法流程图如图3所示，主要用于对迭代式计算结果进行校核检验。反馈校核算法首先根据区域判别算法给定的区域选择相应的吉布斯自由能函数g(p,T)或亥姆霍兹自由能函数f(ρ,T)。进而对g(p,T)或f(ρ,T)进行偏导导出式运算，求得各参数的近似分析解。通过反馈校核算法对偏导导出式计算结果与迭代式计算结果进行比较，当两计算结果差值小于1%时，判定结果可信，输出主程序结果；当其计算结果差值大于等于1%时，反馈校核算法将要求主程序输入参数精度提高三位，再次运行主程序和反馈校核算法，直到计算结果差值小于1%为止。

图3 反馈校核算法流程图Fig.3 Feedback calculation chart

由于IAPWS-IF97在区域3的计算十分繁琐，区域划分多，迭代次数高，程序中并没有沿用其他常用的遍历区域判别算法和多次迭代的参数计算方法。笔者采用数学建模中常用的求最优解的模拟退火算法[9]中的区域判别算法函数，其算法流程图见图4。

图4 区域判别算法流程图Fig.4 Flow chart of regional discriminating algorithm

对区域判别算法函数在区域3的应用进行特殊处理。先通过转换函数将输入参数(ρ,T)转化为(p,T)，将(p,T)视作坐标点，即p为横坐标，T为纵坐标，再计算坐标点(p,T)与参考点(0，273.15)的距离L1，比较参考点到各区域中心点的距离L2与L1的大小，当L2与L1接近时(其差值小于区域长度的一半)，则认为待求参数落在该区域内。将输入参数(p,T)与匹配区域的边界进行比较、校核，满足则确认待求参数落在该区域内，若不满足则将区域移至相邻区域再进行判别，直至找到满足区域或判断其超出计算边界为止。

此外，区域判别算法在待计算参数的全局区域也能使用，经过多组数据的测试，其在全局区域判定中无需遍历5个区域即可准确判断出待求参数所在区域，避免了IAPWS-IF97标准中提供的全区遍历算法须全区遍历的缺点。

为了适应不同输入参数的计算需求，本文程序还附带了不同独立参数之间的转换模块，各模块的导出式也可由麦克斯韦关系式得出，因此能够灵活适配各种输入情况。各模块如表3所示。

表3 输入参数的转换模块Tab.3 Converting module for input variables

2.2 程序主要实现功能

编程语言采用科学计算语言Matlab。根据核电站的实际情况，程序的计算精度设定为64位，其显示精度设定为10位。该程序主要实现以下功能：

(1)计算工质的饱和状态参数，即根据工质的饱和温度求取饱和压力，或根据饱和压力求取饱和温度。此功能采用IAPWS-IF97标准提供的饱和压力psat(T)和饱和温度Tsat(p)函数实现。

(2)根据工质的温度和压力求取工质的物性参数。基于IAPWS-IF97标准提供的g(p,T)函数与自拟合的迭代式计算出区域1、区域2和区域5工质的物性参数。对于区域3的求解，采用国际水和水蒸气物性协会(IAPWS)发布的IF97-S05论文[7]提供的补充方程v(p,T)先行求解工质密度ρ，再代入f(ρ,T)进行计算。

(3)程序根据实际要求，已经预留其与其他相关程序连接的接口。目前，本文程序已经能够实现基于Aduino单片机平台进行水和水蒸气物性参数的计算，即本文程序能够作为其他更大规模的分析计算软件或者嵌入式系统，提供水和水蒸气物性参数的计算数据。

3 程序验证

3.1 程序计算精度验证

选取IAPWS-IF97标准提供的推荐检验数值进行计算，本文程序计算结果与IAPWS-IF97标准对比数据见表4。由表4可知，通过对推荐检验数值的计算可以发现，所有计算结果的误差均低于0.1%，这说明本文程序达到了IAPWS-IF97标准的精度要求。

另外，通过对比本文程序求取的300 K时不同压力下水和水蒸气比体积的结果与2008年版国际水和水蒸气表[10]中相应条件下的比体积检验数据，来验证所开发的程序。以上述数据为例进行误差分析，结果如图5所示，计算结果与标准值之间相对误差低于文献[11]给出的工程标准10-6%。

表4 计算结果与标准值的对比Tab.4 Comparison of calculation results with standard values

图5 300 K时比体积标准值与计算结果的相对误差Fig.5 Relative error of calculated specific volume over standard value at 300 K

根据IAPWS-IF97标准规定，满足工业及科学计算的水和水蒸气物性参数计算结果与标准值相对误差应小于0.96%。根据文献[11]，中广核大亚湾岭澳核电站要求水和水蒸气物性参数计算相对误差不高于0.96%，一回路关键物性参数(温度T、压力p、比体积v和比定压热容cp)计算相对误差不高于0.81%。与2008年版国际水和水蒸气表对比，本文程序计算结果相对误差小于0.96%，关键物性参数计算结果的相对误差小于0.7%。满足一般工业及压水堆核电厂的精度要求。

3.2 程序性能验证

所开发的水和水蒸气物性参数计算程序采用了反馈校核算法，在区域判别上借鉴了模拟退火算法，具有优秀的计算精度以及速度、高度的可移植性和拓展性。

通过与一些水和水蒸气物性参数计算开源程序进行对比来验证本文程序的运算速度。选用程序为哥德堡大学开发的XSteam水和水蒸气物性参数计算程序和加州大学伯克利分校开发的IAPWS-IF97计算程序。在计算精度一致的条件下对35 000组无规律的水和水蒸气(p,T)数据进行计算，对3个程序进行储存所占空间、批量计算时占用系统资源和运算速度进行对比，结果见表5。

表5 3个程序性能指标对比Tab.5 Performance comparison for three programs

本次对比测试采用控制变量法，结合Matlab软件自有的在线监测系统记录其运行状态从而得出相关对比结果。

从表5可以看出，本文程序在计算性能上优于XSteam和IAPWS-IF97计算程序。本文程序储存所占空间比XSteam计算程序小了47.02%，比IAPWS-IF97计算程序小了11.52%；本文程序的运算时长为19.57 s，比XSteam计算程序快43.01%，比IAPWS-IF97计算程序快26.29%；对于计算时占用系统资源，本文程序峰值为18.65%，低于XSteam计算程序的23.41%，更低于IAPWS-IF97计算程序的42.81%。

在程序结构上，本文程序采用了积木式的模块化结构。通过适应所测参数的需要，采用已编写好的专门的参数转换模块，对程序的输入参数进行适当调整，克服了Xsteam计算程序在非常见参数输入下结果偏差较大的问题。因此本文程序能适应不同场合可测参数的变化。

4 结 论

(1)基于科学计算语言Matlab进行数值计算时精度要高于C++和Fortran，同时能够更快、更精确地求取水和水蒸气物性参数，程序求取结果对比2008年版国际水和水蒸气表，其相对误差最大不超过0.96%，符合IAPWS-IF97标准要求。求取精度为64位精度，达到轻水堆核电站热工计算的应用要求。

(2)目前基于IAPWS-IF97标准的物性参数计算软件，其大量机时被消耗在参数的区域判定上。本文程序借鉴模拟退火算法提出了一种区域判别算法，能实现更快、更高效的输入响应，为实现大量数据计算和在线监测的即时反馈提供可能。

(3)在核心计算模块上采用了尽量简单和通用的函数，故程序易于移植，能够在不同平台和设备上运行。同时其所占空间很小，能够移植至嵌入式系统进行运算。

(4)本文程序首次引入反馈校验算法，通过比对分析解与迭代式数值解的方式，相互校核计算结果，能够有效保证计算精度，进一步确保计算结果的可靠性。