第6章: 熱力学と状態方程式

第1章から第5章では超臨界流体を定性的に説明してきました。本章ではそれを定量化する道具立てを提供します。密度と圧縮率を予測する状態方程式、臨界点近傍の異常を支配するスケーリング則、何が何に溶けるかを決める相平衡の条件、そしてこれらすべてを具体的な数値に変換するフガシティの関係式です。8つのPythonコード例により、本章のすべての結果を自分の手で再現できます。

学習目標

本章を修了すると、以下のことができるようになります：

• 理想気体の法則が臨界点近傍で破綻する理由を説明し、圧縮因子 $Z$ でそのずれを定量化できる
• van der Waalsパラメータ $a$、$b$ の物理的意味を理解し、そこから臨界定数と $Z_c = 3/8$ を導出できる
• Peng-Robinson状態方程式で密度と圧縮因子を計算し、偏心因子 $\omega$ の役割を説明できる
• 臨界乳光と $\kappa_T$、$C_P$ の発散を普遍的な臨界指数で記述できる
• 音速が臨界点で極小値を取る理由と、それがプロセス制御に及ぼす影響を説明できる
• 二成分系相図を分類（Type I、II、III）し、Chrastil式で固体溶解度をモデル化できる
• 立方型状態方程式からフガシティ係数を計算し、相平衡条件 $\phi_i^L x_i = \phi_i^V y_i$ を適用できる
• 二成分相互作用パラメータ $k_{ij}$ を含む混合則を多成分系に適用できる

6.1 超臨界流体の状態方程式

理想気体の法則の限界

理想気体の状態方程式

は、分子を体積の無視できる点粒子とみなし、分子間力が存在せず、衝突はすべて弾性的であると仮定しています。しかし超臨界流体が実際に使われる条件では、これらの仮定はいずれも成り立ちません。

1. 高圧下：分子間距離が縮まり、van der Waals力による引力が無視できなくなる。
2. 高密度：分子自身が占める排除体積が全体積に対して無視できない割合になる。
3. 臨界点近傍：気液の区別が消失し、密度ゆらぎが極端に大きくなる。

理想性からのずれは圧縮因子で定量化されます。

ここで $M$ はモル質量です。理想気体では $Z = 1$ ですが、臨界点近傍では $Z$ は概ね 0.2 から 1.2 の範囲で変化します。

van der Waals方程式

van der Waals（1873年）は、分子の大きさと分子間引力の2つの補正を理想気体の式に加えました。

ここで $V_m$ はモル体積（m³/mol）、$a$ は引力パラメータ（Pa·m⁶·mol⁻²）、$b$ は排除体積パラメータ（m³/mol）です。

圧力について解くと、等温線を描くのに便利な形になります。

等温線の形状は温度によって次のように変わります。

• $T > T_c$（超臨界領域）：単調減少し、相転移は現れない
• $T = T_c$（臨界等温線）：水平な変曲点をもつ
• $T < T_c$（二相領域）：S字型のループが現れ、その非物理的な部分はMaxwellの等面積則により水平な共存線に置き換えられる

臨界点の条件。臨界点では等温線が水平な変曲点をもちます。

この2つの条件を連立して解くと、臨界定数が得られます。

これを逆に解けば、文献値の臨界定数から $a$ と $b$ を求められます。

Peng-Robinson方程式

Peng-Robinson方程式（1976年）は、炭化水素、CO₂、超臨界プロセス設計における実用上の標準的な立方型状態方程式です。

ここで各パラメータは次のように与えられます。

偏心因子 $\omega$ は、分子が球対称な単純流体からどれだけ離れているかを表す指標です。

立方形式。数値計算では $Z$ についての3次方程式に書き換えます。

ここで

二相領域の内部では、この3次方程式は3つの実正根をもちます。最大の根が気相、最小の根が液相に対応し、中間の根は熱力学的に不安定で物理的意味をもちません。単相の超臨界領域では実正根は1つだけになります。

状態方程式の比較

6.2 臨界現象

臨界乳光

臨界点に近づくと、透明だった流体が乳白色に濁ります。これが臨界乳光であり、臨界点で何か特異なことが起きていることを最も直接的に示す視覚的証拠です。

1. 密度ゆらぎ：$T_c$ 近傍ではゆらぎを生じるための自由エネルギー的コストがゼロに近づくため、局所密度が大きな振幅で長距離にわたってゆらぐ。
2. 相関長：相関したゆらぎの空間的広がり $\xi$ が次のように発散する。 $$ \xi \sim |T - T_c|^{-\nu} $$ ここで $\nu \approx 0.63$ は相関長の臨界指数。
3. 光散乱：$\xi$ が可視光の波長（$\lambda \sim 400$〜$700$ nm）と同程度になると、Rayleigh散乱が劇的に強まる。

散乱強度は次のように増大します。

これが、$T_c$ から数ケルビン離れれば超臨界流体が完全に透明に見え、$T_c$ の1ケルビン以内では不透明になる理由です。

応答関数の発散

等温圧縮率：

ここで $\gamma \approx 1.24$ です。定圧熱容量：

ここで $\alpha \approx 0.11$ です。

普遍性と臨界指数

驚くべきことに、CO₂、水、キセノンという構造的に何の共通点もない分子が、同じ臨界指数を共有します。これが普遍性です。

流体における秩序変数は、共存する2相の密度差です。

また臨界等温線（$T = T_c$）上では次が成り立ちます。

これらの指数は独立ではなく、スケーリング関係で結ばれています。

くりこみ群理論

くりこみ群理論（Wilson, 1971年）は普遍性を説明します。その中心的な考え方は、臨界点では系が長さスケールの変換に対して不変になり、短距離の詳細が次々と洗い流されていくというものです。最終的に残るのは次の3つだけです。

• 空間次元（流体では $d = 3$）
• 秩序変数の次元（スカラー量である密度では $n = 1$）
• 相互作用の対称性

