Exerciseschevron_rightChapter 12chevron_right12.9
fitness_center

Exercise 12.9

Declining Reservoir Pressure - Stable-Limb Operating Point & Artificial-Lift Trigger

Level 3
Chapter 12: Nodal Analysis
descriptionProblem

A well starts with Pe = 5,000 psi. Over five years, reservoir pressure declines to 3,000 psi in 500-psi increments. For each pressure, calculate the operating point and plot how the production rate changes over time. At what reservoir pressure does the well require artificial lift (natural flow ceases)?

---

We'll work this on an OML-58 oil well. The full Beggs-Brill VLP stack (tubing_performance and its PVT helpers) is embedded for you. Do not modify it or re-derive the physics. You add the nodal layer: the operating point is where the inflow (IPR) and outflow (TPR) curves cross, and as the reservoir depletes that crossing slides down a J-shaped TPR until, at some pressure, it disappears entirely and the well dies (natural flow ceases, artificial lift is required).

A real TPR is J-shaped, so the IPR can cross it in two places: an unstable low-rate point on the liquid-loading limb and a stable high-rate point on the right-hand limb. The well actually flows at the stable (highest-rate) crossing.

Your tasks:

  1. Write stable_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id) returning the

stable-limb operating rate, the highest-rate IPR/TPR intersection:

  • Re-derive the inverse composite IPR inline as ipr_pwf(q) (linear above the

bubble point, Vogel below it): qb = J*(Pe - Pb), qo_max = qb + J*Pb/1.8. For q <= qb return Pe - q/J; otherwise solve Vogel for Pwf via the quadratic 0.8*x**2 + 0.2*x + (frac - 1) = 0 with frac = (q - qb)/(qo_max - qb), x = (-0.2 + sqrt(0.04 - 3.2*(frac - 1)))/1.6, returning max(x*Pb, 0.0) (return 0.0 if the discriminant is negative).

  • Define the residual r(q) = ipr_pwf(q) - tubing_performance(Pwh, depth, tubing_id, max(q, 1)).
  • Scan q over np.linspace(50, qo_max*0.99, 400), find all sign changes

of r(q), run brentq on each bracket, and return the maximum root.

  • If there are no sign changes, return 0.0. The well has died and needs

artificial lift.

  1. Write depletion_profile(Pe_values, Pb, J, Pwh, depth, tubing_id) returning a

numpy array of stable operating rates, one per Pe.

Embedded constants (Exercise 12.9):

PB = 2800.0           # bubble-point pressure, psi
J = 5.0               # productivity index, STB/d/psi
PWH = 150.0           # wellhead pressure, psi
DEPTH = 8000.0        # TVD, ft
TUBING_ID = 2.992     # 3-1/2" tubing ID, inches
PE_VALUES = np.arange(5000.0, 2999.0, -500.0)   # 5000, 4500, 4000, 3500, 3000

Driver (run once, exact output variable names the tests read):

rates = depletion_profile(PE_VALUES, PB, J, PWH, DEPTH, TUBING_ID)
rate_at_5000 = float(rates[0])    # ~8137.3 STB/d
rate_at_3000 = float(rates[-1])   # ~3118.8 STB/d

The starter sets rates, rate_at_5000, and rate_at_3000 to None.

> Think about it: the production rate falls monotonically as Pe drops, but > the well does not coast smoothly to zero. Below about 2,300 psi the IPR no > longer reaches the stable limb of the TPR, so stable_operating_rate returns > 0.0 and the well must go on artificial lift. Why does the highest root > (not the lowest) tell you the rate a well actually flows at, and why does a > J-shaped TPR mean a well can die long before its pressure hits zero?

lightbulbHints (0/3)

Stuck? Reveal hints one at a time — they progress from nudge to near-solution.

codeYour solution
main.py
visibilityReveal reference solutionexpand_more

Try solving it yourself first — the hints walk you through it. The solution below is one valid approach; yours may differ and still be correct.

import numpy as np
from scipy.optimize import brentq


# ── Verified Beggs-Brill VLP stack from ch12 (do not edit) ───────────────
def _api_to_sg(API):
    return 141.5 / (131.5 + API)

def _standing_pb(rs, gas_sg, T_F, API):
    return 18.2 * ((rs / gas_sg) ** 0.83 * 10 ** (0.00091 * T_F - 0.0125 * API) - 1.4)

def _standing_rs(P, gas_sg, T_F, API):
    return gas_sg * ((P / 18.2 + 1.4) * 10 ** (0.0125 * API - 0.00091 * T_F)) ** (1 / 0.83)

def _standing_bo(rs, gas_sg, oil_sg, T_F):
    F = rs * np.sqrt(gas_sg / oil_sg) + 1.25 * T_F
    return 0.972 + 1.47e-4 * F ** 1.175

