Cubic equations of state

Quick refresher

A cubic equation is repulsion, $RT/(V_m - b)$ , minus attraction, $a/(\ldots)$ . We calculate the parameters $a$ and $b$ of each component from its critical temperature and pressure. In the Soave–Redlich–Kwong (SRK) and Peng–Robinson (PR) equations the attraction also depends on temperature, through the acentric factor $\omega$ .
In terms of $Z$ the equation is cubic. With one real root there is nothing to choose. With three real roots, the smallest belongs to the liquid and the largest to the gas, and the middle one has no meaning. The number of roots does not tell us whether a mixture is single-phase.
A mixture uses the same formulas with “average” $a$ and $b$ . The value $b$ is the average by mole fraction. The value $a$ is a sum over all pairs of molecules, with corrections $k_{ij}$ for the attraction between unlike molecules.
Fugacity is a “corrected” pressure. It measures how strongly a component tends to leave a phase. Phases are in equilibrium when the fugacity of each component is the same in all of them.
SRK and PR reproduce vapour pressures well (the pressure at which a liquid and its vapour are in equilibrium). Their liquid densities have errors from a few per cent (PR) to 7–18% (SRK), and larger for heavy substances. Volume translation corrects this error and does not change the phase equilibrium.
For a real fluid we tune the equation to laboratory PVT experiments. We tune first of all the constants of the heavy pseudo-components (groups of heavy components treated as one).

Key relationships

p = \dfrac{RT}{V_m - b} - \dfrac{a_c\,\alpha(T)}{(V_m + \delta_1 b)(V_m + \delta_2 b)}

general form of two-parameter cubic equations; pressure and temperature are absolute

a_c = \Omega_a\dfrac{R^2T_c^2}{p_c},\quad b = \Omega_b\dfrac{RT_c}{p_c}

parameters of a component from the critical constants

A = \dfrac{a\,p}{(RT)^2},\quad B = \dfrac{b\,p}{RT}

dimensionless coefficients of the cubic equation in $Z$

a = \sum_i\sum_j x_ix_j(1 - k_{ij})\sqrt{a_ia_j},\quad b = \sum_i x_ib_i

mixing rules

f_i = \varphi_i\,x_i\,p,\qquad f_i^{\,L} = f_i^{\,V}

fugacity and the condition of phase equilibrium

Common pitfalls

Taking the middle root, or the gas root for the liquid.
Using gauge pressure, or temperature in degrees Celsius.
Mixing the constants of different equations: $m(\omega)$ of SRK in the PR equation, or $k_{ij}$ from another equation.
Trusting the liquid density from an equation without volume translation.
Treating an untuned equation as accurate for a mixture with heavy fractions.

After this lesson you can

explain how a cubic equation describes the repulsion and attraction of molecules and why one formula gives both gas and liquid
write the van der Waals, Soave–Redlich–Kwong and Peng–Robinson equations and calculate their parameters from the critical constants and the acentric factor (a number that shows how far a molecule is from a simple spherical one)
solve the cubic equation for $Z$ and choose the right root (solution) for the liquid and for the gas
calculate the parameters of a mixture with the mixing rules and binary interaction coefficients (corrections for pairs of unlike molecules)
explain what fugacity is, and find the vapour pressure of a pure substance by making the fugacities of its two phases equal
judge where cubic equations are accurate, where they make errors, and how volume translation and tuning to experiments correct them

From the generalised chart to an equation of state

In the lesson Ideal and real gas: the gas deviation factor (Z) we described a real gas by one correction, the gas deviation factor $Z$ , which we found from the generalised chart. For natural gas this is convenient and accurate, but the method has limits. The chart describes only a gas. It says nothing about a liquid. It was built from gases that are mostly methane, so for mixtures with a lot of heavy components, carbon dioxide or hydrogen sulfide it needs corrections. Above all, it cannot tell us whether a mixture will split into gas and liquid, or what the compositions of the phases will be.

An equation of state links the pressure $p$ , temperature $T$ and molar volume $V_m$ in one formula. Suppose the formula correctly reproduces the behaviour of both a dense gas and a liquid. Then the rest follows from it: $Z$ and the density of each phase, and the vapour pressure of a pure substance. For a mixture it gives the bubble point and dew point pressures (the pressures at which the first gas bubbles, or the first drops of liquid, appear) and the compositions of the equilibrium phases (Vapour–liquid equilibrium). That is why equations of state are the core of compositional reservoir simulation models (Types of reservoir simulation model) and of programs that calculate PVT properties. PVT stands for pressure, volume, temperature: how the volume and phase state of a fluid depend on pressure and temperature.

Many equations of state have been proposed. Multi-parameter ones, for example the Benedict–Webb–Rubin equation (1940) with eight constants for each substance, describe pure light gases accurately. But for mixtures each constant needs its own mixing rule, and the calculation is heavy. In petroleum practice the cubic equations have become the standard. They describe a substance with two parameters, calculated from the critical constants, and the equation can be solved in closed form (by a formula, without repeated trials). This is a decisive advantage when we must find the equilibrium of a mixture of 10–15 components in every cell of a model at every time step, which means millions of times in one run.

Generalised chart and cubic equation of state

: What it describes
Generalised $Z$ chart: gas only
Cubic equation: gas, liquid and the equilibrium between them

: Input data
Generalised $Z$ chart: pseudo-critical properties of the mixture
Cubic equation: $T_c$ , $p_c$ , $\omega$ of each component and the coefficients $k_{ij}$

: Phase compositions and saturation pressures
Generalised $Z$ chart: no
Cubic equation: yes

: Accuracy of $Z$ for a dry gas
Generalised $Z$ chart: 1–2%
Cubic equation: usually 1–5% without tuning

: Where it is used
Generalised $Z$ chart: gas in place, $B_g$ , gas pipelines
Cubic equation: compositional models, gas condensate, gas injection, PVT calculations

	Generalised $Z$ chart	Cubic equation
What it describes	gas only	gas, liquid and the equilibrium between them
Input data	pseudo-critical properties of the mixture	$T_c$ , $p_c$ , $\omega$ of each component and the coefficients $k_{ij}$
Phase compositions and saturation pressures	no	yes
Accuracy of $Z$ for a dry gas	1–2%	usually 1–5% without tuning
Where it is used	gas in place, $B_g$ , gas pipelines	compositional models, gas condensate, gas injection, PVT calculations

The accuracy estimates are approximate and refer to a calculation without tuning to laboratory experiments.

The van der Waals equation: repulsion and attraction

An ideal gas is made of point molecules that do not interact. Van der Waals introduced two corrections, and both are still the basis of modern equations. The first is repulsion. The molecules take up room, so the space left for them to move in is not the whole volume but only $V_m - b$ . The quantity $b$ is called the co-volume: it is the volume to which a kilomole of the substance can be compressed at infinitely high pressure. As $V_m$ comes close to $b$ , the pressure grows without limit. This is how the equation reproduces the low compressibility of a liquid.

The second correction is attraction. A molecule at the wall of a vessel is pulled inwards by its neighbours, and its impact on the wall is weaker. The number of attracting pairs is proportional to the square of the density, so the pressure is reduced by $a/V_m^2$ . The larger the molecules and the more strongly they “stick together”, the larger the attraction parameter $a$ .

Van der Waals equation

p = \frac{RT}{V_m - b} - \frac{a}{V_m^2}

where:

$p$ absolute pressure, MPa
$T$ absolute temperature, K
$V_m$ molar volume, m³/kmol
$R$ universal gas constant, 0.008314 MPa·m³/(kmol·K), MPa·m³/(kmol·K)
$a$ attraction parameter, MPa·m⁶/kmol²
$b$ co-volume — the volume of the molecules themselves in a kilomole of the substance, m³/kmol

