5. Virial Equation of State

5.1. Theory

The second order virial equation of state is [GP07]

(1)\[Z = 1 + B\rho = 1 + \frac{BP}{RT}\]

Where the composition dependency of \(B\) is given by the exact mixing rule

(2)\[B = \sum_i \sum_j y_i y_j B_{ij}\]

where \(B_{ij}=B_{ji}\), and \(B_{ii}\) and \(B_{jj}\) are virial coefficients for the pure species

In this package, the useful correlation

\[\frac{BP_\text{c}}{RT_\text{c}} = B^0 + \omega B^1\]

or

(3)\[B = \frac{RT_\text{c}}{P_\text{c}}\left(B^0 + \omega B^1\right)\]

So that, combining Equations (1) and (3), the compressibility can be calculated from dimensionless quantities as

(4)\[Z = 1 + \left(B^0 + \omega B^1\right)\frac{P_\text{r}}{T_\text{r}}\]

where [SVanNessA05]

(5)\[B^0 = 0.083 - \frac{0.422}{T_\text{r}^{1.6}}\]

(6)\[B^1 = 0.139 - \frac{0.172}{T_\text{r}^{4.2}}\]

so that

so that the following derivatives can be computed as

(7)\[\frac{dB^0}{dT_\text{r}} = \frac{0.675}{T_\text{r}^{2.6}}\]

(8)\[\frac{dB^1}{dT_\text{r}} = \frac{0.722}{T_\text{r}^{5.2}}\]

Which allow the \(H^\text{R}\), \(S^\text{R}\), and \(G^\text{R}\) to be readily computed as follows [GP07]

(9)\[\frac{G^\text{R}}{RT} = \left(B^0 + \omega B^1\right)\frac{P_\text{r}}{T_\text{r}}\]

(10)\[\frac{H^\mathrm{R}}{RT} = P_\mathrm{r}\left[ \frac{B^0}{T_\mathrm{r}} - \frac{\mathrm{d} B^0}{\mathrm{d}T_\mathrm{r}} + \omega\left(\frac{B^1}{T_\mathrm{r}} - \frac{\mathrm{d}B^1}{\mathrm{d}T_\mathrm{r}}\right) \right]\]

(11)\[\frac{S^\text{R}}{R} = -P_\text{r}\left(\frac{d B^0}{d T_\text{r}} - \omega\frac{d B^1}{d T_\text{r}}\right)\]

The cross coefficients are calculated as

(12)\[B_{ij} = \frac{RT_{\text{c}ij}}{P_{\text{c}ij}}\left(B^0 + \omega_{ij}B^1\right)\]

so that the cross derivatives can be computed as

\[ \begin{align}\begin{aligned}\frac{dB_{ij}}{dT} = \frac{RT_{\text{c}ij}}{P_{\text{c}ij}}\left(\frac{dB^0}{d T} + \omega_{ij}\frac{dB^1}{dT}\right)\\\frac{dB_{ij}}{dT} = \frac{R}{P_{\text{c}ij}}\left(\frac{dB^0}{d T_{\text{r}ij}} + \omega_{ij}\frac{dB^1}{dT_{\text{r}ij}}\right)\end{aligned}\end{align} \]

or, equivalently,

(13)\[\frac{dB_{ij}}{dT_{\text{r}ij}} = \frac{RT_{\text{c}ij}}{P_{\text{c}ij}}\left(\frac{dB^0}{d T_{\text{r}ij}} + \omega_{ij}\frac{dB^1}{dT_{\text{r}ij}}\right)\]

5.2. Fugacity Coefficients

The Fugacity coefficients are calculated as

\[\ln\hat{\phi}_i=\left(\frac{\partial (nG^\text{R}/R/T)}{\partial n_i}\right)_{P,T,n_j}\]

For the virial equation of state, this becomes [VanNessA82]

(14)\[\ln\hat{\phi}_k = \frac{P}{RT}\left[B_{kk} + \frac{1}{2}\sum_i\sum_jy_iy_j\left(2\delta_{ik} - \delta_{ij}\right)\right]\]

where both \(i\) and \(j\) indices run over all species

(15)\[\delta_{ik} = 2B_{ik} - B_{ii} - B_{kk} = \delta_{ki}\]

and

\[\delta_{ii} = 0\]

5.3. Residual Partial Molar Properties

The partial molar residual free energy of component k is

(16)\[\frac{\bar{G}^\text{R}_k}{RT} = \ln{\hat{\phi}_k} = \frac{P}{RT}\left[B_{kk} + \frac{1}{2}\sum_i\sum_jy_iy_j\left(2\delta_{ik} - \delta_{ij}\right)\right]\]

The partial molar residual enthalpy of component k is

