Zenn
🦔

非理想気体の状態方程式

2023/05/04に公開
物理
化学
tech

概要

仕事で非理想気体の状態方程式を調べる機会があったので、様々な状態方程式をまとめておきます。
以下の英語での情報をもとに書いています。

https://en.wikipedia.org/wiki/Real_gas
https://en.wikipedia.org/wiki/Equation_of_state
https://en.wikipedia.org/wiki/Cubic_equations_of_state
http://www.sklogwiki.org/SklogWiki/index.php/Category:Equations_of_state

ここに書いたのはほんの一部です⋯。また、固体の状態方程式は省略しています。
残念ながら元論文は読めていません。

式の比較

Pが圧力、Vが1モルあたりの体積、Tが温度、Rが気体定数。

\begin{aligned} \text{理想気体} \rule{0pt}{5ex}&& \frac{PV}{RT} &= 1 \\ \text{Rankine} \rule{0pt}{5ex}&& \left(PV + \frac{a}{TV}\right) \frac{1}{RT} &= 1 \\ \text{van der Waals} \rule{0pt}{5ex}&& \left(P + \frac{a}{V^2}\right) \frac{V-b}{RT} &= 1 \\ \text{Clausius} \rule{0pt}{5ex}&& \left\{P+ \frac{a}{T(V+c)^2}\right\}\frac{V-b}{RT} &= 1 \\ \text{Berthelot} \rule{0pt}{5ex}&& \left(P + \frac{a}{T V^2}\right) \frac{V-b}{RT} &= 1 \\ \text{Dieterici} \rule{0pt}{5ex}&& P\:e^{a/VRT}\: \frac{V-b}{RT} &= 1 \\ \text{Wohl} \rule{0pt}{5ex}&& \left\{P + \frac{a}{TV(V-b)} - \frac{c}{T^2V^3}\right\} \frac{V}{RT} &= 1 \\ \text{Boynton-Bramley} \rule{0pt}{5ex}&& \left(P + \frac{a}{V^2}\right)\frac{V-b}{RT} &= \frac{1}{1+\frac{\psi^2}{T^2}} \\ \text{Beattie-Bridgeman} \rule{0pt}{5ex}&& \left(P + \frac{A(V)}{V^2}\right)\frac{V^2}{V+B(V)}\frac{1}{RT} &= 1 - \frac{c}{VT^3} \\ \text{Redlich–Kwong} \rule{0pt}{5ex}&& \left\{P + \frac{a}{T^{1/2}V(V+b)}\right\} \frac{V-b}{RT} &= 1 \\ \text{Redlich–Kwong-Soave} \rule{0pt}{5ex}&& \left\{P+\frac{a(T)}{T^{1/2}V(V+b)}\right\}\frac{V-b}{RT} &= 1 \\ \text{Peng–Robinson} \rule{0pt}{5ex}&&\hspace{1em} \left\{P + \frac{a(T)}{V(V+b)+b(V-b)}\right\} \frac{V-b}{RT} &= 1 \end{aligned}

van der Waals (1873)

\left(P + \frac{a}{V^2}\right) \frac{V-b}{RT} = 1

a, b を臨界圧力、臨界温度を用いて意表すと、

a = \frac{27}{64}\frac{R^2 T_c^2}{P_c}, \quad b = \frac{1}{8}\frac{R T_c}{P_c}

となる。

Berthelot (1899)

van der Waalsの状態方程式に温度依存性を加えたものをClausiusが提唱。それを簡略化した形。

\begin{aligned} \left(P + \frac{a}{T V^2}\right) \frac{V-b}{RT} = 1 \\ \rule{0pt}{6ex} a = \frac{27}{64}\frac{R^2 T_c^3}{P_c}, \quad b = \frac{1}{8}\frac{R T_c}{P_c} \end{aligned}

Dieterici (1899)

指数関数を用いて補正項をいれた方程式。Vが十分大きいとして展開するとvan der Waalsの状態方程式に一致する。

