polykin.thermo.eos

Soave

Soave equation of state.

This EOS is based on the following \(P(v,T)\) relationship:

\[ P = \frac{RT}{v - b_m} -\frac{a_m}{v (v + b_m)} \]

where \(P\) is the pressure, \(T\) is the temperature, \(v\) is the molar volume, \(a_m(T,z)\) and \(b_m(z)\) are the mixture EOS parameters, and \(z\) is the vector of mole fractions.

For a single component, the parameters \(a\) and \(b\) are given by:

\[\begin{aligned} a &= 0.42748 \frac{R^2 T_{c}^2}{P_{c}} [1 + f_\omega(1 - T_{r}^{1/2})]^2 \\ f_\omega &= 0.48 + 1.574\omega - 0.176\omega^2 \\ b &= 0.08664\frac{R T_{c}}{P_{c}} \end{aligned}\]

where \(T_c\) is the critical temperature, \(P_c\) is the critical pressure, and \(T_r = T/T_c\) is the reduced temperature.

References

• RC Reid, JM Prausniz, and BE Poling. The properties of gases & liquids 4th edition, 1986, p. 37, 40, 80, 82.

PARAMETER	DESCRIPTION
Tc	Critical temperatures of all components. Unit = K. TYPE: float | FloatVectorLike
Pc	Critical pressures of all components. Unit = Pa. TYPE: float | FloatVectorLike
w	Acentric factors of all components. TYPE: float | FloatVectorLike
k	Binary interaction parameter matrix. TYPE: FloatSquareMatrix | None DEFAULT: None
name	Name. TYPE: str DEFAULT: ''

Source code in src/polykin/thermo/eos/cubic.py

class Soave(CubicEoS):
    r"""[Soave](https://en.wikipedia.org/wiki/Cubic_equations_of_state#Soave_modification_of_Redlich%E2%80%93Kwong)
    equation of state.

    This EOS is based on the following $P(v,T)$ relationship:

    $$ P = \frac{RT}{v - b_m} -\frac{a_m}{v (v + b_m)} $$

    where $P$ is the pressure, $T$ is the temperature, $v$ is the molar
    volume, $a_m(T,z)$ and $b_m(z)$ are the mixture EOS parameters, and
    $z$ is the vector of mole fractions.

    For a single component, the parameters $a$ and $b$ are given by:

    \begin{aligned}
    a &= 0.42748 \frac{R^2 T_{c}^2}{P_{c}} [1 + f_\omega(1 - T_{r}^{1/2})]^2 \\
    f_\omega &= 0.48 + 1.574\omega - 0.176\omega^2 \\
    b &= 0.08664\frac{R T_{c}}{P_{c}}
    \end{aligned}

    where $T_c$ is the critical temperature, $P_c$ is the critical pressure,
    and $T_r = T/T_c$ is the reduced temperature.

    **References**

    *   RC Reid, JM Prausniz, and BE Poling. The properties of gases &
        liquids 4th edition, 1986, p. 37, 40, 80, 82.

    Parameters
    ----------
    Tc : float | FloatVectorLike
        Critical temperatures of all components. Unit = K.
    Pc : float | FloatVectorLike
        Critical pressures of all components. Unit = Pa.
    w : float | FloatVectorLike
        Acentric factors of all components.
    k : FloatSquareMatrix | None
        Binary interaction parameter matrix.
    name : str
        Name.
    """
    _u = 1.0
    _w = 0.0
    _Ωa = 0.42748
    _Ωb = 0.08664

    def __init__(self,
                 Tc: Union[float, FloatVectorLike],
                 Pc: Union[float, FloatVectorLike],
                 w: Union[float, FloatVectorLike],
                 k: Optional[FloatSquareMatrix] = None,
                 name: str = ''
                 ) -> None:

        super().__init__(Tc, Pc, w, k, name)

    def _alpha(self, T: float) -> FloatVector:
        w = self.w
        Tr = T/self.Tc
        fw = 0.48 + 1.574*w - 0.176*w**2
        return (1.0 + fw*(1.0 - sqrt(Tr)))**2

