余暇の時間で書かれた文章
この文書は余暇の時間に書かれました。現所属組織や、そこで行っている活動等を反映しているものではありません。 また、公開されている文書に基づいて書いています。
NIST REFPROP は、NISTが出している有償の気体物性値計算ソフトウェアである。 学校や会社等、金銭的に余裕がある組織、あるいはすでにREFPROPを所有している組織に属するならREFPROPを使えばよいと思う。 ただ、自宅で突然N2の内部エネルギーや定圧モル比熱を知りたくなった場合、REFPROPは使用できないだろう。 このような場合に有効と思われるのが CoolProp である。
CoolPropはその名が示す通りRefPropの自由な代替として使用されることを意図していると思われる。 公式のドキュメントも非常に充実しているが、日本語の資料はウェブ上に少ない。そこで使用法がわかってきたのでメモ書きを残すことにした。
以下のものについてインストール法を書く。
なお、この文章を書いた時点での最新バージョンは6.4.1だった。
公式ドキュメント の時点でわかりやすい。 pipを使うなら、以下のコマンドだけで終わる。
python3 -m pip install CoolProp
これも公式ドキュメントを参照してほしいが、まずページ一番上の"Installers"をクリックしてSourceForgeのプロジェクトページに移動する。 ソースコードがほしい場合には緑の"Download Latest Version"を選び、インストーラで良い場合はWindows→CoolProp_v6.4.1.0.exeを選ぶ。 リンク。
このインストーラによりMS Excelアドオンがインストールでき、自動的に有効になる。
REFPROPにもExcelアドオンがあるが、githubのページを見ると使い方がずいぶん難しそうだということに気づく。
ページの下のほうには温度密度を求めるにはDensityという関数が必要と書いてある。
では、内部エネルギーは？エンタルピーは？定圧モル比熱は？
別の関数が必要なことは明らかで、ドキュメントもウェブではすぐには見つけられなかった。(これらの英語スペルはすぐわかるだろうか。引数は？)
このgithubページの中にあればよかったのだが…
非英語話者にとっては、状態量ごとに関数名がバラバラでドキュメントもウェブで探し辛いというのは苦痛である。
一方でCoolPropは非常にシンプルであり、ExcelだろうがPythonラッパーだろうが関数名は同じ、PropsSI()・PhaseSI()を使えば良い。 入力は関数名からも明らかなように、SI単位系で与えれば良い。 以下にPythonの例を挙げる。
from CoolProp.CoolProp import PropsSI,PhaseSI
# T=298.15[K], P=101325[Pa]のときの空気の密度
PropsSI('D','T',298.15,'P',101325,'Air')
# T=298.15[K], P=101325[Pa]のときの窒素のモルエンタルピー
PropsSI('HMOLAR','T',298.15,'P',101325,'N2')
# T=100[K],P=3[MPa]のときの窒素の状態
PhaseSI("T",100,"P",3*10**6,"N2")
さらに、使用できる状態量の略称リストは公式ドキュメントに使用例とともに載っている。
まず、CoolPropの対応している気体 はREFPROPよりも少ない。このリストにない気体を計算したいなら商用ソフトウェアを買うしかない。ただ今後追加されるかもしれないし自分で追加する方法も公開されているのでパワーユーザーなら自分で文献を調べて追加することも不可能ではないだろう。
また、入力に制限のあるセットがある。たとえばエンタルピーと温度の組み合わせで計算しようとするとThis pair of inputs [HmolarT_INPUTS] is not yet supported とのメッセージが表示される。 ただ、これはExcelソルバー等を使う労力を惜しまなければ他の物性値から値を出すことはできるだろうと思う。 その手間が惜しいならばやはり商用のソフトウェアを買うべきかもしれない。
一部(あるいは全て?)の流体で凝固点以下の温度の入力を受け付けてくれない。これは一部の固体になりうる流体の計算では問題であるかもしれない。
さらに実例をあげてみる。
ジュールトムソン膨張(以下JT膨張)という物理過程がある。よく言われる例は多孔質から気体を押し出すというものである。 「逆転温度」以上ではJT膨張をすると温度が上がり、それ以下では下がるという。 また、Heの冷凍機はJT効果を使ってHeを冷やすという。(参考-kek中井氏の資料) ではどの程度温度が下がるのか？計算させてみる。
JT膨張が等エンタルピー変化であることを踏まえて、25℃のヘリウムを2MPaから0.1MPaに圧力を下げるときを考える。
from CoolProp.CoolProp import PropsSI
# 25℃, 2MPaのHeのモルエンタルピー
H_1 = PropsSI('HMOLAR','T',273.15+25,'P',2*10**6,'He')
# 0.1MPaのHe、同じモルエンタルピーを持つ温度は？
T = PropsSI('T','HMOLAR',H_1,'P',0.1*10**6,'He')
print(T)
実行すると、299.34 K、つまり26.19℃であり、微妙に温度が上がっていることがわかる。
今度は、30Kのヘリウムを2MPaから0.1MPaに圧力を下げるときを考える。
from CoolProp.CoolProp import PropsSI
# 30K, 2MPaのHeのモルエンタルピー
H_1 = PropsSI('HMOLAR','T',30,'P',5*10**6,'He')
# 0.1MPaのHe、同じモルエンタルピーを持つ温度は？
T = PropsSI('HMOLAR','H',H_1,'P',0.01*10**6,'He')
print(T)
こちらの結果は29.44Kとなり、少し温度が下がることがわかる。
500K、モルエンタルピーが14kJ/molのH2がある。このときの圧力は？ 前述の通りTとHMOLARの組は結果が出ないでエラーになってしまう。 このような場合は解析的に求める方法がある。
from CoolProp.CoolProp import PropsSI
from scipy import optimize
PropsSI("P","T",500,"HMOLAR",14000,"H2") #これはエラーになる
P=lambda x: PropsSI("HMOLAR","T",500,"P",x,"H2") -14000
ans=optimize.newton(P,0.1*10**6)
print(ans)
エンタルピーの圧力依存性は他のパラメータと比べると小さい気がするので例が適切でないかもしれない。 ともあれ上の例のようにニュートン法を用いた方法が可能である。
ここまでわかるとT-S線図やT-H線図を書くことができるようになると思う。
2021.07.06 誤字を修正