Exerciseschevron_rightChapter 12chevron_right12.8

fitness_center

Exercise 12.8

Gas Lift Optimization - Operating Point, Incremental Oil & Economics

Level 3

Chapter 12: Nodal Analysis

descriptionProblem

A well produces 800 STB/day under natural flow. Using gas lift injection rates of 200, 500, 1000, 1500, and 2000 Mscf/d, calculate the incremental oil gain for each injection level. Plot oil gain vs. gas injected. At what injection rate does the incremental benefit become uneconomical if gas costs $3/Mscf and oil sells for$ 75/bbl?

---

We will work this on a tired OML well that just barely flows on its own. The verified Beggs-Brill VLP stack (tubing_performance and its PVT helpers) plus the chapter's tpr_with_gas_lift(Pwh, depth, tubing_id, rate, gas_injection_rate_mscfd=0) are embedded for you. Do not edit them or re-derive the physics. Gas lift lightens the column through lightening_factor = 1/(1 + 0.0003*gl) and a gravity_reduction = 0.433*0.85*depth*(1 - lightening_factor) subtracted from the natural-flow TPR.

The well constants (Exercise 12.8, tuned so natural flow lands near 800 STB/d): PE = 3600.0 psi, PB = 3200.0 psi, J = 0.85 STB/d/psi, PWH = 150.0 psi, DEPTH = 9000.0 ft, TUBING_ID = 1.995". The injection sweep is GL_VALUES = [200.0, 500.0, 1000.0, 1500.0, 2000.0] Mscf/d, with OIL_PRICE = 75.0 ( $/bbl) and$ /Mscf).

Your tasks:

Write gaslift_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id, gl_mscfd):

Build the inverse composite IPR ipr_pwf(q) exactly as the chapter does

(linear above the bubble point, Vogel solved for Pwf below it).

Form the residual r(q) = ipr_pwf(q) - tpr_with_gas_lift(Pwh, depth, tubing_id, max(q, 1), gl_mscfd).
Root-find with brentq over [50, qo_max*0.99] where qo_max = J*(Pe-Pb) + J*Pb/1.8.
Return 0.0 on ValueError (the well will not flow at that injection level).

Write gaslift_economics(Pe, Pb, J, Pwh, depth, tubing_id, gl_values, oil_price, gas_price):

Let q_natural = gaslift_operating_rate(..., 0.0).
For each gl in gl_values, compute q = gaslift_operating_rate(..., gl),

oil_gain = q - q_natural, and net_value = oil_gain*oil_price - gl*gas_price.

Return a list of dicts with keys gl, q, oil_gain, net_value

(one entry per injection level, same order as gl_values).

Driver: compute

q_natural = gaslift_operating_rate(PE, PB, J, PWH, DEPTH, TUBING_ID, 0.0)
econ = gaslift_economics(PE, PB, J, PWH, DEPTH, TUBING_ID, GL_VALUES, OIL_PRICE, GAS_PRICE)
gains = np.array([e["oil_gain"] for e in econ])
q_at_2000 = gaslift_operating_rate(PE, PB, J, PWH, DEPTH, TUBING_ID, 2000.0)
gain_at_2000 = float(q_at_2000 - q_natural)

The tests read these exact names: gaslift_operating_rate, gaslift_economics, tpr_with_gas_lift, q_natural, gains, q_at_2000, gain_at_2000.

> Think about it: the first slug of gas buys a lot of oil; the last slug > buys far less. Each extra Mscf injected lightens an already-lighter column, so > the incremental rate flattens out: diminishing returns. With gas at $3/Mscf > and oil at$ 75/bbl the net value still climbs across this whole sweep, but > watch how the slope of gains keeps falling. Where would the next Mscf finally > stop paying for itself?

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
from scipy.optimize import brentq

# ── Verified PVT helpers (Standing + Papay Z + Beggs-Robinson / Lee) (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)

# ── Verified Beggs-Brill two-phase gradient (do not edit) ────────────────
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."""
    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

# ── Verified tubing performance (full Beggs-Brill VLP march) (do not edit) ─
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."""
    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