P\:e^{a/VRT}\: \frac{V-b}{RT} = 1 \\ \rule{0pt}{6ex} a = \frac{4}{e^2}\frac{R^2 T_c^2}{P_c}, \quad b = \frac{1}{e^2}\frac{R T_c}{P_c}

Redlich–Kwong (1949)

圧力補正項の温度依存性を、温度の平方根の形で入れた式。
この辺りから、定数を臨界圧力・温度で表すための計算が自分には苦しくなってくる⋯。

\left\{P + \frac{a}{T^{1/2}V(V+b)}\right\} \frac{V-b}{RT} = 1 \\ \rule{0pt}{6ex} a = \frac{1}{9 (2^{1/3}-1)}\frac{R^2 T_c^{5/2}}{P_c}, \quad b = \frac{2^{1/3}-1}{3}\frac{R T_c}{P_c}

Redlich–Kwong-Soave (1972)

Redlich–Kwongの式に偏心因子\omegaと温度に依存するパラメーター\alpha(T)をいれて補正したもの。このパラメーターの意味、導出方法は勉強中。

\left\{P + \frac{a\alpha(T)}{T^{1/2}V(V+b)}\right\} \frac{V-b}{RT} = 1 \\ \rule{0pt}{6ex} \alpha(T)= \left\{1 + \Big(0.48508+1.55171\omega-0.15613\omega^2\Big) \left(1-\sqrt{\frac{T}{T_c}}\right)\right\}^2

Peng–Robinson (1976)

定数a, bを臨界圧力・温度で表すための計算は諦めました

\left\{P + \frac{a\alpha(T)}{V(V+b)+b(V-b)}\right\} \frac{V-b}{RT} = 1 \\ \rule{0pt}{6ex} a \simeq 0.45724\: \frac{R^2 T_c^{5/2}}{P_c}, \quad b \simeq 0.07780\: \frac{R T_c}{P_c} \\ \rule{0pt}{6ex} \alpha(T) = \left\{1 + \Big(0.37464+1.54226\omega-0.26992\omega^2\Big) \left(1-\sqrt{\frac{T}{T_c}}\right)\right\}^2

Pythonで図示

これだけだとエンジニアのための新しい情報共有コミュニティの記事としてふさわしくないので
せっかくなのでPythonで図示してみましょう。温度は一定にして、H_2OのP-V図を描きます。
ただし状態方程式は一部のものに絞ります。

自分が書いたコードがこれです。
各パラメーターは臨界圧力・温度を使って計算した値を使っているので、このやり方だと臨界点以外では実験値と一致しないと思います。

import numpy as np
import matplotlib.pyplot as plt
from scipy import constants

molar_mass_air = 28.96 * 10 ** (-3)  # kg mol^{-1}
molar_mass_H2O = 18.02 * 10 ** (-3)  # kg mol^{-1}
Rg = constants.gas_constant  # J K^{-1} mol^{-1}
omega_H2O = 0.344
Tc_H2O = 647.14  # (K)
Pc_H2O = 22.060 * 10 ** 6  # (Pa)

# van der Waals
a_Waals_H2O = (27 * (Rg ** 2) / 64) * (Tc_H2O ** 2) / Pc_H2O  # a m^6 mol^{-2}
b_Waals_H2O = (Rg / 8) * Tc_H2O / Pc_H2O  # m^3 mol^{-1}

# Berthelot
a_Berth_H2O = (27 * (Rg ** 2) / 64) * (Tc_H2O ** 3) / Pc_H2O  # K Pa m^6 mol^{-2}
b_Berth_H2O = (Rg / 8) * Tc_H2O / Pc_H2O  # m^3 mol^{-1}

# Dieterici
a_Diete_H2O = (4 * (Rg ** 2) * (Tc_H2O ** 2)) / (np.exp(2) * Pc_H2O)  # Pa m^6 mol^{-2}
b_Diete_H2O = (Rg / np.exp(2)) * Tc_H2O / Pc_H2O  # mol^{-1}