Bm

Bm(T: float, z: FloatVector) -> float

Calculate the second virial coefficient of the mixture.

\[ B_m = b_m - \frac{a_m}{R T} \]

References

• RC Reid, JM Prausniz, and BE Poling. The properties of gases & liquids 4th edition, 1986, p. 82.

PARAMETER	DESCRIPTION
T	Temperature. Unit = K. TYPE: float
z	Mole fractions of all components. Unit = mol/mol. TYPE: FloatVector

RETURNS	DESCRIPTION
float	Mixture second virial coefficient, \(B_m\). Unit = m³/mol.

Source code in src/polykin/thermo/eos/cubic.py

def Bm(self,
       T: float,
       z: FloatVector
       ) -> float:
    r"""Calculate the second virial coefficient of the mixture.

    $$ B_m = b_m - \frac{a_m}{R T} $$

    **References**

    *   RC Reid, JM Prausniz, and BE Poling. The properties of gases &
        liquids 4th edition, 1986, p. 82.

    Parameters
    ----------
    T : float
        Temperature. Unit = K.
    z : FloatVector
        Mole fractions of all components. Unit = mol/mol.

    Returns
    -------
    float
        Mixture second virial coefficient, $B_m$. Unit = m³/mol.
    """
    return self.bm(z) - self.am(T, z)/(R*T)

DA

DA(T, V, n, v0)

Calculate the departure of Helmholtz energy.

\[ A(T,V,n) - A^{\circ}(T,V,n)\]

PARAMETER	DESCRIPTION
T	Temperature. Unit = K. TYPE: float
V	Volume. Unit = m³. TYPE: float
n	Mole amounts of all components. Unit = mol. TYPE: FloatVector
v0	Molar volume in reference state. Unit = m³/mol. TYPE: float

RETURNS	DESCRIPTION
float	Helmholtz energy departure, \(A - A^{\circ}\). Unit = J.

Source code in src/polykin/thermo/eos/cubic.py

def DA(self, T, V, n, v0):
    nT = n.sum()
    z = n/nT
    u = self._u
    w = self._w
    d = sqrt(u**2 - 4*w)
    am = self.am(T, z)
    bm = self.bm(z)

    return nT*am/(bm*d)*log((2*V + nT*bm*(u - d))/(2*V + nT*bm*(u + d))) \
        - nT*R*T*log((V - nT*bm)/(nT*v0))

K

K(
    T: float, P: float, x: FloatVector, y: FloatVector
) -> FloatVector

Calculate the K-values of all components.

\[ K_i = \frac{y_i}{x_i} = \frac{\hat{\phi}_i^L}{\hat{\phi}_i^V} \]

PARAMETER	DESCRIPTION
T	Temperature. Unit = K. TYPE: float
P	Pressure. Unit = Pa. TYPE: float
x	Liquid mole fractions of all components. Unit = mol/mol. TYPE: FloatVector
y	Vapor mole fractions of all components. Unit = mol/mol. TYPE: FloatVector

RETURNS	DESCRIPTION
FloatVector	K-values of all components.

Source code in src/polykin/thermo/eos/base.py

def K(self,
      T: float,
      P: float,
      x: FloatVector,
      y: FloatVector
      ) -> FloatVector:
    r"""Calculate the K-values of all components.

    $$ K_i = \frac{y_i}{x_i} = \frac{\hat{\phi}_i^L}{\hat{\phi}_i^V} $$

    Parameters
    ----------
    T : float
        Temperature. Unit = K.
    P : float
        Pressure. Unit = Pa.
    x : FloatVector
        Liquid mole fractions of all components. Unit = mol/mol.
    y : FloatVector
        Vapor mole fractions of all components. Unit = mol/mol.

    Returns
    -------
    FloatVector
        K-values of all components.
    """
    phiV = self.phi(T, P, y, 'V')
    phiL = self.phi(T, P, x, 'L')
    return phiL/phiV

N property

N: int

Number of components.

P

