Exerciseschevron_rightChapter 12chevron_right12.7
fitness_center

Exercise 12.7

Wellhead Pressure Sensitivity - When the Well Dies

Level 3
Chapter 12: Nodal Analysis
descriptionProblem

For the well in Exercise 12.6, calculate the operating rate at wellhead pressures ranging from 50 to 600 psi in 50-psi increments. Plot rate vs. Pwh. At what Pwh does the well stop flowing (die)?

---

This is the same OML well from Exercise 12.6 (Pe = 4000, Pb = 2800, J = 7, depth = 8500 ft, tubing ID = 2.441"), but now you raise the wellhead (backpressure) in 50-psi steps and watch the operating point walk along the IPR/TPR intersection, until it doesn't. Past a certain Pwh the tubing-performance curve no longer crosses the inflow curve at all: the well dies. The verified Beggs-Brill VLP stack and find_operating_rate() are embedded for you. Do not modify them or re-derive the physics. When the IPR and TPR no longer intersect, find_operating_rate catches the brentq "no bracket" ValueError and returns 0.0; that 0.0 is the well-dead signal you must detect.

The embedded constants (verbatim from Exercise 12.7) are: PE = 4000.0, PB = 2800.0, J = 7.0, DEPTH = 8500.0, TUBING_ID = 2.441, and the sweep grid PWH_VALUES = np.arange(50.0, 650.0, 50.0) (50 → 600 psi in 50-psi steps).

Your tasks:

  1. Write pwh_sweep(Pe, Pb, J, depth, tubing_id, pwh_values) that returns a

numpy array of operating rates, one entry per Pwh in pwh_values, each obtained from find_operating_rate(Pe, Pb, J, pwh, depth, tubing_id).

  1. Write flowing_limit(rates, pwh_values) that returns the lowest Pwh

at which the rate has dropped to 0.0 (the well dies): min(pwh for pwh, r in zip(pwh_values, rates) if r <= 0.0), or None if the well never dies anywhere in the range.

  1. Run the sweep on the embedded well:

``python rates = pwh_sweep(PE, PB, J, DEPTH, TUBING_ID, PWH_VALUES) kill_pwh = flowing_limit(rates, PWH_VALUES) ` then expose rate_at_50 = float(rates[0]) and rate_at_250 = float(rates[4])`.

The tests will read these names: find_operating_rate, pwh_sweep, flowing_limit, rates, kill_pwh, rate_at_50, rate_at_250.

> Think about it: the rate is nearly flat, even climbing slightly on the > low-backpressure limb, and then drops abruptly to zero once the TPR no > longer reaches the IPR. The well does not "fade away"; it falls off a cliff > between two grid points. Why does raising backpressure eventually kill flow > entirely, and why can no amount of more backpressure ever bring it back?

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 + nodal stack from ch.12 (do not edit) ──────
# PVT helpers (Standing + Papay Z + Beggs-Robinson / Lee viscosity)
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


def find_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id):
    """Find the operating rate for given system parameters."""

    def ipr_pw(q):
        qb = J * (Pe - Pb)
        qo_max = qb + J * Pb / 1.8
        if q <= qb:
            return Pe - q / J
        else:
            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)

    def tpr_pw(q):
        return tubing_performance(Pwh, depth, tubing_id, max(q, 1))

    def residual(q):
        return ipr_pw(q) - tpr_pw(q)

    qb = J * (Pe - Pb)
    qo_max = qb + J * Pb / 1.8

    try:
        return brentq(residual, 50, qo_max * 0.99)
    except ValueError:
        return 0.0


# ── Exercise 12.7 constants (verbatim, same well as Exercise 12.6) ───────
PE = 4000.0          # reservoir (external) pressure, psi
PB = 2800.0          # bubble-point pressure, psi
J = 7.0              # productivity index, STB/day/psi
DEPTH = 8500.0       # tubing depth (TVD), ft
TUBING_ID = 2.441    # 2-7/8" tubing ID, inches
PWH_VALUES = np.arange(50.0, 650.0, 50.0)   # 50 -> 600 psi in 50-psi steps


def pwh_sweep(Pe, Pb, J, depth, tubing_id, pwh_values):
    """Operating rate at each wellhead pressure in pwh_values.

    Returns a numpy array of operating rates (STB/day), one per Pwh, each from
    find_operating_rate. A 0.0 entry means the well is dead at that backpressure
    (the IPR and TPR no longer intersect).
    """
    return np.array([
        find_operating_rate(Pe, Pb, J, pwh, depth, tubing_id)
        for pwh in pwh_values
    ])


def flowing_limit(rates, pwh_values):
    """Lowest wellhead pressure at which the well has died (rate <= 0.0).

    Returns that Pwh, or None if the well never dies anywhere in the range.
    """
    dead = [pwh for pwh, r in zip(pwh_values, rates) if r <= 0.0]
    return min(dead) if dead else None


rates = pwh_sweep(PE, PB, J, DEPTH, TUBING_ID, PWH_VALUES)
kill_pwh = flowing_limit(rates, PWH_VALUES)
rate_at_50 = float(rates[0])
rate_at_250 = float(rates[4])

print("Pwh (psi):", PWH_VALUES.tolist())
print("rate (STB/d):", [round(float(r), 1) for r in rates])
print("rate at Pwh=50 psi :", round(rate_at_50, 1), "STB/d")
print("rate at Pwh=250 psi:", round(rate_at_250, 1), "STB/d")
print("well dies at Pwh =", kill_pwh, "psi")

lockCopying code is a Full Access feature.