分子量や結合角、化学種の違いは、専門用語の意味で「無関係（irrelevant）」になります。実用上の見返りは大きく、1つの流体の臨界挙動を丁寧に測定すれば、すべての流体について知ったことになります。これは対応状態原理の理論的基礎でもあり、単一の換算物性チャートが多くの物質に使える理由でもあります。

6.3 臨界点近傍の熱力学的性質

エンタルピーとエントロピー

理想的な挙動からのずれは残差項の積分としてまとめられます。

$T_c$ 近傍ではこれらの積分項が大きく、かつ強く非線形になります。熱膨張係数 $(\partial V/\partial T)_P$ 自体が発散するためです。

臨界点では共存する2相が同一になるため、蒸発潜熱と蒸発エントロピーはいずれもゼロに収束します。

その結果、Clausius-Clapeyronの式は不定形になります。

これは、蒸気圧曲線が超臨界領域へ延びていくのではなく、臨界点で文字どおり終端することの熱力学的な表現です。

実用上の帰結：超臨界プロセスの熱交換器は、狭い温度範囲で生じる非常に大きなエンタルピー変化を前提に設計しなければなりません。

熱容量の異常

定圧熱容量

は臨界点で鋭いピークを示します。$P \approx P_c$ におけるCO₂の場合：

物理的には、$T_c$ 近傍で加えた熱は温度上昇ではなく分子構造の再編成（クラスターの生成と解裂）に使われます。定積熱容量 $C_V$ も異常を示しますが、$C_P$ よりはるかに緩やかです。

熱容量比は標準的な熱力学関係から得られます。

ここで $\alpha_P = (1/V)(\partial V/\partial T)_P$ です。臨界点では両方の熱容量が発散しますが、その比は有限に留まり1に近づきます。

音速

音速は断熱微分で与えられます。

実在流体では次のように書けます。

臨界点近傍では $\kappa_T$ の発散により $(\partial P/\partial \rho)_T \to 0$ となり、一方で $C_P/C_V \to 1$、断熱圧縮率 $\kappa_S$ は有限に留まります。その結果、音速は $T_c$ で極小値を取ります。$P_c$ におけるCO₂では：

• $T = 290$ K（液相）：$c \approx 900$ m/s
• $T = 304$ K（$T_c$）：$c \approx 180$ m/s — 極小値
• $T = 320$ K（超臨界）：$c \approx 250$ m/s

プロセス設計上の帰結：伝播時間から距離を換算する超音波式の液面計や流量計は、仮定した音速が成り立たなくなるため、臨界点近傍では信頼できません。

プロセス設計への影響

6.4 超臨界流体系における相平衡

気液平衡

亜臨界領域では蒸気圧曲線に沿って2相が共存します。平衡条件は化学ポテンシャルの一致、すなわちフガシティの一致です。

蒸気圧曲線の傾きはClausius-Clapeyronの式で与えられます。

蒸気を理想気体とみなし、蒸発潜熱が温度に依らないと仮定すると、よく知られた積分形 $\ln P = -\Delta H_{vap}/(RT) + C$ が得られます。実務では経験式であるAntoine式が用いられます。

ここで $A$、$B$、$C$ は物質固有の定数で、適用温度範囲と単位系を明記して表にまとめられています。$T \to T_c$ では $\Delta H_{vap} \to 0$、$V_V - V_L \to 0$ となり、曲線は臨界点で終端します。

超臨界流体を含む二成分系相図

二成分系（超臨界溶媒＋溶質）の相挙動は $P$-$T$-$x$ 曲面上に表されます。全体像を整理する鍵は次の2つです。

• 臨界曲線（critical locus）：純成分1の臨界点と純成分2の臨界点を結ぶ曲線。単調である必要はありません。
• 三相線：固体・液体・気体が共存する線で、純成分の三重点から延びています。

古典的なvan Konynenburg-Scott分類はいくつかの型を区別しますが、超臨界プロセスで重要なのは次の3つです。

Type IIIの挙動は厄介な現象ではなく道具です。そこに現れる強い圧力感受性こそが、選択性の高い超臨界分別を可能にしています。

溶解度のモデル化：Chrastil式

Chrastil（1982年）は、超臨界流体中の固体溶解度に対する半経験的な相関式を提案しました。

ここで $S$ は溶解度、$\rho$ は超臨界流体の密度（kg/m³）、$k$ は会合数、$a$ は溶媒化熱と蒸発熱に関係する定数（K）、$b$ は定数です。$S$ は文献によって単位体積あたり（kg/m³）または単位質量あたり（kg/kg）で報告されており、その違いはフィッティングされた $b$ が吸収してしまいます。公表パラメータを再利用する前に必ず単位を確認してください。

一定温度では、溶解度は密度（したがって圧力）とともに急激に増加します。一定密度では、溶媒化が発熱的であるため（$a < 0$）、通常は温度上昇とともに減少します。

逆行凝縮

逆行凝縮とは、一定圧力下で温度を上げると液相が現れるという直観に反する現象です（ガスコンデンセート貯留層では等温減圧でも同様のことが起こります）。

1. 臨界圧力より高い単相の超臨界領域から出発する。
2. 等圧のまま温度を上げる。
3. 相境界を越え、液滴が現れる。
4. さらに加熱すると液滴は再溶解して消える。

分子論的な説明：低温側では超臨界流体の密度、したがって溶媒力が高く、すべてが溶解したままです。温度が上がると、溶質の蒸気圧の増加よりも密度の低下のほうが速く進むため、溶質が析出します。この現象は $T > T_c$、$P \approx P_c$ で、混合物の臨界曲線が負の勾配をもつ領域で起こります。

重要となる場面：天然ガス処理（メタンからのコンデンセート回収）、超臨界逆溶媒析出法（SAS）、CO₂圧入による石油増進回収。

6.5 熱力学計算

フガシティとフガシティ係数

フガシティ $f$ は、化学ポテンシャルの理想気体的表現を厳密に成り立たせる「実効圧力」です。

フガシティ係数は理想性からのずれを表します。

理想気体では $\phi = 1$ です。任意の状態方程式から積分により求められます。