\[\begin{split}\begin{align} \frac{\bar{H}^\mathrm{R}_k}{RT} &= -T \left(\frac{\partial \ln\hat{\phi}_k}{\partial T}\right)_{P,y} \\ &= - \frac{T}{T_c} \left(\frac{\partial \ln\hat{\phi}_k}{\partial T_r}\right)_{P,y} \\ &= - \frac{T}{T_c}\left\{ \frac{P}{RT}\left[ \frac{\partial B_{kk}}{\partial T_r} + \frac{1}{2}\sum_i\sum_j y_iy_j\left( 2\frac{\partial \delta_{ik}}{\partial T_r} - \frac{\partial \delta_{ij}}{\partial T_r} \right) \right] - \frac{T_c\ln\hat{\phi}_i}{T} \right\} \end{align}\end{split}\]

where \(\frac{\partial \delta_{ij}}{\partial T_r}\) is given by (25) so that we obtain

(17)\[\frac{\bar{H}^\mathrm{R}_k}{RT} = - \frac{P}{RT_\text{c}}\left[ \frac{\partial B_{kk}}{\partial T_r} + \frac{1}{2}\sum_i\sum_j y_iy_j\left( 2\frac{\partial \delta_{ik}}{\partial T_r} - \frac{\partial \delta_{ij}}{\partial T_r} \right) \right] + \ln\hat{\phi}_i\]

Since

\[G^\mathrm{R} = H^\mathrm{R} - T S^\mathrm{R}\]

In terms of partial molar properties, then

\[\begin{split}\begin{align} \bar{S}_i^\mathrm{R} &= \frac{\bar{H}_i^\mathrm{R} - \bar{G}_i^\mathrm{R}}{T} \\ \frac{\bar{S}_i^\mathrm{R}}{R} &= \frac{\bar{H}_i^\mathrm{R}}{RT} - \frac{\bar{G}_i^\mathrm{R}}{RT} \\ \end{align}\end{split}\]

By comparing Equation (16) and (17) it is observed that

(18)\[\frac{\bar{S}_k^\mathrm{R}}{R} = - \frac{P}{RT_\text{c}}\left[ \frac{\partial B_{kk}}{\partial T_r} + \frac{1}{2}\sum_i\sum_j y_iy_j\left( 2\frac{\partial \delta_{ik}}{\partial T_r} - \frac{\partial \delta_{ij}}{\partial T_r} \right) \right]\]

where \(\frac{\partial \delta_{ij}}{\partial T_r}\) is given by (25)

The partial molar residual volume of component i is calculated as

\[\begin{split}\begin{align} \frac{\bar{V}_k^\mathrm{R}}{RT} &= \left(\frac{\partial \ln\hat{\phi}_i}{\partial P}\right)_{T,y}\\ &= \frac{\partial}{\partial P}\left\{\frac{P}{RT}\left[B_{kk} + \frac{1}{2}\sum_i\sum_jy_iy_j\left(2\delta_{ik} - \delta_{ij}\right)\right]\right\} \end{align}\end{split}\]

which simplifies to

(19)\[\frac{\bar{V}_k^\mathrm{R}}{RT} = \frac{1}{RT}\left[B_{kk} + \frac{1}{2}\sum_i\sum_jy_iy_j\left(2\delta_{ik} - \delta_{ij}\right)\right]\]

from which we obtain the intuitive result of

\[\bar{V}_k^\mathrm{R} = B_{kk} + \frac{1}{2}\sum_i\sum_jy_iy_j\left(2\delta_{ik} - \delta_{ij}\right)\]

class gasthermo.eos.virial.Virial(pow: callable = <ufunc 'power'>, exp: callable = <ufunc 'exp'>)

Parameters

• pow (callable, optional) – function for computing power, defaults to numpy.power

• exp (callable, optional) – function for computing logarithm, defaults to numpy.exp

B0_expr(T_r)

Parameters

T_r – Reduced temperature

Returns

Equation (5)

B1_expr(T_r)

Parameters

T_r – reduced temperature

Returns

Equation eq:B1_expr

B_expr(T_r, w, T_c, P_c)

Parameters

• T_r – reduced temperature

• w – accentric factor

• T_c – critical temperautre [K]

• P_c – critical pressure [Pa]

Returns

Equation (3)

d_B0_d_Tr_expr(T_r)

Parameters

T_r – reduced temperature

Returns

Equation (7)

d_B1_d_Tr_expr(T_r)

Parameters

T_r – reduced temperature

Returns

Equation (8)

hat_phi_i_expr(*args)

expression for fugacity coefficient :returns: \(\exp\left({\ln{\hat{\phi}_i}}\right)\)

