三次状態方程式

**三次状態方程式**（さんじじょうたいほうていしき、）は、温度と密度との関数として気体の圧力をモデル化する特定のクラスの熱力学モデルであり、モル体積の三次関数として書き換えることができる。

状態方程式は一般に物理化学や化学工学の分野で適用され、特に気液平衡のモデル化や化学工学のプロセス設計に用いられる。

## ファンデルワールスの状態方程式
ファンデルワールスの状態方程式は、以下のように表すことができる。

: \left(p + \frac{a}{V_\text{m}^2}\right)\left(V_\text{m} - b\right) = RT

ここで、Tは絶対温度、pは圧力、V_\text{m}はモル体積、Rは気体定数である。なお、V_\text{m} = V / nであり、ここでVは体積、またn=N/N_\text{A}であり、ここでnはモル数、Nは粒子数、N_\text{A}はアボガドロ定数である。これらの定義は、以下に示すすべての状態方程式に適用される。

1873年に提案されたファンデルワールスの状態方程式は、理想気体の状態方程式よりも大幅に精度の高い方程式の一つだった。この方程式では、通常aを引力パラメータ、bを斥力パラメータ（または有効分子体積）と呼ぶ。ファンデルワールスの状態方程式は、理想気体の状態方程式よりも優れており、液相の形成を予測できるものの、気液平衡に関する実験データとの一致性は限定的である。この方程式は歴史的な理由などから教科書や論文でよく引用されるが、その後、わずかに複雑な他の状態方程式が開発され、それらの多くはファンデルワールス方程式よりもはるかに正確であった。

ファンデルワールスの状態方程式は、理想気体の状態方程式に対し、「非理想的な寄与」を2つ加えて改良したものと考えることができる。この方程式を考察すると、次のようになる。

: p = \frac{RT}{V_\text{m}-b} - \frac{a}{V_\text{m}^2} 

ここで、理想気体の状態方程式と比較すると次のようになる。

: p = \frac{RT}{V_\text{m}} 

ファンデルワールスの状態方程式の形式は、以下のように動機づけられる。

# 分子は点粒子ではなく有限の体積を占めるため、全ての分子が同時に利用できる空間が制限される。これにより、点粒子を仮定した場合よりも圧力がわずかに上昇するため、第一項ではV_\text{m}の代わりに、「有効なモル体積」としてV_\text{m} - bを用いる。
# 理想気体分子は相互作用しないが、実際の分子は十分に接近するとファンデルワールス力による引力を及ぼす。これにより、分子が容器の壁に衝突する速度が遅くなり、結果として圧力が低下する。この影響を受ける衝突の数は密度\rhoに比例するため、圧力の低下量は\rho^2に比例し、すなわちモル体積の2乗に反比例することになる。

物質特有の定数aとbは、臨界特性（臨界圧力p_\text{c}と臨界モル体積V_\text{c}）を用いて次のように計算できる。

: a = 3 p_\text{c} V_\text{c}^2
: b = \frac{V_\text{c}}{3} 

ファンデルワールスの状態方程式の定数(a,b)は、臨界
温度および臨界圧力(T_\text{c},p_\text{c})の関数として表すこともでき、これらの実験的に容易に測定できる値を用いてパラメータ化されることが一般的である。

: a = \frac{27(R T_\text{c})^2}{64p_\text{c}} 
: b = \frac{R T_\text{c}}{8p_\text{c}} 

ファンデルワールスの状態方程式は、以下の減少状態変数（reduced state variables）を導入することで一般化できる。

V_\text{r}=V_\text{m}/V_\text{c}
P_\text{r}=p/p_\text{c}
T_\text{r}=T/T_\text{c}

これより、ファンデルワールス方程式の縮約形を次のように変形することができる。

: \left(P_\text{r} + \frac{3}{V_\text{r}^2}\right)\left(3V_\text{r} - 1\right) = 8T_\text{r}

この形式の利点は、与えられたT_\text{r}とP_\text{r}に対して、液相と気相の縮約体積を直接求めることができる点である。特に、縮約形の三次方程式をカルダノの公式を用いて解くことで、これらの値を求めることができる。

: V_\text{r}^3 - \left(\frac{1}{3} + \frac{8T_\text{r}}{3P_\text{r}}\right)V_\text{r}^2 + \frac{3V_\text{r}}{P_\text{r}} - \frac{1}{P_\text{r}} = 0

P_\text{r}およびT_\text{r}という条件下では、系は気液平衡にあり、縮約形の三次方程式は3つの解を持つ。このとき、最大の解は気相、最小の解は液体減容に対応する。気液平衡状態においては、モル体積に対する圧力の変化をより正確にモデル化するために、マクスウェル構成（Maxwell construction）が用いられることがある。

