6. Cubic Equation of State¶

6.1. Theory¶

The generic cubic equation of state is [GP07]

\[P = \frac{RT}{V - b} - \frac{a(T)}{(V + \epsilon b)(V + \sigma b)}\]

where \(\epsilon\) and \(\sigma\) are pure numbers (the same for all substances), and \(a(T)\) and \(b\) are given by the following equations

(1)¶\[a(T) = \Psi \frac{\alpha(T_r)R^2T_{\text{c}}^2}{P_{\text{c}}}\]

\[b = \Omega\frac{RT_\text{c}}{P_\text{c}}\]

The compressibility factor can be calculated by solving the following equation

(2)¶\[Z - \left(1 + \beta - \frac{q\beta (Z - \beta)}{(Z + \epsilon\beta)(Z + \sigma\beta)}\right) = 0\]

where

(3)¶\[\beta = \Omega\frac{P_\text{r}}{T_\text{r}}\]

and

(4)¶\[q = \frac{\Psi \alpha(T_\text{r})}{\Omega T_\text{r}}\]

An iterative routine to calculate \(Z\) using Equation (2) is implemented, following [GP07], where

(5)¶\[Z^\text{new} = 1 + \beta - \frac{q\beta (Z^\text{old} - \beta)}{(Z^\text{old} + \epsilon\beta)(Z^\text{old} + \sigma\beta)}\]

that is continued until the following is true

(6)¶\[\|\frac{Z^\text{new} - Z^\text{old}}{Z^\text{new} + Z^\text{old}}\times 200\| < \text{tol}\]

6.2. Residual Molar Properties¶

(7)¶\[\frac{H^\text{R}}{RT} = Z - 1 + \left[\frac{\text{d} \ln \alpha(T_\text{r})}{\text{d} \ln T_\text{r}} - 1\right]qI\]

(8)¶\[\frac{G^\text{R}}{RT} = Z - 1 + \ln(Z-\beta) - qI\]

(9)¶\[\frac{S^\text{R}}{R} = \ln(Z-\beta) + \frac{\text{d} \ln \alpha(T_\text{r})}{\text{d} \ln T_\text{r}} qI\]

where the following functions have been defined

(10)¶\[I = \frac{1}{\sigma - \epsilon}\ln{\left( \frac{Z + \sigma \beta}{Z + \epsilon \beta} \right)}\]

(11)¶\[1 + \beta - \frac{q\beta (Z - \beta)}{(Z + \epsilon\beta)(Z + \sigma\beta)}\]

6.3. Fugacity Coefficients¶

The pure-component fugacity coefficient is defined as

\[\frac{G_i^\text{R}}{RT} = \ln\hat{\phi}_i\]

So that, from Equation (8),

(12)¶\[\ln\hat{\phi}_i = Z_i - 1 + \ln(Z_i-\beta_i) - q_iI_i\]

6.4. Mixtures¶

Warning

Not implemented yet!

Todo

Implement mixtures with cubic equations of state

class gasthermo.eos.cubic.Cubic(sigma: float, epsilon: float, Omega: float, Psi: float, dippr_no: str = None, compound_name: str = None, cas_number: str = None, log: callable = <ufunc 'log'>, exp: callable = <ufunc 'exp'>, name: str = 'cubic', **kwargs)[source]¶

Generic Cubic Equation of State

Parameters

sigma – Parameter defined by specific equation of state, \(\sigma\)
epsilon – Parameter defined by specific equation of state, \(\epsilon\)
Omega – Parameter defined by specific equation of state, \(\Omega\)
Psi – Parameter defined by specific equation of state, \(\Psi\)
tol – tolerance for iteration (see Equation (6)), set to 0.01
log – function for computing natural log, defaults to numpy.log
exp – function for computing exponential, defaults to numpy.exp

G_R_RT_expr(P, V, T)[source]¶

Dimensionless residual gibbs

Parameters

P – pressure in Pa
V – molar volume in m**3/mol
T – temperature in K

Returns

