Equilibrium ratios and flash calculations

Sofia, the new facilities engineer, asks Ethan a question he cannot answer from black-oil tables. 'If we drop the separator to 115 psia and run it at 90 °F in summer, what comes out of the gas line and what stays in the oil? I need compositions for the gas sales contract, not just a GOR.' Nadia overhears and smiles. 'Then you need a flash calculation,' she says. 'Ethan, it is time you wrote one.'

Phase equilibrium and the K-value

When a hydrocarbon mixture splits into gas and liquid at a given pressure and temperature, each component distributes itself between the two phases. At equilibrium, the chemical potential, or equivalently the fugacity, of each component is the same in both phases. The practical measure of the distribution is the equilibrium ratio, or K-value, of each component: its mole fraction in the gas, yᵢ, divided by its mole fraction in the liquid, xᵢ. Light components such as methane have K-values much greater than 1, because they prefer the gas. Heavy components have K-values much less than 1. Components with K near 1 split roughly equally.

K-values depend on pressure, temperature and, at higher pressures, on the composition of the mixture. At low pressures, they follow Raoult's law: K is the component's vapor pressure divided by the system pressure, so K rises with temperature and falls with pressure. As pressure increases toward the convergence pressure of the mixture, all K-values approach 1, because the gas and liquid become identical at the critical point. Charts of K-values against pressure for each component, at several temperatures and convergence pressures, were the standard tool before computers, and the GPSA Engineering Data Book still publishes them.

The Wilson correlation

For quick calculations and as starting values for equation-of-state flashes, the Wilson correlation estimates K-values from each component's critical pressure, critical temperature and acentric factor. It is accurate at low pressures, below about 3.5 MPa (500 psia), and becomes increasingly approximate at higher pressures, because it does not know about convergence. It is the standard first guess in almost every flash program.

Temperatures in the Wilson equation must be absolute, in kelvin or degrees Rankine, and the critical and system pressures in the same units. The constant 5.373 comes from fitting the acentric factor definition to vapor pressures; it equals 7/3 × ln 10. For a heavy pseudo-component such as C₇+, the critical properties come from a characterization correlation, discussed in Lesson 12.

Material balance for a flash

Consider one mole of feed with composition zᵢ that splits into V moles of gas and L = 1 − V moles of liquid. The material balance for component i is zᵢ = V yᵢ + L xᵢ. Combining this with yᵢ = Kᵢ xᵢ gives each phase composition in terms of V. Since the liquid and gas mole fractions must each sum to 1, the difference of the two sums must be zero. This gives the Rachford-Rice equation, a single equation in one unknown, V.

The Rachford-Rice function is monotonic in V: it decreases steadily as V increases. That makes it easy and safe to solve. If f(0) is negative, the feed is all liquid, below its bubble point curve (a subcooled liquid). If f(1) is positive, the feed is all gas, above its dew point curve. Otherwise there is one root between 0 and 1, which can be found by bisection, by Newton's method, or by more robust bounded methods for difficult cases with very large or very small K-values.

Bubble point and dew point calculations

The same equations give saturation pressures. At the bubble point, V = 0 and the liquid composition equals the feed. The first bubble of gas has composition yᵢ = Kᵢ zᵢ, so the condition is that the sum of Kᵢ zᵢ equals 1. At the dew point, V = 1 and the first drop of liquid has composition xᵢ = zᵢ/Kᵢ, so the sum of zᵢ/Kᵢ equals 1. Because K-values depend on pressure, these are solved by iteration: guess a pressure, calculate K-values, check the sum, and adjust the pressure. With Wilson K-values, the bubble point equation can be solved directly, because each K is inversely proportional to pressure.

# Wilson K-values and a Rachford-Rice flash of a reservoir oil at separator conditions
import math
SI = False  
comp = {  # z, Tc (R), Pc (psia), acentric factor   (approximate critical data)
    'C1': (0.36, 343.0, 667.8, 0.011), 'C2': (0.10, 549.6, 707.8, 0.099),
    'C3': (0.07, 665.7, 616.3, 0.152), 'nC4': (0.05, 765.3, 550.7, 0.200),
    'nC5': (0.04, 845.4, 488.6, 0.251), 'C6': (0.03, 913.3, 436.9, 0.301),
    'C7+': (0.35, 1180.0, 290.0, 0.520)}
p, T = 115.0, 90 + 459.67                     # separator pressure (psia) and temperature (R)
K = {c: (pc / p) * math.exp(5.373 * (1 + w) * (1 - tc / T)) for c, (z, tc, pc, w) in comp.items()}

def rr(V):                                    # Rachford-Rice objective
    return sum(z * (K[c] - 1) / (1 + V * (K[c] - 1)) for c, (z, *_) in comp.items())

lo, hi = 0.0, 1.0
for _ in range(60):                           # bisection on the gas fraction
    mid = 0.5 * (lo + hi)
    if rr(mid) > 0: lo = mid
    else: hi = mid