Units: metric: MPa, K, m³/kmol (any consistent system)

Holds when

$a$ and $b$ are constants for a substance and do not depend on temperature. This is the main simplification of the equation.
Pressure and temperature are absolute.
If $a = b = 0$ , we get the ideal gas equation $pV_m = RT$ .

The first term is repulsion. It makes the pressure higher than for an ideal gas. The second is attraction. It lowers the pressure. In a rarefied gas both are small, and the gas is nearly ideal. In a liquid both are large and almost cancel each other, so a small change in the volume of the liquid changes the pressure sharply.

In other unit systems (2)

Oilfield units

p = \frac{RT}{V_m - b} - \frac{a}{V_m^2},\qquad R = 10.73\ \frac{\text{psia}\cdot\text{ft}^3}{\text{lb-mol}\cdot{}^{\circ}\text{R}}

$p$ — psia, $T$ — °R, $V_m$ — ft³/lb-mol, $a$ — psia·ft⁶/lb-mol², $b$ — ft³/lb-mol.

The same gas constant converted to oilfield units (ideal gas equation).

Cubic equation in

Z

Z^3 - (1 + B)\,Z^2 + A\,Z - AB = 0,\qquad A = \frac{a\,p}{(RT)^2},\quad B = \frac{b\,p}{RT}

$A$ and $B$ are dimensionless and the same in any system of units.

It comes from substituting $V_m = ZRT/p$ .

Van der Waals tied the parameters $a$ and $b$ to the critical point (critical point). On the isotherm $T = T_c$ , the pressure at the critical point neither falls nor rises as the volume changes, and the curve changes its bend. The first and second derivatives of the pressure with respect to the volume are zero. These two conditions, together with the equation itself, give three relations. From them we find $a$ , $b$ and the critical molar volume.

Parameters of the van der Waals equation from the critical constants

a = \frac{27}{64}\,\frac{R^2T_c^2}{p_c},\qquad b = \frac{1}{8}\,\frac{RT_c}{p_c},\qquad Z_c = \frac{p_cV_c}{RT_c} = \frac{3}{8}

where:

$T_c,\ p_c$ critical temperature and pressure of the substance, K, MPa
$V_c$ critical molar volume; in this equation $V_c = 3b$ , m³/kmol
$Z_c$ critical Z-factor (the value of $Z$ at the critical point), —

Units: any consistent units; with $R = 0.008314$ MPa·m³/(kmol·K), $a$ comes out in MPa·m⁶/kmol² and $b$ in m³/kmol

Holds when

At the critical point the isotherm has a horizontal inflection (it flattens and changes its bend): $\partial p/\partial V_m = \partial^2 p/\partial V_m^2 = 0$ at $T = T_c$ .
The parameters are chosen so that the equation passes exactly through the critical temperature and pressure. The critical volume then comes out as the equation gives it.

The van der Waals equation gives all substances the same critical Z-factor, 0.375. For real hydrocarbons it is from 0.29 for methane to 0.25 for decane (0.23 for water). So near the critical point the equation overestimates the molar volume and underestimates the density.

Below the critical temperature the isotherm of a cubic equation behaves strangely (see the figure). As the volume decreases, the pressure first rises, then falls along a loop, and then rises again. The section where the pressure rises together with the volume is physically impossible, because such a state would split by itself into two phases. In reality, at the vapour pressure the isotherm is horizontal, because gas and liquid coexist. The equation itself fixes the position of the horizontal line. It cuts off two equal areas from the loop, which is the same as saying that the fugacities of the gas and the liquid are equal (see the section on fugacity). The ends of the horizontal line are the molar volumes of the saturated liquid and the saturated gas.

CalculatedPropane isotherms from the Peng–Robinson equation (the shape is the same as for van der Waals). One is below the critical temperature (340 K), one is at the critical temperature (369.9 K) and one is above it (400 K). The dashed line is the boundary of the two-phase region from the same equation, and

C

is the critical point. The horizontal line at 340 K is drawn at the calculated vapour pressure of 2.44 MPa. Points 1, 2 and 3 are the three roots of the equation at this pressure: liquid, a root without physical meaning, and gas. The volume axis is logarithmic.

Redlich–Kwong, Soave, Peng–Robinson

The next equations kept the repulsion $RT/(V_m - b)$ but changed the attraction. Redlich and Kwong replaced $V_m^2$ in the denominator with $V_m(V_m + b)$ , and made the attraction decrease with temperature: $a \propto 1/\sqrt{T}$ . Soave chose the temperature function $\alpha(T)$ so that the equation reproduces the vapour pressure of each substance. For this he added a third constant of the substance, the acentric factor. Peng and Robinson changed the denominator once more, to bring the liquid density and the critical Z-factor closer to reality.

General form of two-parameter cubic equations

p = \frac{RT}{V_m - b} - \frac{a_c\,\alpha(T)}{(V_m + \delta_1 b)(V_m + \delta_2 b)},\qquad a_c = \Omega_a\frac{R^2T_c^2}{p_c},\quad b = \Omega_b\frac{RT_c}{p_c}

where:

$a_c$ attraction parameter at the critical temperature, MPa·m⁶/kmol²
$\alpha(T)$ temperature function of the attraction; $\alpha(T_c) = 1$ , —
$\delta_1,\ \delta_2$ numbers that set the form of the denominator of the attraction term (table below), —
$\Omega_a,\ \Omega_b$ constants of the equation, found from the conditions at the critical point, —
$p,\ T,\ V_m,\ R,\ b$ as in the van der Waals equation, MPa, K, m³/kmol

Units: any consistent system; $\Omega_a$ , $\Omega_b$ , $\delta_1$ , $\delta_2$ are dimensionless

Holds when

A substance is described by three constants: $T_c$ , $p_c$ and (for SRK and PR) the acentric factor $\omega$ .
The equation passes exactly through the critical temperature and pressure of each substance, but not through its critical volume.
$b$ does not depend on temperature. All the dependence on temperature is in $\alpha(T)$ .

All four equations are special cases of one form. So a program solves the same cubic equation and changes only the constants.

Constants of the cubic equations

Equation: Van der Waals
$\Omega_a$: 27/64 = 0.4219
$\Omega_b$: 1/8 = 0.125
$\delta_1;\ \delta_2$: 0; 0
$Z_c$: 0.375

Equation: Redlich–Kwong, SRK
$\Omega_a$: 0.42748
$\Omega_b$: 0.08664
$\delta_1;\ \delta_2$: 1; 0
$Z_c$: 0.333

Equation: Peng–Robinson
$\Omega_a$: 0.45724
$\Omega_b$: 0.07780
$\delta_1;\ \delta_2$: $1 + \sqrt2;\ 1 - \sqrt2$
$Z_c$: 0.307

Equation	$\Omega_a$	$\Omega_b$	$\delta_1;\ \delta_2$	$Z_c$
Van der Waals	27/64 = 0.4219	1/8 = 0.125	0; 0	0.375
Redlich–Kwong, SRK	0.42748	0.08664	1; 0	0.333
Peng–Robinson	0.45724	0.07780	$1 + \sqrt2;\ 1 - \sqrt2$	0.307

$Z_c$ is the critical Z-factor that the equation gives to any substance. For real hydrocarbons it is 0.25–0.29. The function $\alpha(T)$ is $\alpha = 1$ for van der Waals, $\alpha = 1/\sqrt{T_r}$ for Redlich–Kwong, and the formula below for SRK and PR.

Peng–Robinson equation

