Determining Equilibrium Composition using Equation of State Method

Determining Equilibrium Composition using Equation of State Method#

Teng-Jui Lin

Content adapted from UW CHEME 326, Chemical Engineering Thermodynamics, in Winter 2022.

Problems reprinted from Koretsky, M. D. Engineering and Chemical Thermodynamics. (Wiley, 2012).

Using van der Waals EOS to predict VLE behavior#

Koretsky Example 8.15. Using the equation of state method, determine the equilibrium composition in the vapor phase and the system pressure of a mixture of methane (1) and n-pentane (2) with a liquid mole fraction of $x_{1} = 0.3$ at 37.78°C. Use the van der Waals equation of state to determine fugacity for both vapor and liquid phases. Iterate until $| Σ y_{i} - 1 | < 0.001$ .

Property calculations#

All properties are converted to SI units to avoid confusions.

Given temperature

$T = {37.78}^{\circ} C = 310.93 K$

and look up Antoine coefficients (Appendix A1)

Methane (1)
- $A_{1} = 8.6041$
- $B_{1} = 897.84$
- $C_{1} = - 7.16$
n-pentane (2)
- $A_{2} = 9.2131$
- $B_{2} = 2477.07$
- $C_{2} = - 39.94$

we can calculate vapor pressure of methane and n-pentane using Antoine’s equation:

$P_{s a t 1} = \exp (A_{1} - \frac{B_{1}}{T + C_{1}}) = 283.85 bar = 283.85 \times 10^{5} Pa$
$P_{s a t 2} = \exp (A_{2} - \frac{B_{2}}{T + C_{2}}) = 1.07 bar = 1.07 \times 10^{5} Pa$

Given Liquid compositions

$x_{1} = 0.3$
$x_{2} = 1 - x_{1} = 0.7$

we can calculate the pressure as the weighted sum of the saturation pressure:

$P = \sum x_{i} P_{s a t, i} = x_{1} P_{s a t 1} + x_{2} P_{s a t 2} = 85.9 \times 10^{5} Pa$

Using Raoult’s law, we can estimate the vapor compositions:

$y_{1} = \frac{x_{1} P_{s a t 1}}{P} = 0.991$
$y_{2} = 1 - y_{1} = 0.009$

vdW EOS property calculations#

The results from Raoult’s law is a rough estimate, and we want to get more accurate results using the vdW EOS. Recall that the vdW EOS requires two vdW parameters $a$ and $b$ .

Look up critical temperatures and pressures of methane and n-pentane

$T_{c 1} = 190.6 \times 10^{5} K$
$T_{c 2} = 469.6 \times 10^{5} K$
$P_{c 1} = 46 \times 10^{5} Pa$
$P_{c 2} = 33.74 \times 10^{5} Pa$

we can calculate the vdW parameters of pure species

$a_{1} = \frac{27}{64} \frac{R^{2} T_{c 1}^{2}}{P_{c 1}} = 0.2303 J m^{3} / {mol}^{3}$
$a_{2} = \frac{27}{64} \frac{R^{2} T_{c 2}^{2}}{P_{c 2}} = 1.906 J m^{3} / {mol}^{3}$
$b_{1} = \frac{R T_{c 1}}{8 P_{c 1}} = 4.31 \times 10^{- 5} m^{3} / mol$
$b_{2} = \frac{R T_{c 2}}{8 P_{c 2}} = 1.45 \times 10^{- 4} m^{3} / mol$

and interaction term

$a_{12} = \sqrt{a_{1} a_{2}} = 1.462 J m^{3} / {mol}^{3}$

where the ideal gas constant is $R = 8.314 J / mol K$ . We can then calculate vdW parameters of the mixed system liquid and vapor phase using mixing rules (note that we assumed the vdW EOS works well in liquid phase)

Liquid phase (we know for sure, same for every iteration)
- $a_{l} = x_{1}^{2} a_{1} + x_{2}^{2} a_{2} + 2 x_{1} x_{2} a_{12} = 1.23 J m^{3} / {mol}^{3}$
- $b_{l} = x_{1} b_{1} + x_{2} b_{2} = 1.14 \times 10^{- 5} J m^{3} / {mol}^{3}$
Vapor phase (an estimate for first iteration, to be recalculated for each iteration)
- $a_{v} = y_{1}^{2} a_{1} + y_{2}^{2} a_{2} + 2 y_{1} y_{2} a_{12} = 0.238 J m^{3} / {mol}^{3}$
- $b_{v} = y_{1} b_{1} + y_{2} b_{2} = 4.40 \times 10^{- 5} m^{3} / mol$

vdW EOS calculations#

The vdW EOS

P = \frac{R T}{v - b} - \frac{a}{v^{2}}