P(T: float, v: float, z: FloatVector) -> float

Calculate the pressure of the fluid.

PARAMETER	DESCRIPTION
T	Temperature. Unit = K. TYPE: float
v	Molar volume. Unit = m³/mol. TYPE: float
z	Mole fractions of all components. Unit = mol/mol. TYPE: FloatVector

RETURNS	DESCRIPTION
float	Pressure. Unit = Pa.

Source code in src/polykin/thermo/eos/cubic.py

def P(self,
      T: float,
      v: float,
      z: FloatVector
      ) -> float:
    r"""Calculate the pressure of the fluid.

    Parameters
    ----------
    T : float
        Temperature. Unit = K.
    v : float
        Molar volume. Unit = m³/mol.
    z : FloatVector
        Mole fractions of all components. Unit = mol/mol.

    Returns
    -------
    float
        Pressure. Unit = Pa.
    """
    am = self.am(T, z)
    bm = self.bm(z)
    u = self._u
    w = self._w
    return R*T/(v - bm) - am/(v**2 + u*v*bm + w*bm**2)

Psat

Psat(T: float, Psat0: float | None = None) -> float

Calculate the saturation pressure of the fluid.

Note

The saturation pressure is only defined for single-component systems. For multicomponent systems, a specific flash solver must be used: polykin.thermo.flash.

PARAMETER	DESCRIPTION
T	Temperature. Unit = K. TYPE: float
Psat0	Initial guess for the saturation pressure. By default, the value is estimated using the Wilson equation. Unit = Pa. TYPE: float | None DEFAULT: None

RETURNS	DESCRIPTION
float	Saturation pressure. Unit = Pa.

Source code in src/polykin/thermo/eos/cubic.py

def Psat(self,
         T: float,
         Psat0: float | None = None
         ) -> float:
    r"""Calculate the saturation pressure of the fluid.

    !!! note

        The saturation pressure is only defined for single-component
        systems. For multicomponent systems, a specific flash solver
        must be used: [polykin.thermo.flash](../flash/index.md).

    Parameters
    ----------
    T : float
        Temperature. Unit = K.
    Psat0 : float | None
        Initial guess for the saturation pressure. By default, the value
        is estimated using the Wilson equation. Unit = Pa.

    Returns
    -------
    float
        Saturation pressure. Unit = Pa.
    """

    if self.N != 1:
        raise ValueError("Psat is only defined for single-component systems.")

    if T >= self.Tc[0]:
        raise ValueError(f"T >= Tc = {self.Tc[0]} K.")

    if Psat0 is None:
        Psat0 = Psat_guess_Wilson(T, self.Tc[0], self.Pc[0], self.w[0])

    # Solve as fixed-point problem (Newton fails when T is close to Tc)
    def g(x, T=T, z=np.array([1.0])):
        return x*self.K(T, x[0], z, z)

    sol = fixpoint_wegstein(g, np.array([Psat0]))

    if sol.success:
        return sol.x[0]
    else:
        raise ConvergenceError(
            f"Psat failed to converge. Solution: {sol}.")

Z

Z(T: float, P: float, z: FloatVector) -> FloatVector

Calculate the compressibility factors for the possible phases of a fluid.

The calculation is handled by Z_cubic_roots.

PARAMETER	DESCRIPTION
T	Temperature. Unit = K. TYPE: float
P	Pressure. Unit = Pa. TYPE: float
z	Mole fractions of all components. Unit = mol/mol. TYPE: FloatVector

RETURNS	DESCRIPTION
FloatVector	Compressibility factors of the possible phases. If two phases are possible, the first result is the lowest value (liquid).

Source code in src/polykin/thermo/eos/cubic.py