Peng-Robinson方程式についてこの積分を実行すると閉じた表式が得られます。立方型状態方程式が今なお使われ続ける主要な理由の1つです。

$A$ と $B$ は6.1節で定義したものです。液相と気相の間の平衡条件は次のように書けます。

ここで $x_i$、$y_i$ はそれぞれ液相と気相のモル分率です。

超臨界流体相における化学ポテンシャル

混合物中の成分 $i$ について、

液体的な密度領域では、活量係数を用いるほうが便利なことが多いです。

固相と超臨界相における溶質の化学ポテンシャルを等置する条件 $\mu_i^{solid} = \mu_i^{SCF}$ から溶解度が直接得られ、無限希釈の極限では

となります。これにより純成分の物性のみから溶解度を予測できます。最後の項がPoynting補正です。

多成分系の混合則

混合物に対する立方型状態方程式のパラメータは、van der Waals型の一流体混合則を用いて純成分パラメータから構成します。

二成分相互作用パラメータ $k_{ij}$ は、異種分子対に対する幾何平均近似を補正します。実験データからフィッティングされ、$k_{ii} = 0$ を満たし、値は小さいものの決して無視できません。

$k_{ij}$ が正であることは、異種分子間の引力が幾何平均より弱いことを意味します（CO₂と炭化水素の組み合わせでは通常 $k_{ij} \approx 0.1$〜$0.15$）。負の値はまれです。より高度な混合則（Wong-Sandler、MHV2など）は活量係数モデルを混合則に組み込み、強い極性をもつ混合物に対してかなり高い精度を与えます。

6.6 Pythonコード例

例1: van der Waals等温線

van der Waals方程式からCO₂の $P$-$V$ 等温線を、亜臨界・臨界・超臨界の各温度について描き、$Z_c = 3/8$ を数値的に確認します。

import numpy as np
import matplotlib.pyplot as plt

# CO2のvan der Waalsパラメータ
R = 8.314  # J/(mol·K)
Tc = 304.1  # K
Pc = 7.38e6  # Pa

# 臨界定数からaとbを算出
a = 27 * R**2 * Tc**2 / (64 * Pc)  # Pa·m^6/mol^2
b = R * Tc / (8 * Pc)  # m^3/mol

# モル体積の範囲（Vm = b の特異点を避ける）
Vm = np.linspace(1.5*b, 50*b, 1000)

# 温度リスト: 亜臨界、臨界、超臨界
temperatures = [280, 304.1, 320, 350]
colors = ['blue', 'red', 'green', 'purple']

plt.figure(figsize=(10, 6))

for T, color in zip(temperatures, colors):
    # van der Waals式による圧力
    P = R * T / (Vm - b) - a / Vm**2

    # プロット用にMPaへ変換
    P_MPa = P / 1e6

    label = f'T = {T} K'
    if T == Tc:
        label += '（臨界）'

    plt.plot(Vm * 1e6, P_MPa, color=color, linewidth=2, label=label)

# 臨界点をマーク
Vm_c = 3 * b
P_c_vdw = R * Tc / (Vm_c - b) - a / Vm_c**2
plt.plot(Vm_c * 1e6, P_c_vdw / 1e6, 'ro', markersize=10, label='臨界点')

plt.xlabel('モル体積 (cm³/mol)', fontsize=12)
plt.ylabel('圧力 (MPa)', fontsize=12)
plt.title('CO₂のvan der Waals等温線', fontsize=14, fontweight='bold')
plt.legend(fontsize=10)
plt.grid(alpha=0.3)
plt.xlim(0, 500)
plt.ylim(0, 15)
plt.tight_layout()
plt.show()

print("CO2のvan der Waalsパラメータ:")
print(f"a = {a:.4e} Pa·m^6/mol^2")
print(f"b = {b:.4e} m^3/mol")
print(f"臨界圧縮因子: Zc = {Pc * Vm_c / (R * Tc):.3f}")

着目点：亜臨界の等温線（$T < T_c$）にはMaxwellの等面積則で水平な共存線に置き換えられるべき非物理的なS字ループが現れます。臨界等温線は水平な変曲点をもち、超臨界の等温線は単調です。出力された $Z_c = 0.375$ はちょうど $3/8$ です。

例2: 臨界点の決定

等温線の1階微分と2階微分が同時にゼロとなる点として臨界点を求め、解析解を数値的に検証します。

import numpy as np
import matplotlib.pyplot as plt

def van_der_waals_P(Vm, T, a, b, R=8.314):
    """van der Waals式による圧力。"""
    return R * T / (Vm - b) - a / Vm**2

def dP_dVm(Vm, T, a, b, R=8.314):
    """圧力のモル体積に関する1階微分。"""
    return -R * T / (Vm - b)**2 + 2 * a / Vm**3

def d2P_dVm2(Vm, T, a, b, R=8.314):
    """圧力のモル体積に関する2階微分。"""
    return 2 * R * T / (Vm - b)**3 - 6 * a / Vm**4

def find_critical_point(a, b, R=8.314):
    """
    van der Waalsパラメータから臨界点を求める。

    Returns:
    --------
    Tc, Pc, Vm_c : float
        臨界温度、臨界圧力、臨界モル体積
    """
    # van der Waals式の解析解
    Tc = 8 * a / (27 * R * b)
    Vm_c = 3 * b
    Pc = a / (27 * b**2)

    return Tc, Pc, Vm_c

def verify_critical_conditions(Vm_c, Tc, a, b):
    """
    臨界点で1階微分と2階微分が消えることを検証する。

    SI単位ではこれらの微分値が巨大なスケール因子を含むため
    （dP/dVm は Pc/Vm_c ~ 1e11 Pa per m^3/mol のオーダー）、
    絶対値による許容判定は意味をもたない。無次元化した残差で比較する。
    """
    dP = dP_dVm(Vm_c, Tc, a, b)
    d2P = d2P_dVm2(Vm_c, Tc, a, b)

    Pc = a / (27 * b**2)
    res1 = dP * Vm_c / Pc
    res2 = d2P * Vm_c**2 / Pc

    print("臨界点における検証:")
    print(f"  (dP/dVm)   * Vm_c/Pc   = {res1:.2e} （ほぼ0であるべき）")
    print(f"  (d2P/dVm2) * Vm_c^2/Pc = {res2:.2e} （ほぼ0であるべき）")

    return np.abs(res1) < 1e-9 and np.abs(res2) < 1e-9

