6.2 Cubic Equations of State (vdW, RK, SRK, and Peng-Robinson)

6.2 Cubic Equations of State (vdW, RK, SRK, and Peng-Robinson)

In modern chemical engineering process design, cubic equations of state (EOS) are the workhorses of process simulators (such as Aspen Plus, HYSYS, and Pro/II) and flash calculation algorithms. On the NCEES PE Chemical Exam, questions test your conceptual mastery of cubic equation formulations, parameter dependencies on critical constants ($T_c, P_c, \omega$), root interpretation in two-phase regimes, and mixing rules for multi-component fluids.

A cubic equation of state is the lowest-order polynomial equation in molar volume ($v$) or compressibility factor ($Z$) capable of describing both the liquid phase and the vapor phase, as well as the continuity between them through the critical point.

---

1. Mathematical Architecture of Cubic Equations of State

All two-parameter cubic equations of state can be expressed in a unified pressure-explicit form:

$$P = \frac{R T}{v - b} - \frac{a(T)}{(v + \epsilon b)(v + \sigma b)}$$

Where:

• $\frac{R T}{v - b}$ is the repulsive pressure term, representing hard-sphere molecular collisions and molecular co-volume ($b$).
• $\frac{a(T)}{(v + \epsilon b)(v + \sigma b)}$ is the attractive pressure term, representing intermolecular cohesive forces ($a$).
• $b$ = effective molar volume of the molecules themselves ($\text{m}^3/\text{mol}$ or $\text{ft}^3/\text{lbmol}$).
• $a(T)$ = measure of intermolecular attractive force strength ($\text{Pa}\cdot\text{m}^6/\text{mol}^2$ or $\text{psia}\cdot\text{ft}^6/\text{lbmol}^2$).
• $\epsilon, \sigma$ = pure numerical constants specific to each equation family.

The Critical Point Inflection Constraints

To determine the parameters $a$ and $b$ for any cubic equation of state, we apply the mathematical conditions of the critical isotherm at the critical point ($T = T_c, P = P_c, v = v_c$). At this unique point, the liquid and vapor phases become indistinguishable, producing a horizontal inflection point:

$$\left( \frac{\partial P}{\partial v} \right)_{T = T_c} = 0$$ $$\left( \frac{\partial^2 P}{\partial v^2} \right)_{T = T_c} = 0$$

Solving these two simultaneous differential constraints alongside $P(T_c, v_c) = P_c$ uniquely determines the parameter scaling constants $\Omega_a$ and $\Omega_b$, as well as the theoretical critical compressibility factor $Z_c = P_c v_c / (R T_c)$.

---

2. Systematic Derivation: From van der Waals to Peng-Robinson

1. van der Waals (vdW, 1873)

The historical foundation of real fluid thermodynamics:

$$P = \frac{R T}{v - b} - \frac{a}{v^2} \quad (\epsilon = 0, \sigma = 0)$$

Applying the critical inflection conditions yields:

$$a = \frac{27 R^2 T_c^2}{64 P_c} = 0.42188 \frac{R^2 T_c^2}{P_c}$$ $$b = \frac{R T_c}{8 P_c} = 0.12500 \frac{R T_c}{P_c}$$ $$Z_c = \frac{3}{8} = 0.37500$$

• Flaws: The universal critical compressibility of $Z_c = 0.375$ is far too high for real substances (hydrocarbons typically have $Z_c \approx 0.26 - 0.29$; water is $0.229$). Consequently, vdW severely underpredicts liquid densities (overpredicts liquid volumes by $30-50\%$) and cannot accurately match pure-component vapor pressures across wide temperature ranges.

2. Redlich-Kwong (RK, 1949)

Otto Redlich and Joseph N. S. Kwong introduced a temperature-dependent attraction term proportional to $T^{-0.5}$ and modified the volume denominator:

$$P = \frac{R T}{v - b} - \frac{a}{\sqrt{T} v(v + b)} \quad (\epsilon = 0, \sigma = 1)$$ $$a = 0.42748 \frac{R^2 T_c^{2.5}}{P_c}$$ $$b = 0.08664 \frac{R T_c}{P_c}$$ $$Z_c = \frac{1}{3} \approx 0.33333$$

