http://www.paper.edu.cn
水和水蒸气热物理性质的高精度计算
王 黎,郭吉敏
华中科技大学能源与动力工程学院,武汉(430074)
E-mail:caspernsync@sina.com
摘 要:本文基于 IAPWS-IF97 公式计算模型,运用 C 语言编制了水和水蒸气热力学性质的通用计
算软件,以实现水和水蒸气热物理性质的高精度计算。研究表明,IAPWS-IF97 作为一项国际标准
具有计算速度快,精度高等优点。软件的基本功能是计算水和水蒸气的压力、温度、比容、比熵、
比焓、干度并判断水和水蒸气的性质和状态。此外,该软件还兼顾了边界区域的连续性,迭代间
隔和精确度的一致性,如何平衡计算速度与精度的矛盾和提高程序可移植性等关键问题。
关键词:水和水蒸气; 热物理性质; 计算软件; 1997年工业用水和水蒸气性质计算模型
1.引 言
水和水蒸气作为热能动力工程中最常用的一种工质,已被广泛应用于工业领域,了解水和水
蒸气的热力学性质已成为工业生产和科学研究的需要,如何快速而准确地计算它的热力性质十分
重要。为此国际公式化委员会于1997年提出工业用新型水和水蒸气热力学性质计算模型(简称
IAPWS-IF97)。目前国内对IAPWS-IF97的文献报道都偏重介绍其特点和数学模型,而没有全面介
绍关于其数学模型的计算机应用程序[1][2],利用IF97公式制成水和水蒸气热力学性质图表,存在查
图的不准确性和表格的不连续性。因此我将编制一个计算机应用程序用于工业设计和科研工作中。
符号说明见表1。
p 压力,MPa
T
绝对温度,K
h
比焓,kJ/kg
s
比熵,kJ/(kg·K)
v
比容,m3/kg
g
无量纲吉布斯自由能
f
无量纲赫尔姆霍茨自由能
t
温度升高比
d
比密度,kg/m3
上标:
o
理想气体部分
r
剩余部分
'
饱和液态状态
"
饱和气态状态
表 1 符号说明
常量
x
R
气体常数,kJ/(kg·K)
M
摩尔质量,kg/mol
mR
摩尔气体常数,J/(mol·K)
g
吉布斯自由能,J
f
赫尔姆霍茨自由能,J
q
温度降低比
p
压力降低比
r
质量密度
下标:
c
临界点
max 最大量
s
t
饱和状态
三相点
- 1 -
2.IAPWS-IF97 的计算模型
http://www.paper.edu.cn
图1 IAPWS-IF97的区域划分和计算模型
IAPWS-IF97的计算区域被划分为五个区域,如图1所示。它们分别是:1-过冷水区,2-过热蒸
气区,3-临界水和临界蒸气区,4-饱和区,5-高温(过热)蒸气区。其中2区可细分为2a,2b,2c三个
子区域,三区可细分为3a,3b两个子区域[3]。
IAPWS-IF97 的计算公式模型:区域 2 和区域 3 间的边界是由 B23 曲线方程式定义。区域 1
gp T 的基本方程式来分别描述,区域 3 则可由特定的赫尔姆霍茨
sp T 来描述。
)
Tr 的基本方程式来描述,区域 4 也就是饱和曲线是由饱和压力方程式 (
,
gp T 方程式。在图 1 中还显示了这五个基本方程式的逆向方程式。
)
)
和 2 可由特定的吉布斯自由能 (
自由能 (
f
高温区 5 亦适用 (
)
,
,
IAPWS-IF97 公式在适用范围、区域划分和计算模型方面对 IFC67 公式进行了极大的改善和
优化,更适合工业实际应用和编程计算[4]。IAPWS-IF97 公式还给出了水的动力粘度、导热系数、
普朗特数、表面张力、静介电常数、折射率等新的计算公式和表格。对 IAPWS-IF97 和 IFC-67 进
行测试,发现在相同条件下,IAPWS-IF97 的计算结果的精度要比 IFC-67 高出一个数量级[5]。
实验表明:在最为重要的1区、2区和4区,IAPWS-IF97公式计算速度平均比IFC67公式快5.1倍。
在3区IAPWS-IF97公式计算速度比IFC-67公式快3.6倍。IFC-67公式不包括5区,但是沿5区的低温
边界区,它接近于IFC-67的最高温度区,IAPWS-IF97公式计算速度比IFC- 67公式快12.2倍[6][7]。
与IFC-67相比,IAPWS-IF97另一个品质上的飞跃是在子区域边界上计算结果更符合连续性要
求[8]。内部一致性的IAPWS-IF97公式在分区边界上的偏差非常小,远小于IFC-67公式的允许偏差。
例如在2区和3区边界上,IAPWS-IF97公式计算的比焓的最大可能偏差为0.0008%。
- 2 -
http://www.paper.edu.cn
3.计算软件的编制
3.1 程序基本内容及适用范围
本水和水蒸气性质的通用计算程序以水和水蒸气性质计算工业标准IAPWS-IF97公式为依据,
运用 C 语言编程,依靠 Microsoft Visual C++ 6.0开发环境来完成的。软件的输出界面如图 2。
图2 水和水蒸气性质计算软件界面输出
本程序可实现在IAPWS-IF97公式区域范围内已知压力 p ,温度T ,比容 v ,比焓 h ,比熵s
中任意两个值,求出该状态参数下所在的区域,并求出水和水蒸气热物性的其他参数值。也可已
知压力 p ,温度T 中一个值和干度 x 求出其他状态参数,还可判断出水蒸气的性质和状态。
适用范围:IAPWS-IF97把水和水蒸气分为五个区域,描绘了温度在 273.15KT1073.15K
,
压力在 P100MPa
的范围
下水和水蒸气的性质和状态[9]。本软件能够在工程上,尤其在热能动力、石油、化工、核能、航
的范围下及高温区温度在1073.15KT2273.15K
,压力在 P10MPa
天、海洋等设计、制造及运行领域得到广泛应用,还可应用于计算机辅助教学和计算机辅助设计。
3.2 程序软件的功能和特点
⑴ 全区域求解:该程序可实现根据用户输入的热力学参数自动进行区域判断和计算,并能进
行多区域求解,无须确定计算点所在的子区域,这样无论对于熟悉用户还是不熟悉用户均可适用。
⑵ 9种组合:该程序可以在水和水蒸气的整个计算范围内根据压力、温度、比焓、比熵、比
容、湿蒸气干度6个状态参量中两个已知参量求解出其它所有参量。由于已知焓和熵、已知干度和
焓等6种情况下涉及到二元迭代的问题,而且有时非线性方程的根不是唯一的,比如已知温度和焓
或干度和焓时可能会有单解、双解,已知干度和熵时可能会得到单解、双解甚至三解,处理起来
存在较大的困难,基于此目前完成的是15种组合中的9种,并提出另外6种的解决办法。
⑶ 该程序运用了IAPWS-IF97基本公式和多个补充推导公式:T(p,h) (表示由压力、焓求温度,
下类同)、T(p,s)、v(p,h)、v(p,s)等,避免了二次迭代,提高了运行速度和精确度。
- 3 -
£
£
£
£
£
£
⑷ 程序中使用多个限制条件退出不必要的运行,并及时输出结果,最大限度地减少了运
行时间。
⑸ 在该程序中,利用函数的单调性进行迭代,提高精确度。程序精确度的要求见表 2。
http://www.paper.edu.cn
单项区
精确度
饱和区
精确度
比容
±0.05%
焓
/kJkg
0.2
饱和压力
±0.05%
表 2 程序的精确度
熵
/kJkg K
0.0002
饱和温度
±0.02%
吉布斯自由能
0.2
/kJkg
吉布斯自由能
/kJkg
0.2
3.3 区域的判断
3.3.1 饱和区的判断
饱和区判断方法是在611.213Pap22.064MPa
的范围内,计算
出已知状态下饱和水和饱和蒸气的状态参数,判断已知的 h , s , v 是否在两者之间,即在饱和
区。例如,已知压力 p 、熵 s ,由 p 分别求出饱和水和饱和蒸气的 _s a 、 _s b ,若 _
_
就可判断为饱和区,可用饱和区和一、二区的公式求出四区的状态参数。
和 273.15KT647.096K
sass b< <
3.3.2 非饱和区的判断
根据边界方程 23B 和饱和线方程 (
sp T 或 (
)
)
sT p 以及等温线( 623.15
T
=
K
, 1073.15
T
=
K
)
pMPa=
以及等压线( 10
)即可判断水和水蒸气的区域。
如果已知温度T 和压力 p ,可根据温度T 求得边界上的压力 p ,与输入压力 p 相比较,即
可判断工质所在区域。如果已知温度T 与其它热力学参数中的任一个,可以由边界方程 23B 计算
出等温线上的压力 p ,然后用求来的压力 p 和输入的温度T 来计算出相应的热力学参数,通过与
输入参数进行比较,即可判定工质所在区域。例如,如果已知温度T 和熵 s ,首先,通过等温线
与温度T 的比较,得出工质所处的大概区域;其次,将温度T 代入4区的公式,得到对应的饱和
压力 p ;再次,将温度T 和压力 p 分别代入1区和2区的计算公式,计算得到对应的 's (饱和水熵)
和 "s (饱和汽熵);最后,将已知 s 和计算而得的 's 和 "s 进行比较,即可判断工质所在的区域。
该程序对区域的判断有两种基本方法:
⑴ 以等压线 p 为基准,求出在此压力的等压线与各区分界线的交点上的状态参数值,根据
两部分,
求出压力为 p 的等压线与1、3,2、3区域分界线交点上
s> 则判断在2区。如果
s< 则判断在1区,若 1323
s
ss
求出压力为 p 的等压线与饱和线的交点值,若
则判断在1区,否则判断
已知值和所求值的关系进行区域判断。程序将 p 分为 16.529
例如,已知 p 、 s ,如果 16.529
的 13s 、 23s ,若
pMPa
在2区。(饱和区前面已经判断)
s< < 则判断在3区,若
ss b<
_
和 16.529
16.529
pMPa>
pMPa>
pMPa
s
13
23
⑵ 以等温线T 为基准,求出在此温度下的等温线与各区分界线交点上的状态参数值,根据
已知值和所求值的关系判断。程序将T分为273.15623.15
823.151273.15
,1273.152273.15
KT
KT
KT
K
K
K
, 623.15823.15
KT
K
,
四个区域(饱和区前面已经判断)。判断方法
- 4 -
–
–
–
–
£
£
£
£
£
£
£
£
£
£
£
£
£
£
和等压线判断相同。具体的区域判断的方法可由程序图中具体的表现出来。
3.4 计算原理和方法
http://www.paper.edu.cn
)
IAPWS-IF97公式将有效范围分成了5个分区。在1区、2区和5区内使用吉布斯自由能函数
)
gp T ,在3区内使用赫尔姆霍茨自由能 (
(
f
Tr ,在4区内使用饱和压力公式。2区和3区的边界
,
线由公式给出。IAPWS -IF97公式的分区及公式水和水蒸气的性质参数不是互相独立的。如果取压
力 p 和温度T 为自变量,那么比容 v 、比焓 h 和比熵 s 可以利用下面的正则函数的偏微分得到:
(
式中 g ——吉布斯自由能函数;如果取密度 r 和温度T 为自变量,那么压力 p 、比焓
ggp T=
式中 f ——赫尔姆霍茨自由能。
h 和比熵 s 可以利用下面的正则函数的偏微分得到:
Tp h 、 (
),
除基本公式,1区、2区和4区还给出了导出方程,即用于1区和2区的方程 (
Tp s 以
),
hp s 就可
,
求得。
Tp h 和 (
sT p 。应用导出方程可减少迭代求解的复杂性。如:求 (
(
(
hpTp s
,
及用于4区的方程 (
以使用导出方程而不需要迭代求解。对 (
hp s 等应用,可通过关联式
),
(
Tr
,
)
ff
=
)
)
,
)
,
)
,
)
,
)
3.5 公式的迭代
3.5.1 迭代方法
计算模型的通用性指对于水和水蒸气的热物性性质参数,可以由任何2个相互独立的状态参
数,确定其它的热力学参数。为了提高计算模型的通用性,从而设计出适用面更广的计算机程序,
我们要有能以多组相互独立的参数作为自变量来确定其他的状态参数的功能,用迭代的方法来实
现这一功能。因此,迭代逻辑和算法就成了计算模型通用性的关键。
迭代求解问题可大概地分为两类:如果自变量组合为已知温度T 与焓h 、熵 s 、比容v 中的
任一个参数,则其它热物性性质参数的计算为隐式方程。我们称这类至少知道一个基本参数的自
变量组合形式为第一类问题。如果自变量组合形式为已知焓 h 、熵 s 、比容 v 或压力 p (在3区)中
的任2个参数,则其它热力学参数的计算也为隐式方程。我们称这类不含有基本参数的自变量组合
形式为第二类问题。
第一类问题的求解:如果自变量组合形式为已知温度T 与焓h 、熵 s 、比容v 中的任一个参
数,则可根据相应的基本公式进行一维迭代求解。例如,已知温度T 和焓 h 求压力 p 、熵 s 、比
容 v ,可由基本公式 (
计算压力 p ,比容v 和熵 s 可由压力
hp T 建立迭代方程
p 与温度T 根据相应的基本公式计算。
fhhp T=-
=
,
0
(
)
)
,
第二类问题的求解:如果已知自变量组合形式为已知焓 h 、熵 s 、比容 v 或压力 p (在3区)
中的任2个参数,则其它热力学参数的计算需要二维迭代。这类问题中最常见的是已知焓h 、熵s
组合,因为在第一类问题中已经解决了已知温度T 、焓 h 计算熵 s 的问题,则可以构造由焓 h 、
熵 s 迭代温度T ,再根据已知的迭代方法温度T和焓h 计算其它的热力学参数。
3.5.2 利用函数的单调性简化迭代过程
对于水和水蒸气的性质,具有明显的函数单调性,在迭代时如果能够充分地利用水和水蒸气
的单调性,可以简化迭代过程,提高迭代的准确度。例如,温度一定,随着压力 p 的增加,比容
- 5 -
) 0
(
vv p-
)时即可得到结果,避免了对
v 会减少,若利用到迭代过程中,把 p 从小到大(或从大到小)地迭代,当
) 0
= 判断时精确度的取值,结果更准确。对
(
于水和水蒸气中的一些区域性复杂的单调函数,例如 p 、h 的关系,当温度一定时,水的h 是随 p
的增大而增加的,水蒸气的h 是随 p 的增大而减少的,利用这个性质在迭代时可以简化迭代的过
程,同时提高了迭代的准确度。本程序正是利用水和水蒸气的单调性进行迭代的。
对已知 h 、 s 求解的主要方法是:
(
vv p-
http://www.paper.edu.cn
(
vv p-
) 0
⑴ 假设法:假设其所在区域已知,求出结果,再验证计算结果的是否在该区域。(这种方法
计算结果精确,但是计算复杂,计算速度慢。)
⑵ 拟合曲线:由于曲线方程准确度不可靠,程序中拟合的曲线只用来作为辅助判断,以便简
化计算过程,加快计算速度,这种方法误差较大,许多地方要做必要的修正。
⑶ 利用水和水蒸气热力学性质的国际协会在 2001 年 9 月补充发表的关于 1997 年 IAPWS 对
Th s )可计算 1
ph s 和 (
于水和水蒸气热力学性质陈述中压力、焓和熵函数的反推式[10](即 (
区和 2 区的压力 p 和温度T 。从而利用基本方程求出其他参数,如比容 v 等。
),
),
发展反推公式 (
),
ph s 的一个非常重要的目的是:减少从给定的 h 、s 值得到 p 和T 的计算时
),
间。在 IAPWS-IF97 中,时间要消耗在反复的过程计算中。例如,牛顿二元迭代法。用公式 (
ph s
和 IAPWS-IF97 反推式 (
Tp h
联合计算,速度是牛顿二元迭代法的 20~30 倍。
)
97
,
3.5.3 精确度的一致性
为了符合水和水蒸气的精确度的要求,在迭代过程中注意了精确度的一致性要求。迭代时 p
的迭代间隔为0.00001,即每相间0.00001取一个 p ,v 的迭代间隔为0.000005,T 的迭代间隔为0.05。
迭代间隔由IAPWS-IF97公式的精确度要求确定的。
3.6 计算速度与程序精确度的矛盾及解决方案
由于IAPWS-IF97公式精确度要求较高,在进行数据迭代时迭代的间隔很小,这就引起了迭
代的数据的数量很多,有的一个循环要迭代上万个甚至更多的数据,造成计算速度慢。
解决方案:可以先“粗”迭代再“细”迭代,迭代时数据的间隔先取得较大,判断出再要计算值
的取值范围。例如,用压力进行迭代时(必须是该区域单调函数),数据的迭代间隔先取0.1,若
从小到大取值,当迭代到某值 p 时后再从 1p - 开始迭代到 p ,把间隔先取为0.00001,再次满足
条件要求时就计算出了要求的压力。我们计算一下,若一个循环, p 从16.529~100取值迭代,若
都取到则要取80多万个,然而,用这种方法取值的个数为835加10000,即一万多个,明显减少了
迭代数目。如果比一万还多,还可用同样的原理进行三步迭代,直到满意为止。
3.7 边界区域的连续性
⑴ 1-3区和2-5区的边界处理。以已知P、h为例,当判断处于1区还是3区时,需计算1-3边界上
的焓值。由于IAPWS-IF97公式没有1-3边界线的计算式,所以我们只能做近似处理,即在边界上
先用1区公式计算,再用3区公式计算,将两结果与标准水和水蒸气表对比取最接近的值。若两结
- 6 -
‡
£
果都不吻合则取平均值。在处理2、5区边界问题时亦是如此。
⑵ 2-3区的边界处理。以已知T、S为例,B23曲线概略地描述了一条熵值介于
5.047kJkg K-
1
1
http://www.paper.edu.cn
5.261kJkg K-
1
和
再用3区公式迭代,取两组值中与标准值最接近的一组;若两组值均不吻合则取平均值。
之间的等熵线;S在此区域亦无计算其他参数的公式,因此先用2区公式迭代,
1
3.8 计算算例
随机选取九组数据进行测试并计算误差。标准值由水和水蒸气热力性质表给出。结果见表 3。
表 3 水和水蒸气性质计算软件计算结果算例
输入
PT
区域
1
P(Mpa)
80.000000
T(K)
373.150000
V(m3/kg)
0.001008
0.001008
0%
0.008000
—
—
0.001634
0.001634
0%
0.074846
0.074843
0.004008%
0.009743
0.009743
0%
—
—
992.950000
992.893941
0.005646%
711. 688233
711.701618
0.001880%
1302.850000
1302.802594
0.003639%
584.149488
584.149488
0%
1773.150000
0.300000
—
—
873.150000
—
—
673.150000
—
—
473.150000
—
—
—
—
0.014585
0.014589
0.027418%
0.001725
0.001725
0%
0.064189
0.064189
0%
标准值
计算误差
2
PV
标准值
计算误差
3
PH
标准值
计算误差
5
PS
标准值
计算误差
4
PX
标准值
计算误差
5
TV
标准值
计算误差
2
TH
标准值
计算误差
3
TS
标准值
计算误差
4
TX
标准值
计算误差
—
—
50.000000
—
—
90.000000
—
—
8.000000
—
—
10.000000
—
—
2.740611
2.740625
0.000510%
24.330611
24.330014
0.002454%
50.481871
50.483196
0.002624%
1.554672
1.554672
0%
H(kJ/kg)
479.753824
479.753824
0%
3678.842576
3678.823166
0.000528%
2000.000000
—
—
4694.660432
4694.539892
0.002568%
2066.670034
2066.670034
0%
5951.250344
5951.261222
0.000183%
3500.000000
—
—
1872.419707
1872.497603
0.004160%
1822.227316
1822.227316
0%
S(kJ·kg-1·K-1)
1.250120
1.250120
0%
6.283193
6.283173
0.000318%
4.087355
4.087407
0.001272%
8.000000
—
—
4.488090
4.488090
0%
9.316389
9.316199
0.002039%
6.381909
6.382307
0.006236%
4.000000
—
—
4.380549
4.380549
0%
X
—
—
—
—
—
—
—
—
—
—
—
—
0.5
—
—
—
—
—
—
—
—
—
—
—
0.5
—
—
部分数据组与标准值存在一定的计算误差,究其原因,主要是由 IAPWS-IF97 公式的不连续
性和迭代算法的误差引起的。尽管我们对边界问题进行了有效处理,迭代间隔和精确度保持了一
致,但是理论上讲存在精度误差是在所难免的,并可以控制在允许的范围内。图 2 对应显示了软
件部分数据的输出结果。
- 7 -
-
-
4.程序流程图
http://www.paper.edu.cn
输入 P、T
T≤647.096
True
Ps(T)
P=Ps
False
273.15≤T≤2273.15
10≤P≤100
True
T≤623.15
True
P23(T)
P
Ps
True
g1 (P,T)
P