Tubing Size Comparison - TPR Curves Across Tubing IDs

Level 2

Chapter 12: Nodal Analysis

descriptionProblem

For a well producing from 10,000 ft TVD with Pwh = 150 psi, calculate and plot TPR curves for tubing IDs of 1.995", 2.441", 2.992", and 3.958". At what rate does friction begin to dominate over gravity for each size?

---

A bigger tubing isn't always better. Gravity dominates the column at low rates (a wider pipe holds a heavier liquid column), but friction (which scales like v^2 / d^5) explodes at high rates, where the narrow string chokes the well. The crossover is exactly what an engineer sizing tubing on an OML-58 well needs to see, and the vertical lift performance (VLP) family of TPR curves makes it visible.

The full Beggs-Brill (1973) PVT/VLP stack is embedded for you. _api_to_sg, _standing_pb, _standing_rs, _standing_bo, _zfac, _mu_oil, _mu_gas, beggs_brill_gradient, and tubing_performance. Do not modify them or re-derive the physics. The marcher recomputes the in-situ PVT at the local pressure of every segment, so dissolved gas coming out of solution actually lightens the column.

The constants (verbatim from Exercise 12.5) are: DEPTH = 10000.0 ft TVD, PWH = 150.0 psi wellhead pressure, and four candidate tubing IDs TUBING_IDS = [1.995, 2.441, 2.992, 3.958] inches.

Your tasks:

Write tpr_for_size(tubing_id, rates, Pwh, depth):

For each q in rates, call tubing_performance(Pwh, depth, tubing_id, q)

(default PVT args) to get the required flowing bottomhole pressure.

Return a np.array of those flowing BHPs, in the same order as rates.
It must be a faithful wrapper: one tubing_performance call per rate, no

lookup tables, no fudge factors.

Build the TPR curves over RATES = np.array([500.0, 2000.0, 4000.0]):

``python tpr_curves = {tid: tpr_for_size(tid, RATES, PWH, DEPTH) for tid in TUBING_IDS} ``

Pull out these scalar readings the tests look for:

pwf_small_2000 = float(tpr_curves[1.995][1]), smallest tubing at 2000 STB/d
pwf_large_2000 = float(tpr_curves[3.958][1]), largest tubing at 2000 STB/d
pwf_small_4000 = float(tpr_curves[1.995][2]), smallest tubing at 4000 STB/d

> Think about it: at 2000 STB/d the small 1.995" string already needs > more flowing BHP than the wide 3.958" string; friction is winning. Push > the rate to 4000 STB/d and the gap blows open. Why does the largest tubing > need the most BHP at the very lowest rate (500 STB/d). What is holding that > well back when the pipe is wide and the velocity is tiny?

lightbulbHints (0/3)

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

codeYour solution

main.py

import numpy as np

# ── Verified Beggs-Brill VLP (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

# ── Exercise 12.5 constants (do not edit) ────────────────────────────────
DEPTH = 10000.0                              # ft TVD
PWH = 150.0                                  # psi, wellhead pressure
TUBING_IDS = [1.995, 2.441, 2.992, 3.958]   # candidate tubing IDs, inches
RATES = np.array([500.0, 2000.0, 4000.0])   # oil rates to evaluate, STB/d

def tpr_for_size(tubing_id, rates, Pwh, depth):
    """Required flowing BHP (psi) for each rate through one tubing size.

For each q in `rates`, run a full Beggs-Brill VLP via
    tubing_performance(Pwh, depth, tubing_id, q) and collect the flowing
    bottomhole pressures into a np.array (same order as `rates`).
    """
    # TODO: return np.array([
    # TODO:     tubing_performance(Pwh, depth, tubing_id, q) for q in rates
    # TODO: ])
    pass

# TODO: tpr_curves = {tid: tpr_for_size(tid, RATES, PWH, DEPTH) for tid in TUBING_IDS}
tpr_curves = None

# TODO: pwf_small_2000 = float(tpr_curves[1.995][1])   # smallest tubing @ 2000 STB/d
# TODO: pwf_large_2000 = float(tpr_curves[3.958][1])   # largest  tubing @ 2000 STB/d
# TODO: pwf_small_4000 = float(tpr_curves[1.995][2])   # smallest tubing @ 4000 STB/d
pwf_small_2000 = None
pwf_large_2000 = None
pwf_small_4000 = None

print("pwf small @2000:", pwf_small_2000, " large @2000:", pwf_large_2000)
print("pwf small @4000:", pwf_small_4000)

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


# ── Verified Beggs-Brill VLP (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


# ── Exercise 12.5 constants (do not edit) ────────────────────────────────
DEPTH = 10000.0                              # ft TVD
PWH = 150.0                                  # psi, wellhead pressure
TUBING_IDS = [1.995, 2.441, 2.992, 3.958]   # candidate tubing IDs, inches
RATES = np.array([500.0, 2000.0, 4000.0])   # oil rates to evaluate, STB/d


def tpr_for_size(tubing_id, rates, Pwh, depth):
    """Required flowing BHP (psi) for each rate through one tubing size.

    For each q in `rates`, run a full Beggs-Brill VLP via
    tubing_performance(Pwh, depth, tubing_id, q) and collect the flowing
    bottomhole pressures into a np.array (same order as `rates`).
    """
    return np.array([
        tubing_performance(Pwh, depth, tubing_id, q) for q in rates
    ])


tpr_curves = {tid: tpr_for_size(tid, RATES, PWH, DEPTH) for tid in TUBING_IDS}

pwf_small_2000 = float(tpr_curves[1.995][1])   # smallest tubing @ 2000 STB/d
pwf_large_2000 = float(tpr_curves[3.958][1])   # largest  tubing @ 2000 STB/d
pwf_small_4000 = float(tpr_curves[1.995][2])   # smallest tubing @ 4000 STB/d

print("pwf small @2000:", pwf_small_2000, " large @2000:", pwf_large_2000)
print("pwf small @4000:", pwf_small_4000)

lockCopying code is a Full Access feature.

arrow_back

12.4 Fetkovich from Multi-Rate Test - Log-Log Deliverability Fit

12.6 Operating Point Determination - IPR/TPR Intersection

arrow_forward