\(\frac{G^\text{R}}{RT}\), see Equation (8)

H_R_RT_expr(P, V, T)[source]¶

Dimensionless residual enthalpy

Parameters

P – pressure in Pa
V – molar volume in m**3/mol
T – temperature in K

Returns

\(\frac{H^\text{R}}{RT}\), see Equation (7)

I_expr(P, V, T)[source]¶

Parameters

P – pressure in Pa
V – molar volume in m**3/mol
T – temperature in K

Returns

\(I\) (see Equation (10))

S_R_R_expr(P, V, T)[source]¶

Dimensionless residual entropy

Parameters

P – pressure in Pa
V – molar volume in m**3/mol
T – temperature in K

Returns

\(\frac{S^\text{R}}{R}\), see Equation (9)

Z_right_hand_side(Z, beta, q)[source]¶

Estimate of compressibility of vapor [GP07], used for iterative methods

Returns: RHS of residual, see Equation (11)

a_expr(T)[source]¶

Parameters: T – temperautre in K
Returns: \(a(T)\) (see Equation (1))

alpha_expr(T_r)[source]¶

An empirical expression, specific to a particular form of the equation of state

Parameters: T_r – reduced temperature (T/Tc), dimensionless
Returns: \(\alpha (T_r)\) see Table ??

beta_expr(T, P)[source]¶

Parameters

T – temperature in K
P – pressure in Pa

Returns

\(\beta\) (see Equation (3))

cardano_constants(T, P)[source]¶

Parameters

T – temperature [T]
P – pressure [Pa]

Returns

cardano constants p, q

Return type

tuple

coefficients(T, P)[source]¶

Polynomial oefficients for cubic equation of state

\[Z^3 c_0 + Z^2c_1 + Z*c_2 + c_3 = 0\]

Returns: (c_0, c_1, c_2, c_3)

d_ln_alpha_d_ln_Tr(T_r)[source]¶

Parameters: T_r – reduced temperature [dimensionless]
Returns: Expression for \(\frac{\mathrm{d} \ln \alpha(T_\mathrm{r})}{\mathrm{d} \ln T_\mathrm{r}}\)

hat_phi_i_expr(*args)[source]¶: expression for fugacity coefficient :returns: \(\exp\left({\ln{\hat{\phi}_i}}\right)\)

iterate_to_solve_Z(T, P) → float[source]¶

Iterate to compute \(Z\) using Equation (5) util termination condition (Equation (6) is met)

Parameters

T – temperature in K
P – pressure in Pa

Returns

compressibility factor

ln_hat_phi_k_expr(P, V, T)[source]¶

Parameters

P – pressure in Pa
T – temperature in K
V – molar volume in m**3/mol

Returns