# 例: CO2の臨界点を求める
R = 8.314
a = 0.3658  # Pa·m^6/mol^2（実験データにフィッティングした値）
b = 4.267e-5  # m^3/mol

Tc, Pc, Vm_c = find_critical_point(a, b, R)

print("CO2の臨界点の決定")
print("=" * 50)
print("van der Waalsパラメータ:")
print(f"  a = {a:.4f} Pa·m^6/mol^2")
print(f"  b = {b:.2e} m^3/mol")
print("\n計算された臨界点:")
print(f"  Tc = {Tc:.2f} K")
print(f"  Pc = {Pc/1e6:.2f} MPa")
print(f"  Vm,c = {Vm_c*1e6:.2f} cm^3/mol")
print(f"  Zc = {Pc * Vm_c / (R * Tc):.4f}")
print("\n実験値:")
print("  Tc = 304.1 K")
print("  Pc = 7.38 MPa")
print("  Zc = 0.274")

# 条件の検証
is_critical = verify_critical_conditions(Vm_c, Tc, a, b)
print(f"\n臨界点の条件を満たすか: {is_critical}")

# 臨界温度近傍の等温線をプロット
Vm_range = np.linspace(1.5*b, 10*b, 500)
temps = [Tc - 5, Tc, Tc + 5]
colors = ['blue', 'red', 'green']

plt.figure(figsize=(10, 6))

for T, color in zip(temps, colors):
    P_values = [van_der_waals_P(Vm, T, a, b, R) for Vm in Vm_range]
    label = f'T = {T:.1f} K'
    if T == Tc:
        label += '（臨界）'
    plt.plot(Vm_range*1e6, np.array(P_values)/1e6, color=color,
             linewidth=2, label=label)

# 臨界点をマーク
Pc_calc = van_der_waals_P(Vm_c, Tc, a, b, R)
plt.plot(Vm_c*1e6, Pc_calc/1e6, 'ko', markersize=10, label='臨界点')

# 臨界点における接線（水平になるはず）
tangent_Vm = np.linspace(0.9*Vm_c, 1.1*Vm_c, 50)
tangent_P = np.ones_like(tangent_Vm) * Pc_calc / 1e6
plt.plot(tangent_Vm*1e6, tangent_P, 'k--', linewidth=1.5, alpha=0.7,
         label='接線（水平）')

plt.xlabel('モル体積 (cm³/mol)', fontsize=12)
plt.ylabel('圧力 (MPa)', fontsize=12)
plt.title('臨界点: van der Waals等温線の変曲点',
          fontsize=14, fontweight='bold')
plt.legend(fontsize=10)
plt.grid(alpha=0.3)
plt.xlim(0, 400)
plt.ylim(0, 12)
plt.tight_layout()
plt.show()

例3: Peng-Robinson式による密度計算

再利用可能なPeng-Robinsonクラスです。3次方程式を解いて $Z$ を求め、気相根または液相根を選び、質量密度に変換します。

import numpy as np
import matplotlib.pyplot as plt

# 純成分に対するPeng-Robinson状態方程式
class PengRobinson:
    def __init__(self, Tc, Pc, omega):
        """
        PR状態方程式のパラメータを初期化する。

        Parameters:
        -----------
        Tc : float
            臨界温度 (K)
        Pc : float
            臨界圧力 (Pa)
        omega : float
            偏心因子
        """
        self.R = 8.314  # J/(mol·K)
        self.Tc = Tc
        self.Pc = Pc
        self.omega = omega

        # aとbの算出
        self.a = 0.45724 * self.R**2 * Tc**2 / Pc
        self.b = 0.07780 * self.R * Tc / Pc

        # κパラメータ
        self.kappa = 0.37464 + 1.54226 * omega - 0.26992 * omega**2

    def alpha(self, T):
        """引力項の温度補正。"""
        Tr = T / self.Tc
        return (1 + self.kappa * (1 - np.sqrt(Tr)))**2

    def solve_Z(self, T, P):
        """
        圧縮因子Zについて3次方程式を解く。

        Returns:
        --------
        Z_positive : array
            3次方程式の実正根
        """
        A = self.a * self.alpha(T) * P / (self.R * T)**2
        B = self.b * P / (self.R * T)

        # 3次方程式: Z^3 + p*Z^2 + q*Z + r = 0
        p = -(1 - B)
        q = A - 3*B**2 - 2*B
        r = -(A*B - B**2 - B**3)

        # 3次方程式を解く
        coeffs = [1, p, q, r]
        Z_roots = np.roots(coeffs)

        # 実正根のみを返す
        Z_real = Z_roots[np.isreal(Z_roots)].real
        Z_positive = Z_real[Z_real > 0]

        return Z_positive

    def density(self, T, P, phase='vapor'):
        """
        モル密度を計算する。

        Parameters:
        -----------
        phase : str
            'vapor'（最大のZ）または 'liquid'（最小のZ）

        Returns:
        --------
        rho : float
            モル密度 (mol/m^3)
        Z : float
            選択された圧縮因子
        """
        Z_values = self.solve_Z(T, P)

        if len(Z_values) == 0:
            raise ValueError("実正根が見つかりません")

        if phase == 'vapor':
            Z = np.max(Z_values)
        elif phase == 'liquid':
            Z = np.min(Z_values)
        else:
            raise ValueError("phase は 'vapor' または 'liquid' を指定してください")

        Vm = Z * self.R * T / P  # モル体積 (m^3/mol)
        rho = 1 / Vm  # モル密度 (mol/m^3)

        return rho, Z

