Exercise 12.6
Operating Point Determination - IPR/TPR Intersection
Using Pe = 4,000 psi, Pb = 2,800 psi, J = 7 STB/day/psi, Pwh = 250 psi, depth = 8,500 ft, and tubing ID = 2.441", find the operating point. Report the rate, flowing BHP, and drawdown.
---
The operating point is the single rate at which the reservoir and the tubing agree: it is where the inflow performance relationship (IPR) and the tubing performance relationship (TPR) curves intersect. At that rate, the flowing bottomhole pressure the reservoir can deliver exactly equals the pressure the Beggs-Brill column requires to lift that rate to surface. Find that crossing and you have predicted what the well will actually do.
The full verified Beggs-Brill VLP stack (nine PVT/holdup helpers plus tubing_performance) and the chapter's find_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id) are embedded for you. Do not modify them or re-derive the physics. find_operating_rate builds the composite-IPR pressure, wraps tubing_performance as the TPR, forms the residual IPR(q) - TPR(q), and roots it with scipy.optimize.brentq.
Your task: write operating_point(Pe, Pb, J, Pwh, depth, tubing_id) that returns the tuple (q_op, Pwf_op, drawdown):
q_op = find_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id), the
intersection rate, STB/day.
Pwf_op= the composite-IPR flowing bottomhole pressure atq_op,
recomputed from the same inverse-IPR closed form find_operating_rate uses internally:
qb = J * (Pe - Pb)(the rate at the bubble point)qo_max = qb + J * Pb / 1.8(Vogel AOF on the saturated branch)- if
q_op <= qb: straight-line branch,Pwf = Pe - q_op / J - else (saturated branch): with
frac = (q_op - qb) / (qo_max - qb), solve the
Vogel quadratic 0.8 x*2 + 0.2 x + (frac - 1) = 0 via x = (-0.2 + sqrt(0.04 - 3.2 (frac - 1))) / 1.6, then Pwf = max(x * Pb, 0.0).
drawdown = Pe - Pwf_op, the pressure draw the reservoir is seeing, psi.
Embedded constants (verbatim from Exercise 12.6): PE = 4000.0, PB = 2800.0, J = 7.0, PWH = 250.0, DEPTH = 8500.0, TUBING_ID = 2.441.
Driver: call it once on the Exercise 12.6 well and expose the named outputs:
q_op, Pwf_op, drawdown = operating_point(PE, PB, J, PWH, DEPTH, TUBING_ID)
pwf_op = float(Pwf_op) # ~3411.4 psi
drawdown_psi = float(drawdown) # ~588.6 psi
# q_op ~ 4119.9 STB/dayThe names the tests read are: find_operating_rate, tubing_performance, operating_point, q_op, pwf_op, drawdown_psi.
> Think about it: at the operating point the IPR and TPR pressures must > be equal; that is the very definition of the crossing, and no hardcoded > number can force two independent curves to meet at a chosen rate. A higher > productivity index J should push q_op up; a larger tubing ID (less > friction at these rates) should let the well flow at least as fast. Why does > wider tubing not always help? What happens at very low rates on the > liquid-loading limb?
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 stack (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.6 well constants (do not edit) ───────────────────────────
PE = 4000.0 # external (reservoir) pressure, psi
PB = 2800.0 # bubble-point pressure, psi
J = 7.0 # productivity index, STB/day/psi
PWH = 250.0 # wellhead pressure, psi
DEPTH = 8500.0 # tubing depth (TVD), ft
TUBING_ID = 2.441 # 2-7/8" tubing ID, inches
def operating_point(Pe, Pb, J, Pwh, depth, tubing_id):
"""Operating point = IPR/TPR intersection for one well.
Returns (q_op, Pwf_op, drawdown):
q_op = intersection rate from find_operating_rate (STB/day)
Pwf_op = composite-IPR flowing BHP at q_op (psi), from the same
inverse-IPR closed form find_operating_rate uses internally
drawdown = Pe - Pwf_op (psi)
"""
q_op = find_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id)
qb = J * (Pe - Pb)
qo_max = qb + J * Pb / 1.8
if q_op <= qb:
Pwf_op = Pe - q_op / J
else:
frac = (q_op - qb) / (qo_max - qb)
x = (-0.2 + np.sqrt(0.04 - 3.2 * (frac - 1))) / 1.6
Pwf_op = max(x * Pb, 0.0)
drawdown = Pe - Pwf_op
return q_op, Pwf_op, drawdown
q_op, Pwf_op, drawdown = operating_point(PE, PB, J, PWH, DEPTH, TUBING_ID)
pwf_op = float(Pwf_op)
drawdown_psi = float(drawdown)
print(f"Operating rate q_op : {q_op:,.1f} STB/day")
print(f"Flowing BHP pwf_op : {pwf_op:,.1f} psi")
print(f"Drawdown drawdown : {drawdown_psi:,.1f} psi")
lockCopying code is a Full Access feature.