now has all the parameters. Solve for $v_{l}$ using $R, P, T, a_{l}, b_{l}$ and take the smallest root. Solve for $v_{v}$ using $R, P, T, a_{v}, b_{v}$ and take the largest root:

$v_{l} = 1.59 \times 10^{- 4} m^{3} / mol$
$v_{v} = 2.55 \times 10^{- 4} m^{3} / mol$

vdW fugacity coefficient calculations#

Knowing the molar volume of liquid and vapor phases, we can calculate the fugacity coefficient with equation from Koretsky Table 7.1:

\ln {\hat{φ}}_{i} = \frac{b_{i}}{v - b} - \ln (\frac{(v - b) P}{R T}) - \frac{2 \sum_{n = 1}^{m} y_{k} a_{i k}}{R T v}

with appropriate properties in each phase (vapor and liquid) for each species (methane and n-pentane).

${\hat{φ}}_{1}^{v} = \exp [\frac{b_{1}}{v_{v} - b_{v}} - \ln (\frac{(v_{v} - b_{v}) P}{R T}) - \frac{2 (y_{1} a_{1} + y_{2} a_{12})}{R T v_{v}}] = 0.86$
${\hat{φ}}_{2}^{v} = \exp [\frac{b_{2}}{v_{v} - b_{v}} - \ln (\frac{(v_{v} - b_{v}) P}{R T}) - \frac{2 (y_{2} a_{2} + y_{1} a_{12})}{R T v_{v}}] = 0.37$
${\hat{φ}}_{1}^{l} = \exp [\frac{b_{1}}{v_{l} - b_{l}} - \ln (\frac{(v_{l} - b_{l}) P}{R T}) - \frac{2 (x_{1} a_{1} + x_{2} a_{12})}{R T v_{l}}] = 1.30$
${\hat{φ}}_{2}^{l} = \exp [\frac{b_{2}}{v_{l} - b_{l}} - \ln (\frac{(v_{l} - b_{l}) P}{R T}) - \frac{2 (x_{2} a_{2} + x_{1} a_{12})}{R T v_{l}}] = 0.10$

Warning

The Koretsky Book has typos in the latter two equations. The third term should use $x_{1}$ and $x_{2}$ because they’re in liquid phase, but NOT $y_{1}$ and $y_{2}$ .

Vapor composition calculation#

We can calculate the vapor composition using fugacity coefficients:

$y_{1} = \frac{x_{1} {\hat{φ}}_{1}^{l}}{{\hat{φ}}_{1}^{v}} = 0.45$
$y_{2} = \frac{x_{2} {\hat{φ}}_{2}^{l}}{{\hat{φ}}_{2}^{v}} = 0.18$

Warning

The Koretsky Book has a typo in the $y_{2}$ equation: use $x_{2}$ instead of $x_{1}$ .

Note that the compositions does not sum to 1, indicating that further iteration is needed for the result to converge. We iterate by updating the pressure

$P \leftarrow P (y_{1} + y_{2})$

Iteration#

Repeat calculations including and following vdW EOS property calculation until the pressure and vapor compositions converge (no longer changes a lot from last iteration).

Iteration result#

After iteration, we get $P = 34.5 \times 10^{5} Pa$ and $y_{1} = 0.79$ .

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sympy import Symbol
from sympy.solvers import solve

# universal constant
R = 8.314  # J/mol K
# given properties
T = 37.78 + 273.15  # K
x1 = 0.3
x2 = 1 - x1

# formula for vdW params and saturation pressure
pure_a = lambda Tc, Pc, R: 27 / 64 * (R * Tc) ** 2 / Pc
pure_b = lambda Tc, Pc, R: 1 / 8 * (R * Tc) / Pc
antoine = lambda A, B, C, T: np.exp(A - B / (T + C)) * 1e5  # Pa

# physical properties - appendix a1
# methane (1)
A1 = 8.6041
B1 = 897.84
C1 = -7.16
Tc1 = 190.6  # K
Pc1 = 46.00e5  # Pa
a1 = pure_a(Tc1, Pc1, R)
b1 = pure_b(Tc1, Pc1, R)
# n-pentane (2)
A2 = 9.2131
B2 = 2477.07
C2 = -39.94
Tc2 = 469.6  # K
Pc2 = 33.74e5  # Pa
a2 = pure_a(Tc2, Pc2, R)
b2 = pure_b(Tc2, Pc2, R)
# mixture
a12 = np.sqrt(a1 * a2)

# saturation pressures
Psat1 = antoine(A1, B1, C1, T)
Psat2 = antoine(A2, B2, C2, T)
Psat1, Psat2

(28384910.01139723, 107495.95099231375)

# estimate total pressure
P = x1 * Psat1 + x2 * Psat2
P

8590720.169113789

# estimate vapor compositions
y1 = x1 * Psat1 / P
y2 = x2 * Psat2 / P
y1, y2

