By Saad Iqbal
A gas-lift optimization engineer bumps the injection rate on a well, expecting the extra lift gas to lighten the column and pull more oil. Instead, production drops. The VLP curve the team built last quarter used a single average gradient — 0.10 psi/ft, borrowed from a nearby well — because nobody wanted to run the full multiphase calculation by hand. Above a certain gas rate, friction starts winning the fight that elevation used to win alone, and a model that can’t see that turn will send you chasing the wrong lever every time.
The fix isn’t a bigger simulator. It’s the same correlation every commercial nodal-analysis package quietly runs under the hood: Beggs and Brill (1973). It classifies the flow pattern inside the pipe, works out how much of the pipe is actually liquid versus gas at any instant (the holdup), and from that builds a real two-phase pressure gradient instead of a guessed one. This tutorial walks through every step of the method by hand for one tubing segment, then turns it into a Python function so you never have to repeat the arithmetic.
Prerequisites: what you need before you start
- Tubing inner diameter, the segment length, and its inclination from horizontal (90° for a vertical well)
- In-situ (actual, at flowing pressure and temperature) liquid and gas volumetric rates — pull these from a PVT report, or estimate Rs and Bo with Standing’s correlation if you only have surface rates
- Liquid and gas densities, viscosities, and the liquid–gas interfacial tension at the segment’s average pressure and temperature
- Python 3 with NumPy (a calculator works for the single worked example, but you’ll want the script for anything beyond one point)
- The worked example below uses: 2-7/8″ tubing (2.441″ ID), a 1,000 ft vertical segment just above a gas-lift valve, 650 bbl/d in-situ liquid, 56,000 ft³/d in-situ lift gas, ρL = 51.8 lbm/ft³, ρg = 6.9 lbm/ft³, σL = 25 dynes/cm, μL = 3.0 cp, μg = 0.016 cp
That in-situ gas rate is high relative to the liquid on purpose — this segment sits just above a gas-lift injection point, where injected lift gas has pushed the local gas-liquid ratio well above what comes out of solution naturally. That is exactly the regime where guessing a single gradient gets expensive.
Step 1: Convert flow rates to superficial velocities
Beggs-Brill works in superficial velocities — the velocity each phase would have if it alone filled the entire pipe. With pipe area A = π/4 × d²:
vsl = qL / A and vsg = qg / A, with vm = vsl + vsg
For 2.441″ ID tubing, A = 0.03250 ft². Converting 650 bbl/d to 0.04224 ft³/s and 56,000 ft³/d to 0.64815 ft³/s:
vsl = 0.04224 / 0.03250 = 1.30 ft/s • vsg = 0.64815 / 0.03250 = 19.94 ft/s • vm = 21.24 ft/s
Step 2: Calculate the no-slip liquid holdup
The no-slip holdup λL is just the liquid’s share of the combined volumetric flow — it ignores the fact that gas moves faster than liquid and slips past it. It’s not the real in-situ liquid fraction, but every later step is built from it:
λL = vsl / vm = 1.30 / 21.24 = 0.0612
Step 3: Calculate the Froude number and liquid velocity number
Two dimensionless groups decide which flow pattern the pipe is in. The mixture Froude number compares inertial to gravitational forces; the liquid velocity number folds in liquid density and surface tension:
NFr = vm² / (g × d) • NLv = 1.938 × vsl × (ρL/σL)0.25 (field units: ft/s, lbm/ft³, dynes/cm)
NFr = 21.24² / (32.174 × 0.2034) = 68.96 • NLv = 1.938 × 1.30 × (51.8/25)0.25 = 3.02
Step 4: Classify the flow regime
Beggs and Brill defined four boundary curves, all functions of λL alone, that carve the (λL, NFr) plane into segregated, transition, intermittent, and distributed flow:
L1 = 316 * lamL**0.302
L2 = 0.0009252 * lamL**(-2.4684)
L3 = 0.10 * lamL**(-1.4516)
L4 = 0.5 * lamL**(-6.738)With λL = 0.0612: L1 = 135.9, L2 = 0.91, L3 = 5.77, L4 is irrelevant here (it only matters once λL ≥ 0.4). Since 0.01 ≤ λL < 0.4 and L3 (5.77) < NFr (68.96) ≤ L1 (135.9), this segment is in intermittent flow — alternating slugs of liquid and pockets of gas moving up the string, the pattern most gas-lifted and moderate-GLR oil wells actually see.
Step 5: Calculate the horizontal liquid holdup
Next comes the holdup the pipe would show if it were perfectly horizontal, HL(0), using regime-specific constants a, b, c in HL(0) = a × λLb / NFrc (floored at λL if the formula ever returns less):
| Regime | a | b | c |
|---|---|---|---|
| Segregated | 0.98 | 0.4846 | 0.0868 |
| Intermittent | 0.845 | 0.5351 | 0.0173 |
| Distributed | 1.065 | 0.5824 | 0.0609 |
Using the intermittent coefficients: HL(0) = 0.845 × 0.06120.5351 / 68.960.0173 = 0.176.
Step 6: Apply the inclination correction
A vertical well holds up more liquid than a horizontal one at the same rates, because gravity works against the gas instead of across the pipe. The correction factor ψ adjusts the horizontal value for the actual pipe angle θ (measured from horizontal; θ = 90° for vertical uphill flow):
ψ = 1 + C × [sin(1.8θ) − ⅓sin³(1.8θ)], with C = (1−λL) × ln(d × λLe × NLvf × NFrg), floored at zero
For uphill intermittent flow (d=2.96, e=0.305, f=−0.4473, g=0.0978): C = (1−0.0612) × ln(2.96 × 0.06120.305 × 3.02−0.4473 × 68.960.0978) = 0.143. At θ=90°, sin(1.8×90°)=sin(162°)=0.309, so ψ = 1 + 0.143×(0.309−0.0098) = 1.043, giving the actual holdup HL(θ) = 0.176 × 1.043 = 0.184.
Step 7: Calculate the two-phase friction factor
Friction uses the no-slip mixture properties for the Reynolds number, then corrects the resulting friction factor for the real (slip) holdup. No-slip density and viscosity: ρns = ρL×λL + ρg×(1−λL) = 9.65 lbm/ft³; μns = μL×λL + μg×(1−λL) = 0.199 cp. Reynolds number Re = 1488×ρns×vm×d/μns = 312,385 — solidly turbulent, so the no-slip Darcy friction factor fns comes from a standard smooth-pipe correlation (Chen’s explicit approximation to Colebrook, assuming relative roughness ≈ 0.0006 for in-service tubing): fns = 0.01875.
The slip correction uses y = λL/HL(θ)² = 0.0612/0.184² = 1.814, then S = ln(y) / [−0.0523 + 3.182ln(y) − 0.8725ln(y)² + 0.01853ln(y)⁴] = 0.388 (y falls outside the 1.0–1.2 special case, so the general formula applies), giving the two-phase factor ftp = fns × eS = 0.01875 × 1.474 = 0.0276 — about 47% higher than the no-slip value, because the real gas-liquid interface drags more than a smooth mixture would.
Step 8: Combine elevation and friction into the pressure gradient
The total gradient has two terms — a hydrostatic term using the true (slip) mixture density, and a friction term using the no-slip density with the two-phase friction factor:
dP/dz = [ρm×sinθ + ftp×ρns×vm²/(2gcd)] / 144, where ρm = ρL×HL(θ) + ρg×(1−HL(θ)) = 15.15 lbm/ft³
Elevation term: 15.15 × 1.0 / 144 = 0.1052 psi/ft. Friction term: 0.0276 × 9.65 × 21.24² / (2 × 32.174 × 0.2034) / 144 = 0.0638 psi/ft. Total gradient: 0.1690 psi/ft, or 169.0 psi across the 1,000 ft segment — nearly 40% of it is friction, which a single borrowed gradient would have missed entirely.
Step 9: Automate it in Python
Every step above is five lines of code once you’ve done it by hand once. This function reproduces the worked example exactly, and is the building block for a full VLP curve — call it in a loop over trial rates, or march it down a string in short segments, re-evaluating fluid properties at each step’s average pressure (using petropt’s Standing correlation for Rs and Bo at each new pressure, for example):
import math
def beggs_brill_dp(d_in, theta_deg, L_ft, qL_bbl_d, qg_ft3_d,
rho_L, rho_g, sigma_L, mu_L, mu_g, eps_over_d=0.0006):
g, gc = 32.174, 32.174
d_ft = d_in / 12.0
A = math.pi / 4 * d_ft**2
v_sl = (qL_bbl_d * 5.615 / 86400.0) / A
v_sg = (qg_ft3_d / 86400.0) / A
v_m = v_sl + v_sg
lamL = v_sl / v_m
NFr = v_m**2 / (g * d_ft)
NLv = 1.938 * v_sl * (rho_L / sigma_L) ** 0.25
L1 = 316 * lamL**0.302
L2 = 0.0009252 * lamL**(-2.4684)
L3 = 0.10 * lamL**(-1.4516)
L4 = 0.5 * lamL**(-6.738)
if (lamL < 0.01 and NFr < L1) or (lamL >= 0.01 and NFr < L2):
regime = "segregated"
elif lamL >= 0.01 and L2 <= NFr <= L3:
regime = "intermittent" # treat transition as intermittent
elif (0.01 <= lamL < 0.4 and L3 < NFr <= L1) or (lamL >= 0.4 and L3 < NFr <= L4):
regime = "intermittent"
else:
regime = "distributed"
coeff = {"segregated": (0.98, 0.4846, 0.0868),
"intermittent": (0.845, 0.5351, 0.0173),
"distributed": (1.065, 0.5824, 0.0609)}
a, b, c = coeff[regime]
HL0 = max(a * lamL**b / NFr**c, lamL)
inc = {"segregated": (0.011, -3.768, 3.539, -1.614),
"intermittent": (2.96, 0.305, -0.4473, 0.0978)}
if regime in inc:
d_, e_, f_, g_ = inc[regime]
C = max((1 - lamL) * math.log(d_ * lamL**e_ * NLv**f_ * NFr**g_), 0.0)
else:
C = 0.0
th = math.radians(theta_deg)
psi = 1 + C * (math.sin(1.8 * th) - math.sin(1.8 * th) ** 3 / 3)
HLtheta = HL0 * psi
rho_ns = rho_L * lamL + rho_g * (1 - lamL)
rho_m = rho_L * HLtheta + rho_g * (1 - HLtheta)
mu_ns = mu_L * lamL + mu_g * (1 - lamL)
Re = 1488 * rho_ns * v_m * d_ft / mu_ns
Ac = (eps_over_d / 3.7065) - (5.0452 / Re) * math.log10(
eps_over_d**1.1098 / 2.8257 + (7.149 / Re) ** 0.8981)
f_ns = (1.0 / (-2 * math.log10(Ac))) ** 2
y = lamL / HLtheta**2
S = math.log(2.2*y - 1.2) if 1 < y < 1.2 else (
math.log(y) / (-0.0523 + 3.182*math.log(y) - 0.8725*math.log(y)**2 + 0.01853*math.log(y)**4))
f_tp = f_ns * math.exp(S)
dp_elev = rho_m / 144.0 * math.sin(th)
dp_fric = f_tp * rho_ns * v_m**2 / (2 * gc * d_ft) / 144.0
grad = dp_elev + dp_fric
return grad, grad * L_ft, regime
grad, dP, regime = beggs_brill_dp(2.441, 90, 1000, 650, 56000, 51.8, 6.9, 25.0, 3.0, 0.016)
print(f"Regime: {regime}, gradient: {grad:.4f} psi/ft, total dP: {dP:.1f} psi")
# Regime: intermittent, gradient: 0.1690 psi/ft, total dP: 169.0 psiWrap this in a loop that re-runs it segment by segment down the tubing, feeding each segment’s outlet pressure in as the next segment’s inlet, and you have a full VLP curve generator — exactly what you’d plug into the nodal analysis workflow we built for single-phase wells, now extended to the multiphase case. Automate the gas-lift side with the same inputs feeding the gas lift valve spacing calculation, and the two scripts together size and diagnose the whole string.
How to verify it worked
- Zero-gas check: set qg to a tiny value and λL should approach 1.0, HL should approach 1.0, and the gradient should converge on the single-phase hydrostatic-plus-friction result for liquid alone.
- Bounding check: the computed gradient (0.169 psi/ft here) should sit between the no-slip homogeneous gradient (using λL instead of HL(θ) for density, which is always lower) and the fully-segregated gradient — if it falls outside that range, a sign or exponent got mixed up.
- Regime sanity: intermittent and segregated flow should give noticeably higher holdup (and gradient) than distributed flow at the same rates — if distributed gives a higher holdup, check the regime classification logic.
Common pitfalls
- Forgetting to floor HL(0) at λL. At high NFr the raw correlation can return a holdup below the no-slip value, which is physically impossible — always clamp it.
- Using the wrong density in the wrong term. The hydrostatic term needs the true (slip) mixture density ρm; the friction term needs the no-slip density ρns. Swapping them is the single most common implementation bug.
- Treating the transition zone as a hard boundary. Near L2–L3 the published method interpolates between segregated and intermittent holdup; snapping straight to one or the other creates a visible kink in the VLP curve right where engineers are most likely to be reading it.
Hand-calculating one Beggs-Brill point is a useful exercise exactly once. After that, the function above is the artifact worth keeping — drop it into your nodal analysis script, your gas-lift design sheet, or a Python job that re-runs every well in the field overnight, and the VLP curve stops being a borrowed gradient and starts being an actual answer.
