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 \(|\Sigma y_i - 1| < 0.001\).
Property calculations#
All properties are converted to SI units to avoid confusions.
Given temperature
\(T = 37.78 ^\circ\mathrm{C} = 310.93 \ \mathrm{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:
\(\displaystyle P_{sat1}=\exp\left(A_{1}-\frac{B_{1}}{T+C_{1}}\right) = 283.85 \ \mathrm{bar} = 283.85 \times 10^5 \ \mathrm{Pa}\)
\(\displaystyle P_{sat2}=\exp\left(A_{2}-\frac{B_{2}}{T+C_{2}}\right) = 1.07 \ \mathrm{bar} = 1.07 \times 10^5 \ \mathrm{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:
\(\displaystyle P = \sum x_i P_{sat, i} = x_1 P_{sat1} + x_2 P_{sat2} = 85.9 \times 10^5 \ \mathrm{Pa}\)
Using Raoult’s law, we can estimate the vapor compositions:
\(\displaystyle y_{1}=\frac{x_{1}P_{sat1}}{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_{c1}=190.6 \times 10^5 \ \mathrm{K}\)
\(T_{c2}=469.6 \times 10^5 \ \mathrm{K}\)
\(P_{c1}=46 \times 10^5 \ \mathrm{Pa}\)
\(P_{c2}=33.74 \times 10^5 \ \mathrm{Pa}\)
we can calculate the vdW parameters of pure species
\(\displaystyle a_{1}=\frac{27}{64}\frac{R^{2}T_{c1}^{2}}{P_{c1}} = 0.2303 \ \mathrm{J \ m^3/mol^3}\)
\(\displaystyle a_{2}=\frac{27}{64}\frac{R^{2}T_{c2}^{2}}{P_{c2}} = 1.906 \ \mathrm{J \ m^3/mol^3}\)
\(\displaystyle b_{1}=\frac{RT_{c1}}{8P_{c1}} = 4.31 \times 10^{-5} \ \mathrm{m^3/mol}\)
\(\displaystyle b_{2}=\frac{RT_{c2}}{8P_{c2}} = 1.45 \times 10^{-4} \ \mathrm{m^3/mol}\)
and interaction term
\(a_{12} = \sqrt{a_1a_2} = 1.462 \ \mathrm{J \ m^3/mol^3}\)
where the ideal gas constant is \(R = 8.314 \ \mathrm{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)
\(\displaystyle a_{l}=x_{1}^{2}a_{1}+x_{2}^{2}a_{2}+2x_{1}x_{2}a_{12} = 1.23 \ \mathrm{J m^3/mol^3}\)
\(\displaystyle b_{l}=x_{1}b_{1}+x_{2}b_{2} = 1.14 \times 10^{-5} \ \mathrm{J m^3/mol^3}\)
Vapor phase (an estimate for first iteration, to be recalculated for each iteration)
\(\displaystyle a_{v}=y_{1}^{2}a_{1}+y_{2}^{2}a_{2}+2y_{1}y_{2}a_{12} = 0.238 \ \mathrm{J m^3/mol^3}\)
\(\displaystyle b_{v}=y_{1}b_{1}+y_{2}b_{2} = 4.40\times 10^{-5} \ \mathrm{m^3/mol}\)
vdW EOS calculations#
The vdW EOS
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} \ \mathrm{m^3/mol}\)
\(v_v = 2.55 \times 10^{-4} \ \mathrm{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:
with appropriate properties in each phase (vapor and liquid) for each species (methane and n-pentane).
\(\displaystyle \hat{\varphi}_1^v = \exp\left[\frac{b_{1}}{v_v-b_v}-\ln\left(\frac{\left(v_v-b_v\right)P}{RT}\right)-\frac{2\left( y_1 a_1 + y_2 a_{12} \right)}{RTv_v}\right] = 0.86\)
\(\displaystyle \hat{\varphi}_2^v = \exp\left[\frac{b_{2}}{v_v-b_v}-\ln\left(\frac{\left(v_v-b_v\right)P}{RT}\right)-\frac{2\left( y_2 a_2 + y_1 a_{12} \right)}{RTv_v}\right] = 0.37\)
\(\displaystyle \hat{\varphi}_1^l = \exp\left[\frac{b_{1}}{v_l-b_l}-\ln\left(\frac{\left(v_l-b_l\right)P}{RT}\right)-\frac{2\left( x_1 a_1 + x_2 a_{12} \right)}{RTv_l}\right] = 1.30\)
\(\displaystyle \hat{\varphi}_2^l = \exp\left[\frac{b_{2}}{v_l-b_l}-\ln\left(\frac{\left(v_l-b_l\right)P}{RT}\right)-\frac{2\left( x_2 a_2 + x_1 a_{12} \right)}{RTv_l}\right] = 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:
\(\displaystyle y_1 = \dfrac{x_1 \hat{\varphi}_1^l}{\hat{\varphi}_{1}^{v}} = 0.45\)
\(\displaystyle y_2 = \dfrac{x_2 \hat{\varphi}_2^l}{\hat{\varphi}_{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} \ \mathrm{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:
To solve Peng-Robinson EOS numerically, we use its cubic form:
where
\(z = \dfrac{Pv}{RT}\) is the compressibility factor
\(A = \dfrac{a\alpha P}{(RT)^2}\)
\(B = \dfrac{bP}{RT}\)
We calculate fugacity coefficients using
with careful substitution of \(a\alpha, 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{\varphi}_1^v, \hat{\varphi}_2^v, \hat{\varphi}_1^l, \hat{\varphi}_2^l\)
We calculate pure species Peng-Robinson parameters with
\(a = 0.45724 \dfrac{R^2 T_c^{2}}{P_c}\)
\(b = 0.07780 \dfrac{RT_c}{P_c}\)
\(\alpha(T) = [1 + \kappa(1 - \sqrt{T_r})]^2\)
\(\kappa = 0.37464 + 1.54226\omega - 0.26992\omega^2\)
\(\omega\) is the acentric factor of the species found in Appendix A.1
and calculate Peng-Robinson parameters of mixtures using mixing rules of
\(a\alpha = y_1^2 (a\alpha)_1 + 2y_1 y_1 \sqrt{(a\alpha)_1(a\alpha)_2} (1-k_{12}) + y_2^2 (a\alpha)_2\)
\(k\) is binary interaction parameter of the two species found in literature
\(b = y_1b_1 + y_1b_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 \(|\Sigma y_i - 1| < 0.001\).
Note
Code your own answer before checking the answer!
Show 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
Show code cell source
pd.set_option("display.precision", 2)
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.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 |