圧縮率因子Z=PV_\text{m}/RTは、非理想的な挙動を特徴づけるために用いられる。ファンデルワールスの状態方程式を縮約形で表すと、圧縮率因子も以下のように縮約形に変換できる。

: Z = \frac{V_\text{r}}{V_\text{r}-\frac{1}{3}} - \frac{9}{8 V_\text{r} T_\text{r}} 

特に臨界点では、圧縮率因子は次の値を取る。

 Z_\text{c} = 3/8 = 0.375 

## レドリッヒ・クオンの状態方程式
1949年に発表された、レドリッヒ・クオンの状態方程式は、ファンデルワールスの状態方程式に対する顕著な改良として考えられていた。その比較的単純な形から、現在でも注目を集めている。

ファンデルワールスの状態方程式よりも優れた特性を持つが、液相に対する精度が低いため、正確な気液平衡の計算には適さない。しかし、別の液相相関式と組み合わせることで、この目的に使用することは可能である。この方程式は以下のように表され、そのパラメータと臨界定数の関係式も示される。

: \begin{align}
  p &= \frac{R\,T}{V_\text{m} - b} - \frac{a}{\sqrt{T}\,V_\text{m}\left(V_\text{m} + b\right)} \\[3pt]
  a &= \frac{\Omega_a\,R^2 T_\text{c}^\frac{5}{2}}{p_\text{c}} \approx 0.42748\frac{R^2\,T_\text{c}^\frac{5}{2}}{P_\text{c}} \\[3pt]
  b &= \frac{\Omega_b\,R T_\text{c}}{P_\text{c}} \approx 0.08664\frac{R\,T_\text{c}}{p_\text{c}} \\[3pt]
  \Omega_a &= \left[9\left(2^{1/3}-1\right)\right]^{-1} \approx 0.42748 \\[3pt]
  \Omega_b &= \frac{2^{1/3}-1}{3} \approx 0.08664
\end{align}

また、レドリッヒ・クオンの圧縮率因子を表す、別の等価な形も存在する。

: Z=\frac{p V_\text{m}}{RT} = \frac{V_\text{m}}{V_\text{m} - b} - \frac{a}{R T^{3/2} \left(V_\text{m} + b\right)} 

この方程式は、縮約圧力（前節で定義）が縮約温度に対する比の約1/2以下の場合、気相の性質を計算するのに適している。

: P_\text{r} 

さらに、この方程式は対応状態の法則に一致する。縮約形で表現すると、すべての気体に対して同一の方程式が得られる。