# Redlich-Kwong
a_RK_H2O = (1 / (9 * (2 ** (1/3) - 1))) * (Rg ** 2) * (Tc_H2O ** (5/2)) / Pc_H2O  # K^{1/2} Pa m^6 mol^{-2}
b_RK_H2O = ((2 ** (1/3) - 1) / 3) * Rg * Tc_H2O / Pc_H2O  # m^3 mol^{-1}

# Peng-Robinson
a_PR_H2O = 0.45724 * Rg ** 2 * Tc_H2O ** 2 / Pc_H2O  # Pa m^6 mol^{-2}
b_PR_H2O = 0.07780 * Rg * Tc_H2O / Pc_H2O  # m^3 mol^{-1}


def P_ideal(Vm, T):
    return Rg * T / Vm


def P_Waals(Vm, T, a=a_Waals_H2O, b=b_Waals_H2O):
    return Rg * T / (Vm - b) - a / Vm ** 2


def P_Berth(Vm, T, a=a_Berth_H2O, b=b_Berth_H2O):
    return Rg * T / (Vm - b) - a / (T * Vm ** 2)


def P_Diete(Vm, T, a=a_Diete_H2O, b=b_Diete_H2O):
    return (Rg * T / (Vm - b)) * np.exp(- a / (Vm * Rg * T))


def P_RK(Vm, T, a=a_RK_H2O, b=b_RK_H2O):
    return Rg * T / (Vm - b) - a / (T ** (1/2) * Vm * (Vm + b))


def P_PR(Vm, T, a=a_PR_H2O, b=b_PR_H2O, omega=omega_H2O):
    def alpha(T):
        return (1 + (0.37464 + 1.54226 * omega - 0.26992 * omega ** 2) * (1 - np.sqrt(T/647.14))) ** 2
    alpha = alpha(T)
    return Rg * T / (Vm - b) - a * alpha / (Vm * (Vm + b) + b * (Vm - b))


fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
Vm = np.linspace(max(b_Waals_H2O, b_Berth_H2O, b_Diete_H2O, b_RK_H2O, b_PR_H2O) + 10 ** -9, 0.005, 10000)
ax.plot(10 ** 3 * Vm, 10 ** (-6) * P_ideal(Vm, 400), label='ideal gas')
ax.plot(10 ** 3 * Vm, 10 ** (-6) * P_Waals(Vm, 400), label='van der Waals')
ax.plot(10 ** 3 * Vm, 10 ** (-6) * P_Berth(Vm, 400), label='Berth')
ax.plot(10 ** 3 * Vm, 10 ** (-6) * P_Diete(Vm, 400), label='Diete')
ax.plot(10 ** 3 * Vm, 10 ** (-6) * P_RK(Vm, 400), label='Redlich-Kwong')
ax.plot(10 ** 3 * Vm, 10 ** (-6) * P_PR(Vm, 400), label='Peng-Robinson')
ax.legend()
ax.set_title('H$_2$O, 400 K')
ax.set_xlabel('$V_m$ ($×10^{-3}$ m$^{3}$ /mol)')
ax.set_ylabel('$P$ (MPa)')
ax.set_xlim(0, 1)
ax.set_ylim(0, 15)
ax.grid()
plt.show()


fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
Vm = np.linspace(max(b_Waals_H2O, b_Berth_H2O, b_Diete_H2O, b_RK_H2O, b_PR_H2O) + 10 ** -9, 0.005, 10000)
ax.plot(10 ** 3 * Vm, 10 ** (-6) * P_ideal(Vm, 2000), label='ideal gas')
ax.plot(10 ** 3 * Vm, 10 ** (-6) * P_Waals(Vm, 2000), label='van der Waals')
ax.plot(10 ** 3 * Vm, 10 ** (-6) * P_Berth(Vm, 2000), label='Berth')
ax.plot(10 ** 3 * Vm, 10 ** (-6) * P_Diete(Vm, 2000), label='Diete')
ax.plot(10 ** 3 * Vm, 10 ** (-6) * P_RK(Vm, 2000), label='Redlich-Kwong')
ax.plot(10 ** 3 * Vm, 10 ** (-6) * P_PR(Vm, 2000), label='Peng-Robinson')
ax.legend()
ax.set_title('H$_2$O, 2000 K')
ax.set_xlabel('$V_m$ ($×10^{-3}$ m$^{3}$ /mol)')
ax.set_ylabel('$P$ (MPa)')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(10 ** 3 * max(b_Waals_H2O, b_Berth_H2O, b_Diete_H2O, b_RK_H2O, b_PR_H2O), 10 ** -1)
ax.set_ylim(100, 10000)
ax.grid()
plt.show()

で結果がこちら。

低温低圧領域ではこうなりました。どれも定性的には同じく　\left(\frac{\partial P}{\partial V}\right)_T>0　となる領域があります。


高温高圧領域ではこうなりました。
Redlich-Kwong式とPeng-Robinson式は理想気体の式とvan der Waalsの式の中間の値になっています。

雑感

他にも数多の状態方程式があり、数えきれません。
有名なものだとBenedict–Webb–Rubinの状態方程式というのがあるけれど、パラメーターが8個もあるわりには劇的に精度が上がるわけでもないようなので省略しました。
現代でも色々な方々が状態方程式の改良を続けていますが、パラメーターの数が多すぎては現実的には役に立たないので、結局Redlich–KwongやPeng–Robinsonの状態方程式が使われているようです。

他にはどんな状態方程式があるのか、気になります。
詳しい方がいたら教えてもらえると嬉しいです。

参考文献

  1. William John Macquorn Rankine, "On the Geometrical Representation of the Expansive Action of Heat, and the Theory of Thermo-Dynamic Engines", Philosophical Transactions of the Royal Society of London 144 pp. 115-175 (1854).
  2. J. D. van der Waals, "Over de Continuiteit van den Gasen Vloeistoftoestand", doctoral thesis, Leiden, A,W, Sijthoff (1873).
  3. R. Clausius, "Über das Verhalten der Kohlensäure in Bezug auf Druck, Volumen und Temperatur", Annalen der Physik und Chemie 9 pp. 337-357 (1880).
  4. D. J. Berthelot, "Sur Une Méthode Purement Physique Pour La Détermination des Poids Moléculaires des Gaz et des Poids Atomiques de Leurs Éléments", J. Phys., 8 pp. 263-274 (1899).
  5. C. Dieterici, Ann. Phys. Chem. Wiedemanns Ann. 69, 685 (1899).
  6. A. Wohl, "Investigation of the condition equation", Zeitschrift für Physikalische Chemie (Leipzig) 87 pp. 1-39 (1914).
  7. W. P. Boynton and Arthur Bramley, "A Modification of Van Der Waals' Equation", Physical Review 20 pp. 46-50 (1922).
  8. James A. Beattie and Oscar C. Bridgeman, "A New Equation Of State For Flueds. I. Application To Gaseous Ethyl Ether And Carbon Dioxide", Journal of the American Chemical Society 49 pp. 1665 - 1667 (1927).
  9. Otto Redlich and J. N. S. Kwong, "On the Thermodynamics of Solutions. V. An Equation of State. Fugacities of Gaseous Solutions", Chemical Reviews 44 pp. 233-244 (1949).
  10. Giorgio Soave, "Equilibrium constants from a modified Redlich-Kwong equation of state", Chemical Engineering Science 27 pp. 1197-1203 (1972).
  11. Ding-Yu Peng and Donald B. Robinson, "A New Two-Constant Equation of State", Industrial & Engineering Chemistry Fundamentals 15 pp. 59-64 (1976).

Discussion