Vapour–liquid equilibrium

Quick refresher

The equilibrium ratio $K_i = y_i/x_i$ shows how many times richer in a component the gas is than the liquid. Here $y_i$ and $x_i$ are the mole fractions of the component in the gas and in the liquid. Light components have $K_i > 1$ , heavy ones $K_i < 1$ .
At low pressure, for components below their critical temperature, $K_i \approx p_{s,i}/p$ , where $p_{s,i}$ is the vapour pressure of the pure component (Raoult's and Dalton's laws). At reservoir pressures $K_i$ also depends on the compositions of the phases. It is calculated from an equation of state (a formula that links the pressure, volume and temperature of a fluid) as a ratio of fugacity coefficients.
For known $K_i$ , the gas fraction $V$ is the root of the Rachford–Rice equation, the value of $V$ that makes it equal zero. The function is monotonic (it always falls as $V$ grows), and a root between 0 and 1 exists only if the mixture is two-phase.
Bubble point pressure: $\sum z_iK_i = 1$ . Dew point: $\sum z_i/K_i = 1$ . Here $z_i$ is the mole fraction of the component in the whole mixture. Both conditions are the same Rachford–Rice equation, taken at $V = 0$ and at $V = 1$ .
The stability test is the strict way to check whether a mixture is single-phase: a trial phase with $\sum Y_i > 1$ means that the mixture splits into two phases.
The equation-of-state calculation is a loop: $K_i$ → Rachford–Rice → phase compositions → roots $Z$ (compressibility factors) and fugacity coefficients → new $K_i$ . It repeats until the fugacities of each component are equal in the two phases.

Key relationships

K_i = \dfrac{y_i}{x_i} = \dfrac{\varphi_i^{\,L}}{\varphi_i^{\,V}}

equilibrium ratio; the second equality follows from equality of fugacities

\sum_i \dfrac{z_i\,(K_i - 1)}{1 + V(K_i - 1)} = 0

Rachford–Rice equation for the mole fraction of gas $V$

x_i = \dfrac{z_i}{1 + V(K_i - 1)},\quad y_i = K_ix_i

compositions of the equilibrium phases

\sum_i z_iK_i = 1,\qquad \sum_i z_i/K_i = 1

bubble point pressure and dew point pressure

Common pitfalls

Using ideal $K_i = p_{s,i}/p$ at reservoir pressures and for components above their critical temperature.
Solving for the gas fraction from the single condition $\sum x_i = 1$ . It always holds at $V = 0$ , so the search lands on this false root.
Confusing the mole fraction of gas with the volume fraction: in the lesson example, gas making up 32% of the moles occupies 52% of the volume.
Treating a root $V$ outside the interval from 0 to 1 as a program error. For the given $K_i$ it is a sign of a single-phase mixture.
Mistaking a calculation that converges to identical phases ( $K_i = 1$ ) for a two-phase equilibrium.

After this lesson you can

explain what phase equilibrium and the equilibrium ratio $K_i$ are, and why at reservoir pressures it depends not only on pressure and temperature but also on composition
calculate the gas fraction and the phase compositions from known $K_i$ by solving the Rachford–Rice equation, and check the result
write the conditions for the bubble point pressure and the dew point in terms of $K_i$ and explain how they are found by iteration (improving a first guess again and again)
check whether a mixture is single-phase from the signs of the Rachford–Rice function and with the stability test
describe the equation-of-state equilibrium loop, from the initial $K_i$ to equal fugacities (a fugacity measures the tendency of a component to escape from a phase), and explain what slows its convergence

Phase equilibrium and equilibrium ratios

Imagine a cell in which an oil has separated into gas and liquid. Molecules keep crossing the boundary between the phases in both directions. For each component, the flow from the liquid into the gas can become equal to the flow back. Then the fractions and compositions of the phases stop changing, and the phases are in equilibrium. In a reservoir, pressure and temperature change slowly, over days and years. Equilibrium in the pores is reached much faster, so we treat the gas and oil at every point of the reservoir as being in equilibrium. We make the same assumption for a separator (the vessel that splits the well stream into gas and liquid). For a well or a pipeline, we make it for each section.

The equilibrium ratio (K-value) describes how a component is shared between the phases. It is the ratio of the component's mole fractions in the gas and in the liquid. We use the same symbols as in the lesson Composition of reservoir fluids: $z_i$ is the mole fraction of the component in the whole mixture, $x_i$ in the liquid and $y_i$ in the gas. A component with $K_i > 1$ “prefers” the gas; one with $K_i < 1$ prefers the liquid.

Equilibrium ratio

K_i = \frac{y_i}{x_i}

where:

$K_i$ equilibrium ratio of component $i$ , —
$y_i,\ x_i$ mole fractions of the component in the equilibrium gas and liquid, fraction

Units: dimensionless

Holds when

The gas and the liquid are in equilibrium at the same pressure and temperature.
The name “constant” is conventional: $K_i$ changes with pressure and temperature, and at high pressures also with the compositions of the phases.

For methane at reservoir conditions $K$ is usually greater than 1, for heavy hydrocarbons much less than 1. The more the equilibrium ratios of the components differ, the more cleanly the mixture separates into a light gas and a heavy liquid.

