polykin.thermo.eos¤

RedlichKwong ¤

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}} T_{r}^{-1/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`
`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 RedlichKwong(CubicEoS):
    r"""[Redlich-Kwong](https://en.wikipedia.org/wiki/Redlich%E2%80%93Kwong_equation_of_state)
    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}} T_{r}^{-1/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.
    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],
                 k: Optional[FloatSquareMatrix] = None,
                 name: str = ''
                 ) -> None:

        w = np.zeros_like(Tc) if isinstance(Tc, Iterable) else 0.0
        super().__init__(Tc, Pc, w, k, name)

    def _alpha(self, T: float) -> FloatVector:
        return sqrt(self.Tc/T)

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 RedlichKwong
import numpy as np

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

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

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

0.79

polykin.thermo.eos¤

RedlichKwong ¤

Bm ¤

DA ¤

K ¤

N property ¤

P ¤

Psat ¤

Z ¤

a cached ¤

am ¤

b cached property ¤

bm ¤

f ¤

phi ¤

v ¤

Examples¤

N `property` ¤

a `cached` ¤

b `cached` `property` ¤