Fugacity of pure fluids and cubic equations of state

Next door to the fermenters, Lakeshore Gas compresses the plant's carbon dioxide to about 20 bar before drying and liquefying it, the unit Theo studied in Thermodynamics I. Its engineers want to know how much water the compressed gas will carry, and how much ethanol will condense in the coolers between stages. Both questions are phase equilibrium questions at a pressure where the ideal-gas expression \(\mu = G^\circ + RT\ln(P/P^\circ)\) is no longer accurate. Lewis's answer, in the early 1900s, was to keep the form of that equation and replace the pressure with a corrected pressure he called the fugacity.

Defining fugacity

For a pure ideal gas at constant temperature, \(d\mu = RT\,d\ln P\). For any real fluid, the fugacity \(f\) is defined so that the same form holds exactly:

The fugacity has units of pressure. The ratio \(\phi = f/P\) is the fugacity coefficient, which measures how far the fluid departs from ideal-gas behaviour: \(\phi = 1\) for an ideal gas, less than 1 when attractions dominate, as they do for most gases at moderate pressure, and greater than 1 when repulsions dominate at very high pressure. Because \(\mu\) and \(\ln f\) differ only by a function of temperature, equal chemical potentials mean equal fugacities. The criterion of phase equilibrium can therefore be written in a form that has no awkward infinities: for a pure substance at its vapour pressure, \(f^{\text{V}} = f^{\text{L}}\); for a species in a mixture, \(\hat{f}_i^{\,\alpha} = \hat{f}_i^{\,\beta}\) in every phase.

Fugacity from volumetric data

Integrating \(d\ln f = (V/RT)\,dP\) from zero pressure, and subtracting the ideal-gas result, gives the fugacity coefficient from the compressibility factor \(Z = PV/RT\):

So fugacity is not a new kind of measurement. Any reliable source of \(Z\) along an isotherm, whether tables, a virial equation, a generalized correlation or an equation of state, gives \(\phi\). For moderate pressures the truncated virial equation \(Z = 1 + BP/RT\) gives the simple result \(\ln\phi = BP/RT\). For carbon dioxide at 35 °C, \(B\) is about \(-112\) cm\(^3\)/mol, so at 20 bar \(\ln\phi \approx -0.087\) and \(\phi \approx 0.92\). At higher pressures, and for liquids, a cubic equation of state is the practical tool.

Cubic equations of state

Van der Waals proposed in 1873 that real gases differ from ideal ones because molecules have a finite size, represented by a constant \(b\), and attract one another, represented by a constant \(a\): \(P = RT/(V - b) - a/V^2\). It was the first equation to predict both liquid and vapour from one expression, and its descendants are still the workhorses of process simulation. The two most used today are the Soave-Redlich-Kwong equation (1972) and the Peng-Robinson equation (1976):

In the Peng-Robinson equation, \(a = 0.45724\,R^2T_c^2/P_c\), \(b = 0.07780\,RT_c/P_c\), and \(\alpha(T) = [1 + \kappa(1 - \sqrt{T/T_c})]^2\), with \(\kappa = 0.37464 + 1.54226\omega - 0.26992\omega^2\). The only inputs are the critical temperature, the critical pressure and the acentric factor \(\omega\), which describes how far the vapour pressure curve departs from that of a simple spherical molecule. For carbon dioxide these are about 304.1 K, 73.8 bar and 0.225.

Written in terms of \(Z\), with \(A = aP/(RT)^2\) and \(B = bP/RT\), every cubic equation becomes a cubic polynomial in \(Z\). At conditions where only one phase can exist it has one real root. In the two-phase region, below the critical temperature and near the vapour pressure, it has three real roots: the largest belongs to the vapour, the smallest to the liquid, and the middle one has no physical meaning. Choosing the right root is the first thing a careful calculation checks.

The fugacity coefficient from Peng-Robinson

Substituting the equation into the integral for \(\ln\phi\) gives a closed form:

Evaluated with the vapour root it gives the vapour fugacity; with the liquid root, the liquid fugacity. At the vapour pressure the two are equal, which is how simulators compute vapour pressures from an equation of state: they adjust \(P\) until \(\phi^{\text{L}} = \phi^{\text{V}}\).

# Fugacity of pure carbon dioxide and of compressed liquid water
import math
R = 8.314  # J/(mol K)

def pr_params(Tc, Pc, w, T):
    kappa = 0.37464 + 1.54226 * w - 0.26992 * w * w
    alpha = (1 + kappa * (1 - math.sqrt(T / Tc))) ** 2
    a = 0.45724 * R * R * Tc * Tc / Pc * alpha
    b = 0.07780 * R * Tc / Pc
    return a, b

def cubic_roots(c2, c1, c0):
    # real roots of Z^3 + c2 Z^2 + c1 Z + c0 = 0 (Cardano with trigonometric form)
    q = (3 * c1 - c2 * c2) / 9; r = (9 * c2 * c1 - 27 * c0 - 2 * c2 ** 3) / 54
    d = q ** 3 + r * r
    if d > 0:
        s = math.copysign(abs(r + math.sqrt(d)) ** (1 / 3), r + math.sqrt(d))
        t = math.copysign(abs(r - math.sqrt(d)) ** (1 / 3), r - math.sqrt(d))
        return [s + t - c2 / 3]
    th = math.acos(r / math.sqrt(-q ** 3))
    return sorted(2 * math.sqrt(-q) * math.cos((th + 2 * math.pi * k) / 3) - c2 / 3 for k in range(3))

def pr_Z(A, B):
    return cubic_roots(-(1 - B), A - 3 * B * B - 2 * B, -(A * B - B * B - B ** 3))