def Z(self,
      T: float,
      P: float,
      z: FloatVector
      ) -> FloatVector:
    r"""Calculate the compressibility factors for the possible phases of a
    fluid.

    The calculation is handled by [`Z_cubic_roots`](Z_cubic_roots.md).

    Parameters
    ----------
    T : float
        Temperature. Unit = K.
    P : float
        Pressure. Unit = Pa.
    z : FloatVector
        Mole fractions of all components. Unit = mol/mol.

    Returns
    -------
    FloatVector
        Compressibility factors of the possible phases. If two phases
        are possible, the first result is the lowest value (liquid).
    """
    A = self.am(T, z)*P/(R*T)**2
    B = self.bm(z)*P/(R*T)
    return Z_cubic_roots(self._u, self._w, A, B)

a cached

a(T: float) -> FloatVector

Calculate the attractive parameters of the pure-components that make up the mixture.

PARAMETER	DESCRIPTION
T	Temperature. Unit = K. TYPE: float

RETURNS	DESCRIPTION
FloatVector	Attractive parameters of all components, \(a_i\). Unit = J·m³.

Source code in src/polykin/thermo/eos/cubic.py

@functools.cache
def a(self,
      T: float
      ) -> FloatVector:
    r"""Calculate the attractive parameters of the pure-components that
    make up the mixture.

    Parameters
    ----------
    T : float
        Temperature. Unit = K.

    Returns
    -------
    FloatVector
        Attractive parameters of all components, $a_i$. Unit = J·m³.
    """
    return self._Ωa * (R*self.Tc)**2 / self.Pc * self._alpha(T)

am

am(T: float, z: FloatVector) -> float

Calculate the mixture attractive parameter from the corresponding pure-component parameters.

\[ a_m = \sum_i \sum_j z_i z_j (a_i a_j)^{1/2} (1 - \bar{k}_{ij}) \]

References

• RC Reid, JM Prausniz, and BE Poling. The properties of gases & liquids 4th edition, 1986, p. 82.

PARAMETER	DESCRIPTION
T	Temperature. Unit = K. TYPE: float
z	Mole fractions of all components. Unit = mol/mol. TYPE: FloatVector

RETURNS	DESCRIPTION
float	Mixture attractive parameter, \(a_m\). Unit = J·m³.

Source code in src/polykin/thermo/eos/cubic.py

def am(self,
       T: float,
       z: FloatVector
       ) -> float:
    r"""Calculate the mixture attractive parameter from the corresponding
    pure-component parameters.

    $$ a_m = \sum_i \sum_j z_i z_j (a_i a_j)^{1/2} (1 - \bar{k}_{ij}) $$

    **References**

    *   RC Reid, JM Prausniz, and BE Poling. The properties of gases &
        liquids 4th edition, 1986, p. 82.

    Parameters
    ----------
    T : float
        Temperature. Unit = K.
    z : FloatVector
        Mole fractions of all components. Unit = mol/mol.

    Returns
    -------
    float
        Mixture attractive parameter, $a_m$. Unit = J·m³.
    """
    return geometric_interaction_mixing(z, self.a(T), self.k)

b cached property

b: FloatVector

Calculate the repulsive parameters of the pure-components that make up the mixture.

RETURNS	DESCRIPTION
FloatVector	Repulsive parameters of all components, \(b_i\). Unit = m³/mol.

bm

bm(z: FloatVector) -> float

Calculate the mixture repulsive parameter from the corresponding pure-component parameters.

\[ b_m = \sum_i z_i b_i \]

References

• RC Reid, JM Prausniz, and BE Poling. The properties of gases & liquids 4th edition, 1986, p. 82.

PARAMETER	DESCRIPTION
z	Mole fractions of all components. Unit = mol/mol. TYPE: FloatVector

RETURNS	DESCRIPTION
float	Mixture repulsive parameter, \(b_m\). Unit = m³/mol.

Source code in src/polykin/thermo/eos/cubic.py

def bm(self,
       z: FloatVector
       ) -> float:
    r"""Calculate the mixture repulsive parameter from the corresponding
    pure-component parameters.

    $$ b_m = \sum_i z_i b_i $$

    **References**

    *   RC Reid, JM Prausniz, and BE Poling. The properties of gases &
        liquids 4th edition, 1986, p. 82.

    Parameters
    ----------
    z : FloatVector
        Mole fractions of all components. Unit = mol/mol.

    Returns
    -------
    float
        Mixture repulsive parameter, $b_m$. Unit = m³/mol.
    """
    return dot(z, self.b)

