Skip to content

polykin.thermo.eos¤

PengRobinson ¤

Peng-Robinson equation of state.

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

\[ P = \frac{RT}{v - b_m} -\frac{a_m}{v^2 + 2 v b_m - b_m^2} \]

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

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

\[\begin{aligned} a &= 0.45724 \frac{R^2 T_{c}^2}{P_{c}} [1 + f_\omega(1 - T_{r}^{1/2})]^2 \\ f_\omega &= 0.37464 + 1.54226\omega - 0.26992\omega^2 \\ b &= 0.07780\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

Source code in src/polykin/thermo/eos/cubic.py 
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
class PengRobinson(Cubic):
    r"""[Peng-Robinson](https://en.wikipedia.org/wiki/Cubic_equations_of_state#Peng%E2%80%93Robinson_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^2 + 2 v b_m - b_m^2} $$

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

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

    \begin{aligned}
    a &= 0.45724 \frac{R^2 T_{c}^2}{P_{c}} [1 + f_\omega(1 - T_{r}^{1/2})]^2 \\
    f_\omega &= 0.37464 + 1.54226\omega - 0.26992\omega^2 \\
    b &= 0.07780\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.
    """
    _u = 2.
    _w = -1.
    _Ωa = 0.45724
    _Ωb = 0.07780

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

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

    def _alpha(self, T: float) -> FloatVector:
        w = self.w
        Tr = T/self.Tc
        fw = 0.37464 + 1.54226*w - 0.26992*w**2
        return (1. + fw*(1. - sqrt(Tr)))**2

Bm ¤

Bm(T: float, y: 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

y

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 
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
def Bm(self,
       T: float,
       y: 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.
    y : FloatVector
        Mole fractions of all components. Unit = mol/mol.

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

P ¤

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

Calculate the pressure of the fluid.

PARAMETER DESCRIPTION
T

Temperature. Unit = K.

TYPE: float

v

Molar volume. Unit = m³/mol.

TYPE: float

y

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 
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
def P(self,
      T: float,
      v: float,
      y: FloatVector
      ) -> float:
    r"""Calculate the pressure of the fluid.

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

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

Z ¤

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

Calculate the compressibility factors of the coexisting phases a fluid.

The calculation is handled by Z_cubic_roots.

PARAMETER DESCRIPTION
T

Temperature. Unit = K.

TYPE: float

P

Pressure. Unit = Pa.

TYPE: float

y

Mole fractions of all components. Unit = mol/mol.

TYPE: FloatVector

RETURNS DESCRIPTION
FloatVector

Compressibility factor of the vapor and/or liquid phases.
Source code in src/polykin/thermo/eos/cubic.py 
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
def Z(self,
      T: float,
      P: float,
      y: FloatVector
      ) -> FloatVector:
    r"""Calculate the compressibility factors of the coexisting phases a
    fluid.

    The calculation is handled by
    [`Z_cubic_roots`](RedlichKwong.md#polykin.properties.eos.Z_cubic_roots).

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

    Returns
    -------
    FloatVector
        Compressibility factor of the vapor and/or liquid phases.
    """
    A = self.am(T, y)*P/(R*T)**2
    B = self.bm(y)*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 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
@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, y: FloatVector) -> float

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

\[ a_m = \sum_i \sum_j y_i y_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

y

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 
 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
def am(self,
       T: float,
       y: FloatVector
       ) -> float:
    r"""Calculate the mixture attractive parameter from the corresponding
    pure-component parameters.

    $$ a_m = \sum_i \sum_j y_i y_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.
    y : FloatVector
        Mole fractions of all components. Unit = mol/mol.

    Returns
    -------
    float
        Mixture attractive parameter, $a_m$. Unit = J·m³.
    """
    return geometric_interaction_mixing(y, 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(y: FloatVector) -> float

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

\[ b_m = \sum_i y_i b_i \]

References

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

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 
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def bm(self,
       y: FloatVector
       ) -> float:
    r"""Calculate the mixture repulsive parameter from the corresponding
    pure-component parameters.

    $$ b_m = \sum_i y_i b_i $$

    **References**

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

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

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

fV ¤

fV(T: float, P: float, y: FloatVector) -> FloatVector

Calculate the fugacity of all components in the vapor phase.

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

\(\hat{f}_i\) is the fugacity in the vapor phase, \(\hat{\phi}_i(T,P,y)\) is the fugacity coefficient, \(P\) is the pressure, and \(y_i\) is the mole fraction in the vapor phase.

PARAMETER DESCRIPTION
T

Temperature. Unit = K.

TYPE: float

P

Pressure. Unit = Pa.

TYPE: float

y

Mole fractions of all components. Unit = mol/mol.

TYPE: FloatVector

RETURNS DESCRIPTION
FloatVector

Fugacity coefficients of all components.
Source code in src/polykin/thermo/eos/base.py 
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
def fV(self,
        T: float,
        P: float,
        y: FloatVector
       ) -> FloatVector:
    r"""Calculate the fugacity of all components in the vapor phase.

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

    $\hat{f}_i$ is the fugacity in the vapor phase, $\hat{\phi}_i(T,P,y)$
    is the fugacity coefficient, $P$ is the pressure, and $y_i$ is the mole
    fraction in the vapor phase.

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

    Returns
    -------
    FloatVector
        Fugacity coefficients of all components.
    """
    return self.phiV(T, P, y)*y*P

phiV ¤

phiV(T: float, P: float, y: FloatVector) -> FloatVector

Calculate the fugacity coefficients of all components in the vapor phase.

\[\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 y_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

y

Mole fractions of all components. Unit = mol/mol.

TYPE: FloatVector

RETURNS DESCRIPTION
FloatVector

Fugacity coefficients of all components.
Source code in src/polykin/thermo/eos/cubic.py 
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
def phiV(self,
         T: float,
         P: float,
         y: FloatVector
         ) -> FloatVector:
    r"""Calculate the fugacity coefficients of all components in the vapor
    phase.

    \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 y_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.
    y : FloatVector
        Mole fractions of all components. Unit = mol/mol.

    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, y)
    b = self.b
    bm = self.bm(y)
    k = self.k
    A = am*P/(R*T)**2
    B = bm*P/(R*T)
    b_bm = b/bm
    Z = self.Z(T, P, y)

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

    # get only vapor solution !!!
    z = max(Z)
    lnphi = b_bm*(z - 1) - log(z - B) + A/(B*d) * \
        (b_bm - delta)*log((2*z + B*(u + d))/(2*z + B*(u - d)))

    return exp(lnphi)

v ¤

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

Calculate the molar volumes of the coexisting phases a fluid.

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

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

PARAMETER DESCRIPTION
T

Temperature. Unit = K.

TYPE: float

P

Pressure. Unit = Pa.

TYPE: float

y

Mole fractions of all components. Unit = mol/mol.

TYPE: FloatVector

RETURNS DESCRIPTION
FloatVector

Molar volume of the vapor and/or liquid phases. Unit = m³/mol.
Source code in src/polykin/thermo/eos/base.py 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
def v(self,
      T: float,
      P: float,
      y: FloatVector
      ) -> FloatVector:
    r"""Calculate the molar volumes of the coexisting phases a fluid.

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

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

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

    Returns
    -------
    FloatVector
        Molar volume of the vapor and/or liquid phases. Unit = m³/mol.
    """
    return self.Z(T, P, y)*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 PengRobinson
import numpy as np

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

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

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