# ── Verified gas-lifted TPR (do not edit) ────────────────────────────────
def tpr_with_gas_lift(Pwh, depth, tubing_id, rate, gas_injection_rate_mscfd=0):
    """Simplified TPR with gas lift effect.

Gas injection reduces effective fluid density in the tubing,
    lowering the required bottomhole pressure for a given rate.
    """
    base_Pwf = tubing_performance(Pwh, depth, tubing_id, rate)

# Gas lightening effect (simplified)
    # Each Mscf/d of injected gas reduces the effective liquid gradient
    if gas_injection_rate_mscfd > 0:
        lightening_factor = 1 / (1 + 0.0003 * gas_injection_rate_mscfd)
        gravity_reduction = 0.433 * 0.85 * depth * (1 - lightening_factor)
        return base_Pwf - gravity_reduction
    return base_Pwf

# ── OML well constants (Exercise 12.8) (do not edit) ─────────────────────
PE = 3600.0          # reservoir pressure, psi
PB = 3200.0          # bubble point, psi
J = 0.85             # productivity index, STB/d/psi
PWH = 150.0          # wellhead pressure, psi
DEPTH = 9000.0       # TVD, ft
TUBING_ID = 1.995    # 2-3/8" tubing ID, inches
GL_VALUES = [200.0, 500.0, 1000.0, 1500.0, 2000.0]   # gas injection sweep, Mscf/d
OIL_PRICE = 75.0     # $/bbl
GAS_PRICE = 3.0      # $/Mscf

def gaslift_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id, gl_mscfd):
    """Operating rate where the IPR meets the gas-lifted TPR (STB/day).

Build the inverse composite IPR, form the residual against the gas-lifted
    TPR, and root-find the intersection. Return 0.0 if the well will not flow.
    """
    # TODO: define ipr_pwf(q) -- the inverse composite IPR (chapter's closed form):
    #   qb = J * (Pe - Pb); qo_max = qb + J * Pb / 1.8
    #   if q <= qb:  return Pe - q / J                 # linear above the bubble point
    #   else:        solve 0.8*x^2 + 0.2*x + (frac-1) = 0 for x = Pwf/Pb (Vogel),
    #                frac = (q - qb) / (qo_max - qb); return max(x * Pb, 0)
    #
    # TODO: define residual(q):
    #   return ipr_pwf(q) - tpr_with_gas_lift(Pwh, depth, tubing_id, max(q, 1), gl_mscfd)
    #
    # TODO: qb = J * (Pe - Pb); qo_max = qb + J * Pb / 1.8
    # TODO: try:    return brentq(residual, 50, qo_max * 0.99)
    #       except ValueError:  return 0.0
    pass

def gaslift_economics(Pe, Pb, J, Pwh, depth, tubing_id, gl_values, oil_price, gas_price):
    """Incremental oil and net value for each gas injection level.

Return a list of dicts {gl, q, oil_gain, net_value}, one per injection
    level, with oil_gain measured against natural flow (gl = 0).
    """
    # TODO: q_natural = gaslift_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id, 0.0)
    # TODO: for each gl in gl_values:
    #         q = gaslift_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id, gl)
    #         oil_gain = q - q_natural
    #         net_value = oil_gain * oil_price - gl * gas_price
    #         append {"gl": gl, "q": q, "oil_gain": oil_gain, "net_value": net_value}
    # TODO: return the list
    pass

# TODO: q_natural = gaslift_operating_rate(PE, PB, J, PWH, DEPTH, TUBING_ID, 0.0)
# TODO: econ = gaslift_economics(PE, PB, J, PWH, DEPTH, TUBING_ID, GL_VALUES, OIL_PRICE, GAS_PRICE)
# TODO: gains = np.array([e["oil_gain"] for e in econ])
# TODO: q_at_2000 = gaslift_operating_rate(PE, PB, J, PWH, DEPTH, TUBING_ID, 2000.0)
# TODO: gain_at_2000 = float(q_at_2000 - q_natural)
q_natural = None
econ = None
gains = None
q_at_2000 = None
gain_at_2000 = None

print("Natural-flow rate (STB/day):", q_natural)
print("Rate @ 2000 Mscf/d lift:", q_at_2000)
print("Incremental oil @ 2000:", gain_at_2000)

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 PVT helpers (Standing + Papay Z + Beggs-Robinson / Lee) (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)


