理想気体が嘘をつき始める場所 — 圧縮因子 Z と van der Waals 状態方程式

同じ速度では運動量流束は密度比の2乗です。26%の密度誤差は、ρu²では約60%になります。

100気圧、200 Kの窒素を理想気体で置くと、密度は実在より4分の1近く外れます。圧力も温度も同じです。分子は互いに引き合い、自分の体積も占めます。

そのずれは圧縮因子 $Z$ の一つです。van der Waalsは定数二つで気体を真似て、臨界点の下ではS字ループが液体と気体を分けます。

換算温度のスライダーを1より下へ下げ、Maxwellのチェックを入れると、ループが育ち、面積の等しい水平な飽和線が現れます。

reduced temperature T_r = 0.900Maxwell construction

saturation pressure P_r,sat = 0.647
Z(sat. vapor) = 0.633  ·   Z(sat. liquid) = 0.163

0.85付近では飽和液体の $Z$ が0.1を下回ります。スライダーを1に寄せると、液体と気体の $Z$ が一点に集まります。区別が消える場所です。

100気圧の窒素は理想気体より重い#

$Pv = RT$ は分子を点と見て、間の力を捨てます。希薄なら合います。分子が離れているからです。

圧力を上げると分子の間隔が狭くなります。弱い引力(van der Waals力)が体積を縮めます。同時に、分子自身の体積が残ります。

二つの効果は逆向きです。引力は体積を減らし、有限の体積は体積を増やします。どちらが勝つかは温度と圧力が決めます。だから実在気体が理想気体より重いか軽いかは、条件ごとにひっくり返ります。その密度の順位がひっくり返ると、底の軽い流体が対流を点火する閾値も一緒に動きます。

圧縮因子 Z — ずれを一つの数で#

ずれを一つの数に束ねると圧縮因子です。

Z = \frac{Pv}{RT} = \frac{v_\text{actual}}{v_\text{ideal}}

$P$ は圧力、 $v$ はモル体積、 $R$ は気体定数、 $T$ は温度です。 $Z=1$ なら理想気体です。1から遠いほど、ずれは大きくなります。

$Z<1$ は引力優勢です。実在の体積が理想値より小さいです。 $Z>1$ は有限体積が優勢です。分子が押して体積が増えます。窒素の例の $Z \approx 0.79$ は、実在密度が理想値の $1/0.79 \approx 1.27$ 倍だという意味です。

対応状態原理:すべての気体が似てくる座標

気体ごとに $Z$ 曲線は違います。圧力と温度をそれぞれの臨界値で割ると、曲線はほぼ重なります。

P_R = \frac{P}{P_c}, \qquad T_R = \frac{T}{T_c}

$P_c$ 、 $T_c$ は臨界圧力・温度です。対応状態原理(principle of corresponding states)は、換算圧力 $P_R$ と換算温度 $T_R$ が同じなら $Z$ もほぼ同じだという主張です。

窒素でもメタンでも二酸化炭素でも、 $T_R=1.5$ 、 $P_R=2$ では似た $Z$ を持ちます。分子の名前ではなく、臨界点に対する相対位置が振る舞いを決めます。だから図を一枚で多くの気体を覆えます。

コードの窒素は、その例とは別のマスにいます。 $T = 200$ 、 $T_c = 126.2$ なので $T_R = 200/126.2 \approx 1.585$ です。 $P = 10^7$ 、 $P_c = 3.39\times 10^6$ なので $P_R = 10/3.39 \approx 2.95$ です。出力の $Z = 0.7905$ を換算体積へ戻すと $v_R = Z\,(T/T_c)\,(P_c/P)\,(8/3) \approx 1.13$ です。この点を換算等温線へ戻すと $P_R$ は2.95で戻ってきます。印刷された $Z$ は、その超臨界等温線の上の点です。

van der Waals:二つの定数で分子を真似る#

1873年、van der Waalsは理想気体に補正を二つ入れました。

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

$a/v^2$ は引力が削る圧力です。体積が小さいほど大きくなります。 $b$ は分子が占める体積なので、自由空間は $v-b$ です。

二つの定数は臨界点で決まります。そこでは $P$ – $v$ 等温線が変曲点で、傾きと曲率が同時に0です。

\left(\frac{\partial P}{\partial v}\right)_{T_c} = 0, \qquad \left(\frac{\partial^2 P}{\partial v^2}\right)_{T_c} = 0

解くと定数と臨界圧縮因子が出ます。

a = \frac{27 R^2 T_c^2}{64 P_c}, \qquad b = \frac{R T_c}{8 P_c}, \qquad Z_c = \frac{P_c v_c}{R T_c} = \frac{3}{8}

$Z_c = 3/8 = 0.375$ は気体の種類に依りません。実在気体の $Z_c$ は0.27〜0.29です。0.375は0.28付近に比べて約3分の1大きいです。定量的には過大で、ループと臨界点の絵は合っています。

換算変数で書くと $a$ 、 $b$ が消えます。