# CO2の物性
co2 = PengRobinson(Tc=304.1, Pc=7.38e6, omega=0.225)

# 各条件での試算
test_conditions = [
    (310, 8e6, '超臨界'),
    (350, 15e6, '高圧・高温'),
    (280, 5e6, '亜臨界（液相）'),
    (280, 5e6, '亜臨界（気相）')
]

print("CO2に対するPeng-Robinson式の計算結果:")
print("=" * 70)

for T, P, description in test_conditions:
    try:
        if '液相' in description:
            rho, Z = co2.density(T, P, phase='liquid')
        else:
            rho, Z = co2.density(T, P, phase='vapor')

        # kg/m^3へ変換（CO2の分子量 = 44.01 g/mol）
        rho_kg = rho * 44.01 / 1000

        print(f"\n{description}:")
        print(f"  T = {T} K, P = {P/1e6:.1f} MPa")
        print(f"  密度 = {rho_kg:.1f} kg/m^3")
        print(f"  圧縮因子 Z = {Z:.3f}")

    except Exception as e:
        print(f"\n{description}: エラー - {e}")

# 一定温度における密度の圧力依存性をプロット
T_iso = 313  # K（超臨界）
P_range = np.linspace(7.5e6, 30e6, 50)
densities = []

for P in P_range:
    rho, _ = co2.density(T_iso, P, phase='vapor')
    densities.append(rho * 44.01 / 1000)  # kg/m^3へ変換

plt.figure(figsize=(8, 5))
plt.plot(P_range / 1e6, densities, 'b-', linewidth=2)
plt.xlabel('圧力 (MPa)', fontsize=12)
plt.ylabel('密度 (kg/m³)', fontsize=12)
plt.title(f'CO₂の密度の圧力依存性 T = {T_iso} K（Peng-Robinson）',
          fontsize=13, fontweight='bold')
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

着目点：280 K・5 MPaでは3次方程式が3つの実根をもち、同じ呼び出しから液相根（869 kg/m³）と気相根（203 kg/m³）の両方が得られます。二相領域の本質がここに凝縮されています。立方型状態方程式における気液の区別は、最小根か最大根のどちらを選ぶかという問題に帰着します。$T_c$ より上では実根は1つだけになり、密度-圧力曲線は $P_c$ のすぐ上で超臨界CO₂の調整可能性を生む強い非線形性を示します。参照値からのずれは5〜10%程度、$T_c$ の数ケルビン以内では10〜20%まで拡大します。第7章でCoolPropと比較して定量化します。

例4: 一般化圧縮因子チャート

換算座標（$T_r$、$P_r$）を用いると多くの流体の挙動が1枚のチャートに集約されます。この例では古典的な圧縮因子チャートを図から読み取る代わりに、Peng-Robinson式から再構成します。

import numpy as np
import matplotlib.pyplot as plt

def PR_compressibility(Tr, Pr, omega):
    """
    換算座標におけるPeng-Robinson式から圧縮因子を計算する。

    Parameters:
    -----------
    Tr : float
        換算温度 T/Tc
    Pr : float
        換算圧力 P/Pc
    omega : float
        偏心因子

    Returns:
    --------
    Z : float
        圧縮因子（最大の実根）
    """
    # PRパラメータ
    kappa = 0.37464 + 1.54226 * omega - 0.26992 * omega**2
    alpha = (1 + kappa * (1 - np.sqrt(Tr)))**2

    # 換算パラメータ
    a_r = 0.45724 * alpha / Tr**2
    b_r = 0.07780 / Tr

    A = a_r * Pr
    B = b_r * Pr

    # 3次方程式を解く
    coeffs = [1, -(1-B), A - 3*B**2 - 2*B, -(A*B - B**2 - B**3)]
    roots = np.roots(coeffs)

    # 最大の実根（気相的な相）を採る
    real_roots = roots[np.isreal(roots)].real
    Z = np.max(real_roots) if len(real_roots) > 0 else 1.0

    return Z

# 換算物性のメッシュを作成
Tr_range = np.linspace(0.7, 2.0, 100)
Pr_range = np.linspace(0.1, 5.0, 100)
Tr_grid, Pr_grid = np.meshgrid(Tr_range, Pr_range)

# CO2（omega = 0.225）についてZを計算
Z_grid = np.zeros_like(Tr_grid)

for i in range(len(Pr_range)):
    for j in range(len(Tr_range)):
        Z_grid[i, j] = PR_compressibility(Tr_grid[i, j], Pr_grid[i, j],
                                          omega=0.225)

# 等高線マップ
plt.figure(figsize=(12, 8))

contour = plt.contourf(Tr_grid, Pr_grid, Z_grid, levels=20, cmap='viridis')
plt.colorbar(contour, label='圧縮因子 Z')

# 等高線
contour_lines = plt.contour(Tr_grid, Pr_grid, Z_grid, levels=10,
                            colors='white', linewidths=0.5, alpha=0.5)
plt.clabel(contour_lines, inline=True, fontsize=8, fmt='%.2f')

# 臨界点をマーク
plt.plot(1.0, 1.0, 'r*', markersize=20, label='臨界点 (Tr=1, Pr=1)')

plt.xlabel('換算温度 Tr = T/Tc', fontsize=13)
plt.ylabel('換算圧力 Pr = P/Pc', fontsize=13)
plt.title('一般化圧縮因子チャート（Peng-Robinson, ω=0.225）',
          fontsize=14, fontweight='bold')
plt.legend(fontsize=11)
plt.grid(alpha=0.2, color='white')
plt.xlim(0.7, 2.0)
plt.ylim(0.1, 5.0)
plt.tight_layout()
plt.show()

# 各等温線についてZ対Prをプロット
plt.figure(figsize=(10, 6))

Tr_values = [0.9, 1.0, 1.05, 1.1, 1.2, 1.5]
colors = plt.cm.coolwarm(np.linspace(0, 1, len(Tr_values)))