V = 0.5 * (lo + hi)
print(f"separator at {p:.0f} psia, {T - 459.67:.0f} F" + (f" ({p * 6.895:.0f} kPa, {(T - 491.67) / 1.8:.0f} C)" if SI else "") + f": gas mole fraction {V:.4f}")
print("comp      K        x (liquid)   y (gas)")
for c, (z, *_) in comp.items():
    x = z / (1 + V * (K[c] - 1)); y = K[c] * x
    print(f"{c:>4} {K[c]:9.4f} {x:11.4f} {y:11.4f}")
print(f"check: sum x = {sum(z / (1 + V * (K[c] - 1)) for c, (z, *_) in comp.items()):.4f}")

separator at 115 psia, 90 F: gas mole fraction 0.4983
comp      K        x (liquid)   y (gas)
  C1   44.7669      0.0158      0.7065
  C2    6.1594      0.0280      0.1725
  C3    1.4509      0.0572      0.0829
 nC4    0.3817      0.0723      0.0276
 nC5    0.1142      0.0716      0.0082
  C6    0.0373      0.0577      0.0021
 C7+    0.0002      0.6975      0.0002
check: sum x = 1.0000

At 115 psia and 90 °F, the flash puts 49.8 mole percent of the feed into the gas. Methane has a K-value of about 45 and leaves almost entirely in the gas, which is 70.7 percent methane and 17.3 percent ethane. The C₇+ fraction has a K-value of 0.0002 and stays in the liquid. Propane, with K about 1.45, splits between both phases. Sofia now has gas compositions to compare with the sales contract's heating value and hydrocarbon dew point limits.

Multistage separation

A multistage separator train is a series of flashes: the liquid from each stage is the feed to the next. Each flash gives a gas composition and volume and a liquid composition. The gas volumes added together, divided by the final stock-tank liquid volume, give the total GOR. The stock-tank liquid density, from its composition and a density correlation or an equation of state, gives its volume and API gravity. This is how process simulators and compositional reservoir models calculate the separator test results of Lesson 9 for any proposed separator conditions.

Sales gas and the hydrocarbon dew point

Flash calculations answer many facilities questions. A common one is the hydrocarbon dew point of sales gas. Pipelines limit the temperature at which liquid hydrocarbons can condense from the gas at pipeline pressures, because liquids that drop out in the pipeline reduce capacity, cause slugging and can damage compressors. Gas plants meet the specification by chilling the gas to condense heavier components, which are recovered as natural gas liquids. The dew point of the residue gas is found by solving the dew point equation for temperature at a given pressure, using K-values from an equation of state. Because the dew point of a lean gas is extremely sensitive to tiny amounts of heavy components, gas analyzes for dew point work must be extended to C₉ or C₁₀ and beyond; a few parts per million of C₁₀ can move the dew point by several degrees.

Flash calculations also give the heating value of the gas, from the heating value of each component weighted by its mole fraction in the gas phase, and the liquid yield from each separator stage. A compositional calculation is also the way to evaluate changes such as adding a compressor, chilling the inlet gas or changing the separator temperature in winter and summer. A useful feature of the Rachford-Rice method is that it extends to V outside 0 to 1, the so-called negative flash. When the feed is single phase, a solution with V below 0 or above 1 still satisfies the equations and tells the program how far the feed lies from the phase boundary, which helps compositional simulators converge when a grid block crosses from single phase to two phases.

Engineers often check a flash result by hand in a few seconds: every component's liquid and gas mole fractions should sum to 1, the material balance zᵢ = V yᵢ + (1 − V) xᵢ should close for each component, and the K-values should fall smoothly from the lightest to the heaviest component. A result that fails any of these checks points to an input error, such as compositions that do not sum to 1, gauge instead of absolute pressure, or a temperature left in degrees Fahrenheit or Celsius.

Stability and phase splitting

A flash calculation assumes the mixture splits into two phases. Before trusting the result, modern programs run a stability test, usually Michelsen's tangent plane analysis, to check whether a single phase is stable at the given conditions or whether splitting lowers the Gibbs energy. If a single phase is stable, no flash is needed. The stability test also gives good initial K-values for the flash. Some systems, such as reservoir oils mixed with carbon dioxide at low temperatures, can form three phases, two liquids and a gas, and need a three-phase flash.

Equilibrium ratios from equations of state

Wilson and chart K-values are adequate for low-pressure separators. At reservoir pressures, K-values depend strongly on composition and must be calculated from an equation of state. The procedure, called successive substitution, starts from Wilson K-values, solves Rachford-Rice for V and the phase compositions, then uses the equation of state to calculate the fugacity of each component in each phase. New K-values are obtained by multiplying the old ones by the ratio of liquid to gas fugacity. The process repeats until the fugacities are equal. Near the critical point, convergence slows, and programs switch to Newton methods. Lesson 11 develops the equation of state that supplies these fugacities.