• Significance: RK significantly improved vapor phase properties and gas-mixture PVT predictions over vdW. However, it still overpredicts $Z_c$ ($0.333$ vs actual $0.27$) and remains inadequate for liquid-phase equilibria.

3. Soave-Redlich-Kwong (SRK, 1972)

Giorgio Soave replaced the $1/\sqrt{T}$ term in Redlich-Kwong with a generalized temperature-dependent function $\alpha(T_r, \omega)$ calibrated to reproduce experimental pure-component vapor pressures at $T_r = 0.70$:

$$P = \frac{R T}{v - b} - \frac{a(T)}{v(v + b)} \quad (\epsilon = 0, \sigma = 1)$$ $$a(T) = a_c \alpha(T_r, \omega) = 0.42748 \frac{R^2 T_c^2}{P_c} \alpha(T_r, \omega)$$ $$b = 0.08664 \frac{R T_c}{P_c}$$ $$\alpha(T_r, \omega) = \left[ 1 + m \left( 1 - \sqrt{T_r} \right) \right]^2$$ $$m = 0.480 + 1.574\omega - 0.176\omega^2$$ $$Z_c = \frac{1}{3} \approx 0.33333$$

• Significance: SRK marked the birth of modern industrial vapor-liquid equilibrium (VLE) modeling for hydrocarbon processing. It predicts vapor pressures with high fidelity. However, because $Z_c$ remains fixed at $0.333$, liquid densities are still underestimated by $10-15\%$.

4. Peng-Robinson (PR, 1976)

Ding-Yu Peng and Donald B. Robinson modified the attractive term denominator to yield a lower, more realistic critical compressibility factor ($Z_c = 0.307$):

$$P = \frac{R T}{v - b} - \frac{a(T)}{v(v + b) + b(v - b)} = \frac{R T}{v - b} - \frac{a(T)}{v^2 + 2bv - b^2} \quad (\epsilon = 1 - \sqrt{2}, \sigma = 1 + \sqrt{2})$$ $$a(T) = a_c \alpha(T_r, \omega) = 0.45724 \frac{R^2 T_c^2}{P_c} \alpha(T_r, \omega)$$ $$b = 0.07780 \frac{R T_c}{P_c}$$ $$\alpha(T_r, \omega) = \left[ 1 + \kappa \left( 1 - \sqrt{T_r} \right) \right]^2$$ $$\kappa = 0.37464 + 1.54226\omega - 0.26992\omega^2 \quad (\text{for } \omega \le 0.49)$$ $$Z_c = 0.30740$$

• Significance: Peng-Robinson is the most widely applied equation of state in petroleum refining, natural gas processing, and chemical plant simulations. Its predicted liquid molar volumes are substantially closer to experimental values than those from SRK.

---

3. Comparison Table: Performance Metrics of Major Cubic Equations

---

4. Solving the Cubic Polynomial in Compressibility Factor $Z$

To solve for molar volume or compressibility factor on a computer or calculator, the equation of state is rearranged into a dimensionless cubic polynomial in $Z$ ($Z = P v / R T$).

Define the dimensionless equation parameters:

$$A \equiv \frac{a P}{(R T)^2} = \Omega_a \frac{P_r}{T_r^2} \alpha(T_r, \omega)$$ $$B \equiv \frac{b P}{R T} = \Omega_b \frac{P_r}{T_r}$$

Substituting $v = Z R T / P$ into each equation yields the standard cubic polynomials:

Cubic Form for van der Waals:

$$Z^3 - (1 + B)Z^2 + A Z - A B = 0$$

Cubic Form for SRK:

$$Z^3 - Z^2 + (A - B - B^2)Z - A B = 0$$

Cubic Form for Peng-Robinson:

$$Z^3 - (1 - B)Z^2 + (A - 3B^2 - 2B)Z - (A B - B^2 - B^3) = 0$$

---

5. Root Selection & The Three-Root Phenomenon

When solving the cubic polynomial for a pure substance or mixture at given $T$ and $P$:

• Single Real Root Regime: At temperatures above the critical point ($T > T_c$) or at very high pressures, the polynomial yields one real root and two complex conjugate roots. The single real root represents the fluid state (supercritical fluid, dense gas, or compressed liquid).
• Three Real Roots Regime: In the subcritical two-phase region ($T < T_c$) at pressures corresponding to the two-phase envelope, the polynomial yields three real roots:

$$Z_L < Z_M < Z_V$$

1. Largest Root ($Z_V$): Corresponds to the vapor phase (highest molar volume, lowest density). This root is selected for superheated vapor, saturated vapor, or flash vapor streams.
2. Smallest Positive Root ($Z_L$): Corresponds to the liquid phase (lowest molar volume, highest density). This root is selected for subcooled liquid, saturated liquid, or flash liquid streams.
3. Intermediate Middle Root ($Z_M$): Has no physical reality. It corresponds to the mechanically unstable region of the isotherm where $(\partial P / \partial v)_T > 0$. In this unphysical region, an increase in pressure would cause an increase in volume, violating the Second Law of Thermodynamics. It is discarded in engineering calculations.

   P ^
     |         Subcritical Isotherm
     |           /---\ (Z_L)
     |          /     \               /---
 P_sat-------(L)-------(M)-----------(V)-------
     |        /           \         /
     |       /             \-------/ (Z_V)
     +----------------------------------------> v
             v_liquid     v_unstable  v_vapor

---

6. Multicomponent Mixing Rules: van der Waals One-Fluid (vdW1f)

To apply cubic equations of state to multicomponent gas or liquid mixtures, effective mixture parameters $a_{mix}$ and $b_{mix}$ must be calculated from the pure-component parameters using mixing rules.

The standard industrial formulation is the van der Waals One-Fluid (vdW1f) mixing rules:

Co-Volume Parameter ($b_{mix}$)

Because molecular volumes are essentially additive:

$$b_{mix} = \sum_{i=1}^C x_i b_i$$

Attractive Energy Parameter ($a_{mix}$)

Because attractive interactions occur between pairs of colliding molecules ($i-j$ pairs):

$$a_{mix} = \sum_{i=1}^C \sum_{j=1}^C x_i x_j a_{ij}$$

Where $x_i$ is the mole fraction of component $i$ (liquid or vapor), and $a_{ij}$ is the cross-interaction parameter between dissimilar species $i$ and $j$.

Combining Rule & The Binary Interaction Parameter ($k_{ij}$)

The cross-parameter $a_{ij}$ is evaluated using the geometric mean combining rule, modified by an empirical binary interaction parameter ($k_{ij}$):

$$a_{ij} = \sqrt{a_i a_j} \left( 1 - k_{ij} \right)$$

• $k_{ij} = k_{ji}$ (symmetry).
• $k_{ii} = 0$ (interaction of a molecule with itself).
• Physical Meaning: $k_{ij}$ accounts for differences in molecular size, polarity, and chemical nature between dissimilar molecules.
  • For chemically similar hydrocarbon pairs (e.g., methane + ethane, propane + butane): $k_{ij} \approx 0.00$.
  • For hydrocarbon + non-hydrocarbon pairs (e.g., $\text{CH}_4 + \text{CO}_2$, $\text{CH}_4 + \text{H}_2\text{S}$, $\text{CH}_4 + \text{N}_2$): $k_{ij}$ typically ranges from $0.05$ to $0.15$.
  • A positive $k_{ij}$ decreases the cross-attraction $a_{ij}$ below the geometric mean, modeling positive deviations from ideal solution behavior.

---

7. Comprehensive Worked Numerical Example: Peng-Robinson Molar Volume and Vapor Density of Propane

Problem Statement

A chemical reactor feed drum stores pure propane ($\text{C}_3\text{H}_8$) vapor at a temperature of $T = 320.0\text{ K}$ ($46.85^\circ\text{C}$) and an operating pressure of $P = 1.50\text{ MPa}$ ($15.0\text{ bar} = 1.50 \times 10^6\text{ Pa}$).

Physical Properties for Propane ($MW = 44.097\text{ kg/kmol}$):

• Critical temperature: $T_c = 369.8\text{ K}$
• Critical pressure: $P_c = 4.248\text{ MPa} = 4.248 \times 10^6\text{ Pa}$
• Acentric factor: $\omega = 0.152$