class gasthermo.eos.virial.SecondVirial(dippr_no: str = None, compound_name: str = None, cas_number: str = None, pow: callable = <ufunc 'power'>, **kwargs)

Virial equation of state for one component. See [GP07][SVanNessA05]

G_R_RT_expr(P, T)

Parameters

• P – pressure in Pa

• T – Temperautre in K

Returns

Equation (9)

H_R_RT_expr(P, T)

Parameters

• P – pressure in Pa

• T – temperature in K

Returns

Equation (10)

S_R_R_expr(P, T)

Parameters

• P – pressure in Pa

• T – temperature in K

Returns

Equation (11)

calc_Z_from_dimensionless(P_r, T_r)

Parameters

• P_r – reduced pressure, dimensionless

• T_r – reduced temperature, dimensionless

Returns

Equation (4)

calc_Z_from_units(P, T)

Parameters

• P – pressure in Pa

• T – temperature in K

Returns

Equation (1)

ln_hat_phi_k_expr(P, T)

logarithm of fugacity coefficient

Note

single-component version

In this case, Equation (14) simplifies to

\[\ln\hat{\phi}_i = \frac{PB}{RT}\]

Parameters

• P (float) – pressure in Pa

• T (float) – temperature in K

plot_Z_vs_P(T, P_min, P_max, symbol='o', ax=None, **kwargs)

Plot compressibility as a function of pressure

Parameters

• T (float) – temperature [K]

• P_min (float) – minimum pressure for plotting [Pa]

• P_max (float) – maximum pressure for plotting [Pa]

• phase (str) – phase type (liquid or vapor), defaults to vapor

• symbol (str) – marker symbol, defaults to ‘o’

• ax (plt.axis) – matplotlib axes for plotting, defaults to None

• kwargs – keyword arguments for plotting

5.4. Mixtures

class gasthermo.eos.virial.MixingRule(pow: callable = <ufunc 'power'>, exp: callable = <ufunc 'exp'>)

Van der Walls mixing rule

combining rules from [PLdeAzevedo86]

(20)\[\omega_{ij} = \frac{\omega_i + \omega_j}{2}\]

(21)\[T_{\text{c}ij} = \sqrt{T_{\text{c}i}T_{\text{c}j}}(1-k_{ij})\]

(22)\[P_{\text{c}ij} = \frac{Z_{\text{c}ij}RT_{\text{c}ij}}{V_{\text{c}ij}}\]

(23)\[Z_{\text{c}ij} = \frac{Z_{\text{c}i} + Z_{\text{c}j}}{2}\]

(24)\[V_{\text{c}ij} = \left(\frac{V_{\text{c}i}^{1/3} + V_{\text{c}j}^{1/3}}{2}\right)^3\]

P_cij_rule(Z_ci, V_ci, T_ci, Z_cj, V_cj, T_cj, k_ij)

Returns

Equation (22)

T_cij_rule(T_ci, T_cj, k_ij)

Parameters

• T_ci – critical temperature of component i [K]

• T_cj – ** j [K]

• k_ij – k_ij parameter

Returns

Equation (21)

V_cij_rule(V_ci, V_cj)

Parameters

• V_ci – critical molar volume of component i [m**3/mol]

• V_cj – critical molar volume of component j [m**3/mol]

Returns

Equation eq:Vc_combine

Z_cij_rule(Z_ci, Z_cj)

Parameters

• Z_ci – critical compressibility factor of component i

• Z_cj – critical compressibility factor of component j

Returns

Equation (23)

w_ij_rule(w_i, w_j)

Parameters

• w_i – accentric factor of component i

• w_j – accentric factor of component j

Returns

Equation (20)

class gasthermo.eos.virial.SecondVirialMixture(num_components=0, dippr_nos: List[str] = None, compound_names: List[str] = None, cas_numbers: List[str] = None, MWs: List[float] = None, P_cs: List[float] = None, V_cs: List[float] = None, Z_cs: List[float] = None, T_cs: List[float] = None, ws: List[float] = None, k_ij: Union[float, List[float]] = None, pow: callable = <ufunc 'power'>, exp: callable = <ufunc 'exp'>)

Second virial with mixing rule from MixingRule

Note

can only input both custom critical properties or both from DIPPR–cant have mixed at the moment

Parameters

• num_components (int) – number of components

• dippr_nos (typing.List[str]) – dippr numbers of components

• compound_names (typing.List[str]) – names of components

• cas_numbers (typing.List[str]) – cas registry numbers

• MWs (typing.List[float]) – molecular weights of component sin g/mol

• P_cs (typing.List[float]) – critical pressures of componets in Pa

• T_cs (typing.List[float]) – critical temperatures of pure components in K

