Exercise 12.7
Wellhead Pressure Sensitivity - When the Well Dies
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:
- 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).
- Write
flowing_limit(rates, pwh_values)that returns the lowestPwh
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.
- 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?
Stuck? Reveal hints one at a time — they progress from nudge to near-solution.
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.