# ── Verified Beggs-Brill two-phase gradient (do not edit) ────────────────
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."""
    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


# ── Verified tubing performance (full Beggs-Brill VLP march) (do not edit) ─
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."""
    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


# ── Verified gas-lifted TPR (do not edit) ────────────────────────────────
def tpr_with_gas_lift(Pwh, depth, tubing_id, rate, gas_injection_rate_mscfd=0):
    """Simplified TPR with gas lift effect.

    Gas injection reduces effective fluid density in the tubing,
    lowering the required bottomhole pressure for a given rate.
    """
    base_Pwf = tubing_performance(Pwh, depth, tubing_id, rate)

    # Gas lightening effect (simplified)
    # Each Mscf/d of injected gas reduces the effective liquid gradient
    if gas_injection_rate_mscfd > 0:
        lightening_factor = 1 / (1 + 0.0003 * gas_injection_rate_mscfd)
        gravity_reduction = 0.433 * 0.85 * depth * (1 - lightening_factor)
        return base_Pwf - gravity_reduction
    return base_Pwf


# ── OML well constants (Exercise 12.8) (do not edit) ─────────────────────
PE = 3600.0          # reservoir pressure, psi
PB = 3200.0          # bubble point, psi
J = 0.85             # productivity index, STB/d/psi
PWH = 150.0          # wellhead pressure, psi
DEPTH = 9000.0       # TVD, ft
TUBING_ID = 1.995    # 2-3/8" tubing ID, inches
GL_VALUES = [200.0, 500.0, 1000.0, 1500.0, 2000.0]   # gas injection sweep, Mscf/d
OIL_PRICE = 75.0     # $/bbl
GAS_PRICE = 3.0      # $/Mscf


def gaslift_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id, gl_mscfd):
    """Operating rate where the IPR meets the gas-lifted TPR (STB/day).

    Builds the inverse composite IPR, forms the residual against the
    gas-lifted Beggs-Brill TPR, and root-finds the intersection. Returns
    0.0 if the well will not flow at this injection level.
    """
    def ipr_pwf(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 residual(q):
        return ipr_pwf(q) - tpr_with_gas_lift(Pwh, depth, tubing_id, max(q, 1), gl_mscfd)

    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


def gaslift_economics(Pe, Pb, J, Pwh, depth, tubing_id, gl_values, oil_price, gas_price):
    """Incremental oil and net value for each gas injection level.

    Returns a list of dicts {gl, q, oil_gain, net_value}, one per injection
    level, with oil_gain measured against natural flow (gl = 0).
    """
    q_natural = gaslift_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id, 0.0)
    out = []
    for gl in gl_values:
        q = gaslift_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id, gl)
        oil_gain = q - q_natural
        net_value = oil_gain * oil_price - gl * gas_price
        out.append({"gl": gl, "q": q, "oil_gain": oil_gain, "net_value": net_value})
    return out


q_natural = gaslift_operating_rate(PE, PB, J, PWH, DEPTH, TUBING_ID, 0.0)
econ = gaslift_economics(PE, PB, J, PWH, DEPTH, TUBING_ID, GL_VALUES, OIL_PRICE, GAS_PRICE)
gains = np.array([e["oil_gain"] for e in econ])
q_at_2000 = gaslift_operating_rate(PE, PB, J, PWH, DEPTH, TUBING_ID, 2000.0)
gain_at_2000 = float(q_at_2000 - q_natural)

print(f"Natural-flow rate:        {q_natural:,.1f} STB/day")
print(f"Rate @ 2000 Mscf/d lift:  {q_at_2000:,.1f} STB/day")
print(f"Incremental oil @ 2000:   {gain_at_2000:,.1f} STB/day")
print("\n  gl (Mscf/d)   q (STB/d)   oil gain   net value ($/d)")
for e in econ:
    print(f"  {e['gl']:>9,.0f}   {e['q']:>9,.1f}   {e['oil_gain']:>8,.1f}   {e['net_value']:>12,.0f}")

lockCopying code is a Full Access feature.

arrow_back

12.7 Wellhead Pressure Sensitivity - When the Well Dies

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

arrow_forward