For a mixture of similar molecules at low pressure, the equilibrium ratios follow from Raoult's and Dalton's laws (bubble point pressure of an ideal mixture). The partial pressure of a component above the liquid, $x_ip_{s,i}$ , equals its partial pressure in the gas, $y_ip$ . (The partial pressure is the part of the total pressure that this one component supplies.)

Equilibrium ratio of an ideal mixture

K_i = \frac{p_{s,i}(T)}{p}

where:

$p_{s,i}(T)$ vapour pressure of the pure component (the pressure at which it boils) at the temperature of the mixture, MPa
$p$ absolute pressure in the system, MPa

Units: absolute pressures, in the same units

Holds when

An ideal solution (Raoult's law) and an ideal gas (Dalton's law): low pressures, and molecules close in size and structure.
The temperature is below the critical temperature of the component (the highest temperature at which it can be a liquid). Above it there is no vapour pressure.

The ideal equilibrium ratio depends only on pressure and temperature. When the pressure is halved, it doubles. It does not depend on the composition of the mixture at all.

When the vapour pressure of a component is unknown, we estimate it from a straight line on a plot of $\ln p_s$ against $1/T$ . The line is drawn through the critical point and through the point at $T = 0.7\,T_c$ . The position of that point is set by the acentric factor, a number that shows how far the shape of the molecule departs from a sphere (acentric factor). This gives Wilson's formula: the ideal equilibrium ratio written with the critical constants (the critical pressure and temperature).

Initial estimate of equilibrium ratios by Wilson's formula

K_i = \frac{p_{c,i}}{p}\,\exp\left[5.373\,(1 + \omega_i)\left(1 - \frac{T_{c,i}}{T}\right)\right]

where:

$p_{c,i},\ T_{c,i}$ critical pressure and temperature of the component, MPa, K
$\omega_i$ acentric factor of the component, —
$p,\ T$ absolute pressure and temperature of the mixture, MPa, K

Units: pressures in the same units, temperatures in kelvins

Holds when

This is the ideal ratio $p_{s,i}/p$ , with the vapour pressure found from two reference points. The compositions of the phases are not taken into account.
For components above their critical temperature (methane and nitrogen at reservoir temperature), the “vapour pressure” is only the straight line continued beyond the critical point. It has no physical meaning.
It works as a first guess for iterations, and as a rough estimate at low pressures. The more accurate empirical correlations are in the lesson Equilibrium ratios: correlations and equations of state.

The number 5.373 follows from the definition of the acentric factor. The straight line $\ln(p_s/p_c) = k\,(1 - T_c/T)$ passes through the point $T = 0.7\,T_c$ , where $\log_{10}(p_s/p_c) = -(1 + \omega)$ . So $k = 2.3026\cdot(1 + \omega)\cdot 0.7/0.3 = 5.373\,(1 + \omega)$ . In the literature and in software the formula is usually printed with the rounded number 5.37. The ratios then differ by less than one per cent (for the heaviest fractions at low temperature, by 1–2%), which does not matter for a first guess. The examples of the lesson use 5.373.

At reservoir pressures, ideal equilibrium ratios can be wrong by a factor of several. At 90 °C methane is far above its critical temperature, so its “vapour pressure” is a fiction. The liquid holds methane because methane dissolves in it, not because of evaporation. There is a second reason. At high pressure the gas becomes dense and dissolves heavy components, so the equilibrium ratios of all the components move closer together. They tend towards one where the gas and the liquid become identical. The pressure at which the equilibrium ratios of a mixture would all reach one is called the convergence pressure. Empirical tables of equilibrium ratios are built on it (Equilibrium ratios: correlations and equations of state). The strict way is to calculate the ratios with an equation of state, from the equality of the fugacities of each component in the two phases (fugacity): $K_i = \varphi_i^{\,L}/\varphi_i^{\,V}$ . The last section of the lesson covers this.

CalculatedEquilibrium ratios of the lesson's model oil (methane, propane, n-hexane, n-decane) at 90 °C against pressure. Solid lines: calculation with the Peng–Robinson equation up to the bubble point pressure of 15.2 MPa. Dashed lines: Wilson's formula. Both axes are logarithmic.

Flash calculation: the Rachford–Rice equation

The main task in equilibrium work is the flash calculation. A mixture of known composition $z_i$ , at pressure $p$ and temperature $T$ , separates into gas and liquid. We have to find the gas fraction and the compositions of the phases. Suppose the equilibrium ratios at these conditions are known. Take 1 mole of the mixture. $V$ moles pass into the gas, and $L = 1 - V$ moles remain in the liquid. The balance for one component is $z_i = x_iL + y_iV$ : what the mixture contains equals what the two phases contain. Together with $y_i = K_ix_i$ , it gives the phase compositions at once in terms of one unknown, $V$ .

Phase compositions in a flash calculation

x_i = \frac{z_i}{1 + V(K_i - 1)},\qquad y_i = K_i\,x_i

where:

$z_i$ mole fraction of the component in the original mixture, fraction
$V$ mole fraction of the gas phase: how many moles out of each mole of the mixture ended up in the gas (not a volume), fraction
$K_i$ equilibrium ratio at the given $p$ and $T$ , —
$x_i,\ y_i$ mole fractions of the component in the liquid and in the gas, fraction

Units: dimensionless

Holds when

The phases are in equilibrium and the $K_i$ are known.
There are two hydrocarbon phases; formation water does not take part in the hydrocarbon equilibrium and is accounted for separately.
The balances are written per 1 mole of the mixture; the mole fraction of the liquid phase is $L = 1 - V$ .

As the gas fraction $V$ grows, the liquid holds fewer light components (those with $K_i > 1$ ) and more heavy ones (those with $K_i < 1$ ).

We can find $V$ from the condition that the fractions in each phase add up to one. But the condition $\sum x_i = 1$ taken alone is already true at $V = 0$ (then $x_i = z_i$ ), and $\sum y_i = 1$ is already true at $V = 1$ . A numerical search easily lands on these false roots. It is safer to subtract one condition from the other: $\sum y_i - \sum x_i = 0$ . This is the Rachford–Rice equation, proposed in 1952.

Rachford–Rice equation

g(V) = \sum_i \frac{z_i\,(K_i - 1)}{1 + V(K_i - 1)} = 0

where:

$g(V)$ Rachford–Rice function: the difference between the sums of the mole fractions in the gas and in the liquid, —
$V$ mole fraction of the gas phase, —
$z_i,\ K_i$ composition of the mixture and equilibrium ratios, —

Units: dimensionless

Holds when

The equilibrium ratios at the given $p$ and $T$ are known and do not depend on $V$ (in an equation-of-state calculation they are refined in an outer loop).

The function $g(V)$ decreases all the way (monotonically), so it has a single root. Its values at the two ends tell us the state of the mixture: $g(0) = \sum z_iK_i - 1$ and $g(1) = 1 - \sum z_i/K_i$ . If $g(0) > 0$ and $g(1) < 0$ , the root lies between 0 and 1 and the mixture is two-phase. The equation is solved by Newton's method (improving a guess step by step), backed up by bisection (repeatedly halving an interval that contains the root); see the advanced section.

Worked exampleGas fraction and phase compositions of the model oil at 8 MPa

Given (assumed data)

Mixture, mole fractions: methane 0.45; propane 0.15; n-hexane 0.15; n-decane 0.25
Conditions: 8 MPa (absolute), 90 °C
Equilibrium ratios: 3.416; 0.6471; 0.09523; 0.009552, calculated with the Peng–Robinson equation of state (the last section shows how)

Solution

Check the state of the mixture. $g(0) = \sum z_i(K_i - 1) = 0.45\cdot 2.416 - 0.15\cdot 0.3529 - 0.15\cdot 0.9048 - 0.25\cdot 0.9904 = 0.651 > 0$ , so the mixture is at a pressure below its bubble point pressure. $g(1) = \sum z_i(1 - 1/K_i) = -27.1 < 0$ , so it is above the dew point. So there are two phases and $0 < V < 1$ .
Newton's method, starting from $V_0 = 0.5$ : $g(0.5) = -0.3103$ , derivative $g'(0.5) = -\sum z_i(K_i - 1)^2/[1 + 0.5(K_i - 1)]^2 = -1.9383$ ; $V_1 = 0.5 - (-0.3103)/(-1.9383) = 0.3399$ . Two more steps give $V_2 = 0.3200$ and $V_3 = 0.3200$ . The root is found: $V = 0.3200$ .
Liquid: $x_i = z_i/[1 + V(K_i - 1)]$ = 0.2538; 0.1691; 0.2111; 0.3660. Gas: $y_i = K_ix_i$ = 0.8670; 0.1094; 0.0201; 0.0035. Sums: $\sum x_i = 1.0000$ , $\sum y_i = 1.0000$ .
Where each component went: the share $Vy_i/z_i$ passed into the gas. This is 62% of the methane, 23% of the propane, 4.3% of the hexane and 0.45% of the decane.
Molar masses of the phases: liquid $M_L = \sum x_iM_i = 81.8$ g/mol, gas $M_V = 20.96$ g/mol.

Answer: 32% of the moles of the mixture pass into the gas: the gas is 87% methane, while the liquid is 25% dissolved methane (by moles) and contains almost all the decane.

What it means. 62% of the methane goes into the gas, but 38% of it stays dissolved in the liquid. The heavy components stay almost entirely in the liquid. The gas also carries off a noticeable share of the propane and a little hexane. At the surface they partly condense, and the separator calculation takes this into account (Separator calculations and optimum separator pressure). Newton's method converges in three steps because $g(V)$ is smooth and monotonic. The data in this example are assumed for illustration.

CalculatedThe Rachford–Rice function of the model oil at 90 °C for two pressures (equilibrium ratios from the Peng–Robinson equation). The function decreases monotonically; the root is the mole fraction of gas: 0.320 at 8 MPa and 0.163 at 12 MPa. The curves are cut off below −2: as

V \to 1

they plunge to large negative values because of the decane, which has a very small