(0.9912408780389382, 0.008759121961061624)

def solve_v(T, P, a, b, R, v_before, threshold=1e-10):
    """Solve vdW EOS for v."""
    coeff = [-a * b, a, -(R * T + P * b), P]
    # solve for compressibility factor
    v = np.polynomial.polynomial.polyroots(coeff)
    # z is real and positive
    v = np.real(v)
    v = v[v > 0]
    # double check if solution is good in case of algorithmic instability
    error = np.abs(np.polynomial.polynomial.polyval(v, coeff))
    v = v[error < threshold]
    # use the root closest to the last v
    closest_idx = np.argmin(np.abs(v - v_before))
    v = v[closest_idx]
    return v

# formula for mixed vdW params
mixed_a = (
    lambda comp1, comp2, a1, a2: comp1**2 * a1
    + 2 * comp1 * comp2 * np.sqrt(a1 * a2)
    + comp2**2 * a2
)
mixed_b = lambda comp1, comp2, b1, b2: comp1 * b1 + comp2 * b2
# formula for fugacity coefficients (phi_1, for phi_2, switch arguments of 1 and 2 properties)
phi = lambda comp1, comp2, P, v, T, R, a1, a12, b, b1: np.exp(
    b1 / (v - b)
    - np.log((v - b) * P / (R * T))
    - 2 * (comp1 * a1 + comp2 * a12) / (R * T * v)
)

# initialize guess
v_v = 0
v_l = 0
# generate lists to store iteration results
vdw_a_v = []
vdw_b_v = []
vdw_P = []
vdw_y1 = []
vdw_y2 = []
vdw_v_v = []
vdw_v_l = []
vdw_phi_v1 = []
vdw_phi_v2 = []
vdw_phi_l1 = []
vdw_phi_l2 = []

# iteratre through vdW EOS
num_iter = 0
max_iter = 10000
tol = 0.001
tol_criteria = tol * 2
while tol_criteria > tol and num_iter < max_iter:
    # calculate vdW parameters for mixtures in each phase
    a_v = mixed_a(y1, y2, a1, a2)
    b_v = mixed_b(y1, y2, b1, b2)
    a_l = mixed_a(x1, x2, a1, a2)
    b_l = mixed_b(x1, x2, b1, b2)
    # solve vdW EOS for molar volume in each phase
    v_v = solve_v(T, P, a_v, b_v, R, v_before=v_v)
    v_l = solve_v(T, P, a_l, b_l, R, v_before=v_l)
    # calculate fugacity coefficient of each species in each phase
    phi_v1 = phi(y1, y2, P, v_v, T, R, a1, a12, b_v, b1)
    phi_v2 = phi(y2, y1, P, v_v, T, R, a2, a12, b_v, b2)
    phi_l1 = phi(x1, x2, P, v_l, T, R, a1, a12, b_l, b1)
    phi_l2 = phi(x2, x1, P, v_l, T, R, a2, a12, b_l, b2)
    # update vapor composition
    y1 = x1 * phi_l1 / phi_v1
    y2 = x2 * phi_l2 / phi_v2
    # update total pressure
    P = P * (y1 + y2)
    # calculate if calculation converges
    tol_criteria = abs(y1 + y2 - 1)
    num_iter += 1
    # store iteration info for debugging
    vdw_a_v.append(a_v)
    vdw_b_v.append(b_v)
    vdw_v_v.append(v_v)
    vdw_v_l.append(v_l)
    vdw_phi_v1.append(phi_v1)
    vdw_phi_v2.append(phi_v2)
    vdw_phi_l1.append(phi_l1)
    vdw_phi_l2.append(phi_l2)
    vdw_y1.append(y1)
    vdw_y2.append(y2)
    vdw_P.append(P)

pd.set_option("display.precision", 3)
pd.DataFrame({
    'vdw_a_v': vdw_a_v,
    'vdw_b_v': vdw_b_v,
    'vdw_v_v': vdw_v_v,
    'vdw_v_l': vdw_v_l,
    'vdw_phi_v1': vdw_phi_v1,
    'vdw_phi_v2': vdw_phi_v2,
    'vdw_phi_l1': vdw_phi_l1,
    'vdw_phi_l2': vdw_phi_l2,
    'vdw_y1': vdw_y1,
    'vdw_y2': vdw_y2,
    'vdw_P': vdw_P,
})

	vdw_a_v	vdw_b_v	vdw_v_v	vdw_v_l	vdw_phi_v1	vdw_phi_v2	vdw_phi_l1	vdw_phi_l2	vdw_y1	vdw_y2	vdw_P
