Nodal Analysis in Python: Automate IPR-VLP

By Saad Iqbal

Every time production wants to know what a well will actually flow at — after a workover, a choke change, or before picking new tubing — somebody redraws the same two curves: the reservoir’s inflow performance and the tubing’s intake requirement. Whoever’s doing it either fires up commercial software they may not have a license for on that machine, or drags points onto a spreadsheet chart and eyeballs where the lines cross. Nodal analysis in Python turns that fifteen-minute redraw into a function call, and once it’s a function you can rerun it for every tubing size or choke setting in seconds instead of redrawing anything.

This tutorial builds a working IPR-VLP nodal analysis solver from scratch: the governing Vogel inflow equation, a simplified tubing intake curve, the numerical root-finder that locates the operating point where they cross, and a sensitivity sweep across tubing sizes you can adapt to your own well.

Prerequisites

• Python 3.8+ with the SciPy and NumPy packages — get Python from python.org, then pip install scipy numpy.
• An average reservoir pressure and one well test point (flowing bottomhole pressure and the oil rate measured at it).
• A rough tubing/wellhead setup: perforation TVD, flowing wellhead pressure, and an idea of your average flowing gradient (from a hydraulics correlation or offset well data).

Step 1: Build the Reservoir Inflow (IPR) Curve

The Vogel correlation is the standard way to model inflow from a solution-gas-drive reservoir, where the average reservoir pressure is at or below the bubble point and the flow is two-phase throughout:

qo / qo_max = 1 - 0.2 (Pwf / Pbar) - 0.8 (Pwf / Pbar)^2

where qo is oil rate, qo_max is the absolute open flow potential (the theoretical rate at zero flowing pressure), Pwf is flowing bottomhole pressure, and Pbar is average reservoir pressure. Given one well test point you can solve for qo_max, then generate the full curve:

Pbar = 2900.0      # psia, average reservoir pressure
test_pwf = 2200.0  # psia, the flowing pressure during the test
test_q = 380.0     # stb/d, the rate measured at that pressure

x = test_pwf / Pbar
qmax = test_q / (1 - 0.2*x - 0.8*x**2)
print(f"qmax = {qmax:.0f} stb/d")   # qmax = 980 stb/d

def ipr_q(pwf, Pbar=Pbar, qmax=qmax):
    x = pwf / Pbar
    return qmax * (1 - 0.2*x - 0.8*x**2)

With a single well test point of 2,200 psia and 380 stb/d against a 2,900 psia reservoir, that solves to a 980 stb/d absolute open flow — the theoretical ceiling this well could ever produce at zero backpressure.

Step 2: Build the Tubing Intake (VLP) Curve

The vertical lift performance curve is the flip side: how much bottomhole pressure the tubing string requires to lift a given rate to surface against a fixed wellhead backpressure. A full VLP model uses a multiphase correlation like Hagedorn-Brown or Beggs-Brill — that’s what the commercial tools in our nodal analysis software comparison are built around. For this tutorial, use a simplified average-gradient model that captures the same shape and is good enough to teach the workflow:

Pwh = 150.0     # psia, flowing wellhead pressure
TVD = 6500.0    # ft, perforation depth
grad = 0.28     # psi/ft, average flowing gradient (from your hydraulics model)
C = 60 / (500 ** 1.8)  # friction coefficient, calibrated so 500 stb/d costs 60 psi

def vlp_pwf(q, Pwh=Pwh, TVD=TVD, grad=grad, C=C):
    return Pwh + grad * TVD + C * q ** 1.8

The hydrostatic term (Pwh + grad * TVD) barely changes with rate; the friction term grows with rate to the 1.8 power, which is why VLP curves slope gently upward — more rate means more friction loss, which means more pressure is needed at bottomhole to deliver it.

Step 3: Solve for the Operating Point

The well will settle at the one rate where the reservoir can supply exactly what the tubing demands — where the IPR and VLP curves cross. Rather than reading it off a chart, use scipy.optimize.brentq to root-find it directly:

