Skip to content

polykin.thermo.eos¤

CubicEoS ¤

Base class for cubic equations of state.

This abstract class represents a general two-parameter cubic EOS of the form:

\[ P = \frac{R T}{v - b_m} - \frac{a_m}{v^2 + u b_m v + w b_m^2} \]

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, \(z\) is the vector of mole fractions, and \(u\) and \(w\) are numerical constants that determine the specific EOS.

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

\[\begin{aligned} a &= \Omega_a \left( \frac{R^2 T_{c}^2}{P_{c}} \right) \alpha(T) \\ b &= \Omega_b \left( \frac{R T_{c}}{P_{c}} \right) \end{aligned}\]

where \(T_c\) is the critical temperature, \(P_c\) is the critical pressure, \(\Omega_a\) and \(\Omega_b\) are numerical constants that determine the specific EOS, and \(\alpha(T)\) is a dimensionless temperature-dependent function.

To implement a specific cubic EOS, subclasses must:

  • Define the class-level constants _Ωa, _Ωb, _u and _w.
  • Implement the temperature-dependent function _alpha(T).
  • Optionally override the mixing-rule functions am(T, z) and bm(z).

References

  • RC Reid, JM Prausniz, and BE Poling. The properties of gases & liquids 4th edition, 1986, p. 42.