def _zfac(P, T_F, gas_sg):                      # Papay approximation
    Ppc = 756.8 - 131 * gas_sg - 3.6 * gas_sg ** 2
    Tpc = 169.2 + 349.5 * gas_sg - 74 * gas_sg ** 2
    Ppr, Tpr = P / Ppc, (T_F + 460) / Tpc
    return 1 - 3.52 * Ppr / 10 ** (0.9813 * Tpr) + 0.274 * Ppr ** 2 / 10 ** (0.8157 * Tpr)

def _mu_oil(rs, T_F, API):
    Y = 10 ** (3.0324 - 0.02023 * API); mod = 10 ** (Y * T_F ** -1.163) - 1
    A = 10.715 * (rs + 100) ** -0.515; B = 5.44 * (rs + 150) ** -0.338
    return A * mod ** B

def _mu_gas(P, T_F, gas_sg, Z):
    Mg = 28.97 * gas_sg; Tr = T_F + 460; rho = 2.7 * gas_sg * P / (Z * Tr) / 62.428
    K = (9.4 + 0.02 * Mg) * Tr ** 1.5 / (209 + 19 * Mg + Tr)
    X = 3.5 + 986 / Tr + 0.01 * Mg; Yv = 2.4 - 0.2 * X
    return 1e-4 * K * np.exp(X * rho ** Yv)


def beggs_brill_gradient(P, T_F, q_liq_stbd, wc, gor, API, gas_sg, d_in, theta_deg=90.0):
    """Beggs-Brill (1973) two-phase pressure gradient (psi/ft) at one point.

    Determines the flow pattern (segregated / transition / intermittent /
    distributed), the liquid holdup with inclination correction, and the
    two-phase friction factor -- the in-situ gas fraction (and thus the
    column weight) follows from the PVT at the local pressure.
    """
    G = 32.174
    oil_sg = _api_to_sg(API); water_sg = 1.07
    Pb = _standing_pb(gor, gas_sg, T_F, API)
    rs = gor if P >= Pb else max(_standing_rs(P, gas_sg, T_F, API), 0.0)
    rs = min(rs, gor)
    Z = _zfac(P, T_F, gas_sg); Bo = _standing_bo(rs, gas_sg, oil_sg, T_F)
    Bg = 0.0283 * Z * (T_F + 460) / P                       # rcf/scf
    q_oil = q_liq_stbd * (1 - wc); q_water = q_liq_stbd * wc
    qL = (q_oil * Bo + q_water * 1.0) * 5.615 / 86400.0      # in-situ ft3/s
    qG = max(q_oil * (gor - rs), 0.0) * Bg / 86400.0
    A = np.pi * (d_in / 24.0) ** 2; d = d_in / 12.0
    vsl, vsg = qL / A, qG / A; vm = vsl + vsg
    if vm <= 0:
        return 0.0
    lam = min(max(vsl / vm, 1e-9), 1.0)
    rhoO = (oil_sg * 62.4 + rs * gas_sg * 0.0764 / 5.615) / Bo
    rhoL = (q_oil * rhoO + q_water * water_sg * 62.4) / max(q_oil + q_water, 1e-9)
    rhog = 2.7 * gas_sg * P / (Z * (T_F + 460))
    NFR = vm ** 2 / (G * d)
    L1 = 316 * lam ** 0.302; L2 = 0.0009252 * lam ** -2.4684
    L3 = 0.10 * lam ** -1.4516; L4 = 0.5 * lam ** -6.738
    if (lam < 0.01 and NFR < L1) or (lam >= 0.01 and NFR < L2):
        pat = "seg"
    elif lam >= 0.01 and L2 <= NFR <= L3:
        pat = "trans"
    elif (0.01 <= lam < 0.4 and L3 < NFR <= L1) or (lam >= 0.4 and L3 < NFR <= L4):
        pat = "int"
    else:
        pat = "dist"
    def hl0(p):
        a, b, c = {"seg": (0.98, 0.4846, 0.0868), "int": (0.845, 0.5351, 0.0173),
                   "dist": (1.065, 0.5824, 0.0609)}[p]
        return max(a * lam ** b / NFR ** c, lam)
    sigma = 30.0
    NLv = vsl * (rhoL / (G * sigma)) ** 0.25 if rhoL > 0 else 0.0
    def Cc(p):
        if p == "dist":
            return 0.0
        e, f, gg, h = {"seg": (0.011, -3.768, 3.539, -1.614),
                       "int": (2.96, 0.305, -0.4473, 0.0978)}[p]
        return max((1 - lam) * np.log(e * lam ** f * max(NLv, 1e-9) ** gg * NFR ** h), 0.0)
    def psi(p):
        s = np.sin(1.8 * np.radians(theta_deg)); return 1 + Cc(p) * (s - 0.333 * s ** 3)
    if pat == "trans":
        w = (L3 - NFR) / (L3 - L2)
        HL = w * hl0("seg") * psi("seg") + (1 - w) * hl0("int") * psi("int")
    else:
        HL = min(hl0(pat) * psi(pat), 1.0)
    rhom = rhoL * HL + rhog * (1 - HL)
    dpdz_elev = rhom * np.sin(np.radians(theta_deg)) / 144.0
    rhons = rhoL * lam + rhog * (1 - lam)
    muns = _mu_oil(rs, T_F, API) * lam + _mu_gas(P, T_F, gas_sg, Z) * (1 - lam)
    Re = rhons * vm * d / (muns * 6.7197e-4)
    fns = 0.0056 + 0.5 / Re ** 0.32 if Re > 0 else 0.02
    y = lam / HL ** 2
    if 1.0 < y < 1.2:
        S = np.log(2.2 * y - 1.2)
    else:
        ly = np.log(max(y, 1e-9))
        S = ly / (-0.0523 + 3.182 * ly - 0.8725 * ly ** 2 + 0.01853 * ly ** 4)
    ftp = fns * np.exp(S)
    dpdz_fric = ftp * rhons * vm ** 2 / (2 * G * d) / 144.0
    return dpdz_elev + dpdz_fric