K

Bubble point pressure and dew point

The two edges of the two-phase region are special cases of the same calculation. At the bubble point there is no gas yet ( $V = 0$ ). The liquid has the composition of the mixture, $x_i = z_i$ , and the first gas bubble has the composition $y_i = K_iz_i$ . At the dew point ( $V = 1$ ) the gas has the composition of the mixture, and the first drop of liquid has the composition $x_i = z_i/K_i$ . We require that the fractions in the new phase add up to one. This condition becomes an equation for the pressure.

Conditions for the bubble point pressure and the dew point

\sum_i z_iK_i(p_b, T) = 1,\qquad \sum_i \frac{z_i}{K_i(p_d, T)} = 1

where:

$p_b$ bubble point pressure (for an oil, the saturation pressure), MPa
$p_d$ dew point pressure (the pressure at which condensation begins), MPa
$K_i(p, T)$ equilibrium ratios at this pressure and the composition of the phase that appears, —
$z_i$ composition of the mixture, fraction

Units: absolute pressures

Holds when

The $K_i$ are calculated for the pair “mixture and new phase”: at the bubble point pressure, for a liquid of composition $z_i$ and a bubble of composition $K_iz_i$ .
At one temperature a gas condensate can have two dew points: an upper one (retrograde) and a lower one.

These are the conditions $g(0) = 0$ and $g(1) = 0$ : the Rachford–Rice function at its two ends. The $K_i$ depend on the pressure and on the composition of the phase that appears, so the equations are solved by iteration. Choose a pressure, calculate the $K_i$ , check the sum and improve the pressure. For ideal ratios $K_i = p_{s,i}/p$ , the first equation gives the explicit formula $p_b = \sum z_ip_{s,i}$ . This is Raoult's law from the lesson Phase behaviour: pure substances and binary mixtures.

Worked exampleBubble point pressure of the model oil: ideal estimate and equation of state

Given (assumed data)

Mixture and temperature: the model oil of the example above, 90 °C ( $T = 363.15$ K)
Equation of state: Peng–Robinson; $k_{ij}$ (the binary interaction coefficients, which tune how pairs of components interact): methane–hexane 0.03, methane–decane 0.04, the rest 0 (assumed values of a typical order)

Solution

Ideal estimate. By Wilson's formula $K_i = A_i/p$ with $A_i = p_{c,i}\exp[5.373(1 + \omega_i)(1 - T_{c,i}/T)]$ : $A_i$ = 60.77; 3.788; 0.1875; 0.007695 MPa (methane, propane, hexane, decane). The condition $\sum z_iA_i/p_b = 1$ gives $p_b = \sum z_iA_i = 27.9$ MPa.
Equation of state. Adjust the pressure until the fugacities of each component are equal in the liquid of composition $z_i$ and in the gas bubble. This gives $p_b = 15.2$ MPa. The equilibrium ratios at this pressure are 1.946; 0.6078; 0.1625; 0.03452. Check: $\sum z_iK_i = 1.000$ .
Composition of the first bubble $y_i = z_iK_i$ : methane 0.876, propane 0.091, hexane 0.024, decane 0.0086.
Cross-check with the stability test (next section): at 15.0 MPa the oil is unstable ( $\sum Y_i = 1.006 > 1$ ), and at 15.5 MPa it is stable ( $\sum Y_i = 0.993 < 1$ ). The boundary lies between them.
Sensitivity: if all the $k_{ij}$ are set to zero, $p_b = 14.0$ MPa.

Answer: The bubble point pressure of the model oil at 90 °C is about 15.2 MPa. The ideal equilibrium ratios would give 27.9 MPa, that is 84% more.

What it means. At the bubble point pressure the ideal formula gives methane the equilibrium ratio $A_1/p_b = 4.0$ , while the equation of state gives 1.946. In a dense system the liquid holds methane much more strongly than its “vapour pressure” would suggest (at this temperature methane has no real vapour pressure). The interaction coefficients $k_{ij}$ of methane with the heavy components change the bubble point pressure by 1.2 MPa. That is why they are the ones adjusted to match the measured bubble point pressure when the equation is tuned (Cubic equations of state). The data in this example are assumed for illustration.

Is the mixture single-phase: the stability test

Before we split a mixture into phases, we must know whether it splits at all. If the equilibrium ratios are known, the signs of $g(0)$ and $g(1)$ give the answer. If $\sum z_iK_i \le 1$ , the mixture is a liquid (the pressure is not below the bubble point pressure). If $\sum z_i/K_i \le 1$ , it is a gas. Otherwise it is two-phase. But in an equation-of-state calculation the equilibrium ratios depend on the compositions of phases that do not exist yet, so this check goes in a circle.

The stability test gives a strict answer. Its idea: a mixture is stable if no small portion of a new phase, of any composition, can separate out and lower the Gibbs energy (the energy that a mixture at fixed pressure and temperature tends to make as small as possible). We check this numerically. We search for the composition of a trial phase that forms most “easily” from the mixture. In the trial phase the fugacity of each component must be the same as in the original mixture. If this requires a total amount of trial phase, $\sum Y_i$ , greater than one, the mixture is unstable and will split.