Source code in src/polykin/thermo/eos/cubic.py 
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
class CubicEoS(GasLiquidEoS):
    r"""Base class for cubic equations of state.

    This abstract class represents a general two-parameter cubic EOS of the
    form:

    $$ P = \frac{R T}{v - b_m} - \frac{a_m}{v^2 + u b_m v + w b_m^2} $$

    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, $z$ is the
    vector of mole fractions, and $u$ and $w$ are numerical constants that
    determine the specific EOS.

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

    \begin{aligned}
    a &= \Omega_a \left( \frac{R^2 T_{c}^2}{P_{c}} \right) \alpha(T) \\ 
    b &= \Omega_b \left( \frac{R T_{c}}{P_{c}} \right)
    \end{aligned}

    where $T_c$ is the critical temperature, $P_c$ is the critical pressure,
    $\Omega_a$ and $\Omega_b$ are numerical constants that determine the
    specific EOS, and $\alpha(T)$ is a dimensionless temperature-dependent
    function.

    To implement a specific cubic EOS, subclasses must:

    * Define the class-level constants `_Ωa`, `_Ωb`, `_u` and `_w`.
    * Implement the temperature-dependent function `_alpha(T)`.
    * Optionally override the mixing-rule functions `am(T, z)` and `bm(z)`.

    **References**

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

    Tc: FloatVector
    Pc: FloatVector
    w: FloatVector
    k: FloatSquareMatrix | None
    _u: float
    _w: float
    _Ωa: float
    _Ωb: float

    def __init__(self,
                 Tc: float | FloatVectorLike,
                 Pc: float | FloatVectorLike,
                 w: float | FloatVectorLike,
                 k: FloatSquareMatrix | None,
                 name: str = ""
                 ) -> None:

        Tc, Pc, w = convert_FloatOrVectorLike_to_FloatVector([Tc, Pc, w])

        N = len(Tc)
        super().__init__(N, name)
        self.Tc = Tc
        self.Pc = Pc
        self.w = w
        self.k = k

    @abstractmethod
    def _alpha(self, T: float) -> FloatVector:
        pass

    @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 [K].

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

    @functools.cached_property
    def b(self) -> FloatVector:
        r"""Calculate the repulsive parameters of the pure-components that
        make up the mixture.

        Returns
        -------
        FloatVector (N)
            Repulsive parameters of all components, $b_i$ [m³/mol].
        """
        return self._Ωb*R*self.Tc/self.Pc

    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 [K].
        z : FloatVector (N)
            Mole fractions of all components [mol/mol].

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

    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 (N)
            Mole fractions of all components [mol/mol].

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

    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 [K].
        z : FloatVector (N)
            Mole fractions of all components [mol/mol].

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

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

        Parameters
        ----------
        T : float
            Temperature [K].
        v : float
            Molar volume [m³/mol].
        z : FloatVector (N)
            Mole fractions of all components [mol/mol].

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

    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 [K].
        P : float
            Pressure [Pa].
        z : FloatVector (N)
            Mole fractions of all components [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)
        Z = Z_cubic_roots(self._u, self._w, A, B)
        return Z

    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 [K].
        P : float
            Pressure [Pa].
        z : FloatVector (N)
            Mole fractions of all components [mol/mol].
        phase : Literal['L', 'V']
            Phase of the fluid. Only relevant for systems where both liquid
            and vapor phases may exist.

        Returns
        -------
        FloatVector (N)
            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)

    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 [K].
        Psat0 : float | None
            Initial guess for the saturation pressure [Pa]. By default, the
            value is estimated using the Wilson equation.

        Returns
        -------
        float
            Saturation pressure [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 = PL_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}.")

    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))

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 [K].

TYPE: float

z

Mole fractions of all components [mol/mol].

TYPE: FloatVector(N)

RETURNS DESCRIPTION
float

Mixture second virial coefficient, \(B_m\) [m³/mol].
Source code in src/polykin/thermo/eos/cubic.py 
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
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 [K].
    z : FloatVector (N)
        Mole fractions of all components [mol/mol].

    Returns
    -------
    float
        Mixture second virial coefficient, $B_m$ [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 [K].

TYPE: float

V

Volume [m³].

TYPE: float

n

Mole amounts of all components [mol].

TYPE: FloatVector(N)

v0

Molar volume in reference state [m³/mol].

TYPE: float

RETURNS DESCRIPTION
float

Helmholtz energy departure, \(A - A^{\circ}\) [J].
Source code in src/polykin/thermo/eos/cubic.py 
385
386
387
388
389
390
391
392
393
394
395
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 [K].

TYPE: float

P

Pressure [Pa].

TYPE: float

x

Liquid mole fractions of all components [mol/mol].

TYPE: FloatVector(N)

y

Vapor mole fractions of all components [mol/mol].

TYPE: FloatVector(N)

RETURNS DESCRIPTION
FloatVector(N)

K-values of all components.
Source code in src/polykin/thermo/eos/base.py 
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
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 [K].
    P : float
        Pressure [Pa].
    x : FloatVector (N)
        Liquid mole fractions of all components [mol/mol].
    y : FloatVector (N)
        Vapor mole fractions of all components [mol/mol].

    Returns
    -------
    FloatVector (N)
        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 [K].

TYPE: float

v

Molar volume [m³/mol].

TYPE: float

z

Mole fractions of all components [mol/mol].

TYPE: FloatVector(N)

RETURNS DESCRIPTION
float

Pressure [Pa].
Source code in src/polykin/thermo/eos/cubic.py 
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
def P(self,
      T: float,
      v: float,
      z: FloatVector
      ) -> float:
    r"""Calculate the pressure of the fluid.

    Parameters
    ----------
    T : float
        Temperature [K].
    v : float
        Molar volume [m³/mol].
    z : FloatVector (N)
        Mole fractions of all components [mol/mol].

    Returns
    -------
    float
        Pressure [Pa].
    """
    u = self._u
    w = self._w
    am = self.am(T, z)
    bm = self.bm(z)
    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 [K].

TYPE: float

Psat0

Initial guess for the saturation pressure [Pa]. By default, the value is estimated using the Wilson equation.

TYPE: float | None DEFAULT: None

RETURNS DESCRIPTION
float

Saturation pressure [Pa].
Source code in src/polykin/thermo/eos/cubic.py 
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
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 [K].
    Psat0 : float | None
        Initial guess for the saturation pressure [Pa]. By default, the
        value is estimated using the Wilson equation.

    Returns
    -------
    float
        Saturation pressure [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 = PL_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 [K].

TYPE: float

P

Pressure [Pa].

TYPE: float

z

Mole fractions of all components [mol/mol].

TYPE: FloatVector(N)

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 
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
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 [K].
    P : float
        Pressure [Pa].
    z : FloatVector (N)
        Mole fractions of all components [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)
    Z = Z_cubic_roots(self._u, self._w, A, B)
    return Z

a cached ¤

a(T: float) -> FloatVector

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

PARAMETER DESCRIPTION
T

Temperature [K].

TYPE: float

RETURNS DESCRIPTION
FloatVector(N)

Attractive parameters of all components, \(a_i\) [J·m³].
Source code in src/polykin/thermo/eos/cubic.py 
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
@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 [K].

    Returns
    -------
    FloatVector (N)
        Attractive parameters of all components, $a_i$ [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 [K].

TYPE: float

z

Mole fractions of all components [mol/mol].

TYPE: FloatVector(N)

RETURNS DESCRIPTION
float

Mixture attractive parameter, \(a_m\) [J·m³].
Source code in src/polykin/thermo/eos/cubic.py 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
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 [K].
    z : FloatVector (N)
        Mole fractions of all components [mol/mol].

    Returns
    -------
    float
        Mixture attractive parameter, $a_m$ [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(N)

Repulsive parameters of all components, \(b_i\) [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 [mol/mol].

TYPE: FloatVector(N)

RETURNS DESCRIPTION
float

Mixture repulsive parameter, \(b_m\) [m³/mol].
Source code in src/polykin/thermo/eos/cubic.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
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 (N)
        Mole fractions of all components [mol/mol].

    Returns
    -------
    float
        Mixture repulsive parameter, $b_m$ [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 [K].

TYPE: float

P

Pressure [Pa].

TYPE: float

z

Mole fractions of all components [mol/mol].

TYPE: FloatVector(N)

phase

Phase of the fluid. Only relevant for systems where both liquid and vapor phases may exist.

TYPE: Literal['L', 'V']

RETURNS DESCRIPTION
FloatVector(N)

Fugacities of all components [Pa].
Source code in src/polykin/thermo/eos/base.py 
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
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 [K].
    P : float
        Pressure [Pa].
    z : FloatVector (N)
        Mole fractions of all components [mol/mol].
    phase : Literal['L', 'V']
        Phase of the fluid. Only relevant for systems where both liquid
        and vapor phases may exist.

    Returns
    -------
    FloatVector (N)
        Fugacities of all components [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 [K].

TYPE: float

P

Pressure [Pa].

TYPE: float

z

Mole fractions of all components [mol/mol].

TYPE: FloatVector(N)

phase

Phase of the fluid. Only relevant for systems where both liquid and vapor phases may exist.

TYPE: Literal['L', 'V']

RETURNS DESCRIPTION
FloatVector(N)

Fugacity coefficients of all components.
Source code in src/polykin/thermo/eos/cubic.py 
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
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 [K].
    P : float
        Pressure [Pa].
    z : FloatVector (N)
        Mole fractions of all components [mol/mol].
    phase : Literal['L', 'V']
        Phase of the fluid. Only relevant for systems where both liquid
        and vapor phases may exist.

    Returns
    -------
    FloatVector (N)
        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 [K].

TYPE: float

P

Pressure [Pa].

TYPE: float

z

Mole fractions of all components [mol/mol].

TYPE: FloatVector(N)

RETURNS DESCRIPTION
FloatVector

Molar volumes of the possible phases [m³/mol]. If two phases are possible, the first result is the lowest value (liquid).
Source code in src/polykin/thermo/eos/base.py 
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
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 [K].
    P : float
        Pressure [Pa].
    z : FloatVector (N)
        Mole fractions of all components [mol/mol].

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