p = \frac{RT}{V_m - b} - \frac{a_c\,\alpha(T)}{V_m^2 + 2bV_m - b^2},\qquad a_c = 0.45724\,\frac{R^2T_c^2}{p_c},\quad b = 0.07780\,\frac{RT_c}{p_c}

where:

$a_c,\ b$ attraction parameter at $T_c$ and co-volume of the component, MPa·m⁶/kmol², m³/kmol
$\alpha(T)$ temperature function (next formula), —
$p,\ T,\ V_m,\ R$ absolute pressure and temperature, molar volume, gas constant, MPa, K, m³/kmol

Units: metric: MPa, K, m³/kmol (the constants 0.45724 and 0.07780 are dimensionless and suit any system)

Holds when

The denominator is $V_m^2 + 2bV_m - b^2 = (V_m + 2.414\,b)(V_m - 0.414\,b)$ , that is, $\delta_1 = 1 + \sqrt2$ and $\delta_2 = 1 - \sqrt2$ .
For a mixture, take $a$ and $b$ from the mixing rules (section on mixtures).

Compared with SRK, the co-volume is smaller ( $\Omega_b = 0.0778$ instead of $0.0866$ ). The new denominator of the attraction term also brings the critical Z-factor closer to the real one: 0.307 instead of 0.333. So the liquid comes out denser than with SRK.

In other unit systems (1)

Soave–Redlich–Kwong equation (SRK)

p = \frac{RT}{V_m - b} - \frac{a_c\,\alpha(T)}{V_m(V_m + b)},\qquad a_c = 0.42748\,\frac{R^2T_c^2}{p_c},\quad b = 0.08664\,\frac{RT_c}{p_c}

the same units; in the original Redlich–Kwong equation $\alpha = 1/\sqrt{T/T_c}$ .

Temperature function of the attraction in the SRK and PR equations

\alpha(T) = \left[1 + m\left(1 - \sqrt{T/T_c}\right)\right]^2

where:

$m$ coefficient that depends on the acentric factor (forms below), —
$T/T_c$ reduced temperature $T_r$ , —

Units: dimensionless; temperatures absolute

Holds when

The dependence $m(\omega)$ was fitted to the vapour pressures of hydrocarbons, so each equation has its own.
For heavy substances with $\omega > 0.49$ , the PR equation uses a refined dependence (advanced section).

At $T < T_c$ the function is greater than 1, so the attraction grows on cooling. At $T > T_c$ it is less than 1. The larger $\omega$ , the steeper this dependence.

In other unit systems (2)

Soave–Redlich–Kwong (SRK)

m = 0.480 + 1.574\,\omega - 0.176\,\omega^2

dimensionless

Peng–Robinson (PR)

m = 0.37464 + 1.54226\,\omega - 0.26992\,\omega^2

dimensionless

Acentric factor

\omega = -\lg\left.\frac{p_s}{p_c}\right|_{T = 0.7\,T_c} - 1

where:

$\omega$ acentric factor of the substance, —
$p_s$ vapour pressure at the temperature $0.7\,T_c$ , MPa
$\lg$ base-10 logarithm, —

Units: dimensionless; pressures absolute, in the same units

Holds when

It is defined from the vapour pressure curve of the liquid. For carbon dioxide the temperature $0.7\,T_c$ is below the triple point (where solid, liquid and gas coexist), so the curve is extended beyond it.
For C7+ fractions and pseudo-components, $\omega$ is estimated from correlations.

For substances with small, almost spherical molecules (argon, krypton), the vapour pressure at $T = 0.7\,T_c$ is almost exactly $0.1\,p_c$ , and $\omega \approx 0$ . The longer and less symmetric the molecule, the lower its vapour pressure at the same reduced temperature, and the larger $\omega$ is: 0.011 for methane, 0.152 for propane, 0.49 for n-decane. Through $\omega$ the equation “knows” not only the end of the vapour pressure curve (the critical point) but also its slope.

The calculation needs three constants for each component. For individual substances they are reference values. For the C7+ fraction and pseudo-components we estimate them with correlations from the molar mass and density (C7+ fraction). This estimate becomes the main source of uncertainty in the calculation.

Component constants used in the examples

Component: Nitrogen $\mathrm{N_2}$
$T_c$ , K: 126.2
$p_c$ , MPa: 3.40
$\omega$: 0.037
$M$ , g/mol: 28.01

Component: Carbon dioxide $\mathrm{CO_2}$
$T_c$ , K: 304.2
$p_c$ , MPa: 7.38
$\omega$: 0.225
$M$ , g/mol: 44.01

Component: Hydrogen sulfide $\mathrm{H_2S}$
$T_c$ , K: 373.2
$p_c$ , MPa: 9.00
$\omega$: 0.094
$M$ , g/mol: 34.08

Component: Methane $\mathrm{C_1}$
$T_c$ , K: 190.6
$p_c$ , MPa: 4.60
$\omega$: 0.011
$M$ , g/mol: 16.04

Component: Ethane $\mathrm{C_2}$
$T_c$ , K: 305.4
$p_c$ , MPa: 4.87
$\omega$: 0.099
$M$ , g/mol: 30.07

Component: Propane $\mathrm{C_3}$
$T_c$ , K: 369.9
$p_c$ , MPa: 4.25
$\omega$: 0.152
$M$ , g/mol: 44.10

Component: Isobutane
$T_c$ , K: 407.9
$p_c$ , MPa: 3.64
$\omega$: 0.184
$M$ , g/mol: 58.12

Component: n-Butane
$T_c$ , K: 425.2
$p_c$ , MPa: 3.80
$\omega$: 0.200
$M$ , g/mol: 58.12

Component: n-Pentane
$T_c$ , K: 469.8
$p_c$ , MPa: 3.37
$\omega$: 0.252
$M$ , g/mol: 72.15

Component: n-Hexane
$T_c$ , K: 507.7
$p_c$ , MPa: 3.03
$\omega$: 0.301
$M$ , g/mol: 86.18

Component: n-Heptane
$T_c$ , K: 540.3
$p_c$ , MPa: 2.74
$\omega$: 0.350
$M$ , g/mol: 100.20

Component: n-Decane
$T_c$ , K: 617.8
$p_c$ , MPa: 2.11
$\omega$: 0.490
$M$ , g/mol: 142.28

Component	$T_c$ , K	$p_c$ , MPa	$\omega$	$M$ , g/mol
Nitrogen $\mathrm{N_2}$	126.2	3.40	0.037	28.01
Carbon dioxide $\mathrm{CO_2}$	304.2	7.38	0.225	44.01
Hydrogen sulfide $\mathrm{H_2S}$	373.2	9.00	0.094	34.08
Methane $\mathrm{C_1}$	190.6	4.60	0.011	16.04
Ethane $\mathrm{C_2}$	305.4	4.87	0.099	30.07
Propane $\mathrm{C_3}$	369.9	4.25	0.152	44.10
Isobutane	407.9	3.64	0.184	58.12
n-Butane	425.2	3.80	0.200	58.12
n-Pentane	469.8	3.37	0.252	72.15
n-Hexane	507.7	3.03	0.301	86.18
n-Heptane	540.3	2.74	0.350	100.20
n-Decane	617.8	2.11	0.490	142.28

Reference values. The critical temperatures and pressures, in kelvins, are as in the table of the lesson on the composition of fluids.

The cubic equation in Z: one root or three

Usually we know the pressure and temperature and need the molar volume or density. If we substitute $V_m = ZRT/p$ , we get a cubic equation in $Z$ with dimensionless coefficients $A$ and $B$ . They are the same in any system of units.

Peng–Robinson equation in terms of the Z-factor