0	0.238	4.395e-05	2.550e-04	1.595e-04	0.860	0.367	1.296	0.095	0.452	0.182	5.449e+06
1	0.219	4.580e-05	4.379e-04	1.653e-04	0.908	0.558	1.780	0.120	0.588	0.151	4.027e+06
2	0.241	4.715e-05	5.969e-04	1.689e-04	0.931	0.631	2.245	0.147	0.724	0.163	3.571e+06
3	0.328	5.476e-05	6.492e-04	1.702e-04	0.944	0.606	2.472	0.161	0.786	0.186	3.469e+06
4	0.401	6.070e-05	6.429e-04	1.705e-04	0.956	0.573	2.530	0.164	0.794	0.201	3.451e+06
5	0.433	6.323e-05	6.341e-04	1.706e-04	0.962	0.558	2.541	0.165	0.792	0.207	3.447e+06
6	0.443	6.403e-05	6.309e-04	1.706e-04	0.965	0.553	2.543	0.165	0.791	0.209	3.446e+06

print(
    f"Converged in {num_iter} iterations, with convergence {round(tol_criteria, 4)} < {tol}"
)
print(f"phi_v1 = {round(phi_v1, 2)}")
print(f"phi_v2 = {round(phi_v2, 2)}")
print(f"phi_l1 = {round(phi_l1, 2)}")
print(f"phi_l2 = {round(phi_l2, 2)}")
print(f"y1 = {round(y1, 2)}")
print(f"y2 = {round(y2, 2)}")
print(f"P = {round(P)} Pa")

Converged in 7 iterations, with convergence 0.0003 < 0.001
phi_v1 = 0.96
phi_v2 = 0.55
phi_l1 = 2.54
phi_l2 = 0.17
y1 = 0.79
y2 = 0.21
P = 3445993 Pa

Using Peng-Robinson EOS to predict VLE behavior#

Koretsky Example 8.16. Repeat Example 8.15 using the Peng–Robinson equation of state. Compare your answer to the reported measured value of P = 69.1 bar and $y_{1} = 0.95$ . (Use $k_{12} = 0$ )

The approach is the same with using vdW EOS, but instead we use Peng-Robinson EOS:

P = \frac{R T}{v - b} - \frac{a α}{v^{2} + 2 b v - b^{2}}

To solve Peng-Robinson EOS numerically, we use its cubic form:

z^{3} - (1 - B) z^{2} + (A - 2 B - 3 B^{2}) z - (A B - B^{2} - B^{3}) = 0

where

$z = \frac{P v}{R T}$ is the compressibility factor
$A = \frac{a α P}{(R T)^{2}}$
$B = \frac{b P}{R T}$

We calculate fugacity coefficients using

{\hat{φ}}_{i} = \exp [\frac{b_{i}}{b} (z - 1) - \ln (\frac{(v - b) P}{R T}) + \frac{a α}{2 \sqrt{2} b R T} [\frac{b_{i}}{b} - \frac{2}{a α} \sum_{k = 1}^{m} y_{k} (a α)_{i k}] \ln [\frac{v + (1 + \sqrt{2}) b}{v + (1 - \sqrt{2}) b}]]

with careful substitution of $a α, b$ noting the species and $v$ noting the phase. Use $x_{k}$ in the sum if calculating fugacity coefficient in the liquid phase.

Try it yourself: write out ${\hat{φ}}_{1}^{v}, {\hat{φ}}_{2}^{v}, {\hat{φ}}_{1}^{l}, {\hat{φ}}_{2}^{l}$

{\hat{φ}}_{1}^{v} = \exp [\frac{b_{1}}{b_{v}} (z - 1) - \ln (\frac{(v_{v} - b_{v}) P}{R T}) + \frac{(a α)_{v}}{2 \sqrt{2} b_{v} R T} [\frac{b_{1}}{b_{v}} - \frac{2}{(a α)_{v}} (y_{1} (a α)_{1} + y_{2} (a α)_{12})] \ln [\frac{v_{v} + (1 + \sqrt{2}) b_{v}}{v_{v} + (1 - \sqrt{2}) b_{v}}]]

{\hat{φ}}_{2}^{v} = \exp [\frac{b_{2}}{b_{v}} (z - 1) - \ln (\frac{(v_{v} - b_{v}) P}{R T}) + \frac{(a α)_{v}}{2 \sqrt{2} b_{v} R T} [\frac{b_{2}}{b_{v}} - \frac{2}{(a α)_{v}} (y_{2} (a α)_{2} + y_{1} (a α)_{12})] \ln [\frac{v_{v} + (1 + \sqrt{2}) b_{v}}{v_{v} + (1 - \sqrt{2}) b_{v}}]]