for Tr, color in zip(Tr_values, colors):
    Z_values = [PR_compressibility(Tr, Pr, 0.225) for Pr in Pr_range]
    label = f'Tr = {Tr:.2f}'
    if Tr == 1.0:
        label += '（臨界）'
    plt.plot(Pr_range, Z_values, color=color, linewidth=2.5, label=label)

# 理想気体の線
plt.plot(Pr_range, np.ones_like(Pr_range), 'k--', linewidth=1.5,
         label='理想気体 (Z=1)')

plt.xlabel('換算圧力 Pr = P/Pc', fontsize=12)
plt.ylabel('圧縮因子 Z', fontsize=12)
plt.title('CO₂の圧縮因子の換算圧力依存性',
          fontsize=13, fontweight='bold')
plt.legend(fontsize=10, loc='best')
plt.grid(alpha=0.3)
plt.xlim(0, 5)
plt.ylim(0, 1.5)
plt.tight_layout()
plt.show()

# 代表的な条件での値を出力
print("\n代表的な条件における圧縮因子:")
print("=" * 60)
test_cases = [
    (1.0, 1.0, "臨界点"),
    (1.05, 1.5, "超臨界（典型的な抽出条件）"),
    (1.2, 2.0, "超臨界（高密度）"),
    (0.9, 0.5, "亜臨界の気相"),
]

for Tr, Pr, description in test_cases:
    Z = PR_compressibility(Tr, Pr, 0.225)
    print(f"{description:28s}: Tr={Tr:.2f}, Pr={Pr:.2f} -> Z={Z:.3f}")

着目点：$T_r$ 一定のもとで $P_r$ が増すと $Z$ は低下し（引力が支配的）、$P_r$ が小さい領域では1に戻ります（理想気体極限）。そして $T_r = 1$、$P_r \sim 1$〜$2$ の付近で最小値に近づきます。これはまさに抽出に用いられる領域です。$T_r = P_r = 1$ における値はPeng-Robinson式の臨界圧縮因子で、解析的には0.307ですが、その付近で3次方程式が極端に平坦になるため数値解は0.321を返します。

例5: 臨界指数と秩序変数

べき乗則 $\rho_L - \rho_V \sim |T - T_c|^{\beta}$ は両対数プロット上で傾き $\beta$ の直線になるため、そこで検証するのが最も確実です。

import numpy as np
import matplotlib.pyplot as plt

def density_difference(T, Tc, rho_c, beta=0.326):
    """
    臨界点近傍の気液密度差。

    rho_L - rho_V ~ |T - Tc|^beta
    """
    epsilon = np.abs((T - Tc) / Tc)
    # 振幅B0は物質固有。ここでは 2*rho_c を妥当なスケールとして用いる
    B0 = 2.0 * rho_c
    return B0 * epsilon**beta

# CO2のパラメータ
Tc = 304.1  # K
rho_c = 467.6  # kg/m^3
beta = 0.326  # 臨界指数

# 臨界温度以下の温度範囲
T = np.linspace(280, Tc - 0.01, 100)

# 密度差
delta_rho = density_difference(T, Tc, rho_c, beta)

# 換算温度差
epsilon = (Tc - T) / Tc

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 線形プロット
ax1.plot(T, delta_rho, 'b-', linewidth=2)
ax1.set_xlabel('温度 (K)', fontsize=12)
ax1.set_ylabel('密度差 ρ_L − ρ_V (kg/m³)', fontsize=12)
ax1.set_title('臨界点近傍の気液密度差', fontsize=13)
ax1.axvline(Tc, color='red', linestyle='--', label=f'Tc = {Tc} K')
ax1.legend()
ax1.grid(alpha=0.3)

# 両対数プロット: べき乗則は傾きβの直線になる
ax2.loglog(epsilon, delta_rho, 'bo-', linewidth=2, markersize=3, label='データ')
fit_line = delta_rho[0] * (epsilon / epsilon[0])**beta
ax2.loglog(epsilon, fit_line, 'r--', linewidth=2,
           label=f'べき乗則 (β={beta})')
ax2.set_xlabel('換算温度差 ε = (Tc − T) / Tc', fontsize=12)
ax2.set_ylabel('密度差 (kg/m³)', fontsize=12)
ax2.set_title('べき乗則の検証（両対数）', fontsize=13)
ax2.legend()
ax2.grid(alpha=0.3, which='both')

plt.tight_layout()
plt.show()

# 生成データから線形回帰で指数を復元する
slope, intercept = np.polyfit(np.log(epsilon), np.log(delta_rho), 1)
print(f"仮定した臨界指数 β = {beta}")
print(f"両対数フィットから復元した傾き = {slope:.3f}")

着目点：線形プロットでは $T \to T_c$ に向かって共存曲線が平坦化していく様子が見え、両対数プロットではそれが直線になり、その傾きがそのまま臨界指数になります。これは実験の共存曲線データを解析する標準的な方法であり、同じ手順を実際のCO₂のデータに適用すると0.32〜0.33が得られます。van der Waals式が予測する平均場的な値0.5ではありません。

例6: Peng-Robinson式によるフガシティ係数

フガシティ係数は、状態方程式を相平衡計算へと変換する要です。ここではCO₂について広い圧力範囲で閉じた表式から評価します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

def fugacity_coefficient_PR(T, P, Tc, Pc, omega):
    """
    Peng-Robinson状態方程式からフガシティ係数を計算する。

    Returns
    -------
    phi : float
        フガシティ係数 f/P
    Z : float
        圧縮因子
    """
    R = 8.314  # J/(mol·K)

    # PRパラメータ
    a = 0.45724 * R**2 * Tc**2 / Pc
    b = 0.07780 * R * Tc / Pc
    kappa = 0.37464 + 1.54226*omega - 0.26992*omega**2
    alpha = (1 + kappa*(1 - np.sqrt(T/Tc)))**2

    A = a * alpha * P / (R**2 * T**2)
    B = b * P / (R * T)

    # 立方形式から圧縮因子を求める
    def equation(Z):
        return Z**3 - (1-B)*Z**2 + (A - 3*B**2 - 2*B)*Z - (A*B - B**2 - B**3)

    Z = fsolve(equation, 1.0)[0]

    # フガシティ係数の閉じた表式
    sqrt2 = np.sqrt(2)
    ln_phi = ((Z - 1) - np.log(Z - B)
              - A/(2*sqrt2*B) * np.log((Z + (1+sqrt2)*B) / (Z + (1-sqrt2)*B)))

    return np.exp(ln_phi), Z