Z^3 - (1 - B)\,Z^2 + (A - 3B^2 - 2B)\,Z - (AB - B^2 - B^3) = 0

where:

$Z$ Z-factor of the phase, $Z = pV_m/(RT)$ , —
$A$ $a\,p/(RT)^2$ — dimensionless attraction, —
$B$ $b\,p/(RT)$ — dimensionless co-volume, —
$a,\ b$ parameters of the substance at temperature $T$ (for a mixture, from the mixing rules), MPa·m⁶/kmol², m³/kmol

Units: dimensionless form; $p$ and $T$ absolute

Holds when

$A$ and $B$ are calculated at the same absolute pressure and temperature for which $Z$ is sought.
For a mixture, $a$ and $b$ are calculated from the composition of the phase whose properties we want.

Each root $Z$ immediately gives the molar volume $V_m = ZRT/p$ and the density $\rho = pM/(ZRT)$ . Only roots with $Z > B$ , that is $V_m > b$ , have physical meaning.

In other unit systems (2)

SRK equation

Z^3 - Z^2 + (A - B - B^2)\,Z - AB = 0

$A$ and $B$ are calculated with the SRK constants.

Van der Waals equation

Z^3 - (1 + B)\,Z^2 + A\,Z - AB = 0

$A$ and $B$ are calculated with the van der Waals constants.

A cubic equation has one or three real roots (real-number solutions). For a pure substance, one root means that at the given pressure and temperature it is single-phase, so there is nothing to choose. Three roots appear where the isotherm forms a loop (figure above). The smallest root belongs to the liquid and the largest to the gas. The middle one lies on the section where the pressure rises together with the volume, and it has no meaning. For a mixture, the number of roots says nothing about stability. For the model oil of the topic Vapour–liquid equilibrium at 8 MPa, the equation for the composition of the whole mixture has one root, although the mixture splits into gas and liquid. A stability test shows whether a mixture will split into phases (stability check).

Which root to take

1
The phase is a liquid
If you know the phase is a liquid (for example, the equilibrium liquid in a phase-equilibrium calculation), take the smallest root.
2
The phase is a gas
If the phase is a gas, take the largest root.
3
The phase is unknown
For a single-phase mixture at given $p$ and $T$ , calculate both possibilities. Take the one with the lower Gibbs energy (a measure of how stable a state is), which means the one with the smaller sum $\sum_i x_i\ln\varphi_i$ (section on fugacity). For a pure substance this means that above the vapour pressure the liquid is stable, and below it the gas is.

We solve the cubic equation with Cardano's formulas (advanced section) or with Newton's method. Programs usually start from the solution given by the formulas and refine it with one or two Newton steps.

Worked examplePropane at 20 °C: three roots of the Peng–Robinson equation

Given (assumed data)

Substance: propane: $T_c = 369.9$ K, $p_c = 4.25$ MPa, $\omega = 0.152$ , $M = 44.10$ g/mol
Temperature: 20 °C ( $T = 293.15$ K)
Pressure: 0.84 MPa (absolute), about the vapour pressure

Solution

Reduced temperature: $T_r = 293.15/369.9 = 0.7925$ . Coefficient: $m = 0.37464 + 1.54226\cdot 0.152 - 0.26992\cdot 0.152^2 = 0.6028$ . Then $\alpha = [1 + 0.6028(1 - \sqrt{0.7925})]^2 = 1.13672$ .
$RT_c = 0.00831446\cdot 369.9 = 3.07552$ MPa·m³/kmol (here and below we use $R$ with enough digits for the calculation). $a_c = 0.45724\cdot 3.07552^2/4.25 = 1.01764$ MPa·m⁶/kmol², and $a = a_c\alpha = 1.1568$ . $b = 0.07780\cdot 3.07552/4.25 = 0.05630$ m³/kmol.
$RT = 0.00831446\cdot 293.15 = 2.4374$ . Then $A = a\,p/(RT)^2 = 0.16356$ and $B = b\,p/(RT) = 0.01940$ .
The equation is $Z^3 - 0.98060\,Z^2 + 0.12363\,Z - 0.00279 = 0$ . It has three real roots: $Z = 0.0291$ , $0.1147$ and $0.837$ . All are greater than $B$ .
Molar volumes $V_m = ZRT/p$ : liquid 0.0844 m³/kmol, gas 2.43 m³/kmol. Densities $\rho = M/V_m$ : 523 and 18 kg/m³. The middle root would give 133 kg/m³, which is neither a liquid nor a gas.

Answer: The liquid has

Z = 0.0291

(about 523 kg/m³), and the gas has

Z = 0.837

(about 18 kg/m³). We discard the middle root, 0.1147.

What it means. At a pressure close to the vapour pressure, the equation allows both phases. The fugacities decide which one is stable (next sections). The molar volume of the liquid is only 1.50 times the co-volume $b$ . The liquid is squeezed almost to the limit, so its density is very sensitive to the constants of the equation. The measured density of liquid propane at 20 °C is about 500 kg/m³, and the calculation overestimates it by 4.6%. The section on volume translation shows how to correct this. The properties of propane are reference values.

Mixtures: mixing rules and binary interaction coefficients

We describe a mixture with the same equation, as if it were a single substance with “average” parameters $a$ and $b$ . The co-volume of the mixture is simply the average by mole fraction, because the volumes of the molecules themselves add up. The attraction is more complicated. It acts between pairs of molecules, so the parameter $a$ of the mixture is a sum over all pairs, including pairs of unlike molecules.

Van der Waals mixing rules

a = \sum_i\sum_j x_ix_j\,(1 - k_{ij})\sqrt{a_ia_j},\qquad b = \sum_i x_ib_i

where:

$x_i$ mole fraction of the component in the phase: $x_i$ in the liquid, $y_i$ in the gas, $z_i$ in the mixture as a whole, fraction
$a_i,\ b_i$ parameters of the component at temperature $T$ : $a_i = a_{c,i}\,\alpha_i(T)$ , MPa·m⁶/kmol², m³/kmol
$k_{ij}$ binary interaction coefficient; $k_{ii} = 0$ , $k_{ij} = k_{ji}$ , —

Units: as for the parameters of the components

Holds when

The mixture is treated as one “average” substance (the van der Waals mixing rules).
The $k_{ij}$ do not depend on pressure, and in the simplest version not on temperature either.
Use the composition of the phase whose properties you are calculating. The gas and the liquid of the same mixture have different $a$ and $b$ .

The attraction between unlike molecules is estimated by the geometric mean $\sqrt{a_ia_j}$ . This is a guess, and it is good for similar molecules. The factor $(1 - k_{ij})$ corrects it. If $k_{ij} > 0$ , unlike molecules attract each other more weakly than the mean assumes, and the components separate more readily into different phases.

Binary interaction coefficients are fitted from experiments with binary mixtures (two components) or from PVT experiments with a particular fluid. They are empirical corrections to the mixing rule (they come from experiments, not from theory), and each belongs to one particular equation. For hydrocarbons of similar size they are usually taken as zero. Noticeable values are needed for pairs of methane with heavy hydrocarbons, and especially for pairs with non-hydrocarbon components.

Approximate binary interaction coefficients for the Peng–Robinson equation

Pair of components: Hydrocarbons of similar size ( $\mathrm{C_1}$ — $\mathrm{C_2}$ , $\mathrm{C_2}$ — $\mathrm{C_3}$ and so on)
Order of $k_{ij}$: 0–0.01

Pair of components: Methane — heavy hydrocarbons ( $\mathrm{C_7}$ and heavier)
Order of $k_{ij}$: 0.02–0.06