f

f(
    T: float,
    P: float,
    z: FloatVector,
    phase: Literal["L", "V"],
) -> FloatVector

Calculate the fugacity of all components in a given phase.

For each component, the fugacity is given by:

\[ \hat{f}_i = \hat{\phi}_i z_i P \]

where \(\hat{\phi}_i(T,P,y)\) is the fugacity coefficient, \(P\) is the pressure, and \(z_i\) is the mole fraction.

PARAMETER	DESCRIPTION
T	Temperature. Unit = K. TYPE: float
P	Pressure. Unit = Pa. TYPE: float
z	Mole fractions of all components. Unit = mol/mol. TYPE: FloatVector
phase	Phase of the fluid. Only relevant for systems where both liquid and vapor phases may exist. TYPE: Literal['L', 'V']

RETURNS	DESCRIPTION
FloatVector	Fugacities of all components. Unit = Pa.

Source code in src/polykin/thermo/eos/base.py

def f(self,
      T: float,
      P: float,
      z: FloatVector,
      phase: Literal['L', 'V']
      ) -> FloatVector:
    r"""Calculate the fugacity of all components in a given phase.

    For each component, the fugacity is given by:

    $$ \hat{f}_i = \hat{\phi}_i z_i P $$

    where $\hat{\phi}_i(T,P,y)$ is the fugacity coefficient, $P$ is the 
    pressure, and $z_i$ is the mole fraction.

    Parameters
    ----------
    T : float
        Temperature. Unit = K.
    P : float
        Pressure. Unit = Pa.
    z : FloatVector
        Mole fractions of all components. Unit = mol/mol.
    phase : Literal['L', 'V']
        Phase of the fluid. Only relevant for systems where both liquid
        and vapor phases may exist.

    Returns
    -------
    FloatVector
        Fugacities of all components. Unit = Pa.
    """
    return self.phi(T, P, z, phase)*z*P

phi

phi(
    T: float,
    P: float,
    z: FloatVector,
    phase: Literal["L", "V"],
) -> FloatVector

Calculate the fugacity coefficients of all components in a given phase.

For each component, the fugacity coefficient is given by:

\[\begin{aligned} \ln \hat{\phi}_i &= \frac{b_i}{b_m}(Z-1)-\ln(Z-B^*) +\frac{A^*}{B^*\sqrt{u^2-4w}}\left(\frac{b_i}{b_m}-\delta_i\right)\ln{\frac{2Z+B^*(u+\sqrt{u^2-4w})}{2Z+B^*(u-\sqrt{u^2-4w})}} \\ \delta_i &= \frac{2a_i^{1/2}}{a_m}\sum_j z_j a_j^{1/2}(1-\bar{k}_{ij}) \end{aligned}\]

References

• RC Reid, JM Prausniz, and BE Poling. The properties of gases & liquids 4th edition, 1986, p. 145.

PARAMETER	DESCRIPTION
T	Temperature. Unit = K. TYPE: float
P	Pressure. Unit = Pa. TYPE: float
z	Mole fractions of all components. Unit = mol/mol. TYPE: FloatVector
phase	Phase of the fluid. Only relevant for systems where both liquid and vapor phases may exist. TYPE: Literal['L', 'V']

RETURNS	DESCRIPTION
FloatVector	Fugacity coefficients of all components.

Source code in src/polykin/thermo/eos/cubic.py