Stability test: iterating the trial phase

\ln Y_i = \ln z_i + \ln\varphi_i(\mathbf{z}) - \ln\varphi_i(\mathbf{y}),\qquad y_i = \frac{Y_i}{\sum_j Y_j}

where:

$Y_i$ “amounts” of the components of the trial phase, not yet scaled to add up to one; the initial values are $z_iK_i$ (gas-like trial) or $z_i/K_i$ (liquid-like) with $K_i$ from Wilson's formula, —
$y_i$ composition of the trial phase, fraction
$\varphi_i(\mathbf{z}),\ \varphi_i(\mathbf{y})$ fugacity coefficients of the component in the original mixture and in the trial phase (from the equation of state), —

Units: dimensionless

Holds when

The iterations are repeated until $Y_i$ stops changing. Both trials are checked: the gas-like one and the liquid-like one.
If the trial converges to the composition of the mixture itself ( $y_i = z_i$ , the trivial solution), it tells us nothing.
The mixture is unstable (two-phase) if $\sum Y_i > 1$ for at least one non-trivial trial.

The composition $y_i$ found for an unstable mixture is a good first guess for a flash calculation. At a saturation pressure it is exactly the composition of the first bubble or the first drop. So the stability test is also a way to find saturation pressures: we change the pressure until $\sum Y_i$ equals one.

Worked exampleDew points of a gas condensate from the stability test

Given (assumed data)

Mixture, mole fractions: methane 0.87; propane 0.07; n-hexane 0.03; n-decane 0.03
Temperature: 90 °C; equation and $k_{ij}$ as in the previous example

Solution

At 30 MPa both trial phases converge to the composition of the mixture itself. This is the trivial solution: the mixture is single-phase, a gas.
At 20 MPa the gas-like trial is trivial, while the liquid-like one converges to a composition different from that of the mixture, with $\sum Y_i = 1.078 > 1$ . The mixture is unstable: liquid drops out of the gas.
Halve the interval 20–30 MPa again and again, keeping a stable mixture at one end and an unstable one at the other. After several tens of halvings (bisections) we get the upper dew point $p_d = 24.7$ MPa.
Composition of the first drop (the trial phase at $p_d$ , scaled to add up to one): methane 0.6603, propane 0.0954, hexane 0.0785, decane 0.1658. The equilibrium ratios $K_i = z_i/x_i$ are 1.318; 0.734; 0.382; 0.181, and the check $\sum z_i/K_i = \sum x_i = 1$ holds.
Repeat the same steps at low pressures: at 0.25 MPa $\sum Y_i = 1.115$ , at 0.20 MPa 0.907. The lower dew point is $p_d = 0.22$ MPa. The first drop there is 95% decane.
Between the dew points the mixture is two-phase. A flash calculation with the same model shows that the liquid occupies up to 6.3% of the volume near 18 MPa.

Answer: At 90 °C the gas condensate has two dew points: the upper one at 24.7 MPa and the lower one at 0.22 MPa.

What it means. Below the upper dew point, retrograde condensation begins: liquid appears as the pressure falls, not as it rises (lesson on phase behaviour). The first drop is still 66% methane. It is a light liquid, close to the gas in composition. That is exactly why the calculation is hard at the upper dew point: the trial phase easily slides into the trivial solution. The lower dew point is an ordinary one: the last liquid is almost pure decane. If a gas-condensate reservoir is produced with the pressure falling below 24.7 MPa, the heavy components will drop out as liquid in the pores (Types of reservoir fluids). The data in this example are assumed for illustration.

Equilibrium calculation with an equation of state

At reservoir pressures the equilibrium ratios do not come from tables. They come from the equilibrium condition: the fugacity of each component in the liquid equals its fugacity in the gas (lesson on equations of state). From $x_i\varphi_i^{\,L}p = y_i\varphi_i^{\,V}p$ it follows that $K_i = \varphi_i^{\,L}/\varphi_i^{\,V}$ . But the fugacity coefficients depend on the compositions of the phases, and the compositions depend on the equilibrium ratios. We break this circle by iteration, in the successive substitution method: guess the ratios, calculate the new ones, put them back in, and repeat.

The loop of an equation-of-state flash calculation

1
Initial equilibrium ratios
Take them from Wilson's formula (or from the previous time step, or from the stability test).
2
Gas fraction and compositions
Solve the Rachford–Rice equation and find $x_i$ and $y_i$ .
3
Equation of state for each phase
From the liquid composition, calculate $A$ , $B$ and the smallest root $Z_L$ (the compressibility factor of the liquid). From the gas composition, calculate the largest root $Z_V$ (the compressibility factor of the gas) (choice of root).
4
Fugacities
Calculate the fugacity coefficients $\varphi_i^{\,L}$ and $\varphi_i^{\,V}$ and the fugacities $f_i^{\,L} = x_i\varphi_i^{\,L}p$ , $f_i^{\,V} = y_i\varphi_i^{\,V}p$ .
5
Check
If the fugacities agree (see the criterion below), the calculation is finished. If not, calculate new ratios $K_i = \varphi_i^{\,L}/\varphi_i^{\,V}$ and go back to the second step.