Pair of components: Nitrogen — hydrocarbons
Order of $k_{ij}$: 0.02–0.12

Pair of components: Carbon dioxide — hydrocarbons
Order of $k_{ij}$: 0.09–0.15

Pair of components: Hydrogen sulfide — hydrocarbons
Order of $k_{ij}$: 0.05–0.09

Pair of components	Order of $k_{ij}$
Hydrocarbons of similar size ( $\mathrm{C_1}$ — $\mathrm{C_2}$ , $\mathrm{C_2}$ — $\mathrm{C_3}$ and so on)	0–0.01
Methane — heavy hydrocarbons ( $\mathrm{C_7}$ and heavier)	0.02–0.06
Nitrogen — hydrocarbons	0.02–0.12
Carbon dioxide — hydrocarbons	0.09–0.15
Hydrogen sulfide — hydrocarbons	0.05–0.09

These are orders of magnitude for orientation, not data for calculation. The values depend on the equation, on how the heavy fractions are characterised and on tuning to experiments on a particular fluid. Pairs with water need special rules.

Worked exampleGas deviation factor of a natural gas from the Peng–Robinson equation

Given (assumed data)

Gas composition, mole fractions: methane 0.88; ethane 0.07; propane 0.03; n-butane 0.02
Pressure and temperature: 25 MPa (absolute), 85 °C ( $T = 358.15$ K)
Binary interaction coefficients: all $k_{ij} = 0$ (hydrocarbons of similar size)
Component constants: from the table above

Solution

$RT = 0.00831446\cdot 358.15 = 2.9778$ MPa·m³/kmol. For each component we find $a_i = a_{c,i}\alpha_i(T)$ and $b_i$ as in the propane example. We get $a_i$ = 0.1824; 0.5538; 1.0374; 1.6748 MPa·m⁶/kmol² and $b_i$ = 0.02680; 0.04057; 0.05630; 0.07238 m³/kmol (methane, ethane, propane, butane). Methane here is above its critical temperature, so its $\alpha = 0.731 < 1$ .
With all $k_{ij} = 0$ , the double sum collapses into a square: $a = \left(\sum_i y_i\sqrt{a_i}\right)^2$ . The square roots are $\sqrt{a_i}$ = 0.4271; 0.7442; 1.0185; 1.2941. So $\sum_i y_i\sqrt{a_i} = 0.4844$ and $a = 0.2346$ MPa·m⁶/kmol². The co-volume of the mixture is $b = \sum_i y_ib_i = 0.02956$ m³/kmol.
$A = a\,p/(RT)^2 = 0.6614$ and $B = b\,p/(RT) = 0.2482$ .
The equation $Z^3 - 0.7518\,Z^2 - 0.0198\,Z - 0.0873 = 0$ has one real root, $Z = 0.885$ , so there is nothing to choose. We know in advance that a dry gas of this composition is single-phase at 85 °C. This is far above its cricondentherm, the highest temperature at which two phases are possible for it (Phase behaviour: pure substances and binary mixtures).
The molar mass of the gas is $M = 18.71$ g/mol, and its density is $\rho = pM/(ZRT) = 177$ kg/m³.
For comparison, here is the generalised chart (Ideal and real gas: the gas deviation factor (Z)). By Kay's rule, $T_{pc} = 208.7$ K and $p_{pc} = 4.592$ MPa. Then $T_{pr} = 1.716$ , $p_{pr} = 5.444$ and $Z = 0.903$ .

Answer: The Peng–Robinson equation gives

Z = 0.885

and a gas density of 177 kg/m³. The generalised chart gives

Z = 0.903

, which is 2.0% more.

What it means. A difference of a couple of per cent is typical of an untuned equation. The SRK equation for the same gas gives $Z = 0.943$ , which is 4.4% above the chart, so it errs in the other direction. For a dry gas of this composition the generalised chart is built directly from measurements, and it is usually more accurate. The value of the equation lies elsewhere: with the same parameters it describes the condensate that drops out of a richer gas, oil, and the equilibrium between them. The data in this example are assumed for illustration; the component constants are reference values.