{\hat{φ}}_{1}^{l} = \exp [\frac{b_{1}}{b_{l}} (z - 1) - \ln (\frac{(v_{l} - b_{l}) P}{R T}) + \frac{(a α)_{l}}{2 \sqrt{2} b_{l} R T} [\frac{b_{1}}{b_{l}} - \frac{2}{(a α)_{l}} (x_{1} (a α)_{1} + x_{2} (a α)_{12})] \ln [\frac{v_{l} + (1 + \sqrt{2}) b_{l}}{v_{l} + (1 - \sqrt{2}) b_{l}}]]

{\hat{φ}}_{2}^{l} = \exp [\frac{b_{2}}{b_{l}} (z - 1) - \ln (\frac{(v_{l} - b_{l}) P}{R T}) + \frac{(a α)_{l}}{2 \sqrt{2} b_{l} R T} [\frac{b_{2}}{b_{l}} - \frac{2}{(a α)_{l}} (x_{2} (a α)_{2} + x_{1} (a α)_{12})] \ln [\frac{v_{l} + (1 + \sqrt{2}) b_{l}}{v_{l} + (1 - \sqrt{2}) b_{l}}]]

We calculate pure species Peng-Robinson parameters with

$a = 0.45724 \frac{R^{2} T_{c}^{2}}{P_{c}}$
$b = 0.07780 \frac{R T_{c}}{P_{c}}$
$α (T) = [1 + κ (1 - \sqrt{T_{r}})]^{2}$
- $κ = 0.37464 + 1.54226 ω - 0.26992 ω^{2}$
  - $ω$ is the acentric factor of the species found in Appendix A.1

and calculate Peng-Robinson parameters of mixtures using mixing rules of

$a α = y_{1}^{2} (a α)_{1} + 2 y_{1} y_{1} \sqrt{(a α)_{1} (a α)_{2}} (1 - k_{12}) + y_{2}^{2} (a α)_{2}$
- $k$ is binary interaction parameter of the two species found in literature
$b = y_{1} b_{1} + y_{1} b_{2}$

with appropriate use of $x$ with liquid phase and $y$ with vapor phase.

# universal constant
R = 8.314  # J/mol K
# given properties
T = 37.78 + 273.15  # K
x1 = 0.3
x2 = 1 - x1
k12 = 0

# formula for Peng-Robinson params and saturation pressure
pure_a = lambda Tc, Pc, R: 0.45724 * (R * Tc) ** 2 / Pc
pure_b = lambda Tc, Pc, R: 0.07780 * (R * Tc) / Pc
kappa = lambda omega: 0.37464 + 1.54226 * omega - 0.26992 * omega**2
alpha = lambda Tr, kappa: (1 + kappa * (1 - np.sqrt(Tr))) ** 2
antoine = lambda A, B, C, T: np.exp(A - B / (T + C)) * 1e5  # Pa

# physical properties - appendix a1
# methane (1)
A1 = 8.6041
B1 = 897.84
C1 = -7.16
omega1 = 0.008
Tc1 = 190.6  # K
Pc1 = 46.00e5  # Pa
Tr1 = T / Tc1
a1 = pure_a(Tc1, Pc1, R)
b1 = pure_b(Tc1, Pc1, R)
kappa1 = kappa(omega1)
alpha1 = alpha(Tr1, kappa1)
a_alpha_1 = a1 * alpha1
# n-pentane (2)
A2 = 9.2131
B2 = 2477.07
C2 = -39.94
omega2 = 0.251
Tc2 = 469.6  # K
Pc2 = 33.74e5  # Pa
Tr2 = T / Tc2
a2 = pure_a(Tc2, Pc2, R)
b2 = pure_b(Tc2, Pc2, R)
kappa2 = kappa(omega2)
alpha2 = alpha(Tr2, kappa2)
a_alpha_2 = a2 * alpha2
# mixture
a12 = np.sqrt(a1 * a2)
a_alpha_12 = np.sqrt(a_alpha_1 * a_alpha_2)

# saturation pressures
Psat1 = antoine(A1, B1, C1, T)
Psat2 = antoine(A2, B2, C2, T)
Psat1, Psat2

(28384910.01139723, 107495.95099231375)

# estimate total pressure
P = x1 * Psat1 + x2 * Psat2
P

8590720.169113789

# estimate vapor compositions
y1 = x1 * Psat1 / P
y2 = x2 * Psat2 / P
y1, y2

(0.9912408780389382, 0.008759121961061624)

def solve_v(T, P, a_alpha, b, R, v_before, threshold=1e-10):
    """Solve Peng-Robinson EOS for v."""
    A = a_alpha * P / (R * T) ** 2
    B = b * P / (R * T)
    coeff = [-(A * B - B**2 - B**3), A - 2 * B - 3 * B**2, -(1 - B), 1]
    # solve for compressibility factor
    z = np.polynomial.polynomial.polyroots(coeff)
    # z is real and positive
    z = np.real(z)
    z = z[z > 0]
    # double check if solution is good in case of algorithmic instability
    error = np.abs(np.polynomial.polynomial.polyval(z, coeff))
    z = z[error < threshold]
    # convert to v
    v = z * R * T / P
    v = v[v > b]
    # use the root closest to the last v
    closest_idx = np.argmin(np.abs(v - v_before))
    v = v[closest_idx]
    return v