Convergence criterion of the equilibrium calculation

\varepsilon = \sum_i\left(1 - \frac{f_i^{\,L}}{f_i^{\,V}}\right)^2 < \varepsilon_0

where:

$f_i^{\,L},\ f_i^{\,V}$ fugacities of the component in the liquid and in the gas at the current iteration, MPa
$\varepsilon_0$ tolerance (the allowed error); $10^{-12}$ in the lesson example, —

Units: dimensionless

Holds when

The fugacities are calculated with one equation of state for both phases at the same $p$ and $T$ .

The ratio $f_i^{\,L}/f_i^{\,V}$ equals the ratio of the new equilibrium ratio to the old one, so a small $\varepsilon$ means that the equilibrium ratios have stopped changing.

Worked exampleFlash calculation of the model oil with the Peng–Robinson equation

Given (assumed data)

Mixture and conditions: model oil, 8 MPa (absolute), 90 °C
Equation: Peng–Robinson; component constants from the lesson Cubic equations of state; $k_{ij}$ as in the bubble point pressure example

Solution

Initial ratios from Wilson's formula ( $K_i = A_i/8$ ): 7.60; 0.474; 0.0234; 0.000962.
Rachford–Rice with these ratios: $V^{(1)} = 0.423$ ; liquid $x_i$ = 0.119; 0.193; 0.255; 0.433; gas $y_i$ = 0.902; 0.091; 0.0060; 0.0004.
Equation of state: for the liquid $Z_L = 0.414$ , for the gas $Z_V = 0.898$ . Fugacity coefficients in the liquid: 3.4468; 0.39466; 0.031813; 0.0014841; in the gas: 0.92945; 0.61783; 0.37537; 0.19494.
New ratios $K_i = \varphi_i^{\,L}/\varphi_i^{\,V}$ : 3.708; 0.6388; 0.08475; 0.007613. The criterion is $\varepsilon = 55$ : the fugacities are still far apart, especially for decane.
Further iterations give $\varepsilon$ = 0.064; 0.00093; 0.0000156; … Each step reduces it tens of times. After 9 updates, $\varepsilon < 10^{-12}$ . The ratios are 3.416; 0.6471; 0.09523; 0.009552, and the gas fraction is $V = 0.3200$ . These are the very values from which the flash calculation example started.
Properties of the phases: $Z_L = 0.382$ , $Z_V = 0.874$ ; molar masses 81.8 and 20.96 g/mol; densities $\rho = pM/(ZRT)$ : liquid 567 kg/m³, gas 63.5 kg/m³ (without volume translation, a correction that improves the calculated liquid density).
Molar volumes $ZRT/p$ : liquid 0.1442, gas 0.3299 m³/kmol. The volume of the gas is $0.3200\cdot0.3299 = 0.106$ , of the liquid $0.6800\cdot0.1442 = 0.098$ m³ per kilomole of the mixture: the gas occupies 52% of the volume.

Answer: Equilibrium is reached in 9 iterations: 32% of the moles are in the gas, but 52% of the volume; the density of the liquid is about 567 kg/m³, of the gas about 63.5 kg/m³.

What it means. The first iteration from Wilson's ratios gave the right order of magnitude: 42% of the moles in the gas instead of 32%. But for the heavy components it was wrong several-fold: the equilibrium ratio of decane came out 10 times smaller than the final one. The fugacity coefficient of methane in an oil is greater than one: the liquid “pushes out” methane, which is why its $K > 1$ . In a reservoir simulation model the first ratios come from the previous time step, and a few iterations are usually enough. Convergence is worst near the critical point of the mixture, where the two phases are almost identical. There, successive substitution needs hundreds of iterations, and programs switch to Newton's method. The density of the liquid without volume translation can be off by a few per cent (volume translation). The data in this example are assumed for illustration.

CalculatedThe gas-phase fraction in the model oil at 90 °C below the bubble point pressure of 15.2 MPa, by moles and by volume (flash calculation with the Peng–Robinson equation, without volume translation). The dots show the example at 8 MPa.

What the results give and where the limits of the calculation lie

Volumes and saturations. The equation of state gives the densities and molar volumes of the phases from their compositions. From these we get the volume fractions, which in a model cell become the gas and oil saturations (the shares of the pore space filled with gas and with oil).
Properties for flow. Separate correlations give the viscosity and interfacial tension from the compositions and densities of the phases. Near the critical point the relative permeabilities (how easily each phase flows when both are present) also depend on interfacial tension.
Surface production. The reservoir gas and liquid pass through several stages of separation, and the same flash calculation is done at each stage. The gas–oil ratio and the oil formation volume factor are built up from these calculations (Separator calculations and optimum separator pressure).
Quick estimates. For preliminary calculations we use empirical correlations of equilibrium ratios that take the convergence pressure into account (Equilibrium ratios: correlations and equations of state).

Terms and symbols