# CO2のパラメータ
Tc = 304.1  # K
Pc = 7.38e6  # Pa
omega = 0.225

temperatures = [300, 320, 350]  # K
pressures = np.linspace(0.1e6, 20e6, 100)  # Pa

plt.figure(figsize=(10, 6))
for T in temperatures:
    phi_values = []
    for P in pressures:
        phi, Z = fugacity_coefficient_PR(T, P, Tc, Pc, omega)
        phi_values.append(phi)

    label = f'{T} K'
    label += '（亜臨界）' if T < Tc else '（超臨界）'
    plt.plot(pressures/1e6, phi_values, linewidth=2, label=label)

plt.axhline(1.0, color='gray', linestyle='--', label='理想気体 (φ=1)')
plt.axvline(Pc/1e6, color='red', linestyle='--', alpha=0.5, label='臨界圧力')
plt.xlabel('圧力 (MPa)', fontsize=12)
plt.ylabel('フガシティ係数 φ = f/P', fontsize=12)
plt.title('CO₂のフガシティ係数（Peng-Robinson）', fontsize=14)
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# 特定条件での値
T_target = 320  # K
P_target = 10e6  # Pa
phi, Z = fugacity_coefficient_PR(T_target, P_target, Tc, Pc, omega)
f = phi * P_target
print(f"\nT = {T_target} K, P = {P_target/1e6:.1f} MPa:")
print(f"圧縮因子 Z = {Z:.3f}")
print(f"フガシティ係数 φ = {phi:.3f}")
print(f"フガシティ f = {f/1e6:.2f} MPa")
print(f"理想性からのずれ: {abs(1-phi)*100:.1f}%")

着目点：$P \to 0$ では $\phi \to 1$ となり、圧力が上がるにつれて引力の効果で1を大きく下回ります。10 MPa・320 Kではフガシティはわずか6.0 MPa、圧力より40%低い値です。プロセス条件下で超臨界CO₂を理想気体として扱うのは小さな近似ではなく、理想気体を仮定した溶解度計算は同程度の倍率で誤ります。

例7: Chrastil式による溶解度モデルのフィッティング

Chrastil式は $\ln\rho$ と $1/T$ について線形なので、対数変換したデータに対する通常最小二乗法で3つのパラメータを一度に求められます。非線形ソルバーも初期値も不要です。

import numpy as np
import matplotlib.pyplot as plt

# この例ではCO2密度の算出にコード例3のPengRobinsonクラスを再利用する。
# 先にそのブロックを実行するか、同じスクリプト内に配置すること。

def chrastil_solubility(rho, T, k, a, b):
    """
    Chrastilの溶解度モデル。

    ln(S) = k * ln(rho) + a/T + b

    Parameters:
    -----------
    rho : float or array
        超臨界流体の密度 (kg/m^3)
    T : float or array
        温度 (K)
    k, a, b : float
        Chrastilパラメータ

    Returns:
    --------
    S : float or array
        溶解度 (kg溶質 / m^3 超臨界流体)
    """
    ln_S = k * np.log(rho) + a / T + b
    return np.exp(ln_S)

# 実験データ: 超臨界CO2中のβ-カロテン（文献の代表的な値）
data = {
    'T': np.array([313, 313, 313, 333, 333, 333, 353, 353, 353]),  # K
    'P': np.array([15, 20, 25, 15, 20, 25, 15, 20, 25]) * 1e6,  # Pa
    'S_exp': np.array([0.08, 0.15, 0.22, 0.12, 0.20, 0.30,
                       0.18, 0.28, 0.40])  # kg/m^3
}

co2 = PengRobinson(Tc=304.1, Pc=7.38e6, omega=0.225)
M_CO2 = 44.01e-3  # kg/mol

def co2_density(T, P):
    """Peng-Robinson式によるCO2の質量密度 (kg/m^3)。"""
    rho_molar, _ = co2.density(T, P, phase='vapor')
    return rho_molar * M_CO2

data['rho'] = np.array([co2_density(T, P)
                        for T, P in zip(data['T'], data['P'])])
print("Peng-Robinson式によるCO2密度 (kg/m^3):")
print(np.round(data['rho'], 1))

# 対数変換したモデルに対する線形最小二乗法:
# ln(S) = k * ln(rho) + a * (1/T) + b
ln_S_exp = np.log(data['S_exp'])
ln_rho = np.log(data['rho'])
T_inv = 1 / data['T']

X = np.column_stack([ln_rho, T_inv, np.ones(len(ln_rho))])
y = ln_S_exp

params = np.linalg.lstsq(X, y, rcond=None)[0]
k, a, b = params

print("超臨界CO2中のβ-カロテンに対するChrastilモデルのフィッティング")
print("=" * 60)
print("フィッティングされたパラメータ:")
print(f"  k（会合数） = {k:.3f}")
print(f"  a（熱項, K） = {a:.1f}")
print(f"  b（定数） = {b:.3f}")

# 予測値と適合度
S_pred = chrastil_solubility(data['rho'], data['T'], k, a, b)

SS_res = np.sum((S_pred - data['S_exp'])**2)
SS_tot = np.sum((data['S_exp'] - np.mean(data['S_exp']))**2)
R2 = 1 - SS_res / SS_tot
RMSE = np.sqrt(np.mean((S_pred - data['S_exp'])**2))

print("\n適合度:")
print(f"  R^2 = {R2:.4f}")
print(f"  RMSE = {RMSE:.4f} kg/m^3")

# フィッティング結果の可視化
plt.figure(figsize=(12, 5))