def pr_lnphi_pure(Z, A, B):
    s2 = math.sqrt(2)
    return Z - 1 - math.log(Z - B) - A / (2 * s2 * B) * math.log((Z + (1 + s2) * B) / (Z + (1 - s2) * B))
Tc, Pc, w = 304.13, 73.77e5, 0.225          # CO2 (approximate)
T = 308.15                                  # 35 C
print(" P/bar    Z      phi    f/bar")
for Pbar in (1, 5, 10, 20, 40, 60):
    P = Pbar * 1e5
    a, b = pr_params(Tc, Pc, w, T)
    A = a * P / (R * T) ** 2; B = b * P / (R * T)
    Z = max(pr_Z(A, B))                     # gas root
    phi = math.exp(pr_lnphi_pure(Z, A, B))
    print(f"{Pbar:5.0f}  {Z:6.4f}  {phi:6.4f}  {phi * Pbar:6.2f}")

# Compare three cubic equations for the gas Z at 20 bar, 35 C
P = 20e5
for name, Om_a, Om_b, use_alpha in (("van der Waals", 27 / 64, 1 / 8, False), ("Redlich-Kwong", 0.42748, 0.08664, True)):
    a = Om_a * R * R * Tc * Tc / Pc * ((Tc / T) ** 0.5 if use_alpha else 1.0)
    b = Om_b * R * Tc / Pc
    A = a * P / (R * T) ** 2; B = b * P / (R * T)
    if use_alpha:
        Z = max(cubic_roots(-1, A - B - B * B, -A * B))
    else:
        Z = max(cubic_roots(-(1 + B), A, -A * B))
    print(f"{name:14s}: Z = {Z:.4f}")
a, b = pr_params(Tc, Pc, w, T); A = a * P / (R * T) ** 2; B = b * P / (R * T)
print(f"{'Peng-Robinson':14s}: Z = {max(pr_Z(A, B)):.4f}")

# Liquid water at 80 C compressed to 20 bar: Poynting correction
Psat = 47.39e3; Vl = 18.5e-6; Tw = 353.15     # Pa, m3/mol, K
for Pbar in (1, 20, 100):
    poy = math.exp(Vl * (Pbar * 1e5 - Psat) / (R * Tw))
    print(f"liquid water at {Pbar:3d} bar: Poynting factor = {poy:.4f}, f = {Psat * poy / 1000:.2f} kPa")

 P/bar    Z      phi    f/bar
    1  0.9950  0.9951    1.00
    5  0.9750  0.9755    4.88
   10  0.9494  0.9513    9.51
   20  0.8962  0.9037   18.07
   40  0.7780  0.8119   32.48
   60  0.6313  0.7228   43.37
van der Waals : Z = 0.9163
Redlich-Kwong : Z = 0.9055
Peng-Robinson : Z = 0.8962
liquid water at   1 bar: Poynting factor = 1.0003, f = 47.41 kPa
liquid water at  20 bar: Poynting factor = 1.0124, f = 47.98 kPa
liquid water at 100 bar: Poynting factor = 1.0647, f = 50.46 kPa

Reading the output

At 20 bar and 35 °C, the Peng-Robinson equation gives \(Z = 0.896\), the value Theo found in Thermodynamics I, and \(\phi = 0.904\), so the fugacity is 18.1 bar rather than 20. At 60 bar the departure is much larger: \(\phi = 0.72\), a fugacity of 43 bar. The virial estimate of 0.92 at 20 bar is within 2 % of the cubic result, which is typical: the virial equation is good where \(Z\) is above about 0.9 and unreliable beyond. The three cubic equations agree on the trend but not the number: van der Waals gives \(Z = 0.916\), Redlich-Kwong 0.906 and Peng-Robinson 0.896. Measured data for carbon dioxide at these conditions are closest to Peng-Robinson, which is one reason it is the default for gas processing.

Fugacity of liquids and the Poynting factor

For a pure liquid below its critical point, the fugacity can be built in two steps. At its vapour pressure, the liquid has the same fugacity as the saturated vapour, \(f^{\text{sat}} = \phi^{\text{sat}}P^{\text{sat}}\). Compressing the liquid to a higher pressure raises its fugacity by the integral of \(V/RT\), which for a nearly incompressible liquid gives

The exponential is the Poynting factor. Because liquid molar volumes are small, it is close to 1 unless the pressure is very high: for water at 80 °C it is 1.0003 at 1 bar, 1.012 at 20 bar and 1.065 at 100 bar. At low pressure, \(\phi^{\text{sat}}\) is also near 1, so the fugacity of a pure liquid is simply its vapour pressure. That is the approximation hidden inside Raoult's law, which lesson 5 introduces. At the 20 bar of the Lakeshore knockout drum the Poynting factor adds about 1.4 % to the fugacity of liquid water, a correction lesson 9 uses when it calculates how much water the compressed gas can hold.

Choosing an equation of state

Cubic equations are fast and robust and represent hydrocarbons, light gases and their mixtures well. They are weaker for liquid densities, which they can miss by 5 to 15 % unless a volume translation is added, and for strongly polar or hydrogen-bonding substances such as water, ethanol and acetic acid, for which no simple \(a\) and \(b\) can capture the association. More elaborate equations exist: the Benedict-Webb-Rubin family and multiparameter reference equations for pure fluids, which reproduce measured data almost exactly, and the statistical association fluid theory (SAFT) family for associating molecules. In practice, simulators use reference equations for steam and refrigerants, cubic equations for gas processing, and activity coefficient models for polar liquids at low pressure. Lesson 12 returns to this choice.