Equilibrium ratio $K_i$: The ratio of the mole fractions of a component in the equilibrium gas and liquid, $y_i/x_i$ . It depends on pressure, temperature and the compositions of the phases.
Flash calculation: The calculation of the gas fraction and the compositions of the equilibrium phases of a mixture of known composition, at a given pressure and temperature.
Mole fraction of the gas phase $V$: The fraction of the moles of the mixture that ended up in the gas. Do not confuse it with the volume fraction of gas.
Rachford–Rice equation: The equation $\sum z_i(K_i - 1)/[1 + V(K_i - 1)] = 0$ for the gas fraction. Its left-hand side always decreases as $V$ grows, so the equation is solved reliably.
Stability test: A check of whether a mixture can separate out a new phase and lower its Gibbs energy. Instability means that the mixture splits into two phases.
Trivial solution: An iteration result in which the trial phase, or the “second” phase, has the same composition as the original mixture. It is a sign of a single-phase state, or of difficulties near the critical point.
Successive substitution: Iterations in which the new equilibrium ratios are taken as the ratio of the fugacity coefficients calculated with the old ones.
Convergence pressure: The pressure at which the equilibrium ratios of all the components of a mixture, at a given temperature, would reach one. It is a parameter of empirical correlations of equilibrium ratios.

Common mistakes and misconceptions

Wrong

“The ratios $K = p_s/p$ are good for a reservoir too.”

In fact

Ideal ratios are valid at low pressure and for components below their critical temperature. For the oil of the lesson they overestimated the bubble point pressure by 84%. At reservoir pressures the ratios are calculated from an equation of state.

Wrong

“I'll find $V$ by trial from the condition $\sum x_i = 1$ .”

In fact

This condition holds at $V = 0$ for any ratios, so a search finds a false root. The reliable form is the Rachford–Rice one: the difference $\sum y_i - \sum x_i$ , which is monotonic in $V$ .

Wrong

“I got $V = 1.3$ — a bug in the program.”

In fact

A root outside the interval from 0 to 1, for the given ratios, means a single-phase mixture: $V > 1$ is a gas, $V < 0$ a liquid. The stability test gives the final check.

Wrong

“A third of the moles went into the gas, so the gas occupies a third of the volume.”

In fact

A mole of gas takes more volume than a mole of liquid. In the lesson example at 8 MPa the gas, 32% of the moles, occupies 52% of the volume, and at low pressure the difference is even greater.

Wrong

“The calculation has converged, so the answer is right.”

In fact

Iterations can converge to a trivial solution with identical phases. If the ratios come out close to one, check the mixture with the stability test.

In depth: advanced detail

Why the Rachford–Rice equation is solved reliably

The derivative of the Rachford–Rice function (its slope) is a sum of negative terms, so the function always falls between its poles, $V = 1/(1 - K_{\max})$ and $V = 1/(1 - K_{\min})$ , where the denominators become zero. The root between the poles is unique, and Newton's method backed up by bisection always finds it.

Derivative of the Rachford–Rice function and the Newton step

g'(V) = -\sum_i \frac{z_i\,(K_i - 1)^2}{[1 + V(K_i - 1)]^2} < 0,\qquad V_{n+1} = V_n - \frac{g(V_n)}{g'(V_n)}

where:

$g'(V)$ derivative of the Rachford–Rice function with respect to $V$ , —
$V_n$ approximation of the gas fraction at step $n$ , —

Units: dimensionless

Holds when

At least one equilibrium ratio is greater than one and at least one is less than one. Otherwise there are no poles on both sides, and the mixture is single-phase for any $V$ .
If a Newton step goes outside the current interval in which $g$ changes sign, take a bisection step instead.

The poles lie outside the interval from 0 to 1 (the left one at negative $V$ , the right one at $V > 1$ ). If the root is outside the range from 0 to 1 but still between the poles, the calculation is called a “negative flash”. It shows on which side of the two-phase region the mixture lies, and in some programs it replaces the stability test.

The stability test and Gibbs energy

At constant pressure and temperature, a mixture moves towards the minimum of its Gibbs energy. Suppose a small amount of a new phase of composition $\mathbf{y}$ separates from the mixture. The change in energy is proportional to the distance from the tangent plane to the Gibbs energy surface at the point of the original composition. (The tangent plane touches the surface at that point.) If this distance is negative for some composition, separating such a phase is favourable, because it lowers the energy, and the mixture is unstable.

Tangent plane distance in the form for a trial phase

tm(\mathbf{Y}) = 1 + \sum_i Y_i\left[\ln Y_i + \ln\varphi_i(\mathbf{y}) - \ln z_i - \ln\varphi_i(\mathbf{z}) - 1\right]

where:

$tm$ dimensionless tangent plane distance (up to a factor $RT$ ), —
$Y_i,\ y_i$ amounts (not scaled to add up to one) and composition of the trial phase, —
$z_i$ composition of the original mixture, —

Units: dimensionless

Holds when

The fugacity coefficients come from the equation of state. For the original mixture, take the root with the lower Gibbs energy.

At the stationary points of this function (where its slope is zero), the iteration equation of the stability test holds, and there $tm = 1 - \sum Y_i$ . So $\sum Y_i > 1$ at a non-trivial stationary point means a negative distance, that is, instability.