\(\ln{\hat{\phi}_k\), see Equation \(ln_hat_phi_i\)

num_roots(T, P)[source]¶

Find number of roots

See [ML12][Dei02]

Parameters

T – temperature in K
P – pressure in Pa

Returns

number of roots

plot_Z_vs_P(T: float, P_min: float, P_max: float, symbol='o', ax: matplotlib.pyplot.axis = None, fig: matplotlib.pyplot.figure = None, **kwargs)[source]¶

Plot compressibility as a function of pressure

Parameters

T – temperature [K]
P_min – minimum pressure for plotting [Pa]
P_max – maximum pressure for plotting [Pa]
symbol – marker symbol, defaults to ‘o’
ax – matplotlib axes for plotting, defaults to None
kwargs – keyword arguments for plotting

print_roots(T, P)[source]¶: Check to see if all conditions have one root

q_expr(T)[source]¶

Parameters: T – temperautre in K
Returns: \(q\) (see Equation (4))

residual(P, V, T)[source]¶

Parameters

P – pressure in Pa
V – volume in [mol/m**3]
T – temperature in K

Returns

residual for cubic equation of state (Equation (2))

class gasthermo.eos.cubic.RedlichKwong(**kwargs)[source]¶

Redlich-Kwong Equation of State [RK49]

This Equation of state has the following parameters [SVanNessA05]

Symbol	Value
\(\alpha(T_\text{r})\)	\(1/\sqrt{T_\text{r}}\)
\(\sigma\)	1
\(\epsilon\)	0
\(\Omega\)	0.08664
\(\Psi\)	0.42748

>>> from gasthermo.eos.cubic import RedlichKwong
>>> model = RedlichKwong(compound_name='Propane')
>>> model.sigma
1
>>> model.epsilon
0
>>> model.Omega
0.08664
>>> model.Psi
0.42748

alpha_expr(T_r)[source]¶

An empirical expression, specific to a particular form of the equation of state

Parameters: T_r – reduced temperature (T/Tc), dimensionless
Returns: \(\alpha (T_r)\) see Table ??

d_ln_alpha_d_ln_Tr(T_r)[source]¶

Parameters: T_r – reduced temperature [dimensionless]
Returns: Expression for \(\frac{\mathrm{d} \ln \alpha(T_\mathrm{r})}{\mathrm{d} \ln T_\mathrm{r}}\)

class gasthermo.eos.cubic.SoaveRedlichKwong(**kwargs)[source]¶

Soave Redlich-Kwong Equation of State [Soa72]

This equation of state has the following parameters

Symbol	Value
\(\alpha(T_\text{r})\)	\(\left[1 + f_w(1-\sqrt{T_\text{r}})\right]^2\)
\(\sigma\)	\(1\)
\(\epsilon\)	\(0\)
\(\Omega\)	\(0.08664\)
\(\Psi\)	\(0.42748\)

where

(13)¶\[f_w = 0.480 + 1.574\omega - 0.176\omega^2\]

>>> from gasthermo.eos.cubic import SoaveRedlichKwong
>>> model = SoaveRedlichKwong(compound_name='Water')
>>> model.Omega
0.08664
>>> model.sigma
1
>>> model.epsilon
0
>>> model.Psi
0.42748
>>> model.f_w_rule(0.)
0.48
>>> model.f_w_rule(1.)
1.878

Parameters: f_w (float, derived) – empirical expression used in \(\\alpha\) [dimensionless], see Equation (13)

alpha_expr(T_r)[source]¶

An empirical expression, specific to a particular form of the equation of state

Parameters: T_r – reduced temperature (T/Tc), dimensionless
Returns: \(\alpha (T_r)\) see Table ??

d_ln_alpha_d_ln_Tr(T_r)[source]¶

Parameters: T_r – reduced temperature [dimensionless]
Returns: Expression for \(\frac{\mathrm{d} \ln \alpha(T_\mathrm{r})}{\mathrm{d} \ln T_\mathrm{r}}\)

f_w_rule(w)[source]¶

Parameters: w – accentric factor
Returns: \(f_w\), see Equation (13)

class gasthermo.eos.cubic.PengRobinson(**kwargs)[source]¶

Peng-Robinson Equation of State [PR76]

This equation of state has the following parameters

Symbol	Value
\(\alpha(T_\text{r})\)	\(\left[1 + f_w(1-\sqrt{T_\text{r}})\right]^2\)
\(\sigma\)	\(1 + \sqrt{2}\)
\(\epsilon\)	\(1 - \sqrt{2}\)
\(\Omega\)	\(0.07780\)
\(\Psi\)	\(0.45724\)

where

(14)¶\[f_w = 0.480 + 1.574\omega - 0.176\omega^2\]

>>> from gasthermo.eos.cubic import PengRobinson
>>> model = PengRobinson(compound_name='Water')
>>> model.Omega
0.0778
>>> model.sigma
2.4142135
>>> model.epsilon
-0.414213
>>> model.Psi
0.45724
>>> model.f_w_rule(0.)
0.37464
>>> model.f_w_rule(1.)
1.64698

Parameters: f_w (float, derived) – empirical expression used in \(\\alpha\) [dimensionless?]

f_w_rule(w)[source]¶

Parameters: w – accentric factor
Returns: \(f_w\), see Equation (14)