: P_\text{r} = \frac{3 T_\text{r}}{V_\text{r} - b'} - \frac{1}{b' \sqrt{T_\text{r}} V_\text{r} \left(V_\text{r}+b'\right)} 

ここで、b'は以下のように定義される。

: b' = 2^{1/3}-1 \approx 0.25992

さらに、臨界点における圧縮因子は、すべての物質に対して同じ値をとる。

: Z_\text{c}=\frac{p_\text{c} V_\text{c}}{R T_\text{c}}=1/3 \approx 0.33333

これはZ_\text{c} = 3/8 = 0.375である、ファンデルワールスの状態方程式による臨界圧縮率因子の予測よりも改善されている。典型的な実験値は以下の通りである。

* 二酸化炭素（CO₂）: Z_\text{c} = 0.274

* 水（H₂O）: Z_\text{c} = 0.235

* 窒素（N₂）: Z_\text{c} = 0.29

## ソアベによるレドリッヒ・クオンの状態方程式の変形
レドリッヒ・クオンの状態方程式を修正した形がソアベ（Soave）によって以下のように提案された。

: p = \frac{R\,T}{V_\text{m}-b} - \frac{a \alpha}{V_\text{m}\left(V_\text{m}+b\right)}
: a = \frac{\Omega_a\,R^2 T_\text{c}^2}{P_\text{c}} = \frac{0.42748\,R^2 T_\text{c}^2}{P_\text{c}}
: b = \frac{\Omega_b\,R T_\text{c}}{P_\text{c}} = \frac{0.08664\,R T_\text{c}}{P_\text{c}}
: \alpha = \left(1 + \left(0.48508 + 1.55171\,\omega - 0.15613\,\omega^2\right) \left(1-T_\text{r}^{0.5}\right)\right)^2
: T_\text{r} = \frac{T}{T_\text{c}}
: \Omega_a = \left[9\left(2^{1/3}-1\right)\right]^{-1} \approx 0.42748
: \Omega_b = \frac{2^{1/3}-1}{3} \approx 0.08664

ここで、*ω*はその物質の偏心因子（acentric factor）である。上記の\alphaの式は、実際にはグラボスキー（Graboski）とダウバート（Daubert）によるものである。ソアベによる元の式は以下の通りである。

: \alpha = \left(1 + \left(0.480 + 1.574\,\omega - 0.176\,\omega^2\right) \left(1-T_\text{r}^{0.5}\right)\right)^2

水素の場合では、以下の通りとなる。

: \alpha = 1.202 \exp\left(-0.30288\,T_\text{r}\right) 

ここで、以下の縮約形の変数と臨界点での圧縮率因子を代入することで、

: \{p_\text{r}=p/P_\text{c}, T_\text{r}=T/T_\text{c}, V_\text{r}=V_\text{m}/V_\text{c}, Z_\text{c}=\frac{P_\text{c} V_\text{c}}{R T_\text{c}}\}

次の式が得られる。

: p_\text{r} P_\text{c} = \frac{R\,T_\text{r} T_\text{c}}{V_\text{r} V_\text{c}-b} - \frac{a \alpha\left(\omega, T_\text{r}\right)}{V_\text{r} V_\text{c}\left(V_\text{r} V_\text{c+}b\right)} 
 = \frac{R\,T_\text{r} T_\text{c}}{V_\text{r} V_\text{c}-\frac{\Omega_b\,R T_\text{c}}{P_\text{c}}} - \frac{\frac{\Omega_a\,R^2 T_\text{c}^2}{P_\text{c}} \alpha\left(\omega, T_\text{r}\right)}{V_\text{r} V_\text{c}\left(V_\text{r} V_\text{c}+\frac{\Omega_b\,R T_\text{c}}{P_\text{c}}\right)} =
: = \frac{R\,T_\text{r} T_\text{c}}{V_\text{c}\left(V_\text{r}-\frac{\Omega_b\,R T_\text{c}}{P_\text{c} V_\text{c}}\right)} - \frac{\frac{\Omega_a\,R^2 T_\text{c}^2}{P_\text{c}} \alpha\left(\omega, T_\text{r}\right)}{V_\text{r} V_\text{c}^2\left(V_\text{r}+\frac{\Omega_b\,R T_\text{c}}{P_\text{c} V_\text{c}}\right)}
 = \frac{R\,T_\text{r} T_\text{c}}{V_\text{c}\left(V_\text{r}-\frac{\Omega_b}{Z_\text{c}}\right)} - \frac{\frac{\Omega_a\,R^2 T_\text{c}^2}{P_\text{c}} \alpha\left(\omega, T_\text{r}\right)}{V_\text{r} V_\text{c}^2\left(V_\text{r}+\frac{\Omega_b}{Z_\text{c}}\right)}

これにより、以下の結果が導かれる。

: p_\text{r} = \frac{R\,T_\text{r} T_\text{c}}{P_\text{c} V_\text{c}\left(V_\text{r}-\frac{\Omega_b}{Z_\text{c}}\right)} - \frac{\frac{\Omega_a\,R^2 T_\text{c}^2}{P_\text{c}^2} \alpha\left(\omega, T_\text{r}\right)}{V_\text{r} V_\text{c}^2\left(V_\text{r}+\frac{\Omega_b}{Z_\text{c}}\right)} = 
 \frac{T_\text{r}}{Z_\text{c}\left(V_\text{r}-\frac{\Omega_b}{Z_\text{c}}\right)} - \frac{\frac{\Omega_a}{Z_\text{c}^2} \alpha\left(\omega, T_\text{r}\right)}{V_\text{r} \left(V_\text{r}+\frac{\Omega_b}{Z_\text{c}}\right)} 

したがって、ソアベ・レドリッヒ・クオンの状態方程式（SRK式）の縮約形は、物質の偏心因子*ω*とZ_\text{c}のみに依存する。これは、すべての物質に対して統一された縮約形を持つ対応状態の法則と一致するVdW方程式やRK方程式とは異なる。

: p_\text{r} = \frac{T_\text{r}}{Z_\text{c}\left(V_\text{r}-\frac{\Omega_b}{Z_\text{c}}\right)} - \frac{\frac{\Omega_a}{Z_\text{c}^2} \alpha\left(\omega, T_\text{r}\right)}{V_\text{r} \left(V_\text{r}+\frac{\Omega_b}{Z_\text{c}}\right)} 

また、次のように多項式の形で表すこともできる。

: A = \frac{a \alpha P}{R^2 T^2}
: B = \frac{bP}{RT}

圧縮率因子で表すと、次の式が得られる。

: 0 = Z^3-Z^2+Z\left(A-B-B^2\right) - AB

この方程式は最大3つの解を持つ可能性があり、一般に最大の解は蒸気状態、最小の解は液体状態に対応する。したがって、気液平衡の計算などで三次方程式を使用する際には注意が必要である。

1972年、G.ソアベは、レドリッヒ・クオンの状態方程式の\frac{1}{\sqrt{T}}の項を、温度と偏心因子を含む関数α*α*(*T*,*ω*)に置き換えた。この修正された方程式はソアベ・レドリッヒ・クオンの状態方程式と呼ばれる。この関数*α*は、炭化水素の蒸気圧データに適合するよう設計されており、炭化水素に対しては比較的高い精度を持つ。

特に、この置き換えにより、パラメータ*a*の定義が若干変更される点に注意が必要である。これは、臨界温度T_\text{c}が二乗の形で関係するためである。

## ペネルーらによる体積補正 (1982)
ソアベ・レドリッヒ・クオンの状態方程式は次のように表すことができる。

: p = \frac{R\,T}{V_{m,\text{SRK}} - b} - \frac{a}{V_{m,\text{SRK}} \left(V_{m,\text{SRK}} + b\right)}

ここで、

: \begin{align}
    a &= a_\text{c}\, \alpha \\
  a_\text{c} &\approx 0.42747\frac{R^2\,T_\text{c}^2}{P_\text{c}} \\
    b &\approx 0.08664\frac{R\,T_\text{c}}{P_\text{c}}
\end{align}

である。また、\alphaおよびその他の要素については、ソアベ・レドリッヒ・クオンの状態方程式の節で定義される。

ソアベ・レドリッヒ・クオンの状態方程式やその他の三次体状態方程式の欠点の一つは、液体のモル体積が気体のモル体積と比べて著しく精度が低いことである。この問題に対し、ペネルー（Peneloux）らは、体積補正（volume translation） を導入することで単純な修正を提案した。

: V_{\text{m},\text{SRK}} = V_\text{m} + c

ここで、cは追加の流体成分パラメータであり、モル体積をわずかに変換する役割を持つ。液体状態では、モル体積の小さな変化が圧力の大きな変化に対応する。気体状態では、モル体積の小さな変化が圧力のごくわずかな変化に対応する。そのため、気体のモル体積に対する補正の影響は小さい。しかしながら、この体積補正の方法には科学界や産業界で2つのバージョンが存在する。

第1のバージョンでは、V_{\text{m},\text{SRK}}のみが補正され、状態方程式は次のようになる。

: p = \frac{R\,T}{V_\text{m} + c - b} - \frac{a}{\left(V_\text{m} + c\right) \left(V_\text{m} + c + b\right)}

第2のバージョンでは、V_{\text{m},\text{SRK}}とb_\text{SRK}の両方が補正されるか、V_{\text{m},\text{SRK}}の補正の後に合成パラメータの名称変更が行われる。これにより、次のような形になる。

: \begin{align}
  b_\text{SRK} &= b + c \quad \text{or} \quad b - c \curvearrowright b \\
             p &= \frac{R\,T}{V_\text{m} - b} - \frac{a}{\left(V_\text{m} + c\right) \left(V_\text{m} + 2c + b\right)}
\end{align}

流体混合物のパラメータ*c*は次の式で計算される。

: c = \sum_{i=1}^n z_i c_i

石油ガスやオイル中の各流体成分のパラメータ*c*は、以下の相関式を用いて推定できる。

: c_i \approx 0.40768\ \frac{RT_{ci}}{P_{ci}} \left(0.29441 - Z_{\text{RA},i}\right) 

ここで、ランケット圧縮率因子 Z_{\text{RA},i} は、次の方法で推定できる。

: Z_{\text{RA},i} \approx 0.29056 - 0.08775\ \omega_i

ペネルーらの体積補正法の優れた点は、気液平衡計算に影響を与えないことである。この体積補正法は、他の三次状態方程式にも適用可能であり、その場合は、選択した状態方程式に適合するようにパラメータ*c*の相関式を調整すればよい。

## ペン・ロビンソンの状態方程式
ペン・ロビンソンの状態方程式（PR式） は、1976年にアルバータ大学で（Ding-Yu Peng）とドナルド・ロビンソン（Donald Robinson）により開発され、以下の目標を満たすように設計された。

# パラメータは臨界特性と偏心因子の関数として表現できること。
# 臨界点付近で適切な精度を提供し、特に圧縮率因子や液密度の計算に対して正確であること。
# 混合則（mixing rules）は、温度・圧力・組成に依存しない単一の二成分相互作用パラメータのみを使用すること。
# 天然ガスプロセスにおける全ての流体特性の計算に適用可能であること。

この状態方程式は、次の形で表される。

p = \frac{R\,T}{V_\text{m} - b} - \frac{a\,\alpha}{V_\text{m}^2 + 2bV_\text{m} - b^2}

a = \Omega_a \frac{R^2\,T_\text{c}^2}{p_\text{c}};  \Omega_a = \frac{8+40\eta_c}{49-37\eta_c} \approx 0.45724 


b = \Omega_b \frac{R\,T_\text{c}}{p_\text{c}}; \Omega_b = \frac{\eta_c}{3+\eta_c} \approx 0.07780 

\eta_c = [1+(4-\sqrt{8})^{1/3}+(4+\sqrt{8})^{1/3}]^{-1}

\alpha = \left(1 + \kappa \left(1 - \sqrt{T_\text{r}}\right)\right)^2; T_\text{r} = \frac{T}{T_\text{c}}

\kappa \approx 0.37464 + 1.54226\omega - 0.26992\omega^2

また、多項式形式は次のようになる。

: A = \frac{\alpha a p}{R^2\,T^2}
: B = \frac{bp}{RT}
: Z^3 - (1 - B)Z^2 + \left(A - 2B - 3B^2\right)Z - \left(AB - B^2 - B^3\right) = 0

概ね、ペン・ロビンソンの状態方程式はソアベの状態方程式と同等の性能を示すが、特に無極性物質の液密度の予測においては一般的に優れている。ペン・ロビンソンの状態方程式の詳細な性能については、密度、熱特性、相平衡に関する報告がなされている。簡単に言えば、オリジナルの形式では蒸気圧や相平衡の誤差が改良版の約3倍であることが分かっている。また、ペン・ロビンソンの状態方程式の偏倚関数（Departure functions）については、別の記事で詳細に説明されている。

特性定数の解析値は以下の通りとなっている。

: Z_\text{c} = \frac{1}{32} \left( 11 - 2\sqrt{7} \sinh\left(\frac{1}{3} \operatorname{arsinh}\left(\frac{13}{7 \sqrt{7}}\right)\right) \right) \approx 0.307401
: b' = \frac{b}{V_{\text{m},\text{c}}} = \frac{1}{3} \left( \sqrt{8} \sinh\left(\frac{1}{3} \operatorname{arsinh}\left(\sqrt{8}\right) \right) - 1 \right) \approx 0.253077 \approx \frac{0.07780}{Z_\text{c}} 
:  \frac{P_\text{c} V_{\text{m},\text{c}}^2}{a\,b'} = \frac{3}{8} \left( 1 + \cosh\left(\frac{1}{3} \operatorname{arcosh}(3) \right) \right) \approx 0.816619 \approx \frac{Z_\text{c}^2}{0.45724 \, b'} 

## ペン・ロビンソン・ストリエク・ベラの状態方程式
### PRSV1
1986年、ストリエク（Stryjek）とベラ（Vera）によるペン・ロビンソンの状態方程式の引力項の修正は、調整可能な純成分パラメータの導入と偏心因子の多項式フィットの修正によって、モデルの精度を大幅に向上させた。

修正式は以下の通りである。

: \begin{align}
    \kappa &= \kappa_0 + \kappa_1 \left(1 + T_\text{r}^\frac{1}{2}\right) \left(0.7 - T_\text{r}\right) \\
  \kappa_0 &= 0.378893+1.4897153\,\omega - 0.17131848\,\omega^2 + 0.0196554\,\omega^3
\end{align}

ここで、\kappa_1は調整可能な純成分パラメータである。ストリエクとベラは、工業的に重要な多くの化合物の純成分パラメータを論文で発表している。

還元温度が0.7を超える場合、彼らは\kappa_1 = 0 とし、\kappa = \kappa_0 をそのまま使用することを推奨している。アルコールや水の場合、 \kappa_1 は臨界温度まで使用可能であり、それ以上の温度では0に設定する。

### PRSV2
1986年に発表された別の修正式では、引力項の修正に2つの純成分パラメータを追加することで、モデルの精度がさらに向上した。

その修正式は以下の通りである。

: \begin{align}
    \kappa &= \kappa_0 + \left[\kappa_1 + \kappa_2\left(\kappa_3 - T_\text{r}\right)\left(1 - T_\text{r}^\frac{1}{2}\right)\right]\left(1 + T_\text{r}^\frac{1}{2}\right) \left(0.7 - T_\text{r}\right) \\
  \kappa_0 &= 0.378893 + 1.4897153\,\omega - 0.17131848\,\omega^2 + 0.0196554\,\omega^3
\end{align}

ここで、\kappa_1、\kappa_2、\kappa_3はすべて調整可能な純成分パラメータである。

PRSV2は、特に気液平衡の計算に有利である。PRSV1は、ペン・ロビンソンモデル に比べ、熱力学的挙動をより正確に記述できるが、相平衡計算には十分な精度を持たない。相平衡計算の高度に非線形な挙動は、通常なら許容範囲の小さな誤差を大きく増幅させる傾向があるため、設計に適用する際はPRSV2の使用が推奨される。ただし、一度平衡状態が決まれば、平衡状態における相ごとの熱力学的特性値は、より単純なモデルを用いて十分な精度で求めることができる。

ペン・ロビンソン・ストリエク・ベラの状態方程式では、パラメータのフィットは特定の温度範囲（通常は臨界温度以下）で行われる。そのため、臨界温度を超えるとPRSVのアルファ関数は発散し、任意に大きくなる傾向がある。これは特に、水素を含む系（しばしば臨界点を大きく超えた温度で存在する）において問題となる。この問題を解決するため、臨界点以上では異なる定式化もいくつか提案されている。 代表的な修正モデルとして、トゥ（Twu）らの方法や、マティアス（Mathias）とコープマン（Copeman）の方法が挙げられる。また、ジョベール（Jaubert）らによりトゥの方法を用いた1,700種類以上の化合物に関する詳細な研究が報告されている。ジョベールらにより更新されたペン・ロビンソンの状態方程式の詳細な性能については、密度、熱特性、相平衡が報告されている。簡単に言うと、蒸気圧や相平衡の誤差が、オリジナルの実装の約1/3に減少している。

## ペン・ロビンソン・ババロラ・ススの状態方程式
ババロラ（Babalola）とスス（Susu）によるペン・ロビンソンの状態方程式の修正は次の通りである。

P =\left ( \frac{RT}{V_\mathrm{m}-b} \right ) -\left [ \frac{(a_1P+a_2)\alpha}{V_\mathrm{m}(V_\mathrm{m}+b)+b(V_\mathrm{m}-b)} \right ]

ペン・ロビンソン状態方程式では、引力項*a*は圧力に対して一定とされていた。ババロラとススは、多成分・多相の高密度貯留層システムにおける物性予測の精度向上のために、*a*を圧力に依存する変数して扱うとする修正を加えた。この変化は線形方程式で表され、パラメータ*a*の値を圧力に対してプロットしたときの直線の傾き*a*1と切片*a*2で表現される。

この修正により、特に30 MPa以上の高圧条件下における重質流体の精度が向上し、従来のペン・ロビンソンの状態方程式のチューニングが不要となった。チューニングは、この修正の過程で本質的に組み込まれている。

ペン・ロビンソン・ババロラ・ススの状態方程式は、2005年に開発され、それ以来約20年にわたり、さまざまな温度・圧力条件下で貯留層データに適用されてきた。その結果、特に超深部の高温高圧条件下における貯留層流体の物性予測において、高精度な状態方程式の一つとして評価されている。これらの研究成果は、権威ある学術誌に掲載されている。

1976年に発表された従来のペン・ロビンソンの状態方程式は、27 MPa（約4,000 psi）までの通常の貯留層流体の物性を高精度で予測できるが、それ以上の圧力では誤差が増大する。一方、新しいペン・ロビンソン・ババロラ・ススの状態方程式は、120MPa（約17,500 psi）という極めて高い圧力範囲においても、PVT挙動を正確にモデル化できる。

## エリオット・スレーシュ・ドナヒューの状態方程式
エリオット・スレーシュ・ドナヒューの状態方程式（ESD式）は、1990年に提案された。この方程式は、ペン・ロビンソンの状態方程式にも適用されているファンデルワールスの反発項の不正確さを修正している。引力項には、スクエアウェル球の第二ビリアル係数に関連する寄与が含まれており、また、トゥの温度依存性の特徴も一部取り入れられている。この状態方程式は、分子の形状の影響を考慮しており、臨界特性を用いる代わりに溶解度パラメータと液体体積によって特徴付けられる分子パラメータを用いることで、直接ポリマーへ拡張可能である。エリオット・スレーシュ・ドナヒューの状態方程式自体は、コンピュータシミュレーションとの比較を通じて開発されており、サイズ・形状・水素結合の基本的な物理的特性を捉えるよう設計されている。特に、直鎖分子（n-アルカンなど）を基にした推論に基づいて構築されている。

: \frac{p V_\text{m}}{RT}=Z=1 + Z^{\rm{rep}} + Z^{\rm{att}}

ここで、

: Z^{\rm{rep}} = \frac{4 c \eta}{1-1.9 \eta}

: Z^{\rm{att}} = -\frac{z_\text{m} q \eta Y}{1+ k_1 \eta Y}

である。この方程式では、形状因子cを用いて分子の形状の影響を表現している。球状分子の場合は、c=1である。

また、非球状分子の場合、形状因子cは、不規則性因子と以下の関係式で結び付けられる。

: c=1+3.535\omega+0.533\omega^2

縮約数密度\etaは、以下の式で定義される。

\eta=b \rho

ここで、bは分子の特徴的なサイズパラメータ[cm3/mol]であり、\rho = \frac{1}{V_\text{m}}= N/(N_\text{A}V)はモル密度 [mol/cm3]である。これらは、形状因子cと次の関係を持つ。

: b=\frac{RT_\text{c}}{P_\text{c}}\Phi

ここで、

: \Phi=\frac{Z_\text{c}^2}{2A_q}{[-B_q+\sqrt{B_q^2+4A_qC_q }] }
: 3Z_\text{c}=([(-0.173/\sqrt{c}+0.217)/\sqrt{c}-0.186]/\sqrt{c}+0.115)/\sqrt{c}+1

: A_q=[1.9(9.5q-k_1)+4ck_1](4c-1.9)

: B_q=1.9k_1Z_\text{c}+3A_q/(4c-1.9)

: C_q=(9.5q-k_1)/Z_\text{c}


となる。引力項に現れる形状パラメータqおよび項Yは、以下の式で与えられる。

: q=1+k_3(c-1)
: Y=\exp\left(\frac{\epsilon}{kT}\right) - k_2

このため、球状分子の場合はq=1となる。また、
項\epsilonは井戸型ポテンシャルの深さで、以下の式で定義される。

: Y_\text{c} =(\frac{R T_\text{c}}{b P_\text{c}})^2 \frac{Z_\text{c}^3}{A_q}

ここで、エリオット・スレーシュ・ドナヒューの状態方程式には、以下の定数が含まれる。

: z_\text{m}、k_1、k_2、k_3

: z_\text{m} = 9.5、 k_1 = 1.7745、 k_2 = 1.0617、 k_3 = 1.90476

このモデルは、会合性成分（水やアルコールなど）および、非会合性成分を含む混合系にも拡張可能である。詳細については、J.R. エリオット Jr. らの論文（1990年）に記載されている。

4(k_3-1)/k_3 = 1.900であることに注目すると、Z^\text{rep}はSAFTの形式で次のように書き換えられる。

:Z^{\rm{rep}} = 4q \eta g-(q-1)\frac{\eta}{g}\frac{dg}{d\eta}= \frac{4q\eta}{1-1.9\eta}-\frac{(q-1)1.9\eta}{1-1.9\eta};g=\frac{1}{1-1.9\eta}

好みに応じて、qをSAFT記法のmに置き換えることで、エリオット・スレーシュ・ドナヒューの状態方程式は以下のように記述できる。

:Z=1 + m(\frac{4\eta}{1-1.9\eta} - \frac{9.5Y\eta}{1+k_1Y\eta})-\frac{(m-1)1.9\eta}{1-1.9\eta}


この形では SAFTのセグメント的視点が明確になり、マイケル・ヴェルトハイム（Michael Wertheim）の理論が直接適用可能となる。SAFTにおいては、各分子は*m*個の球状セグメントから構成され、それぞれが独立して球状相互作用を持つと考えられる。しかし、それらのセグメントは連結されているため、修正項(*m* − 1)を加えて補正する。*m*が整数でない場合、有効な接触球セグメントの数として解釈される。

ヴェルトハイムの理論を解くのは複雑だが、単純化することで実装が容易になる。密度と温度を入力としてZ^{\rm{assoc}}を求めるには、いくつかの追加ステップが必要である。例えば、水素結合の供与体数と受容体数が等しい場合、エリオット・スレーシュ・ドナヒューの状態方程式は以下のように簡略化される。

:\frac{p V_\text{m}}{RT}=Z=1 + Z^{\rm{rep}} + Z^{\rm{att}}+ Z^{\rm{assoc}}

ここで、

:Z^{\rm{assoc}} = -gN^\text{AD}(1-X^\text{AD});X^\text{AD}=2/[1+\sqrt{1+4N^\text{AD}\alpha^\text{AD}}];\alpha^\text{AD}=\rho N_\text{A}K^\text{AD}[\exp{(\epsilon^\text{AD}/kT)-1]}

である。また、N_\text{A}はアボガドロ定数、K^\text{AD}および\epsilon^\text{AD}は、水素結合の体積K^\text{AD} = \mathrm{0.001\ nm^3}、およびエネルギーのパラメータ（保存された入力値）\epsilon^\text{AD}/k_\text{B}=\mathrm{2000\ K}が格納されている。N^\text{AD}は、水素結合の受容体数（供与体数と等しい）である。例えば、メタノールやエタノールなどのアルコール類の場合は、N^\text{AD} = 1であり、水の場合はN^\text{AD} = 2となる。の場合は、N^\text{AD} = 重合度となる。計算手順としては、密度と温度を用いて\alpha^\text{AD}を計算し、\alpha^\text{AD}を利用して他の物理量を求める。技術的には、会合項を含めると エリオット・スレーシュ・ドナヒューの状態方程式はもはや三次方程式ではなくなるが、不自然な挙動（アーティファクト）は生じず、密度に関しては依然として3つの解のみが存在する。また、自己会合（Self-associating）、相互会合（Cross-associating）、非会合（Non-associating）化合物の混合系を含め、任意の数の電子受容体（酸）および電子供与体（塩基）を効率的に扱うための拡張が示されている。エリオット・スレーシュ・ドナヒューの状態方程式の密度・熱物性・相平衡の詳細な評価が報告されている。簡単にまとめると、ESD方程式は蒸気圧と気液平衡において、ジョベールらが改良したペン・ロビンソンの状態方程式のおよそ2倍の偏差を示すが、液-液平衡における偏差はおよそ40％低減した。

## キュービック・プラス会合
キュービック・プラス会合（CPA）方程式は、ソアベ・レドリッヒ・クオンの状態方程式に自己会合を考慮する項を組み込んだ状態方程式である。この会合項は、マイケル・ヴェルトハイムによる会合分子の理論を基に、チャップマンによって拡張・簡略化されたの概念を取り入れている。この方程式の開発は、シェルの資金提供を受けた研究プロジェクトとして1995年に始まり、1996年に公開された。

: p = \frac{RT}{(V_\mathrm{m} - b)} - \frac{a}{V_\mathrm{m} (V_\mathrm{m} + b)} + \frac{RT}{V_\mathrm{m}} \rho \sum_{A} \left[ \frac{1}{X^\text{A}} - \frac{1}{2} \right] \frac{\partial X^\text{A}}{\partial \rho}

会合項において、X^\text{A}は部位Aで結合していない分子のモル分率を表す。

## キュービック・プラス鎖状態方程式
キュービック・プラス鎖（CPC）状態方程式は、古典的な三次状態方程式との鎖項を組み合わせたハイブリッド状態方程式である。

鎖項を追加することで、短鎖および長鎖の非会合成分（アルカンから高分子まで）の物理特性を適切に捉えることが可能となる。また、CPCの単量体（モノマー）項は特定の立方体型状態方程式に限定されず、さまざまな形態を同じ枠組みで適用できる。CPC 状態方程式は、縮約残余ヘルムホルツの自由エネルギー（F^\mathrm{CPC}）の形式で表される。

F^\mathrm{CPC} = \frac{A^\mathrm{R} (T,V,{\textbf{n}})}{RT}= m(F^\mathrm{rep}+F^\mathrm{att})+F^\mathrm{chain} 

ここで、A^\mathrm{R}は残余ヘルムホルツの自由エネルギー、mは鎖長、"rep"および"att"は、それぞれ三次状態方程式の反発項と引力項、"chain"項は、SAFTに基づくモノマービーズの結合寄与を示す。レドリッヒ・クオンの状態方程式をモノマー項として使用した場合、CPCは次のように表される。

p = \frac{nRT}{V} \Biggl(1+\frac{\bar{m}^2B}{V -\bar{m}B}\Biggr) - \frac{\bar{m}^2A}{V(V +\bar{m}B )} - \frac{ nRT }{V} \left[\sum_{i}  n_i(m_i-1)\beta \frac{g'(\beta)}{g(\beta) }  \right] 

ここで、Aは分子間相互作用エネルギーパラメータ、Bは共体積パラメータ、 \bar{m} はモル平均鎖長、*g(β)*は接触時の動径分布関数（RDF）、*β*は縮約体積を表す。

CPCモデルは、従来の高分子モデルに比べてシンプルかつ高速に動作する。シスコ（Sisco）らは、キュービック・プラス鎖状態方程式を用いてさまざまな高分子混合系をモデル化し、高圧・高温条件、溶媒種類、分子量分布（ポリ分散性） などの影響を分析した。CPCモデルは、実験データと比較した際の適合性が高いことが示された。

また、アラジミ（Alajmi）らは短距離のソフト反発（short-range soft repulsion）を、CPCの枠組みに組み込み、蒸気圧や液体密度の予測精度を向上 させた。彼らは、n-アルカン、アルケン、分岐アルカン、シクロアルカン、ベンゼン誘導体、気体など50種類以上の成分を含むデータベースを提供した。このCPCバージョンでは、摂動（perturbation theory）に基づいた温度依存の共体積パラメータを導入し、分子間の短距離ソフト反発を考慮している。

## 脚注