Connections to other topics

Phase behaviour: pure substances and binary mixtures — Raoult's and Dalton's laws, bubble point pressure and dew point, retrograde condensation.
Cubic equations of state — cubic equations, the choice of root and the fugacity coefficients from which the equilibrium ratios are calculated.
Equilibrium ratios: correlations and equations of state — empirical correlations of equilibrium ratios and the convergence pressure; what composition data they need.
Separator calculations and optimum separator pressure — multistage separation: gas–oil ratio and formation volume factor from successive flash calculations.
Types of reservoir fluids — oils and gas condensates, dew points and condensate dropout in the reservoir.
Types of reservoir simulation model — the compositional model, which runs the stability test and the equilibrium calculation in each cell.

Check yourself

1
A mixture of three components: $z$ = 0.5; 0.3; 0.2; the equilibrium ratios at the given conditions are $K$ = 2.5; 0.8; 0.1. What is the state of the mixture?
1. A.Gas only
2. B.Liquid only
3. C.It cannot be said from these data
4. D.Gas and liquid
Show the answer
Correct answer: D. Gas and liquid
$\sum z_iK_i = 1.25 + 0.24 + 0.02 = 1.51 > 1$ , so the pressure is below the bubble point pressure. $\sum z_i/K_i = 0.2 + 0.375 + 2.0 = 2.575 > 1$ , so it is above the dew point. So $g(0) > 0$ and $g(1) < 0$ , and the root of the Rachford–Rice equation lies between 0 and 1. The mixture is two-phase (the gas fraction is about 0.54).
2
A binary mixture: $z_1 = 0.6$ , $z_2 = 0.4$ , $K_1 = 2.0$ , $K_2 = 0.25$ . Find the mole fraction of gas and the compositions of the phases.
Show the answer
The Rachford–Rice equation: $0.6\cdot 1/(1 + V) - 0.4\cdot 0.75/(1 - 0.75V) = 0$ . This gives $0.6(1 - 0.75V) = 0.3(1 + V)$ and $V = 0.4$ . Liquid: $x_1 = 0.6/1.4 = 0.4286$ , $x_2 = 0.4/0.7 = 0.5714$ ; gas: $y_1 = 2\cdot 0.4286 = 0.857$ , $y_2 = 0.25\cdot 0.5714 = 0.143$ . The sums of the fractions in both phases equal one.
3
How do the gas-phase fraction and the equilibrium ratio of methane change when the pressure of an oil is reduced from the bubble point pressure at constant temperature?
1. A.Both increase
2. B.The gas fraction increases, the methane ratio falls
3. C.The gas fraction increases, the methane ratio does not change
4. D.Both fall
Show the answer
Correct answer: A. Both increase
Below the bubble point pressure more and more gas is released from the oil, and $V$ grows from zero. The equilibrium ratio of methane grows too: in the figure of the lesson, from 1.95 at 15.2 MPa to 3.4 at 8 MPa and 25 at 1 MPa. It would stay unchanged only over a narrow range, where the change is too small to see. For heavy components the dependence may even be non-monotonic (it can fall first and then rise).
4
Why, for the lesson's model oil at 0.2 MPa and 90 °C, is the ideal equilibrium ratio of methane (304) 2.4 times greater than the one calculated from the equation of state (127), while for n-hexane the two estimates almost coincide?
Show the answer
The ideal equilibrium ratio is the vapour pressure divided by the pressure. Hexane at 90 °C is below its critical temperature, so it has a vapour pressure, and at low pressure Raoult's and Dalton's laws work for it. Methane at 90 °C is far above its critical temperature. Its “vapour pressure” is only a straight line continued beyond the critical point. The liquid holds methane because methane dissolves in it, and the ideal formula does not describe that.
5
A stability test of an oil at 16 MPa gave $\sum Y_i = 0.98$ for the gas-like trial, while the liquid-like trial converged to the composition of the oil itself. What does this mean?
1. A.The oil is single-phase: the pressure is above the bubble point pressure
2. B.The oil has turned into a gas
3. C.The calculation did not converge
4. D.The oil is two-phase: gas is released
Show the answer
Correct answer: A. The oil is single-phase: the pressure is above the bubble point pressure
A non-trivial trial with $\sum Y_i < 1$ means that releasing gas would not lower the energy, and a trivial trial tells us nothing. So the oil is stable and single-phase: the pressure is above the bubble point pressure (in the lesson example it is 15.2 MPa).
6
At 10 MPa the lesson's model oil splits as follows: the mole fraction of gas is 0.25, but the gas occupies 38% of the volume. Explain the difference and say what the ratio of the volume and mole fractions depends on.
Show the answer
The volume of a phase is its number of moles times the molar volume $ZRT/p$ . A mole of gas takes more volume than a mole of liquid (at 10 MPa about 1.9 times more: the ratio $Z_V/Z_L$ ). So the gas takes a larger share of the volume than of the moles. The lower the pressure, the greater this ratio: at 8 MPa 32% of the moles occupy 52% of the volume, and near atmospheric pressure the gas occupies almost the whole volume.

Updated: 28 October 2026