from scipy.optimize import brentq

def ipr_pwf(q):
    """Invert the IPR curve: given a rate, find the Pwf that produces it."""
    return brentq(lambda p: ipr_q(p) - q, 0, Pbar)

def g(q):
    return vlp_pwf(q) - ipr_pwf(q)

q_op = brentq(g, 1, qmax - 1)
pwf_op = ipr_pwf(q_op)
print(f"Operating point: {q_op:.0f} stb/d at {pwf_op:.0f} psia")
# Operating point: 462 stb/d at 2022 psia

That’s the whole nodal analysis: this well, with this tubing and this reservoir, settles at roughly 462 stb/d — well under the 980 stb/d absolute open flow, because the tubing string is the limiting factor at this rate, not the reservoir.

Step 4: Rerun It Across Tubing Sizes

This is where automation actually pays for itself. Instead of redrawing the chart for every tubing size a completions engineer wants evaluated, loop the same solver over a few friction coefficients — smaller tubing means more friction, larger tubing means less, up to a point where gas slippage starts working against you in real wells:

scenarios = {
    '2 3/8" tubing': C * 2.6,
    '2 7/8" tubing': C * 1.0,
    '3 1/2" tubing': C * 0.42,
}

for name, C_i in scenarios.items():
    q_i = brentq(lambda q: vlp_pwf(q, C=C_i) - ipr_pwf(q), 1, qmax - 1)
    pwf_i = ipr_pwf(q_i)
    print(f"{name:16s} -> {q_i:6.1f} stb/d @ {pwf_i:7.1f} psia")

# 2 3/8" tubing    ->  431.7 stb/d @  2089.7 psia
# 2 7/8" tubing    ->  462.2 stb/d @  2022.1 psia
# 3 1/2" tubing    ->  475.0 stb/d @  1993.0 psia

Stepping up from 2 3/8″ to 2 7/8″ tubing buys roughly 30 stb/d by cutting friction losses; jumping again to 3 1/2″ only adds another 13 stb/d — diminishing returns you’d otherwise have to redraw three separate charts to see. If artificial lift is on the table instead of a tubing change, the same operating-point logic is what drives the artificial lift optimization tools covered in our roundup, and it’s the same curve that shifts underneath you as an ESP degrades — see how AI predicts ESP failures before they show up as a falling operating point.

Expected Result and How to Verify It

Running the Step 3 script against the numbers in this tutorial should return an operating point of roughly 462 stb/d at 2,022 psia. As a sanity check, confirm two things: the operating rate should always be comfortably below qmax (here, 980 stb/d) — if your solver returns a rate above that, something’s wrong with your IPR inputs — and the operating pressure from both curves should match to within rounding, since that agreement is what “operating point” means. If brentq raises a sign-change error, your two curves don’t cross in the bracket you gave it, which usually means your VLP curve is priced too high (too much friction or too deep a TVD) for the well to flow at all.

Common Pitfalls

• The VLP model here is simplified. A constant average gradient ignores how gas breaks out of solution as pressure drops up the tubing, which real multiphase correlations account for. Treat this tutorial’s numbers as directionally correct for learning the workflow, not as a substitute for a licensed multiphase VLP tool on an actual well.
• Wrong bubble point assumption. The pure Vogel form used here assumes the reservoir is already at or below bubble point. For an undersaturated reservoir with Pbar above the bubble point, you need the composite IPR form instead, or your qmax will be wrong.
• Stale well test data. qmax is only as good as the well test point you calibrated it from. Reservoir pressure declines over the life of the well, so re-anchor the IPR curve against a recent test before trusting the operating point for a completions or workover decision.

Once this solver is in your toolkit, the same root-finding pattern extends to choke performance curves, gas lift response curves, or swapping in a real multiphase correlation from a library like petropt instead of the simplified VLP model here. Read the related tutorial below to keep building out your production-engineering automation.