Full Nodal Analysis Report - Operating Point, Sensitivity & Recommendation

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 composite IPR (linear above Pb, Vogel below) (do not edit) ────
def composite_ipr(Pe, Pb, J, Pwf_array):
    """Composite IPR: linear above bubble point, Vogel below."""
    Pwf = np.asarray(Pwf_array)
    q = np.zeros_like(Pwf, dtype=float)
    qb = J * (Pe - Pb)                  # rate at the bubble point
    qo_max = qb + J * Pb / 1.8          # absolute open flow
    for i, pw in enumerate(Pwf):
        if pw >= Pb:
            q[i] = J * (Pe - pw)
        else:
            q[i] = qb + (qo_max - qb) * (1 - 0.2*(pw/Pb) - 0.8*(pw/Pb)**2)
    return q


# ── Verified operating-rate solver (IPR ∩ VLP via brentq) (do not edit) ────
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


# ── Composite-IPR closed-form inverse: Pwf(q) (do not edit) ────────────────
# Same branch find_operating_rate uses internally: linear above Pb, Vogel below.
def composite_pwf(Pe, Pb, J, q):
    """Flowing bottomhole pressure (psi) at oil rate q from the composite IPR."""
    qb = J * (Pe - Pb)
    qo_max = qb + J * Pb / 1.8
    if q <= qb:
        return Pe - q / J
    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.0)


# ── Embedded OML well (Exercise 12.6) + tubing-ID ladder (do not edit) ─────
BASE_PARAMS = {'Pe': 4000.0, 'Pb': 2800.0, 'J': 7.0,
               'Pwh': 250.0, 'depth': 8500.0, 'tubing_id': 2.441}
TUBING_SIZES = [1.995, 2.441, 2.992, 3.476, 3.958]   # standard tubing IDs (in)


def nodal_report(params):
    """Full nodal-analysis report for one well.

    Returns a dict with keys:
      'q_op'           operating rate (STB/d) from find_operating_rate
      'pwf_op'         flowing bottomhole pressure at q_op (composite_pwf)
      'drawdown'       Pe - pwf_op (psi)
      'aof'            composite absolute open flow = J*(Pe-Pb) + J*Pb/1.8
      'sensitivity'    {param: (low_rate, base_rate, high_rate)} for Pe, J, Pwh, tubing_id
      'tornado'        {param: max(triple) - min(triple)}
      'recommendation' 'upsize tubing' | 'stimulate' | 'optimize wellhead/choke'
    """
    Pe = params['Pe']; Pb = params['Pb']; J = params['J']
    Pwh = params['Pwh']; depth = params['depth']; tubing_id = params['tubing_id']

    # Base operating point and its IPR pressure.
    q_op = find_operating_rate(Pe, Pb, J, Pwh, depth, tubing_id)
    pwf_op = composite_pwf(Pe, Pb, J, q_op)
    drawdown = Pe - pwf_op
    aof = J * (Pe - Pb) + J * Pb / 1.8

    # Next-smaller / next-larger tubing size on the standard ladder.
    idx = min(range(len(TUBING_SIZES)),
              key=lambda i: abs(TUBING_SIZES[i] - tubing_id))
    tid_low = TUBING_SIZES[max(idx - 1, 0)]
    tid_high = TUBING_SIZES[min(idx + 1, len(TUBING_SIZES) - 1)]

    # One-at-a-time low/high variations.
    variations = {
        'Pe': (Pe * 0.85, Pe * 1.15),
        'J': (J * 0.70, J * 1.30),
        'Pwh': (max(Pwh - 100.0, 50.0), Pwh + 100.0),
        'tubing_id': (tid_low, tid_high),
    }

    sensitivity = {}
    for key, (lo, hi) in variations.items():
        p_lo = dict(params); p_lo[key] = lo
        p_hi = dict(params); p_hi[key] = hi
        r_lo = find_operating_rate(p_lo['Pe'], p_lo['Pb'], p_lo['J'],
                                   p_lo['Pwh'], p_lo['depth'], p_lo['tubing_id'])
        r_hi = find_operating_rate(p_hi['Pe'], p_hi['Pb'], p_hi['J'],
                                   p_hi['Pwh'], p_hi['depth'], p_hi['tubing_id'])
        sensitivity[key] = (r_lo, q_op, r_hi)   # middle is the unperturbed base

    tornado = {k: max(v) - min(v) for k, v in sensitivity.items()}

    # The lever with the biggest swing drives the recommendation.
    max_param = max(tornado, key=tornado.get)
    if max_param == 'tubing_id':
        recommendation = 'upsize tubing'
    elif max_param == 'J':
        recommendation = 'stimulate'
    else:
        recommendation = 'optimize wellhead/choke'

    return {
        'q_op': q_op,
        'pwf_op': pwf_op,
        'drawdown': drawdown,
        'aof': aof,
        'sensitivity': sensitivity,
        'tornado': tornado,
        'recommendation': recommendation,
    }


report = nodal_report(BASE_PARAMS)

q_op = float(report['q_op'])
aof = float(report['aof'])
recommendation = report['recommendation']

print("Operating rate q_op (STB/d):", q_op)
print("Absolute open flow  aof  (STB/d):", aof)
print("Recommendation:", recommendation)

lockCopying code is a Full Access feature.