CalculatedGas deviation factor of the gas in the example at 85 °C, from three equations of state and from the generalised chart (Kay's rule and the approximation from the lesson on real gas behaviour). The dot is the Peng–Robinson calculation at 25 MPa.

Fugacity: the “corrected” pressure

How can we tell, from an equation of state, at what pressure propane boils at 20 °C? With three roots the equation gives both the liquid and the gas. But it does not say which phase is stable, or at what pressure the two coexist. For that we need a condition of equilibrium. For an ideal gas it is simple: a component moves to where its partial pressure (its own share of the total pressure) is lower, until the partial pressures become equal. For real phases, the role of the partial pressure is played by the fugacity (from the Latin fuga, “flight”). It is a “corrected” pressure that measures how strongly the molecules of a component tend to leave the phase.

Fugacity and the condition of phase equilibrium

f_i = \varphi_i\,x_i\,p,\qquad f_i^{\,L} = f_i^{\,V}\ \ \text{for each component}

where:

$f_i$ fugacity of component $i$ in the phase, MPa
$\varphi_i$ fugacity coefficient of the component in the phase (not to be confused with porosity); it tends to 1 as $p \to 0$ , —
$x_i$ mole fraction of the component in the phase (in the gas it is denoted $y_i$ ), fraction
$p$ absolute pressure, MPa
$L,\ V$ labels for the liquid and the gas phase, —

Units: pressure in any units; $\varphi_i$ is dimensionless

Holds when

The phases are at the same pressure and temperature.
The fugacities of each component are equal in all phases. This is the condition of equilibrium, and it is equivalent to equal chemical potentials (the quantities that drive a component to move from one phase to another).
For an ideal gas $\varphi_i = 1$ , and the fugacity equals the partial pressure $y_ip$ .

A component moves from the phase where its fugacity is higher to the phase where it is lower, until the two become equal. If $\varphi_i < 1$ , the attraction holds the molecules more strongly than in an ideal gas. If $\varphi_i > 1$ , the component is “pushed out” of the phase. At equilibrium the ratio of the mole fractions is the equilibrium ratio, $K_i = y_i/x_i = \varphi_i^{\,L}/\varphi_i^{\,V}$ (Vapour–liquid equilibrium).

We calculate the fugacity coefficient from the same equation of state. It shows how far the isotherm of the real phase departs from the isotherm of an ideal gas when we compress from zero pressure (advanced section). For cubic equations the integral can be written in closed form. For the Peng–Robinson equation:

Fugacity coefficient of a mixture component from the Peng–Robinson equation

\ln\varphi_i = \frac{b_i}{b}(Z - 1) - \ln(Z - B) - \frac{A}{2\sqrt2\,B}\left(\frac{2\sum_j x_ja_{ij}}{a} - \frac{b_i}{b}\right)\ln\frac{Z + (1 + \sqrt2)B}{Z + (1 - \sqrt2)B}

where:

$a_{ij}$ $(1 - k_{ij})\sqrt{a_ia_j}$ — the attraction of a pair of components, MPa·m⁶/kmol²
$a,\ b,\ A,\ B$ parameters of the phase from the mixing rules and their dimensionless forms, MPa·m⁶/kmol², m³/kmol; $A$ and $B$ are dimensionless
$Z$ the root chosen for this phase: the smallest for the liquid, the largest for the gas, —
$b_i$ co-volume of the component, m³/kmol

Units: dimensionless; natural logarithms

Holds when

Van der Waals mixing rules with the coefficients $k_{ij}$ .
$Z$ , $A$ and $B$ are calculated for the composition and the root of this particular phase.

The first two terms are the contribution of repulsion, and the last term is that of attraction. The sum $\sum_j x_ja_{ij}$ shows how strongly the component is attracted to its surroundings in the given phase. That is why the same methane has different fugacity coefficients in oil and in gas.

In other unit systems (2)

Pure substance

\ln\varphi = Z - 1 - \ln(Z - B) - \frac{A}{2\sqrt2\,B}\ln\frac{Z + (1 + \sqrt2)B}{Z + (1 - \sqrt2)B}

the same quantities for one component ( $x = 1$ , $b_i = b$ , $\sum_j x_ja_{ij} = a$ ).

SRK equation

\ln\varphi_i = \frac{b_i}{b}(Z - 1) - \ln(Z - B) - \frac{A}{B}\left(\frac{2\sum_j x_ja_{ij}}{a} - \frac{b_i}{b}\right)\ln\left(1 + \frac{B}{Z}\right)

the same notation, $A$ and $B$ with the SRK constants.

Worked exampleVapour pressure of propane from the equality of fugacities

Given (assumed data)

Substance and temperature: propane, 20 °C, as in the example above
Initial approximation: $p^{(0)} = 0.84$ MPa, with the three roots of $Z$ at this pressure from the previous example
Measured vapour pressure: 0.8364 MPa

Solution

For the liquid root, $Z_L = 0.0291$ , the formula for a pure substance gives $\ln\varphi_L = -0.157155$ . For the gas root, $Z_V = 0.837$ , it gives $\ln\varphi_V = -0.152665$ . The fugacities are $f_L = 0.7178$ MPa and $f_V = 0.7211$ MPa.
$f_L < f_V$ : the liquid holds the molecules better than the gas, so the pressure is above the vapour pressure. The correction is $p^{(1)} = p^{(0)}\,\varphi_L/\varphi_V = 0.84\cdot \e^{-0.00449} = 0.83624$ MPa.
At $p^{(1)}$ we find the roots and fugacity coefficients again. This gives $p^{(2)} = 0.83552$ MPa and $p^{(3)} = 0.83538$ MPa. Each time the correction shrinks by a factor of about 5. After 6 steps the pressure stops changing in the sixth decimal place: $p_s = 0.8353$ MPa.
At $p_s$ : $Z_L = 0.0289$ , $Z_V = 0.8379$ and $\varphi_L = \varphi_V = 0.859$ . The density of the saturated liquid is 523 kg/m³, and that of the saturated gas is 18.0 kg/m³.
The same calculations with other equations give 0.844 MPa for SRK and 1.560 MPa for van der Waals.

Answer: The Peng–Robinson equation gives a vapour pressure of propane at 20 °C of 0.8353 MPa, which is 0.1% below the measured value. SRK gives 0.9% more than the measured value, and the van der Waals equation gives 86.5% more.

What it means. The SRK and PR equations reproduce the vapour pressure well because their temperature function is fitted, through the acentric factor, to exactly this property. The van der Waals equation has no such tuning. The density of the saturated liquid is still overestimated by a few per cent, because equal fugacities do not correct it. The same procedure, refining the pressure by the ratio of the fugacities, is the basis for calculating the bubble point and dew point pressures of mixtures (Vapour–liquid equilibrium). The properties of propane are reference values.

Volume translation: how the liquid density is corrected

A two-parameter equation gives all substances the same critical Z-factor, so it reproduces the density of a liquid with a systematic error (an error that always goes the same way). We can see its size by comparing the calculation with the reference densities of liquid hydrocarbons.

Density of liquid n-alkanes at 15.6 °C: reference and from the equations, kg/m³

Substance: Propane
Reference: 506.5
SRK (error): 468.6 (−7.5%)
PR (error): 531.6 (+5.0%)

Substance: n-Butane
Reference: 583.4
SRK (error): 543.3 (−6.9%)
PR (error): 614.9 (+5.4%)

Substance: n-Pentane
Reference: 630.4
SRK (error): 575.9 (−8.6%)
PR (error): 650.2 (+3.1%)

Substance: n-Hexane
Reference: 663.3
SRK (error): 594.5 (−10.4%)
PR (error): 670.0 (+1.0%)

Substance: n-Heptane
Reference: 687.3
SRK (error): 603.0 (−12.3%)
PR (error): 678.6 (−1.3%)

Substance: n-Decane
Reference: 733.3
SRK (error): 603.9 (−17.6%)
PR (error): 677.6 (−7.6%)

Substance	Reference	SRK (error)	PR (error)
Propane	506.5	468.6 (−7.5%)	531.6 (+5.0%)
n-Butane	583.4	543.3 (−6.9%)	614.9 (+5.4%)
n-Pentane	630.4	575.9 (−8.6%)	650.2 (+3.1%)
n-Hexane	663.3	594.5 (−10.4%)	670.0 (+1.0%)
n-Heptane	687.3	603.0 (−12.3%)	678.6 (−1.3%)
n-Decane	733.3	603.9 (−17.6%)	677.6 (−7.6%)

The reference densities come from the specific gravities at 60 °F (15.6 °C), which are 0.507, 0.584, 0.631, 0.664, 0.688 and 0.734, with a water density of 999.0 kg/m³. The calculation uses the equations of the lesson at the same temperature. Propane and butane are on the saturation line (at their own calculated vapour pressure), and the rest are at atmospheric pressure.

SRK underestimates the density of all liquids, and the more so the heavier the hydrocarbon. The Peng–Robinson equation overestimates the density for light liquids and underestimates it for heavy ones. In an oil the errors of the different components partly cancel each other, but you cannot rely on this.

In 1982 Péneloux and co-authors proposed a simple correction, volume translation. We subtract a constant correction, with its own value for each component, from the calculated molar volume.

Volume translation (Péneloux correction)

V_m = V_m^{\text{EOS}} - c,\qquad c = \sum_i x_ic_i

where:

$V_m^{\text{EOS}}$ molar volume of the phase from the equation of state, m³/kmol
$c_i$ translation of the component; it can also be negative, m³/kmol
$c$ translation of the phase — the average over its composition, m³/kmol
$x_i$ mole fraction of the component in the phase, fraction

Units: m³/kmol

Holds when

The $c_i$ do not depend on pressure, and in the simplest version not on temperature either.
The $c_i$ are fitted to the measured density of a liquid (the pure component, stock-tank oil, which is oil at surface conditions, or reservoir oil) or estimated from correlations. For $c > 0$ the volume decreases and the density rises.

The translation changes only volumes and densities. The fugacity coefficient of each component is multiplied by the same factor, $\e^{-c_ip/(RT)}$ , in both phases. So the ratios $\varphi_i^{\,L}/\varphi_i^{\,V}$ do not change, and neither do the equilibrium ratios, saturation pressures and phase compositions. For a gas, whose molar volume is many times larger than $c$ , the correction is small. For a liquid it is significant.

Worked exampleVolume translation for propane

Given (assumed data)

Reference density of liquid propane: at 15.6 °C: $0.507\cdot 999.0$ kg/m³; at 20 °C: about 500 kg/m³
Equation: Peng–Robinson, propane constants from the table

Solution

At 15.6 °C ( $T = 288.71$ K) the calculated vapour pressure is 0.741 MPa. The equation gives a molar volume of the saturated liquid of $V_m^{\text{EOS}} = 0.08296$ m³/kmol, and a density of 531.6 kg/m³.
The reference density is $0.507\cdot 999.0 = 506.5$ kg/m³, so the molar volume is $44.10/506.5 = 0.08707$ m³/kmol.
The translation is $c = V_m^{\text{EOS}} - V_m = 0.08296 - 0.08707 = -0.00411$ m³/kmol. It is negative, because the equation gives too small a volume. As a fraction of the co-volume, $c/b = -0.073$ .
At 20 °C the equation gives $V_m^{\text{EOS}} = 0.08436$ m³/kmol (522.8 kg/m³). With the translation, $V_m = V_m^{\text{EOS}} - c$ , so the density is $\rho = 44.10/(0.08436 + 0.00411) = 498.5$ kg/m³.
For the saturated gas at 20 °C, $V_m^{\text{EOS}} = 2.445$ m³/kmol. The same translation changes its volume by only 0.17%. The vapour pressure, 0.835 MPa, does not change.

Answer: The translation

c \approx -0.00411

m³/kmol, found at 15.6 °C, reduces the error of the density of liquid propane at 20 °C from 4.6% to −0.3%.

What it means. A constant translation works well near the temperature at which it was fitted. Closer to the critical temperature the error of the equation changes, and a constant correction fixes it less well. So programs also use translations that depend on temperature. For C7+ fractions we fit the translations so that the equation reproduces the measured densities of stock-tank and reservoir oil. The reference densities are rounded.

Where cubic equations are accurate and where they err

Accuracy of an untuned cubic equation (approximate)

Property: Vapour pressure of pure hydrocarbons
Typical error: from fractions of a per cent to a few per cent
Reason: $\alpha(T)$ is fitted to the vapour pressure through $\omega$

Property: $Z$ and gas density
Typical error: 1–5%
Reason: two parameters are not enough to describe a dense gas accurately, and at very high pressures the error grows

Property: Liquid density
Typical error: PR a few per cent, SRK 7–18%; larger for heavy substances
Reason: one $Z_c$ for all substances; volume translation corrects it

Property: Properties near the critical point of a mixture
Typical error: large
Reason: a cubic equation does not reproduce the shape of the isotherms near the critical point

Property: Mixtures with a C7+ fraction
Typical error: depends on how well the C7+ fraction is described
Reason: $T_c$ , $p_c$ and $\omega$ of the pseudo-components are estimated from correlations

Property: Water, methanol, glycols
Typical error: poor without special rules
Reason: hydrogen bonds (strong links between molecules) cannot be described by two parameters

Property	Typical error	Reason
Vapour pressure of pure hydrocarbons	from fractions of a per cent to a few per cent	$\alpha(T)$ is fitted to the vapour pressure through $\omega$
$Z$ and gas density	1–5%	two parameters are not enough to describe a dense gas accurately, and at very high pressures the error grows
Liquid density	PR a few per cent, SRK 7–18%; larger for heavy substances	one $Z_c$ for all substances; volume translation corrects it
Properties near the critical point of a mixture	large	a cubic equation does not reproduce the shape of the isotherms near the critical point
Mixtures with a C7+ fraction	depends on how well the C7+ fraction is described	$T_c$ , $p_c$ and $\omega$ of the pseudo-components are estimated from correlations
Water, methanol, glycols	poor without special rules	hydrogen bonds (strong links between molecules) cannot be described by two parameters

The equation does not give viscosity or interfacial tension (the tension at the boundary between two fluids). They are calculated with separate correlations, which use the densities and phase compositions that the equation gives.

That is why we tune an equation of state for a particular fluid. We fit the critical constants and acentric factors of the heavy pseudo-components, the $k_{ij}$ of methane with the heavy fractions, and the volume translations. We fit them to the results of laboratory PVT experiments, which are tests on a fluid sample at different pressures and temperatures. These are the saturation pressure test, constant-composition expansion, constant-volume depletion, differential liberation and separator tests (Laboratory PVT tests). A tuned equation is a model of this particular fluid. It cannot be carried over to the fluid of another reservoir.

Terms and symbols

Equation of state: A formula that links the pressure, temperature and molar volume of a substance or mixture. Cubic equations describe both gas and liquid.
Co-volume $b$: The volume of the molecules themselves in a kilomole of the substance. It is the limiting molar volume at infinite pressure.
Attraction parameter $a$: A quantity that describes the attraction between molecules. In SRK and PR it depends on temperature through the function $\alpha(T)$ .
Acentric factor $\omega$: It is $-\lg(p_s/p_c) - 1$ at $T = 0.7\,T_c$ . It measures the slope of the vapour pressure curve, and it grows as the molecule gets larger and less symmetric.
Critical Z-factor $Z_c$: It equals $p_cV_c/(RT_c)$ . For real hydrocarbons it is 0.25–0.29, while a cubic equation gives one value to all substances (0.307 for PR).
Mixing rules: Formulas that calculate the parameters of a mixture from the parameters of the components and the composition of the phase.
Binary interaction coefficient $k_{ij}$: An empirical correction to the attraction between unlike molecules. It is specific to each equation.
Fugacity $f_i$ and fugacity coefficient $\varphi_i$: The fugacity is the “corrected” partial pressure of a component, and the fugacity coefficient is its ratio to $x_ip$ . At equilibrium the fugacities of each component are equal in all the phases.
Volume translation $c$: A correction to the molar volume from the equation. It fixes the liquid density without changing the phase equilibrium.

Common mistakes and misconceptions

Wrong

“The equation gave three roots — I'll take the middle one, it is closer to the truth.”

In fact

The middle root lies on the unstable part of the loop and has no meaning. For a liquid take the smallest root, and for a gas take the largest. For a single-phase mixture of unknown nature, take the root with the lower Gibbs energy.

Wrong

“All programs use the Peng–Robinson equation, so it is accurate.”

In fact

Without tuning, the liquid density is wrong by a few per cent, and by more for heavy components. Properties near the critical point are wrong by more still. For mixtures with C7+ everything depends on the constants of the pseudo-components. For a real fluid we tune the equation to PVT experiments.

Wrong

“Volume translation corrected the oil density, so the saturation pressure will shift too.”

In fact

Volume translation changes neither the equilibrium ratios nor the saturation pressures, because the fugacity coefficients of each component change in the same way in both phases.

Wrong

“The $k_{ij}$ are physical constants of a pair of substances, they can be taken from any table.”

In fact

They are empirical corrections fitted for one particular equation, and often for one particular way of describing the heavy fractions. Coefficients from another equation make the calculation worse.

Wrong

“I'll substitute 85 °C and the pressure from the gauge.”

In fact

Use only absolute values in equations of state: 358.15 K, and the pressure with the atmospheric pressure added.

In depth: advanced detail

How a cubic equation is solved

Write the equation in $Z$ in the general form $Z^3 + c_2Z^2 + c_1Z + c_0 = 0$ (for PR $c_2 = B - 1$ , $c_1 = A - 3B^2 - 2B$ , $c_0 = -(AB - B^2 - B^3)$ ). The substitution $Z = t - c_2/3$ removes the $Z^2$ term:

Depressed cubic equation

t^3 + P\,t + Q = 0,\qquad P = c_1 - \frac{c_2^2}{3},\quad Q = \frac{2c_2^3}{27} - \frac{c_2c_1}{3} + c_0

where:

$t$ the new unknown, $t = Z + c_2/3$ , —
$P,\ Q$ coefficients of the depressed equation, —
$c_2,\ c_1,\ c_0$ coefficients of the original equation in $Z$ , —

Units: dimensionless

Holds when

The coefficients are real numbers. The sign of $D = Q^2/4 + P^3/27$ decides how many real roots there are.

This is a standard trick of algebra. The roots of the depressed equation (a cubic with no squared term) can be written with explicit formulas.

In other unit systems (2)

One real root (

D > 0

)

t = \sqrt[3]{-\frac{Q}{2} + \sqrt{D}} + \sqrt[3]{-\frac{Q}{2} - \sqrt{D}}

Cardano's formula

Three real roots (

D \le 0

)

t_k = 2\sqrt{-\frac{P}{3}}\cos\left(\frac{\theta}{3} - \frac{2\pi k}{3}\right),\quad \cos\theta = \frac{3Q}{2P}\sqrt{-\frac{3}{P}},\quad k = 0,\,1,\,2

trigonometric form

When roots are close together, the formulas lose accuracy, because they subtract numbers that are almost equal. So we refine the roots with Newton steps, $Z \leftarrow Z - F(Z)/F'(Z)$ , where $F(Z)$ is the left-hand side of the equation. Then we discard the roots with $Z \le B$ and choose the one we need by the rules of the main section.

Where the formula for the fugacity coefficient comes from

Take a phase of $n$ moles of a mixture in a volume $V$ . The fugacity coefficient of a component is an integral of the equation of state, taken at constant temperature from the given volume to infinity, where any gas is ideal:

Fugacity coefficient from an equation of state

\ln\varphi_i = \int_V^{\infty}\left[\frac{1}{RT}\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_{j\ne i}} - \frac{1}{V}\right]\dd V - \ln Z

