用於水和水蒸汽物性計算的Python模塊——iapws

在進行熱力循環分析、流動傳熱計算時，需獲得水和水蒸汽的物性參數。網上主流的水蒸汽物性計算程序是上海成套所的楊宇教授開發的，有Fortran、C、C#、VB等多個語言版本，還有桌面版本，被本專業學生和研究人員廣泛使用。可以說楊教授為同行的便利做了很大貢獻，本想貼一下他的個人博客的，但發現他的博客鏈接掛了。

最近打算使用混合編程，寫個Python版本的水蒸汽物性計算的接口程序，搞個在線查詢的web，然而google后發現已經有了用於水和水蒸汽物性計算的Python模塊——iapws，不得不安利給大家。其實工科領域用Python的不多，要是有人用的話希望別重復造輪子。

iapws官網：https://github.com/jjgomera/iapws/tree/master/docs

1、簡介

iapws是IAPWS標准的Python實現，包含一下幾個模塊：

• IAPWS-IF97——水蒸汽
• IAPWS-95——水蒸汽
• IAPWS-06——冰
• IAPWS-08——海水
• IAPWS-17—— 重水
• ......

iapws依賴於numpy-scipy科學計算模塊。本文主要介紹IAPWS-IF97模塊的使用。IAPWS-IF97實現了5個區域的基本方程(下圖)。可以看出壓力、溫度的范圍是很寬的，足夠滿足工程需要。

2、使用

直接在控制台執行：pip install iapws，安裝iapws時會自動安裝numpy scipy這兩個模塊。

2.1 IAPWS97類

可以使用該IAPWS97類創建一個特定熱力學狀態的對象，該類的構造函數的關鍵字參數包括：

• T：溫度[K]
• P：壓力[MPa]
• h：比焓[kJ/kg]
• s：比熵[kJ/kgK]
• x：干度[-]

有效的參數組合有：

• T, P: 對兩相無效
• P, h
• P, s
• h, s
• T, x: 僅適用於兩相
• P, x: 僅適用於兩相

計算的物性參數如下表所示。我將自己認為常用的參數列在前面了。許多參數不知道什么意思，翻譯也不知道有沒有問題。可參考：http://iapws.readthedocs.io/en/latest/iapws.iapws97.html

P：壓力[MPa]	a：亥姆霍茲自由能[kJ / kg]	joule：焦耳 - 湯姆森系數[K / MPa]
T：溫度[K]	Z：壓縮系數[ - ]	deltat：等溫節流系數[kJ / kg·MPa]
v：比容量[m³/ kg]	fi：逸度系數[ - ]	region：地區
rho：密度[kg /m³]	f：逸度[MPa]	v0：理想比容[m³/ kg]
h：比焓[kJ / kg]	γ：等熵指數[ - ]	u0：理想的內能[kJ / kg]
u：特定內能[kJ / kg]	alfav：等壓膨脹系數[1 / K]	h0：理想比焓[kJ / kg]
s：比熵[kJ / kg·K]	xkappa：等溫壓縮率[1 / MPa]	s0：理想比熵[kJ / kg·K]
cp：定壓比熱[kJ / kg·K]	kappas：絕熱可壓縮率[1 / MPa]	a0：理想亥姆霍茲自由能[kJ / kg]
cv：定容比熱[kJ / kg·K]	alfap：相對壓力系數[1 / K]	g0：理想比吉布斯自由能[kJ / kg]
g：比Gibbs自由能[kJ / kg]	betap：等溫應力系數[kg /m³]	cp0：理想定壓比熱[kJ / kg·K]
n：折射率[ - ]	Pr：折算壓力[ - ]	cv0：理想定容比熱[kJ / kg·K]
Prandt：普朗特數[ - ]	Tr：折算溫度[ - ]	Svap：蒸發熵[kJ / kg·K]
μ：動態粘度[Pa·s]	w0：理想音速[m / s]	gamma0：理想等熵指數[ - ]
nu：運動粘度[m²/ s]	k：導熱系數[W / m·K]	epsilon：介電常數[ - ]
w：音速[m / s]	alfa：熱擴散系數[m²/ s]
Hvap：汽化熱[kJ / kg]	sigma：表面張力[N / m]

查看源碼發現，使用IAPWS97類進行物性查詢時，先執行calculable方法判斷輸入條件是否可計算，若可以計算則執行calculo方法，判斷輸入參數確定的物性狀態所處的區域，隨后計算各物性。下面是示例：