P_R = \frac{8 T_R}{3 v_R - 1} - \frac{3}{v_R^2}

この一行が対応状態原理の式です。臨界点は $(v_R, P_R, T_R) = (1, 1, 1)$ に固定されます。分母 $3v_R - 1$ が0になる場所は $v_R = 1/3$ です。その左側は分子体積 $b$ より小さく、等温線は描きません。

臨界点の下の S 字曲線と Maxwell の作図#

$T_R > 1$ なら等温線は単調に減少します。圧力を上げると体積は減ります。

$T_R < 1$ ならS字ループができます。真ん中は $(\partial P/\partial v)_T > 0$ です。圧力を上げたのに体積が増えるという意味で、力学的に不安定です。実在の物質はその区間を飛び越えます。

代わりに水平線の一本が液体と気体を結びます。高さが飽和圧力です。Maxwellの等面積作図(equal-area rule)は、ループから切り取られる上下の面積が等しくなる高さを選びます。液体と気体のGibbs自由エネルギーが等しいという条件と同値です。二相が共存できる圧力は、その一つです。

Python で Z を解く#

van der Waalsを $Z$ で整理すると3次です。臨界定数があれば、 $P$ 、 $T$ で $Z$ が出ます。

import numpy as np
 
R = 8.314  # J/(mol·K)
 
def vdw_constants(Tc, Pc):
    """臨界定数から van der Waals a, b を求める。"""
    a = 27.0 * R**2 * Tc**2 / (64.0 * Pc)
    b = R * Tc / (8.0 * Pc)
    return a, b
 
def vdw_compressibility(P, T, Tc, Pc):
    """Z^3 - (1+B)Z^2 + A Z - A B = 0 の物理的な根。"""
    a, b = vdw_constants(Tc, Pc)
    A = a * P / (R * T)**2          # 引力項の無次元の大きさ
    B = b * P / (R * T)            # 有限体積項の無次元の大きさ
    coeffs = [1.0, -(1.0 + B), A, -A * B]
    roots = np.roots(coeffs)
    real = roots[np.abs(roots.imag) < 1e-8].real
    real = real[real > B]          # v > b  ⇒  Z > B だけが物理的
    return real.max(), real.min(), real.size
 
# 窒素: Tc = 126.2 K, Pc = 3.39 MPa, M = 28 g/mol
Tc, Pc, M = 126.2, 3.39e6, 0.028
P, T = 1.0e7, 200.0                # 100 bar, 200 K
Z_gas, Z_liq, n = vdw_compressibility(P, T, Tc, Pc)
 
rho_ideal = P * M / (R * T)        # 理想気体の密度
rho_real = rho_ideal / Z_gas       # 実在気体の密度
print(f"根の数 = {n},  Z = {Z_gas:.4f}")
print(f"理想気体 密度 = {rho_ideal:6.1f} kg/m^3")
print(f"実在   密度 = {rho_real:6.1f} kg/m^3  ({100*(rho_real/rho_ideal-1):+.0f}%)")

実行結果は3行です。

根の数 = 1,  Z = 0.7905
理想気体 密度 =  168.4 kg/m^3
実在   密度 =  213.0 kg/m^3  (+27%)

$213.0/168.4 = 1.265$ です。増加分は26.5%で、出力の $+27\%$ はその値を整数へ上げたものです。 $T_R \approx 1.59$ なので超臨界です。実根が一つなので液体・気体の区別はありません。それでも理想気体の密度は、実在より約4分の1軽いです。

超臨界インジェクションと CFD — なぜ状態方程式が格子を揺らすのか#

液体ロケットの燃焼室は、ふつう推進剤の臨界圧力より上で動きます。LOX(液体酸素)や極低温窒素がその室へ入ると、見慣れた液滴は生まれません。

表面張力が消えるからです。臨界点を越えると液体と気体の境が崩れます。ジェットは滴へ割れず、密度勾配が櫛の歯のようにほどけて周囲と混ざります。超臨界(または遷臨界)インジェクションです。

状態方程式が格子を揺らすのは、密度が圧力・温度と非線形に結ばれているからです。理想気体を使うとジェット密度が数十%外れ、運動量流束 $\rho u^2$ と混合長さが一緒にずれます。実際のコードは cubic EOS(Peng–Robinson、SRK)や NASG(Noble–Abel Stiffened Gas)を使います。 $Z$ は、そのモデルがセルごとに解く量です。

密度誤差を運動量でも27%と書くと失敗する#

同じ速度では $\rho u^2$ は密度比の2乗です。 $1.265^2 \approx 1.60$ なので、26.5%の密度誤差は運動量流束では約60%になります。密度が27%外れたから流束も27%だと書くと、2乗を落としています。

その圧力を力へそのまま入れると、Maxwellが選んだ水平線の上の界面は静かでしょうか。圧力を直接入れたら静止した界面が震えた非理想格子ボルツマンは、その力が自由エネルギーへ書き換わらないと静まらないのか、まだ開いています。