Calculate:

1. The reduced temperature $T_r$ and reduced pressure $P_r$.
2. The Peng-Robinson parameters $\kappa$, $\alpha(T_r, \omega)$, $a(T)$, and $b$.
3. The dimensionless coefficients $A$ and $B$.
4. Formulate the cubic polynomial in $Z$ and solve for the vapor compressibility factor $Z_V$.
5. Calculate the molar volume $v$ ($\text{m}^3/\text{mol}$ and $\text{L/mol}$) and mass density $\rho$ ($\text{kg/m}^3$).
6. Compare the PR result with the ideal gas law prediction.

---

Step 1: Reduced Coordinates

$$T_r = \frac{T}{T_c} = \frac{320.0\text{ K}}{369.8\text{ K}} = \mathbf{0.8653}$$ $$P_r = \frac{P}{P_c} = \frac{1.50\text{ MPa}}{4.248\text{ MPa}} = \mathbf{0.3531}$$

---

Step 2: Peng-Robinson Parameters

Evaluate $\kappa$:

$$\kappa = 0.37464 + 1.54226\omega - 0.26992\omega^2$$ $$\kappa = 0.37464 + 1.54226(0.152) - 0.26992(0.152)^2 = 0.37464 + 0.23442 - 0.00624 = \mathbf{0.6028}$$

Evaluate $\alpha(T_r, \omega)$:

$$\alpha = \left[ 1 + \kappa \left( 1 - \sqrt{T_r} \right) \right]^2 = \left[ 1 + 0.6028 \left( 1 - \sqrt{0.8653} \right) \right]^2$$ $$\sqrt{0.8653} = 0.93022 \implies 1 - \sqrt{T_r} = 0.06978$$ $$\alpha = \left[ 1 + 0.6028(0.06978) \right]^2 = \left[ 1 + 0.04206 \right]^2 = (1.04206)^2 = \mathbf{1.0859}$$

Evaluate pure component attraction parameter $a(T)$:

$$a_c = 0.45724 \frac{R^2 T_c^2}{P_c} = 0.45724 \frac{(8.31447)^2 (369.8)^2}{4.248 \times 10^6} = 0.45724 \frac{(69.1304)(136,752)}{4.248 \times 10^6} = \mathbf{1.0175\text{ Pa}\cdot\text{m}^6/\text{mol}^2}$$ $$a(T) = a_c \cdot \alpha = 1.0175 \times 1.0859 = \mathbf{1.1049\text{ Pa}\cdot\text{m}^6/\text{mol}^2}$$

Evaluate pure component co-volume parameter $b$:

$$b = 0.07780 \frac{R T_c}{P_c} = 0.07780 \frac{(8.31447)(369.8)}{4.248 \times 10^6} = \mathbf{5.633 \times 10^{-5}\text{ m}^3/\text{mol}}$$

---

Step 3: Dimensionless Coefficients $A$ and $B$

$$A = \frac{a(T) P}{(R T)^2} = \frac{(1.1049) \times (1.50 \times 10^6)}{[(8.31447)(320.0)]^2} = \frac{1.65735 \times 10^6}{(2660.63)^2} = \frac{1.65735 \times 10^6}{7.07895 \times 10^6} = \mathbf{0.23412}$$ $$B = \frac{b P}{R T} = \frac{(5.633 \times 10^{-5}) \times (1.50 \times 10^6)}{(8.31447)(320.0)} = \frac{84.495}{2660.63} = \mathbf{0.031758}$$

---

Step 4: Formulation and Solution of the PR Cubic Polynomial

The Peng-Robinson cubic polynomial is:

$$Z^3 - (1 - B)Z^2 + (A - 3B^2 - 2B)Z - (A B - B^2 - B^3) = 0$$

Calculate coefficients:

• Coefficient of $Z^2$: $-(1 - B) = -(1 - 0.031758) = \mathbf{-0.96824}$
• Coefficient of $Z$: $A - 3B^2 - 2B = 0.23412 - 3(0.031758)^2 - 2(0.031758) = 0.23412 - 0.00303 - 0.06352 = \mathbf{+0.16757}$
• Constant term: $-(A B - B^2 - B^3) = -[(0.23412)(0.031758) - (0.031758)^2 - (0.031758)^3] = -[0.007435 - 0.001009 - 0.000032] = \mathbf{-0.006394}$