def tubing_performance(Pwh, depth_ft, tubing_id_in, oil_rate_stbd,
                       oil_sg=0.85, gas_sg=0.65, gor_scf_stb=500,
                       wc=0.1, oil_visc_cp=2.0,
                       temp_surface_F=100.0, temp_bottom_F=210.0, n_seg=60):
    """Flowing bottomhole pressure (psi) from a full Beggs-Brill VLP.

    Marches the two-phase pressure gradient down the tubing, recomputing the
    PVT (solution gas, FVF, densities) at the local pressure of each segment --
    so dissolved gas coming out of solution actually lightens the column. The
    total liquid rate is oil_rate / (1 - wc).
    """
    API = 141.5 / oil_sg - 131.5
    q_liq = oil_rate_stbd / (1.0 - wc)
    P = Pwh
    dz = depth_ft / n_seg
    for i in range(n_seg):
        T = temp_surface_F + (temp_bottom_F - temp_surface_F) * (i + 0.5) / n_seg
        P += beggs_brill_gradient(P, T, q_liq, wc, gor_scf_stb, API, gas_sg,
                                  tubing_id_in) * dz
    return P


# ── Exercise 12.9 constants (do not edit) ────────────────────────────────
PB = 2800.0           # bubble-point pressure, psi
J = 5.0               # productivity index, STB/d/psi
PWH = 150.0           # wellhead pressure, psi
DEPTH = 8000.0        # TVD, ft
TUBING_ID = 2.992     # 3-1/2" tubing ID, inches
PE_VALUES = np.arange(5000.0, 2999.0, -500.0)   # 5000, 4500, 4000, 3500, 3000


def stable_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id):
    """Stable-limb operating rate (STB/d): the HIGHEST-rate IPR/TPR crossing.

    A J-shaped TPR can be crossed twice by the IPR -- an unstable low-rate
    point on the liquid-loading limb and a stable high-rate point on the
    right-hand limb. The well actually flows at the stable (highest) crossing.
    Returns 0.0 when the IPR no longer reaches the TPR (the well has died and
    requires artificial lift).
    """
    def ipr_pwf(q):
        qb = J * (Pe - Pb)
        qo_max = qb + J * Pb / 1.8
        if q <= qb:
            return Pe - q / J
        frac = (q - qb) / (qo_max - qb)
        a, b, c = 0.8, 0.2, (frac - 1)
        disc = b ** 2 - 4 * a * c
        if disc < 0:
            return 0.0
        x = (-b + np.sqrt(disc)) / (2 * a)
        return max(x * Pb, 0.0)

    def residual(q):
        return ipr_pwf(q) - tubing_performance(Pwh, depth, tubing_id, max(q, 1))

    qb = J * (Pe - Pb)
    qo_max = qb + J * Pb / 1.8
    qs = np.linspace(50.0, qo_max * 0.99, 400)
    rs = np.array([residual(q) for q in qs])
    roots = []
    for i in range(len(qs) - 1):
        if rs[i] == 0.0:
            roots.append(float(qs[i]))
        if rs[i] * rs[i + 1] < 0:
            try:
                roots.append(float(brentq(residual, qs[i], qs[i + 1])))
            except ValueError:
                pass
    if not roots:
        return 0.0
    return max(roots)


def depletion_profile(Pe_values, Pb, J, Pwh, depth, tubing_id):
    """Stable operating rate at each reservoir pressure -> numpy array."""
    return np.array([stable_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id)
                     for Pe in Pe_values])


rates = depletion_profile(PE_VALUES, PB, J, PWH, DEPTH, TUBING_ID)
rate_at_5000 = float(rates[0])
rate_at_3000 = float(rates[-1])

print("Pe (psi)   stable rate (STB/d)")
for pe, q in zip(PE_VALUES, rates):
    print(f"  {pe:7.0f}   {q:10.1f}")
print("rate_at_5000:", rate_at_5000)
print("rate_at_3000:", rate_at_3000)

lockCopying code is a Full Access feature.

arrow_back

Previous

12.8 Gas Lift Optimization - Operating Point, Incremental Oil & Economics

Next

13.10 Full Field Optimization Report - Constrained Oil & Profit Allocation

arrow_forward