# formula for mixed Peng-Robinson params
mixed_a_alpha = (
    lambda comp1, comp2, a_alpha_1, a_alpha_2, k12: comp1**2 * a_alpha_1
    + 2 * comp1 * comp2 * np.sqrt(a_alpha_1 * a_alpha_2) * (1 - k12)
    + comp2**2 * a_alpha_2
)
mixed_b = lambda comp1, comp2, b1, b2: comp1 * b1 + comp2 * b2
# formula for fugacity coefficients (phi_1, for phi_2, switch arguments of 1 and 2 properties)
phi = lambda comp1, comp2, P, v, T, R, a_alpha, a_alpha_1, a_alpha_12, b, b1: np.exp(
    b1 / b * (P * v / (R * T) - 1)
    - np.log((v - b) * P / (R * T))
    + a_alpha
    / (2 * np.sqrt(2) * b * R * T)
    * (b1 / b - 2 / a_alpha * (comp1 * a_alpha_1 + comp2 * a_alpha_12))
    * np.log((v + (1 + np.sqrt(2)) * b) / (v + (1 - np.sqrt(2)) * b))
)

# initialize guess
v_v = 0
v_l = 0
# generate lists to store iteration results
pr_a_alpha_v = []
pr_b_v = []
pr_P = []
pr_y1 = []
pr_y2 = []
pr_v_v = []
pr_v_l = []
pr_phi_v1 = []
pr_phi_v2 = []
pr_phi_l1 = []
pr_phi_l2 = []
pr_tol_criteria = []

# iteratre through Peng-Robinson EOS
num_iter = 0
max_iter = 10000
tol = 0.001
tol_criteria = tol * 2
while tol_criteria > tol and num_iter < max_iter:
    # calculate Peng-Robinson parameters for mixtures in each phase
    a_alpha_v = mixed_a_alpha(y1, y2, a_alpha_1, a_alpha_2, k12)
    b_v = mixed_b(y1, y2, b1, b2)
    a_alpha_l = mixed_a_alpha(x1, x2, a_alpha_1, a_alpha_2, k12)
    b_l = mixed_b(x1, x2, b1, b2)
    # solve Peng-Robinson EOS for molar volume in each phase
    v_v = solve_v(T, P, a_alpha_v, b_v, R, v_before=v_v)
    v_l = solve_v(T, P, a_alpha_l, b_l, R, v_before=v_l)
    # v_v = solve_v(T, P, a_alpha_v, b_v, R, phase='v')
    # v_l = solve_v(T, P, a_alpha_l, b_l, R, phase='l')
    # calculate fugacity coefficient of each species in each phase
    phi_v1 = phi(y1, y2, P, v_v, T, R, a_alpha_v, a_alpha_1, a_alpha_12, b_v, b1)
    phi_v2 = phi(y2, y1, P, v_v, T, R, a_alpha_v, a_alpha_2, a_alpha_12, b_v, b2)
    phi_l1 = phi(x1, x2, P, v_l, T, R, a_alpha_l, a_alpha_1, a_alpha_12, b_l, b1)
    phi_l2 = phi(x2, x1, P, v_l, T, R, a_alpha_l, a_alpha_2, a_alpha_12, b_l, b2)
    # update vapor composition
    y1 = x1 * phi_l1 / phi_v1
    y2 = x2 * phi_l2 / phi_v2
    # update total pressure
    P = P * (y1 + y2)
    # calculate if calculation converges
    tol_criteria = abs(y1 + y2 - 1)
    num_iter += 1
    # store iteration info for debugging
    pr_a_alpha_v.append(a_alpha_v)
    pr_b_v.append(b_v)
    pr_v_v.append(v_v)
    pr_v_l.append(v_l)
    pr_phi_v1.append(phi_v1)
    pr_phi_v2.append(phi_v2)
    pr_phi_l1.append(phi_l1)
    pr_phi_l2.append(phi_l2)
    pr_y1.append(y1)
    pr_y2.append(y2)
    pr_P.append(P)
    pr_tol_criteria.append(tol_criteria)

pd.set_option("display.precision", 2)
pd.DataFrame({
    'pr_a_alpha_v': pr_a_alpha_v,
    'pr_b_v': pr_b_v,
    'pr_v_v': pr_v_v,
    'pr_v_l': pr_v_l,
    'pr_phi_v1': pr_phi_v1,
    'pr_phi_v2': pr_phi_v2,
    'pr_phi_l1': pr_phi_l1,
    'pr_phi_l2': pr_phi_l2,
    'pr_y1': pr_y1,
    'pr_y2': pr_y2,
    'pr_P': pr_P,
    'pr_tol_criteria': pr_tol_criteria,
})

	pr_a_alpha_v	pr_b_v	pr_v_v	pr_v_l	pr_phi_v1	pr_phi_v2	pr_phi_l1	pr_phi_l2	pr_y1	pr_y2	pr_P	pr_tol_criteria
