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

圧縮因子・対応状態原理・van der Waals ループで見る実在気体

100気圧・200 K の窒素を理想気体の法則で計算すると、密度を実際より4分の1以上も外します。同じ圧力・温度なのに、分子は互いを引き合い、自分の体積ぶんの空間を占めます。この記事では、そのずれを一つの数(圧縮因子 $Z$ )でとらえる方法、van der Waals がたった二つの定数で実在気体を真似た仕組み、そして臨界点の下に現れる S 字曲線が液体と気体をどう分けるかを扱います。最後に、この話がロケットインジェクターの超臨界流れと CFD 格子にまで直接届く理由を見ます。

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

理想気体の法則 $Pv = RT$ は分子を点とみなし、分子間の力を無視します。希薄な気体ではよく合います。分子がまばらに離れているからです。

圧力を上げると崩れます。分子が密集します。弱い引力(van der Waals 力)が分子を引き寄せ、体積を縮めます。同時に、分子自身の体積が無視できなくなります。

二つの効果は逆向きです。引力は体積を縮め、有限の大きさは体積を膨らませます。どちらが勝つかは温度と圧力が決めます。だから「実在気体は理想気体より重い、あるいは軽い」は条件ごとにひっくり返ります。

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

このずれを一つの数にまとめたものが圧縮因子です。

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

$P$ は圧力、 $v$ はモル体積、 $R$ は気体定数、 $T$ は温度です。 $Z=1$ なら理想気体。 $Z$ が 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$ を持ちます。分子の種類ではなく、臨界点に対する相対位置が振る舞いを決めます。この普遍性のおかげで、たった一枚の一般化チャートで数多くの気体を扱えます。

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 なので、van der Waals は定量的には過大評価しますが、定性的な描像は正確に当てます。

換算変数で書き直すと、定数 $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)$ に固定されます。

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

$T_R > 1$ なら等温線は単調に減少します。圧力を上げると体積が縮む、直感に合う曲線です。

$T_R < 1$ では曲線がうねります。S 字のループが現れます。真ん中の区間は $(\partial P/\partial v)_T > 0$ です。圧力を上げたのに体積が大きくなるという意味で、力学的に不安定です。実在の物質はこの区間を飛ばします。

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

下のシミュレーションで直接動かしてみましょう。換算温度のスライダーを 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

$T_R$ を 0.85 まで下げると、飽和気体と飽和液体の $Z$ の差が大きく開きます。液体側の $Z$ は 0.1 を下回ります。分子がほぼ接して引力が圧倒するからです。 $T_R$ を 1 に近づけると二つの $Z$ が一点に集まります。臨界点で液体と気体の区別が消える瞬間です。

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}%)")

実行結果は次の通りです。

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

$T_R = 200/126.2 \approx 1.59$ なので超臨界領域です。実根が一つだけで、液体・気体の区別はありません。それでも理想気体で求めた密度は実際より 27% 小さい。冒頭の「4分の1以上」の誤差がここから出ます。

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

この話は黒板の上だけのものではありません。液体ロケットエンジンの燃焼室は、通常、推進剤の臨界圧力より高い圧力で動きます。LOX(液体酸素)や極低温窒素のジェットがそんなチャンバーに噴射されると、おなじみの液滴はできません。

表面張力が消えるからです。臨界点を越えると、液体と気体の境界が崩れます。ジェットは液滴に割れず、密度勾配が櫛のようにほどけて周囲と混ざります。これが超臨界(あるいは遷臨界)インジェクションで、次世代エンジン設計の中心的な流れ現象です。

ここで状態方程式が格子を揺らします。密度が圧力・温度と非線形に結びついているからです。CFD ソルバーが理想気体を使うと、ジェット密度を数十 % 外し、それに伴って運動量フラックス $\rho u^2$ と混合長がまるごとずれます。だから実際のコードは実在気体モデル — cubic EOS(Peng–Robinson、SRK)や NASG(Noble–Abel Stiffened Gas)— を使います。圧縮因子 $Z$ は、そのモデルが格子の一マスごとに解く量そのものです。

覚えておくこと#

圧縮因子 $Z = Pv/RT$ は、実在気体が理想気体からどれだけ外れているかを一つの数でとらえます。 $Z<1$ は引力優勢、 $Z>1$ は有限体積優勢。
van der Waals 式は、引力( $a/v^2$ )と分子体積( $b$ )という二つの定数で、臨界点・ループ・相転移を定性的にすべて再現します。換算変数では $Z_c = 3/8$ に普遍化されます。
臨界点の向こうでは超臨界流れが表面張力なしで混ざり、それを CFD で解くには理想気体ではなく実在の状態方程式が格子の一マスごとに必要です。