The polynomial is:

$$f(Z) = Z^3 - 0.96824 Z^2 + 0.16757 Z - 0.006394 = 0$$

Because the propane is superheated vapor at $T = 320\text{ K}$ and $P = 1.5\text{ MPa}$ (propane saturation pressure at $320\text{ K}$ is $P^{sat} \approx 1.68\text{ MPa} > 1.50\text{ MPa}$), we seek the largest root $Z_V$.

Iterating via Newton-Raphson or substitution:

• For $Z = 0.76$: $$f(0.76) = (0.76)^3 - 0.96824(0.76)^2 + 0.16757(0.76) - 0.006394$$ $$f(0.76) = 0.43898 - 0.55925 + 0.12735 - 0.00639 = +0.00069$$
• For $Z = 0.758$: $$f(0.758) = 0.43552 - 0.55631 + 0.12702 - 0.00639 = -0.00016$$

Interpolating yields the vapor compressibility factor:

$$Z_V = \mathbf{0.7583} \approx \mathbf{0.758}$$

---

Step 5: Molar Volume and Mass Density

Actual molar volume:

$$v = \frac{Z_V R T}{P} = \frac{(0.7583) \times (8.31447\text{ J/(mol}\cdot\text{K)}) \times (300.0\text{ K} \times 320.0/300.0)}{1.50 \times 10^6\text{ Pa}}$$ $$v = \frac{(0.7583) \times (2660.63)}{1.50 \times 10^6} = \mathbf{1.345 \times 10^{-3}\text{ m}^3/\text{mol}} = \mathbf{1.345\text{ L/mol}}$$

Mass density:

$$\rho = \frac{MW}{v} = \frac{0.044097\text{ kg/mol}}{0.001345\text{ m}^3/\text{mol}} = \mathbf{32.79\text{ kg/m}^3}$$

---

Step 6: Comparison with Ideal Gas Prediction

Ideal gas molar volume:

$$v^{ig} = \frac{R T}{P} = \frac{2660.63}{1.50 \times 10^6} = 1.774 \times 10^{-3}\text{ m}^3/\text{mol} = 1.774\text{ L/mol}$$ $$\rho^{ig} = \frac{MW}{v^{ig}} = \frac{0.044097}{0.001774} = 24.86\text{ kg/m}^3$$ $$\text{Ideal Gas Density Error} = \frac{24.86 - 32.79}{32.79} \times 100\% = \mathbf{-24.18\%}$$

Assuming an ideal gas underpredicts the true mass of propane in the drum by over $24\%$, creating a severe risk of vessel overfilling or mass balance accounting failure.

---

8. Critical PE Exam Traps & Pitfalls

Trap 1: Root Inversion (Selecting the Liquid Root for Vapor Flow)
When a cubic equation solver produces three roots, always double-check the physical phase. The largest root ($Z_V$) represents vapor, while the smallest root ($Z_L$) represents liquid. If a problem states "superheated vapor" or asks for the volume of gas exiting a reboiler, selecting $Z_L$ (which is typically around $0.03-0.08$) will produce an answer off by a factor of 10 to 30.

Trap 2: Mixing Up Equation-Specific Parameters
Never mix parameters across equations! Using $\Omega_a = 0.42748$ (from SRK) inside the Peng-Robinson polynomial, or using Peng-Robinson's $\kappa$ formula with the SRK denominator, will produce completely corrupted results. Ensure you match the numerical constants exactly as given in the reference handbook.

Trap 3: Binary Interaction Parameter Arithmetic
In the cross-term $a_{ij} = \sqrt{a_i a_j}(1 - k_{ij})$, remember that $k_{ij}$ reduces the cross-attraction if $k_{ij} > 0$. A common exam trap tests your ability to expand $a_{mix} = y_1^2 a_1 + 2 y_1 y_2 a_{12} + y_2^2 a_2$. Forgetting the factor of $2$ on the cross term $y_1 y_2 a_{12}$ is a frequent algebra error.