# パリティプロット
plt.subplot(1, 2, 1)
plt.scatter(data['S_exp'], S_pred, c=data['T'], cmap='coolwarm', s=100,
            edgecolor='black')
plt.plot([0, max(data['S_exp'])], [0, max(data['S_exp'])], 'k--',
         label='完全一致')
plt.xlabel('実験溶解度 (kg/m³)', fontsize=12)
plt.ylabel('予測溶解度 (kg/m³)', fontsize=12)
plt.title('パリティプロット: Chrastilモデル', fontsize=13, fontweight='bold')
plt.colorbar(label='温度 (K)')
plt.legend()
plt.grid(alpha=0.3)

# 各温度における溶解度の圧力依存性
plt.subplot(1, 2, 2)

for T_iso in [313, 333, 353]:
    P_smooth = np.linspace(15e6, 25e6, 50)
    rho_smooth = np.array([co2_density(T_iso, P) for P in P_smooth])
    S_smooth = chrastil_solubility(rho_smooth, T_iso, k, a, b)

    plt.plot(P_smooth / 1e6, S_smooth, linewidth=2, label=f'T = {T_iso} K（モデル）')

    mask = data['T'] == T_iso
    plt.scatter(data['P'][mask] / 1e6, data['S_exp'][mask], s=100,
                edgecolor='black')

plt.xlabel('圧力 (MPa)', fontsize=12)
plt.ylabel('溶解度 (kg/m³)', fontsize=12)
plt.title('超臨界CO₂中のβ-カロテンの溶解度', fontsize=13, fontweight='bold')
plt.legend(fontsize=10)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# 新しい条件での予測
T_new, P_new = 323, 22e6  # K, Pa
rho_new = co2_density(T_new, P_new)
S_new = chrastil_solubility(rho_new, T_new, k, a, b)

print("\n新しい条件における予測:")
print(f"  T = {T_new} K, P = {P_new/1e6:.1f} MPa")
print(f"  推定CO2密度 = {rho_new:.1f} kg/m^3")
print(f"  予測溶解度 = {S_new:.3f} kg/m^3")

例8: 二成分系の混合則

最後に、純成分の状態方程式を混合物へ拡張する混合則と、混合物パラメータが二成分相互作用パラメータ $k_{12}$ にどれだけ敏感かを見ます。

import numpy as np
import matplotlib.pyplot as plt

def mixing_rule(x1, a1, a2, b1, b2, k12=0):
    """
    van der Waals型の一流体混合則。

    Parameters:
    -----------
    x1 : float
        成分1のモル分率
    a1, a2 : float
        各純成分のaパラメータ
    b1, b2 : float
        各純成分のbパラメータ
    k12 : float
        二成分相互作用パラメータ

    Returns:
    --------
    a_mix, b_mix : float
        混合系のパラメータ
    """
    x2 = 1 - x1

    # クロスパラメータ
    a12 = np.sqrt(a1 * a2) * (1 - k12)

    # 混合則
    a_mix = x1**2 * a1 + 2*x1*x2*a12 + x2**2 * a2
    b_mix = x1 * b1 + x2 * b2

    return a_mix, b_mix

R = 8.314  # J/(mol·K)

# CO2
Tc1, Pc1 = 304.1, 7.38e6  # K, Pa
a1 = 27 * R**2 * Tc1**2 / (64 * Pc1)
b1 = R * Tc1 / (8 * Pc1)

# エタノール
Tc2, Pc2 = 513.9, 6.14e6  # K, Pa
a2 = 27 * R**2 * Tc2**2 / (64 * Pc2)
b2 = R * Tc2 / (8 * Pc2)

k12_values = [0.0, 0.05, 0.10, 0.15]
x_CO2 = np.linspace(0, 1, 100)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

for k12 in k12_values:
    a_mix_values = []
    b_mix_values = []

    for x1 in x_CO2:
        a_mix, b_mix = mixing_rule(x1, a1, a2, b1, b2, k12)
        a_mix_values.append(a_mix)
        b_mix_values.append(b_mix)

    ax1.plot(x_CO2, a_mix_values, linewidth=2, label=f'k₁₂ = {k12}')
    ax2.plot(x_CO2, np.array(b_mix_values)*1e6, linewidth=2,
             label=f'k₁₂ = {k12}')

ax1.set_xlabel('CO₂モル分率', fontsize=12)
ax1.set_ylabel('a_mix (Pa·m⁶/mol²)', fontsize=12)
ax1.set_title('混合系の引力パラメータ', fontsize=13)
ax1.legend()
ax1.grid(alpha=0.3)

ax2.set_xlabel('CO₂モル分率', fontsize=12)
ax2.set_ylabel('b_mix (cm³/mol)', fontsize=12)
ax2.set_title('混合系の排除体積パラメータ', fontsize=13)
ax2.legend()
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.show()

# 等モル組成における相互作用パラメータの影響
print("x_CO2 = 0.5 における k₁₂ の影響:")
for k12 in k12_values:
    a_mix, b_mix = mixing_rule(0.5, a1, a2, b1, b2, k12)
    print(f"  k12 = {k12:.2f}: a_mix = {a_mix:.4f} Pa·m^6/mol^2, "
          f"b_mix = {b_mix*1e6:.2f} cm^3/mol")

着目点：$b_m$ は組成に対して厳密に線形で、$k_{12}$ には全く影響されません。一方 $a_m$ は $k_{12}$ が大きくなるにつれて下に凸に垂れ下がります。補正を担うのは引力項だけです。だからこそペアごとに $k_{12}$ を1つフィッティングすれば二成分の気液平衡データを再現できることが多く、またエタノールのような共溶媒がモル分率に不釣り合いなほど溶媒力を変える理由にもなっています。

まとめ

実務への示唆

• 流体の種類と、その判断に本当に必要な精度に応じて状態方程式を選ぶ。
• 臨界点近傍での急激な物性変化に対してマージンを確保して設計する。
• 一般化相関と妥当性の確認には換算物性を用いる。
• 多成分の超臨界系で理想混合を仮定してはならない。