where:

$V$ total volume of the phase, m³
$n_i$ amount of component $i$ in the phase, kmol
$p(T, V, n_1, n_2, \ldots)$ the equation of state written for $n$ kilomoles of the mixture: it contains $n^2a = \sum_i\sum_j n_in_ja_{ij}$ and $nb = \sum_i n_ib_i$ , MPa

Units: dimensionless; $R$ , $T$ , $V$ , $p$ in consistent units

Holds when

We integrate at constant temperature and constant amounts of the other components.
At infinite volume the phase is an ideal gas.

If we put in a cubic equation and differentiate with respect to $n_i$ using the mixing rules, we get the explicit formulas of the main section. They can be checked by numerical integration. For a pure substance, the equality $f^{\,L} = f^{\,V}$ is equivalent to Maxwell's equal-area rule: $\int_{V_L}^{V_V} p\,\dd V_m = p_s\,(V_V - V_L)$ . In words, the horizontal line cuts off two equal areas from the loop of the isotherm.

Other equations and refinements

The virial equation $Z = 1 + \beta_2/V_m + \beta_3/V_m^2 + \ldots$ has $\beta_2$ and $\beta_3$ as its second and third virial coefficients. Molecular theory justifies it rigorously, but with two or three coefficients it suits only a gas of moderate density.
Multi-parameter equations (for example, Benedict–Webb–Rubin) are accurate for pure light gases, but for mixtures they need a mixing rule for each constant.
Equations for associating substances (substances whose molecules link up), for example CPA (cubic plus association), add a term for hydrogen bonds to the cubic equation. They are used when water, methanol or glycols matter in the calculation. Methanol and glycols are hydrate inhibitors (chemicals that stop gas hydrates, ice-like solids, from forming).

