Skip to content

polykin.thermo.eos¤

This module implements equations of state (EOS) for gas and liquid mixtures.

Virial ¤

Virial equation of state truncated to the second coefficient.

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

\[ P = \frac{R T}{v - B_m} \]

where \(P\) is the pressure, \(T\) is the temperature, \(v\) is the molar volume, \(B_m\) is the mixture second virial coefficient. The parameter \(B_m=B_m(T,y)\) is estimated based on the critical properties of the pure components and the mixture composition \(y\).

Important

This equation is an improvement over the ideal gas model, but it should only be used up to moderate pressures such that \(v/v_c > 2\).

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

Zc

Critical compressibility factors of all components.

TYPE: float | FloatVectorLike

w

Acentric factors of all components.

TYPE: float | FloatVectorLike

Source code in src/polykin/thermo/eos/virial.py 
 26
 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
 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
class Virial(GasEoS):
    r"""Virial equation of state truncated to the second coefficient.

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

    $$ P = \frac{R T}{v - B_m} $$

    where $P$ is the pressure, $T$ is the temperature, $v$ is the molar
    volume, $B_m$ is the mixture second virial coefficient. The parameter
    $B_m=B_m(T,y)$ is estimated based on the critical properties of the pure
    components and the mixture composition $y$.

    !!! important

        This equation is an improvement over the ideal gas model, but it
        should only be used up to moderate pressures such that $v/v_c > 2$.

    **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.
    Zc : float | FloatVectorLike
        Critical compressibility factors of all components.
    w : float | FloatVectorLike
        Acentric factors of all components.
    """

    Tc: FloatVector
    Pc: FloatVector
    Zc: FloatVector
    w: FloatVector

    def __init__(self,
                 Tc: Union[float, FloatVectorLike],
                 Pc: Union[float, FloatVectorLike],
                 Zc: Union[float, FloatVectorLike],
                 w: Union[float, FloatVectorLike]
                 ) -> None:

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

        self.Tc = Tc
        self.Pc = Pc
        self.Zc = Zc
        self.w = w

    def Z(self,
          T: float,
          P: float,
          y: FloatVector
          ) -> float:
        r"""Calculate the compressibility factor of the fluid.

        $$ Z = 1 + \frac{B_m P}{R T} $$

        where $Z$ is the compressibility factor, $P$ is the pressure, $T$ is
        the temperature, and $B_m=B_m(T,y)$ is the mixture second virial
        coefficient.

        **References**

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

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

        Returns
        -------
        float
            Compressibility factor of the fluid.
        """
        Bm = self.Bm(T, y)
        return 1. + Bm*P/(R*T)

    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.
        """
        Bm = self.Bm(T, y)
        return R*T/(v - Bm)

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

        $$ B_m = \sum_i \sum_j y_i y_j B_{ij} $$

        **References**

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

        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 quadratic_mixing(y, self.Bij(T))

    @functools.cache
    def Bij(self,
            T: float,
            ) -> FloatSquareMatrix:
        r"""Calculate the matrix of interaction virial coefficients.

        The calculation is handled by
        [`B_mixture`](.#polykin.properties.eos.B_mixture).

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

        Returns
        -------
        FloatSquareMatrix
            Matrix of interaction virial coefficients, $B_{ij}$.
            Unit = m³/mol.
        """
        return B_mixture(T, self.Tc, self.Pc, self.Zc, self.w)

    def phiV(self,
             T: float,
             P: float,
             y: FloatVector
             ) -> FloatVector:
        r"""Calculate the fugacity coefficients of all components in the vapor
        phase.

        $$
        \ln \hat{\phi}_i = \left(2\sum_j {y_jB_{ij}} -B_m \right)\frac{P}{RT}
        $$

        where $\hat{\phi}_i$ is the fugacity coefficient, $P$ is the pressure,
        $T$ is the temperature, $B_{ij}$ is the matrix of interaction virial
        coefficients, $B_m$ is the second virial coefficient of the mixture,
        and $y_i$ is the mole fraction in the vapor phase.

        **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.
        """
        B = self.Bij(T)
        Bm = self.Bm(T, y)
        return exp((2*dot(B, y) - Bm)*P/(R*T))

    def DA(self, T, V, n, v0):
        nT = n.sum()
        y = n/nT
        Bm = self.Bm(T, y)
        return -nT*R*T*log((V - nT*Bm)/(nT*v0))

Bij cached ¤

Bij(T: float) -> FloatSquareMatrix

Calculate the matrix of interaction virial coefficients.

The calculation is handled by B_mixture.

PARAMETER DESCRIPTION
T

Temperature. Unit = K.

TYPE: float

RETURNS DESCRIPTION
FloatSquareMatrix

Matrix of interaction virial coefficients, \(B_{ij}\). Unit = m³/mol.
Source code in src/polykin/thermo/eos/virial.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
@functools.cache
def Bij(self,
        T: float,
        ) -> FloatSquareMatrix:
    r"""Calculate the matrix of interaction virial coefficients.

    The calculation is handled by
    [`B_mixture`](.#polykin.properties.eos.B_mixture).

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

    Returns
    -------
    FloatSquareMatrix
        Matrix of interaction virial coefficients, $B_{ij}$.
        Unit = m³/mol.
    """
    return B_mixture(T, self.Tc, self.Pc, self.Zc, self.w)

Bm ¤

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

Calculate the second virial coefficient of the mixture.

\[ B_m = \sum_i \sum_j y_i y_j B_{ij} \]

References

  • RC Reid, JM Prausniz, and BE Poling. The properties of gases & liquids 4th edition, 1986, p. 79.
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/virial.py 
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
def Bm(self,
       T: float,
       y: FloatVector
       ) -> float:
    r"""Calculate the second virial coefficient of the mixture.

    $$ B_m = \sum_i \sum_j y_i y_j B_{ij} $$

    **References**

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

    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 quadratic_mixing(y, self.Bij(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/virial.py 
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
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.
    """
    Bm = self.Bm(T, y)
    return R*T/(v - Bm)

Z ¤

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

Calculate the compressibility factor of the fluid.

\[ Z = 1 + \frac{B_m P}{R T} \]

where \(Z\) is the compressibility factor, \(P\) is the pressure, \(T\) is the temperature, and \(B_m=B_m(T,y)\) is the mixture second virial coefficient.

References

  • RC Reid, JM Prausniz, and BE Poling. The properties of gases & liquids 4th edition, 1986, p. 37.
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
float

Compressibility factor of the fluid.
Source code in src/polykin/thermo/eos/virial.py 
 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
def Z(self,
      T: float,
      P: float,
      y: FloatVector
      ) -> float:
    r"""Calculate the compressibility factor of the fluid.

    $$ Z = 1 + \frac{B_m P}{R T} $$

    where $Z$ is the compressibility factor, $P$ is the pressure, $T$ is
    the temperature, and $B_m=B_m(T,y)$ is the mixture second virial
    coefficient.

    **References**

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

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

    Returns
    -------
    float
        Compressibility factor of the fluid.
    """
    Bm = self.Bm(T, y)
    return 1. + Bm*P/(R*T)

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.

\[ \ln \hat{\phi}_i = \left(2\sum_j {y_jB_{ij}} -B_m \right)\frac{P}{RT} \]

where \(\hat{\phi}_i\) is the fugacity coefficient, \(P\) is the pressure, \(T\) is the temperature, \(B_{ij}\) is the matrix of interaction virial coefficients, \(B_m\) is the second virial coefficient of the mixture, and \(y_i\) is the mole fraction in the vapor phase.

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/virial.py 
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
def phiV(self,
         T: float,
         P: float,
         y: FloatVector
         ) -> FloatVector:
    r"""Calculate the fugacity coefficients of all components in the vapor
    phase.

    $$
    \ln \hat{\phi}_i = \left(2\sum_j {y_jB_{ij}} -B_m \right)\frac{P}{RT}
    $$

    where $\hat{\phi}_i$ is the fugacity coefficient, $P$ is the pressure,
    $T$ is the temperature, $B_{ij}$ is the matrix of interaction virial
    coefficients, $B_m$ is the second virial coefficient of the mixture,
    and $y_i$ is the mole fraction in the vapor phase.

    **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.
    """
    B = self.Bij(T)
    Bm = self.Bm(T, y)
    return exp((2*dot(B, y) - Bm)*P/(R*T))

v ¤

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

Calculate the molar volume the fluid.

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

where \(v\) is the molar volue, \(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
float

Molar volume of the fluid. Unit = m³/mol.
Source code in src/polykin/thermo/eos/base.py 
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
def v(self,
      T: float,
      P: float,
      y: FloatVector
      ) -> float:
    r"""Calculate the molar volume the fluid.

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

    where $v$ is the molar volue, $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
    -------
    float
        Molar volume of the fluid. Unit = m³/mol.
    """
    return self.Z(T, P, y)*R*T/P

B_pure ¤

B_pure(
    T: Union[float, FloatArray],
    Tc: float,
    Pc: float,
    w: float,
) -> Union[float, FloatArray]

Estimate the second virial coefficient of a nonpolar or slightly polar gas.

\[ \frac{B P_c}{R T_c} = B^{(0)}(T_r) + \omega B^{(1)}(T_r) \]

where \(B\) is the second virial coefficient, \(P_c\) is the critical pressure, \(T_c\) is the critical temperature, \(T_r=T/T_c\) is the reduced temperature, and \(\omega\) is the acentric factor.

References

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

Temperature. Unit = K.

TYPE: float | FloatArray

Tc

Critical temperature. Unit = K.

TYPE: float

Pc

Critical pressure. Unit = Pa.

TYPE: float

w

Acentric factor.

TYPE: float

RETURNS DESCRIPTION
float | FloatArray

Second virial coefficient, \(B\). Unit = m³/mol.
Source code in src/polykin/thermo/eos/virial.py 
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
def B_pure(T: Union[float, FloatArray],
           Tc: float,
           Pc: float,
           w: float
           ) -> Union[float, FloatArray]:
    r"""Estimate the second virial coefficient of a nonpolar or slightly polar
    gas.

    $$ \frac{B P_c}{R T_c} = B^{(0)}(T_r) + \omega B^{(1)}(T_r) $$

    where $B$ is the second virial coefficient, $P_c$ is the critical pressure,
    $T_c$ is the critical temperature, $T_r=T/T_c$ is the reduced temperature,
    and $\omega$ is the acentric factor.

    **References**

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

    Parameters
    ----------
    T : float | FloatArray
        Temperature. Unit = K.
    Tc : float
        Critical temperature. Unit = K.
    Pc : float
        Critical pressure. Unit = Pa.
    w : float
        Acentric factor.

    Returns
    -------
    float | FloatArray
        Second virial coefficient, $B$. Unit = m³/mol.
    """
    Tr = T/Tc
    B0 = 0.083 - 0.422/Tr**1.6
    B1 = 0.139 - 0.172/Tr**4.2
    return R*Tc/Pc*(B0 + w*B1)

B_mixture ¤

B_mixture(
    T: float,
    Tc: FloatVector,
    Pc: FloatVector,
    Zc: FloatVector,
    w: FloatVector,
) -> FloatSquareMatrix

Calculate the matrix of interaction virial coefficients using the mixing rules of Prausnitz.

\[\begin{aligned} B_{ij} &= B(T,T_{cij},P_{cij},\omega_{ij}) \\ v_{cij} &= \frac{(v_{ci}^{1/3}+v_{cj}^{1/3})^3}{8} \\ k_{ij} &= 1 -\frac{\sqrt{v_{ci}v_{cj}}}{v_{cij}} \\ T_{cij} &= \sqrt{T_{ci}T_{cj}}(1-k_{ij}) \\ Z_{cij} &= \frac{Z_{ci}+Z_{cj}}{2} \\ \omega_{ij} &= \frac{\omega_{i}+\omega_{j}}{2} \\ P_{cij} &= \frac{Z_{cij} R T_{cij}}{v_{cij}} \end{aligned}\]

The calculation of the individual coefficients is handled by B_pure.

References

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

Temperature. Unit = K.

TYPE: float

Tc

Critical temperatures of all components. Unit = K.

TYPE: FloatVector

Pc

Critical pressures of all components. Unit = Pa.

TYPE: FloatVector

Zc

Critical compressibility factors of all components.

TYPE: FloatVector

w

Acentric factors of all components.

TYPE: FloatVector

RETURNS DESCRIPTION
FloatSquareMatrix

Matrix of interaction virial coefficients \(B_{ij}\). Unit = m³/mol.
Source code in src/polykin/thermo/eos/virial.py 
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
def B_mixture(T: float,
              Tc: FloatVector,
              Pc: FloatVector,
              Zc: FloatVector,
              w: FloatVector,
              ) -> FloatSquareMatrix:
    r"""Calculate the matrix of interaction virial coefficients using the
    mixing rules of Prausnitz.

    \begin{aligned}
        B_{ij} &= B(T,T_{cij},P_{cij},\omega_{ij}) \\
        v_{cij} &= \frac{(v_{ci}^{1/3}+v_{cj}^{1/3})^3}{8} \\
        k_{ij} &= 1 -\frac{\sqrt{v_{ci}v_{cj}}}{v_{cij}} \\
        T_{cij} &= \sqrt{T_{ci}T_{cj}}(1-k_{ij}) \\
        Z_{cij} &= \frac{Z_{ci}+Z_{cj}}{2} \\
        \omega_{ij} &= \frac{\omega_{i}+\omega_{j}}{2} \\
        P_{cij} &= \frac{Z_{cij} R T_{cij}}{v_{cij}}
    \end{aligned}

    The calculation of the individual coefficients is handled by
    [`B_pure`](.#polykin.properties.eos.B_pure).

    **References**

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

    Parameters
    ----------
    T : float
        Temperature. Unit = K.
    Tc : FloatVector
        Critical temperatures of all components. Unit = K.
    Pc : FloatVector
        Critical pressures of all components. Unit = Pa.
    Zc : FloatVector
        Critical compressibility factors of all components.
    w : FloatVector
        Acentric factors of all components.

    Returns
    -------
    FloatSquareMatrix
        Matrix of interaction virial coefficients $B_{ij}$. Unit = m³/mol.
    """
    vc = Zc*R*Tc/Pc
    N = Tc.size
    B = np.empty((N, N), dtype=np.float64)
    for i in range(N):
        for j in range(i, N):
            if i == j:
                B[i, j] = B_pure(T, Tc[i],  Pc[i], w[i])
            else:
                vcm = (vc[i]**(1/3) + vc[j]**(1/3))**3 / 8
                km = 1. - sqrt(vc[i]*vc[j])/vcm
                Tcm = sqrt(Tc[i]*Tc[j])*(1. - km)
                Zcm = (Zc[i] + Zc[j])/2
                wm = (w[i] + w[j])/2
                Pcm = Zcm*R*Tcm/vcm
                B[i, j] = B_pure(T, Tcm, Pcm, wm)
                B[j, i] = B[i, j]
    return B

Examples¤

Estimate the molar volume of a 50 mol% ethylene/nitrogen at 350 K and 10 bar.

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

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

eos = Virial(Tc, Pc, Zc, w)
v = eos.v(T=350., P=10e5, y=np.array([0.5, 0.5]))

print(f"{v:.2e} m³/mol")
2.87e-03 m³/mol