• V_cs (typing.List[float]) – critical molar volumes of pure components in m**3/mol

• ws (typing.List[float]) – accentric factors of components

• k_ij (typing.List[typing.List[float]]) – equation of state mixing rule in calculation of critical temperautre, see Equation (21). When \(i=j\) and for chemical similar species, \(k_{ij}=0\). Otherwise, it is a small (usually) positive number evaluated from minimal \(PVT\) data or, in the absence of data, set equal to zero.

B_ij_expr(i: int, j: int, T)

Parameters

• i – index of first component

• j – index of second component

• T – temperature [K]

Returns

Equation (12)

B_mix_expr(y_k: List[Union[float, Any]], T)

Parameters

• y_k – mole fractions of each component \(k\)

• T – temperature in K

Returns

Equation (2)

G_R_RT(*args)

Residual free energy of mixture \(G^\mathrm{R}/R/T\)

H_R_RT(*args)

Residual enthalpy of mixture \(H^\mathrm{R}/R/T\)

M_R_dimensionless(method: callable, ys: List[Union[float, Any]], P: float, T: float)

Residual property of \(X\) for mixture.

Similar to Equation (3) but in dimensionless form

Parameters

method (callable) – function to compute partial molar property of compound

S_R_R(*args)

Residual entropy of mixture \(S^\mathrm{R}/R\)

bar_GiR_RT(cas_k: str, ys: List[Union[float, Any]], P, T)

Dimensionless residual partial molar free energy of component \(i\)

Parameters

• cas_k – cas number of component k

• ys – mole fractions

• P – pressure in Pa

• T – temperature in K

Returns

Equation (16)

bar_HiR_RT(cas_k: str, ys: List[Union[float, Any]], P, T)

Dimensionless residual partial molar enthalpy of component \(i\)

Parameters

• cas_k – cas number for component of interest

• ys – mole fractions

• P – pressure in Pa

• T – temperature in K

bar_SiR_R(cas_k: str, ys: List[Union[float, Any]], P, T)

Dimensionless residual partial molar entropy of component \(i\)

Parameters

• cas_k – cas number for component of interest

• ys – mole fractions

• P – pressure in Pa

• T – temperature in K

Returns

Equation (18)

bar_ViR_RT(cas_k: str, ys: List[Union[float, Any]], P, T)

residual Partial molar volume for component i

Note

This expression does not depend on \(P\)

Parameters

• cas_k – cas number for component of interest

• ys – mole fractions

• P – pressure in Pa

• T – temperature in K

Returns

Equation (19)

calc_Z_from_units(y_k: List, P, T)

Parameters

• y_k – mole fractions of each component \(k\)

• P – pressure in Pa

• T – temperature in K

Returns

Equation (1)

d_Bij_d_Trij(i: int, j: int, T)

Parameters

• i – index for component i

• j – index for componen t j

• T – temperature in K

Returns

Equation (13)

d_dij_d_Tr(i, j, T)

(25)\[\frac{\partial \delta_{ij}}{\partial T_{\text{r}ij}} = 2\frac{\partial B_{ij}}{\partial T_{\text{r}ij}}-\frac{\partial B_{ii}}{\partial T_{\text{r}ij}}-\frac{\partial B_{jj}}{\partial T_{\text{r}ij}}\]

Todo

test this with symbolic differentiation of d_ik expression

Parameters

• i – index for component i

• j – index for componen t j

• T – temperature [K]

d_ik_expr(i: int, k: int, T)

Parameters

• i – index of component i

• k – index of component k

• T – temperature [K]

Returns

Equation (15)

fugacity_i_expr(cas_i: str, ys: List[Union[float, Any]], P, T)

Fugacity of component i in mixture \(f_i=\hat{\phi}_i y_i P\)

Parameters

• cas_i – cas number for component of interest

• ys – mole fractions

• P – pressure in Pa

• T – temperature in K

get_w_Tc_Pc(i: int, j=None)

Returns critical constants for calculation based off of whetner i = j or not

Returns

(\(w\), \(T_c\), \(P_c\))

Return type

tuple

ln_hat_phi_k_expr(k: int, ys: List[Union[float, Any]], P, T)

logarithm of fugacity coefficient

Parameters

• k – index of component k

• ys – mole fractions

• P – pressure in Pa

• T – temperature in K

Returns

Equation (14)

plot_residual_HSG(P, T, ax=None, fig=None) → Tuple[matplotlib.pyplot.figure, matplotlib.pyplot.subplot]

Plot dimensionless residual properties as a function of mole fraction

Parameters

• P – pressure in Pa

• T – Temperature in K

• ax – matplotlib ax, defaults to None

• fig – matplotlib figure, defautls to None