0	0.21	2.74e-05	2.58e-04	9.74e-05	0.86	0.27	2.21	0.02	0.77	0.05	7.03e+06	1.82e-01
1	0.18	2.48e-05	3.32e-04	9.82e-05	0.88	0.38	2.59	0.02	0.88	0.04	6.46e+06	8.00e-02
2	0.21	2.71e-05	3.54e-04	9.85e-05	0.89	0.37	2.77	0.02	0.94	0.04	6.31e+06	2.30e-02
3	0.24	2.88e-05	3.55e-04	9.86e-05	0.89	0.35	2.83	0.02	0.95	0.04	6.28e+06	5.74e-03
4	0.25	2.95e-05	3.53e-04	9.86e-05	0.90	0.34	2.84	0.02	0.95	0.05	6.27e+06	1.62e-03
5	0.25	2.97e-05	3.53e-04	9.86e-05	0.90	0.34	2.85	0.02	0.95	0.05	6.27e+06	5.04e-04

print(
    f"Converged in {num_iter} iterations, with convergence {round(tol_criteria, 4)} < {tol}"
)
print(f"phi_v1 = {round(phi_v1, 2)}")
print(f"phi_v2 = {round(phi_v2, 2)}")
print(f"phi_l1 = {round(phi_l1, 2)}")
print(f"phi_l2 = {round(phi_l2, 2)}")
print(f"y1 = {round(y1, 2)}")
print(f"y2 = {round(y2, 2)}")
print(f"P = {round(P)} Pa")

Converged in 6 iterations, with convergence 0.0005 < 0.001
phi_v1 = 0.9
phi_v2 = 0.34
phi_l1 = 2.85
phi_l2 = 0.02
y1 = 0.95
y2 = 0.05
P = 6265035 Pa

Practice Problem#

Koretsky Problem 8.61. Using the equation of state method, determine the equilibrium composition in the vapor phase and the system pressure of a mixture of methane (1) and n-pentane (2) with a liquid mole fraction of $x_{1} = 0.1$ at 60°C. Use the van der Waals equation of state to determine fugacity for both vapor and liquid phases. Iterate until $| Σ y_{i} - 1 | < 0.001$ .

Note

Code your own answer before checking the answer!

Show code cell source Hide code cell source

# universal constant
R = 8.314  # J/mol K
# given properties
T = 60 + 273.15  # K
x1 = 0.2
x2 = 1 - x1

# formula for vdW params and saturation pressure
pure_a = lambda Tc, Pc, R: 27 / 64 * (R * Tc) ** 2 / Pc
pure_b = lambda Tc, Pc, R: 1 / 8 * (R * Tc) / Pc
antoine = lambda A, B, C, T: np.exp(A - B / (T + C)) * 1e5  # Pa

# physical properties - appendix a1
# methane (1)
A1 = 8.6041
B1 = 897.84
C1 = -7.16
Tc1 = 190.6  # K
Pc1 = 46.00e5  # Pa
a1 = pure_a(Tc1, Pc1, R)
b1 = pure_b(Tc1, Pc1, R)
# n-pentane (2)
A2 = 9.2131
B2 = 2477.07
C2 = -39.94
Tc2 = 469.6  # K
Pc2 = 33.74e5  # Pa
a2 = pure_a(Tc2, Pc2, R)
b2 = pure_b(Tc2, Pc2, R)
# mixture
a12 = np.sqrt(a1 * a2)

# saturation pressures
Psat1 = antoine(A1, B1, C1, T)
Psat2 = antoine(A2, B2, C2, T)
# estimate total pressure
P = x1 * Psat1 + x2 * Psat2
# estimate vapor compositions
y1 = x1 * Psat1 / P
y2 = x2 * Psat2 / P


def solve_v(T, P, a, b, R, v_before, threshold=1e-10):
    """Solve vdW EOS for v."""
    coeff = [-a * b, a, -(R * T + P * b), P]
    # solve for compressibility factor
    v = np.polynomial.polynomial.polyroots(coeff)
    # z is real and positive
    v = np.real(v)
    v = v[v > 0]
    # double check if solution is good in case of algorithmic instability
    error = np.abs(np.polynomial.polynomial.polyval(v, coeff))
    v = v[error < threshold]
    # use the root closest to the last v
    closest_idx = np.argmin(np.abs(v - v_before))
    v = v[closest_idx]
    return v