Refined dependence m(ω) of the Peng–Robinson equation for heavy components

m = 0.379642 + 1.48503\,\omega - 0.164423\,\omega^2 + 0.016666\,\omega^3

where:

$m$ coefficient of the temperature function $\alpha(T)$ , —
$\omega$ acentric factor, —

Units: dimensionless

Holds when

The authors of the equation proposed it in 1978 for substances with $\omega > 0.49$ . For smaller $\omega$ the formula of 1976 is used.

At $\omega = 0.49$ the two formulas almost agree (1.066 and 1.070). For heavy pseudo-components with $\omega$ about 1 the difference is already noticeable. Check the documentation of a program to see which formula it uses.

Connections to other topics

Ideal and real gas: the gas deviation factor (Z) — the gas deviation factor, the generalised chart and the law of corresponding states, with which the equation is compared.
Phase behaviour: pure substances and binary mixtures — the vapour pressure, critical point and phase envelope that the equation reproduces.
Composition of reservoir fluids — the components, the C7+ fraction and the pseudo-components, for which the constants of the equation are needed.
Vapour–liquid equilibrium — how to calculate the compositions and fractions of the equilibrium phases of a mixture from fugacities.
Equilibrium ratios: correlations and equations of state — equilibrium ratios from correlations and from an equation of state, requirements for compositional data.
Types of reservoir simulation model — the compositional model, which solves an equation of state in every cell.

Check yourself

1
For a certain phase a cubic equation has given three roots: $Z$ = 0.045; 0.19 and 0.71. Which root describes the gas phase?
1. A.Their average
2. B.0.19
3. C.0.71
4. D.0.045
Show the answer
Correct answer: C. 0.71
The largest root belongs to the gas, and the smallest to the liquid. The middle root lies on the section of the loop where the pressure rises together with the volume, and it has no physical meaning. If you do not know which of the two phases is stable, compare their Gibbs energies (or fugacities).
2
At $T = 0.7\,T_c$ the vapour pressure of a substance is $0.05\,p_c$ . Find the acentric factor. Which hydrocarbon in the table of the lesson is this substance similar to?
Show the answer
$\omega = -\lg 0.05 - 1 = 1.301 - 1 = 0.301$ , the same as for n-hexane. The vapour pressure is lower than for a “simple” substance ( $0.1\,p_c$ ), because the molecule is long and asymmetric.
3
A volume translation fitted to the density of stock-tank oil was added to a calculation with the Peng–Robinson equation. What will change?
1. A.Both the saturation pressure and the oil density
2. B.Nothing
3. C.Only the densities and volumes of the phases
4. D.Only the saturation pressure
Show the answer
Correct answer: C. Only the densities and volumes of the phases
The translation multiplies the fugacity coefficient of each component in both phases by the same factor, $\e^{-c_ip/(RT)}$ . The ratios $\varphi_i^{\,L}/\varphi_i^{\,V}$ , the equilibrium ratios, saturation pressures and phase compositions stay as they were. Only the molar volumes and densities change.
4
Why can the van der Waals equation not give the density of liquid n-butane, 584 kg/m³, at 15.6 °C? For n-butane $T_c$ = 425.2 K, $p_c$ = 3.80 MPa, $M$ = 58.12 g/mol.
Show the answer
The co-volume is $b = RT_c/(8p_c) = 0.008314\cdot 425.2/30.4 = 0.116$ m³/kmol. In this equation the molar volume of any phase is greater than $b$ , so the density is less than $M/b = 58.12/0.116 \approx 500$ kg/m³. The equation cannot reach the measured 584 kg/m³ at any pressure, because its $Z_c = 0.375$ is too large and its co-volume is overestimated.
5
For which pair of components is the binary interaction coefficient in the Peng–Robinson equation usually the largest?
1. A.Methane — ethane
2. B.Propane — n-butane
3. C.n-Hexane — n-heptane
4. D.Carbon dioxide — methane
Show the answer
Correct answer: D. Carbon dioxide — methane
For hydrocarbons of similar size the $k_{ij}$ are almost zero, because the geometric mean describes the attraction between similar molecules well. Carbon dioxide interacts with hydrocarbons differently from the way they interact with one another. For pairs with it, $k_{ij}$ is about 0.1.
6
For a dry gas the Peng–Robinson equation gave $Z = 0.88$ and the generalised chart gave 0.90. How much do the estimates of the amount of gas in the reservoir differ, and which estimate would you trust more?
Show the answer
The amount of gas in a given pore volume is inversely proportional to $Z$ : $0.90/0.88 = 1.023$ , so the equation gives 2.3% more gas. For a dry “sweet” gas, the chart is built from measurements on such gases and is usually more accurate than an untuned equation. Best of all is to measure $Z$ in the laboratory on a sample.

Updated: 7 October 2026