>>> from iapws import IAPWS97
#>>>> 常壓常溫水 <<<<
>>> water=IAPWS97(T=24+273.15,P=0.013)
# 焓值
>>> water.h
100.66509664191254
# 密度
>>> water.rho
997.2595184928771
#>>>> 高壓蒸汽 <<<<
>>> vapor = IAPWS97(P=15.5, x=1.0)
# 溫度
>>> vapor.T
617.9415516035506
# 飽和汽焓
>>> vapor.h
2596.2167214338015
# 飽和水焓
>>> IAPWS97(P=15.5, x=0).h
1629.8502994294881

當輸入無效的參數組合時，實例會正常產生但不進行物性計算，實例的calculable方法返回空字符串；當輸入有效的參數組合時，若可查詢則calculable方法返回查詢組合，若超出查詢范圍則拋出異常。

#>>>> 無效的輸入 <<<<
>>> test1 = IAPWS97(h=2000.0, T=300.0)
# 未綁定屬性h
>>> test1.h
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'IAPWS97' object has no attribute 'h'
# calculable返回空字符串
>>> test1.calculable
''
#>>>> 有效的輸入 <<<<
>>> vapor = IAPWS97(P=15.5, x=1.0)
# calculable返回輸入的組合
>>> vapor.calculable
'Px'
#>>>> 超出查詢范圍的輸入 <<<<
>>> test2 = IAPWS97(P=100.0, T=45.0)
Traceback (most recent call last):
...
NotImplementedError: Incoming out of bound

IAPWS97類提供了豐富的物性計算，若僅是為了查詢物性，直接使用該類就行了。但在瞬態計算或其他需要頻繁計算物性的場合下，通過IAPWS97類計算物性是不夠明智的。因為一般情況下僅僅需要輸出一個物性參數就足夠了，如根據壓力、焓值計算溫度，而IAPWS97類計算了許多我們不需要的參數，這會拖慢程序的運行效率。本來程序具有超實時的運行能力，由於物性查詢方法的錯誤使用，程序可能變得巨慢。物性查詢成為程序的性能瓶頸，這顯然是不合理的。

iapws.iapws97模塊提供了一系列的基本方程（fundamental equations）和向后方程（backward equation），用於特定場合下的物性計算。

2.2 基本方程

iapws.iapws97模塊提供了如下基本方程，用於不同區域內的物性計算。形如_Regionr()的函數，根據輸入的參數，確定區域r內的狀態，返回一個字典；_TSat_P、_PSat_T這兩個函數計算飽和線上的物性。

• _Region1(T, P)
• _Region2(T, P)
• _Region3(rho, T)
• _Region4(P, x)
• _Region5(s)
• _TSat_P(P)
• _PSat_T(s)

_Region1和_Region2的輸入均是T、P，那么在使用時要預先判斷T、P所指示的區域：

from iapws import iapws97
# >>>> 區域1計算示例 <<<<
# 確定區域
>>> iapws97._Bound_TP(300,3)
1
# 獲得比容
>>> iapws97._Region1(300,3)['v']
0.0010021516796866943
# 獲得比焓
>>> iapws97._Region1(300,3)['h']
115.3312730214384
# >>>> 區域2計算示例 <<<<
# 確定區域
>>> iapws97._Bound_TP(700,3)
2
# 獲得比容
>>> iapws97._Region1(300,3)['v']
0.0010021516796866943
# 獲得比熵
>>> iapws97._Region1(300,3)['s']
0.39229479240262427

飽和線上的物性計算：

# 計算飽和溫度
>>> iapws97._TSat_P(15.5)
617.9415516035506
# 計算飽和壓力
>>> iapws97._PSat_T(100+273.15)
0.10141797792131013
>>>

2.3 向后方程

向后方程指形如_Backwardr_x_yz的方程，其中r為物性區域，yz為輸入參數，x為返回參數。

• _Backward1_T_Ph
• _Backward1_T_Ps
• _Backward1_P_hs
• _Backward2_T_Ph
• ...

需要注意，和基本方程一樣，向后方程也不進行區域判斷。當選擇錯誤的函數時，將輸出離奇的計算結果：

>>> from iapws import iapws97
# 錯誤的使用示例
>>> iapws97._Backward2_T_Ph(3,500)
4.1313215739117547e+21
# 錯誤的使用示例
>>> iapws97._Backward3_T_Ph(3,500)
-1637746.3600011615
# 判斷區域
>>> iapws97._Bound_Ph(3,500)
1
# 正確的使用示例
>>> iapws97._Backward1_T_Ph(3,500)
391.7985087624256