蒸気表にある物性値をPythonで記述してみる - 口下手エンジニアの悪あがき

空気中の蒸気の物性値を何かの計算で使おうとすると、蒸気表なるものを持ち出して温度や圧力のパラメータで表を読み解いて目的の比容積や飽和水蒸気量なんかを見つけないといけないんですが、どう読めばいいか思い出したり、とびとびの値だからほしい値は近似しなければならなかったり、面倒この上ない。。。こんな時代にいちいち蒸気表を読んでられません、というわけで、Pythonにあるライブラリでそういうものがあるのか調べてみました。

iapws.iapws97モジュールとは
計算してみる
任意の温度における蒸気量

iapws.iapws97モジュールとは

ありました。「iapws」というライブラリの中の「iapws.iapws97」というモジュール。

もともとは「国際水・蒸気性質協会（The International Association for the Properties of Water and Steam）」で、水の状態を表す状態方程式がさまざまな研究者によって提唱されている。

1997年に制定されたWagnerらの手によって「実用国際蒸気状態式（IAPWS Industrial Formulation 1997（略称：IAPWS-IF97）」で定められた蒸気物性値の計算方法を利用することができ、それを活用したモジュールになります。

どのようなモジュールかDocumentationを見てみる。まず解説のグラフが描いてあり、そこに蒸気の状態図にて「IAPWS-IF97」の適用可能領域(implementation)が示してある。Region1（①）からRegion5（⑤）まで分けて与えられており、蒸気曲線としては、Region4（④）であり、「ｃ」点は「Critical point」と呼び、臨界点（647.096Ｋ/22.064MPa）を指す。

Region4（④）の範囲としては、
温度：273.15K $\le T \le$ 647.096K
圧力：0.000611MPa $\le P \le$ 22.064MPa
の範囲で計算が可能なようである。

図１．IAPWS-IF97 standard implementation

※「iapws.iapws97 module Documentation」より引用 https://iapws.readthedocs.io/en/latest/iapws.iapws97.html

計算してみる

早速使ってみる。とりあえず、先の記事であげたような空気中の水蒸気量を計算する基本的な計算を考えてみる。

Documentationでは、Region4（④）では以下の関数しか与えられていないので、「Vapor auality（ $x$ ）」を「１」として、圧力 $P$ の範囲の水準を振ってグラフを描いてみよう。

ganbaruengineer.hatenablog.com

まずはモジュールをpipでインストールしておく。

pip install iapws

次に、簡単にコードを書いていきましょう。最低限しか書いてませんが・・・

import iapws
import matplotlib.pyplot as plt
import numpy as np

# 実行可能範囲で圧力変数を定義（今回は0.001MPa～0.01MPaを100分割）
p = np.linspace(0.001,0.01,100) 
# 圧力に対応する温度を計算
t = np.array([iapws.iapws97._Region4(pi,1)["T"] for pi in p])  
# グラフが見やすいように絶対温度を摂氏温度に換算
tc = np.array([ti-273.15 for ti in t]) 
# 各圧力での比容積（m^3/kg）を計算
v = np.array([iapws.iapws97._Region4(pi,1)["v"] for pi in p])  
# 各圧力での水蒸気量（kg/m^3）を計算
w = np.array([1/(iapws.iapws97._Region4(pi,1)["v"]) for pi in p]) 

plt.plot(tc,w*1000)
plt.xlabel("Temperature(degree celsius)")
plt.ylabel("amount of water vapor(g)")
plt.legend()
plt.grid()

まー、思っていることはできたんじゃないでしょうか。

任意の温度における蒸気量

IAPWS-IF97モジュールの蒸気曲線である「Region4」を計算しようとすると

iapws.iapws97._Region4(p,x)

を用いないといけないため、温度が何度の時の質量を直接求めることができないみたい。
なので、上のコードではまず、 $P$ を決めてから対応する温度を出して、その温度に対応する $P$ から比容積や蒸気量を求めるようにしていました。

ただ、そのやり方では回りくどいので直接任意の温度から蒸気量を求めたい、となったときはどうするかを考えてみると、以下のようになります。

import iapws
import scipy.optimize as scop
import numpy as np

def f(p):
    t = 30+273.15  # 任意の温度を入力
    return iapws.iapws97._Region4(p,1)["T"] - t

p = scop.newton(f,0.001)
w = 1/iapws.iapws97._Region4(p,1)["v"]*1000

print("{:.2f}(g)".format(w))