# formula for mixed vdW params
mixed_a = (
    lambda comp1, comp2, a1, a2: comp1**2 * a1
    + 2 * comp1 * comp2 * np.sqrt(a1 * a2)
    + comp2**2 * a2
)
mixed_b = lambda comp1, comp2, b1, b2: comp1 * b1 + comp2 * b2
# formula for fugacity coefficients (phi_1, for phi_2, switch arguments of 1 and 2 properties)
phi = lambda comp1, comp2, P, v, T, R, a1, a12, b, b1: np.exp(
    b1 / (v - b)
    - np.log((v - b) * P / (R * T))
    - 2 * (comp1 * a1 + comp2 * a12) / (R * T * v)
)

# initialize guess
v_v = 0
v_l = 0

# generate lists to store iteration results
vdw_a_v = []
vdw_b_v = []
vdw_P = []
vdw_y1 = []
vdw_y2 = []
vdw_v_v = []
vdw_v_l = []
vdw_phi_v1 = []
vdw_phi_v2 = []
vdw_phi_l1 = []
vdw_phi_l2 = []

# iteratre through vdW EOS
num_iter = 0
max_iter = 10000
tol = 0.001
tol_criteria = tol * 2
while tol_criteria > tol and num_iter < max_iter:
    # calculate vdW parameters for mixtures in each phase
    a_v = mixed_a(y1, y2, a1, a2)
    b_v = mixed_b(y1, y2, b1, b2)
    a_l = mixed_a(x1, x2, a1, a2)
    b_l = mixed_b(x1, x2, b1, b2)
    # solve vdW EOS for molar volume in each phase
    v_v = solve_v(T, P, a_v, b_v, R, v_before=v_v)
    v_l = solve_v(T, P, a_l, b_l, R, v_before=v_l)
    # calculate fugacity coefficient of each species in each phase
    phi_v1 = phi(y1, y2, P, v_v, T, R, a1, a12, b_v, b1)
    phi_v2 = phi(y2, y1, P, v_v, T, R, a2, a12, b_v, b2)
    phi_l1 = phi(x1, x2, P, v_l, T, R, a1, a12, b_l, b1)
    phi_l2 = phi(x2, x1, P, v_l, T, R, a2, a12, b_l, b2)
    # update vapor composition
    y1 = x1 * phi_l1 / phi_v1
    y2 = x2 * phi_l2 / phi_v2
    # update total pressure
    P = P * (y1 + y2)
    # calculate if calculation converges
    tol_criteria = abs(y1 + y2 - 1)
    num_iter += 1
    # store iteration info for debugging
    vdw_a_v.append(a_v)
    vdw_b_v.append(b_v)
    vdw_v_v.append(v_v)
    vdw_v_l.append(v_l)
    vdw_phi_v1.append(phi_v1)
    vdw_phi_v2.append(phi_v2)
    vdw_phi_l1.append(phi_l1)
    vdw_phi_l2.append(phi_l2)
    vdw_y1.append(y1)
    vdw_y2.append(y2)
    vdw_P.append(P)

# report results
print(
    f"Converged in {num_iter} iterations, with convergence {round(tol_criteria, 4)} < {tol}"
)
print(f"phi_v1 = {round(phi_v1, 2)}")
print(f"phi_v2 = {round(phi_v2, 2)}")
print(f"phi_l1 = {round(phi_l1, 2)}")
print(f"phi_l2 = {round(phi_l2, 2)}")
print(f"y1 = {round(y1, 2)}")
print(f"y2 = {round(y2, 2)}")
print(f"P = {round(P)} Pa")

Converged in 5 iterations, with convergence 0.0 < 0.001
phi_v1 = 1.0
phi_v2 = 0.64
phi_l1 = 3.35
phi_l2 = 0.26
y1 = 0.67
y2 = 0.33
P = 2802522 Pa

	vdw_a_v	vdw_b_v	vdw_v_v	vdw_v_l	vdw_phi_v1	vdw_phi_v2	vdw_phi_l1	vdw_phi_l2	vdw_y1	vdw_y2	vdw_P
0	0.25	4.55e-05	3.46e-04	1.76e-04	0.90	0.49	1.63	0.14	0.36	0.22	4.16e+06
1	0.23	4.79e-05	6.31e-04	1.83e-04	0.94	0.68	2.43	0.19	0.52	0.23	3.09e+06
2	0.31	5.49e-05	8.38e-04	1.87e-04	0.96	0.71	3.09	0.24	0.64	0.27	2.83e+06
3	0.47	6.72e-05	8.68e-04	1.88e-04	0.98	0.67	3.32	0.26	0.68	0.31	2.80e+06
4	0.57	7.41e-05	8.43e-04	1.88e-04	1.00	0.64	3.35	0.26	0.67	0.33	2.80e+06