def phi(self,
        T: float,
        P: float,
        z: FloatVector,
        phase: Literal['L', 'V']
        ) -> FloatVector:
    r"""Calculate the fugacity coefficients of all components in a given
    phase.

    For each component, the fugacity coefficient is given by:

    \begin{aligned}
    \ln \hat{\phi}_i &= \frac{b_i}{b_m}(Z-1)-\ln(Z-B^*)
    +\frac{A^*}{B^*\sqrt{u^2-4w}}\left(\frac{b_i}{b_m}-\delta_i\right)\ln{\frac{2Z+B^*(u+\sqrt{u^2-4w})}{2Z+B^*(u-\sqrt{u^2-4w})}} \\
    \delta_i &= \frac{2a_i^{1/2}}{a_m}\sum_j z_j a_j^{1/2}(1-\bar{k}_{ij})
    \end{aligned}

    **References**

    *   RC Reid, JM Prausniz, and BE Poling. The properties of gases &
        liquids 4th edition, 1986, p. 145.

    Parameters
    ----------
    T : float
        Temperature. Unit = K.
    P : float
        Pressure. Unit = Pa.
    z : FloatVector
        Mole fractions of all components. Unit = mol/mol.
    phase : Literal['L', 'V']
        Phase of the fluid. Only relevant for systems where both liquid
        and vapor phases may exist.

    Returns
    -------
    FloatVector
        Fugacity coefficients of all components.
    """
    u = self._u
    w = self._w
    d = sqrt(u**2 - 4*w)
    a = self.a(T)
    am = self.am(T, z)
    b = self.b
    bm = self.bm(z)
    k = self.k
    A = am*P/(R*T)**2
    B = bm*P/(R*T)
    b_bm = b/bm
    Z = self.Z(T, P, z)

    if k is None:
        δ = 2*sqrt(a/am)
    else:
        δ = np.sum(z * sqrt(a) * (1 - k), axis=1)

    if phase == 'L':
        Zi = Z[0]
    elif phase == 'V':
        Zi = Z[-1]
    else:
        raise ValueError(f"Invalid phase: {phase}.")

    ln_phi = b_bm*(Zi - 1) - log(Zi - B) + A/(B*d) * \
        (b_bm - δ)*log((2*Zi + B*(u + d))/(2*Zi + B*(u - d)))

    return exp(ln_phi)

v

v(T: float, P: float, z: FloatVector) -> FloatVector

Calculate the molar volumes of the possible phases a fluid.

\[ v = \frac{Z R T}{P} \]

where \(v\) is the molar volume, \(Z(T, P, z)\) is the compressibility factor, \(T\) is the temperature, and \(P\) is the pressure, and \(z\) is the mole fraction vector.

PARAMETER	DESCRIPTION
T	Temperature. Unit = K. TYPE: float
P	Pressure. Unit = Pa. TYPE: float
z	Mole fractions of all components. Unit = mol/mol. TYPE: FloatVector

RETURNS	DESCRIPTION
FloatVector	Molar volumes of the possible phases. If two phases are possible, the first result is the lowest value (liquid). Unit = m³/mol.

Source code in src/polykin/thermo/eos/base.py

def v(self,
      T: float,
      P: float,
      z: FloatVector
      ) -> FloatVector:
    r"""Calculate the molar volumes of the possible phases a fluid.

    $$ v = \frac{Z R T}{P} $$

    where $v$ is the molar volume, $Z(T, P, z)$ is the compressibility
    factor, $T$ is the temperature, and $P$ is the pressure, and $z$ is the
    mole fraction vector.

    Parameters
    ----------
    T : float
        Temperature. Unit = K.
    P : float
        Pressure. Unit = Pa.
    z : FloatVector
        Mole fractions of all components. Unit = mol/mol.

    Returns
    -------
    FloatVector
        Molar volumes of the possible phases. If two phases are possible,
        the first result is the lowest value (liquid). Unit = m³/mol.
    """
    return self.Z(T, P, z)*R*T/P

Examples

Estimate the compressibility factor of a 50 mol% ethylene/nitrogen gas mixture at 300 K and 100 bar.

from polykin.thermo.eos import Soave
import numpy as np

Tc = [282.4, 126.2]    # K
Pc = [50.4e5, 33.9e5]  # Pa
w = [0.089, 0.039]

eos = Soave(Tc, Pc, w)
Z = eos.Z(T=300., P=100e5, z=np.array([0.5, 0.5]))

print(f"{Z[0]:.2f}")

0.83