""" Topic 2.4 — Limitations of Linear PI & Non-Linear IPR Models Reference implementation for the Karama KRM-4 worked examples and the chalk gas well back-pressure / LIT analysis. Well Productivity Programme · Course 01 · Module 02 · Topic 2.4 (c) 2026 Abalt Solutions Ltd — training use only. Models implemented ------------------ - vogel_q / vogel_ipr Vogel (1968) dimensionless IPR, P_bar <= Pb - vogel_qmax_single_point Qmax from a single stabilised test point - standing_q Standing modification for damaged wells (EF = FE) - fetkovich_qmax_no_test Qmax estimate when no below-Pb test exists - composite_q / composite_ipr Linear-above-Pb + Vogel-below-Pb composite - vogel_pwf_for_target / composite_pwf_for_target Required Pwf for a target rate (exact quadratic) - rawlins_schellhardt_fit Back-pressure (C, n) from multi-rate gas test - rawlins_schellhardt_q Qg = C (P_bar^2 - Pwf^2)^n - backpressure_aofp Gas AOFP at Pwf = 14.7 psia - lit_fit / lit_q / lit_aofp LIT (Jones) Δp²/q = A + B·q linearisation - linear_aofp Straight-line AOFP = J · P_bar (for comparison) `python3 ipr_models.py` reproduces the Topic 2.4 worked-example numbers: Vogel Qmax = 1,845 STB/d; composite Qb = 330 / Qmax = 1,547 / Q@3,100 = 638 STB/d; back-pressure n = 0.800, C ≈ 0.0190 (page prints 0.01906), AOFP ≈ 10.2 MMscf/d. """ import math PSC = 14.7 # psia, atmospheric back-pressure for AOFP definitions # ---------------------------------------------------------------------- # Vogel IPR (saturated reservoir, P_bar <= Pb) # ---------------------------------------------------------------------- def vogel_factor(pwf, p_bar): """Vogel dimensionless rate: 1 - 0.2 x - 0.8 x^2, x = Pwf/P_bar.""" x = pwf / p_bar return 1.0 - 0.2 * x - 0.8 * x * x def vogel_q(qmax, pwf, p_bar): """Rate (STB/d) from the Vogel equation.""" return qmax * vogel_factor(pwf, p_bar) def vogel_qmax_single_point(q_test, pwf_test, p_bar): """Qmax (AOFP) from one stabilised test point below the bubble point.""" return q_test / vogel_factor(pwf_test, p_bar) def vogel_ipr(qmax, p_bar, pwf_grid=None): """(Pwf, Q) pairs over a Pwf grid (default: P_bar down to 0 in 500-psi steps).""" if pwf_grid is None: pwf_grid = [p_bar] + list(range(int(p_bar) // 500 * 500, -1, -500)) pwf_grid = sorted(set(pwf_grid), reverse=True) return [(p, vogel_q(qmax, p, p_bar)) for p in pwf_grid] def standing_q(qmax, pwf, p_bar, ef=1.0): """Standing-modified Vogel for damaged/stimulated wells (EF = flow efficiency).""" x = pwf / p_bar return qmax * (1.0 - 0.2 * ef * x - 0.8 * (ef ** 2) * x * x) def fetkovich_qmax_no_test(j, pb): """No-test Qmax estimate at P_bar = Pb: Qmax = J·Pb/1.8.""" return j * pb / 1.8 def vogel_pwf_for_target(q_target, qmax, p_bar): """Exact quadratic solve of the Vogel equation for the Pwf giving q_target.""" # 0.8 x^2 + 0.2 x - (1 - q_target/qmax) = 0 c = -(1.0 - q_target / qmax) x = (-0.2 + math.sqrt(0.04 - 4 * 0.8 * c)) / (2 * 0.8) return x * p_bar # ---------------------------------------------------------------------- # Composite IPR (partially saturated: P_bar > Pb, Pwf may fall below Pb) # ---------------------------------------------------------------------- def composite_anchors(j, p_bar, pb): """(Qb, Qv, Qmax): bubble-point rate, Vogel add-on, total AOFP.""" qb = j * (p_bar - pb) qv = j * pb / 1.8 return qb, qv, qb + qv def composite_q(j, p_bar, pb, pwf): """Composite rate: linear for Pwf >= Pb, Qb + Vogel segment below Pb.""" if pwf >= pb: return j * (p_bar - pwf) qb, qv, _ = composite_anchors(j, p_bar, pb) return qb + qv * vogel_factor(pwf, pb) def composite_ipr(j, p_bar, pb, pwf_grid): """(Pwf, Q, segment) tuples across the supplied Pwf grid.""" return [(p, composite_q(j, p_bar, pb, p), 'linear' if p >= pb else 'vogel') for p in pwf_grid] def composite_pwf_for_target(q_target, j, p_bar, pb): """Required Pwf for a target rate on the composite IPR (exact solve).""" qb, qv, _ = composite_anchors(j, p_bar, pb) if q_target <= qb: return p_bar - q_target / j # still on the linear segment c = -(1.0 - (q_target - qb) / qv) x = (-0.2 + math.sqrt(0.04 - 4 * 0.8 * c)) / (2 * 0.8) return x * pb def linear_aofp(j, p_bar): """Straight-line AOFP = J · P_bar (only valid if the whole IPR stays above Pb).""" return j * p_bar # ---------------------------------------------------------------------- # Gas wells — Rawlins-Schellhardt back-pressure equation # ---------------------------------------------------------------------- def rawlins_schellhardt_fit(p_bar, points): """ Fit Qg = C (P_bar^2 - Pwf^2)^n by least squares on log-log coordinates. points: iterable of (Qg_Mscfd, Pwf_psia). Returns (C, n). """ lx = [math.log10(p_bar ** 2 - pwf ** 2) for _, pwf in points] ly = [math.log10(q) for q, _ in points] nn = len(points) mx, my = sum(lx) / nn, sum(ly) / nn n = sum((a - mx) * (b - my) for a, b in zip(lx, ly)) / \ sum((a - mx) ** 2 for a in lx) c = 10 ** (my - n * mx) return c, n def rawlins_schellhardt_q(c, n, p_bar, pwf): """Gas rate (Mscf/d) from the back-pressure equation.""" return c * (p_bar ** 2 - pwf ** 2) ** n def backpressure_aofp(c, n, p_bar, psc=PSC): """Gas AOFP (Mscf/d) at atmospheric flowing pressure.""" return rawlins_schellhardt_q(c, n, p_bar, psc) # ---------------------------------------------------------------------- # Gas wells — LIT / Jones analysis (Δp²/q = A + B·q) # ---------------------------------------------------------------------- def lit_fit(p_bar, points): """ Fit (P_bar^2 - Pwf^2)/Qg = A + B·Qg by least squares. points: iterable of (Qg_Mscfd, Pwf_psia). Returns (A, B). A — laminar (Darcy + mechanical skin) coefficient, psia²/(Mscf/d). B — turbulent (rate-dependent) coefficient, psia²/(Mscf/d)². """ qs = [q for q, _ in points] ys = [(p_bar ** 2 - pwf ** 2) / q for q, pwf in points] nn = len(points) mq, my = sum(qs) / nn, sum(ys) / nn b = sum((q - mq) * (y - my) for q, y in zip(qs, ys)) / \ sum((q - mq) ** 2 for q in qs) a = my - b * mq return a, b def lit_q(a, b, p_bar, pwf): """Gas rate (Mscf/d) from the LIT deliverability quadratic.""" dp2 = p_bar ** 2 - pwf ** 2 return (-a + math.sqrt(a * a + 4 * b * dp2)) / (2 * b) def lit_aofp(a, b, p_bar, psc=PSC): """LIT AOFP (Mscf/d) at atmospheric flowing pressure.""" return lit_q(a, b, p_bar, psc) def lit_turbulent_share(a, b, q): """Fraction of total Δp² consumed by the turbulent (B·q²) term at rate q.""" return b * q / (a + b * q) # ---------------------------------------------------------------------- # __main__ — reproduce the Topic 2.4 worked-example numbers # ---------------------------------------------------------------------- if __name__ == '__main__': print('=' * 72) print('Topic 2.4 reference calculations — KRM-4 & chalk gas well') print('=' * 72) # --- WE1: KRM-4 Vogel IPR at P_bar = Pb = 3,650 psia ----------------- PB = 3650.0 qmax = vogel_qmax_single_point(900.0, 2500.0, PB) print('\n[WE1] Vogel, single-point test 900 STB/d @ 2,500 psia, P_bar = Pb = 3,650:') print(f' Vogel factor at test point = {vogel_factor(2500, PB):.4f} (page: 0.4877)') print(f' Qmax = {qmax:,.0f} STB/d (page: 1,845)') print(f' Q at Pwf = 3,100 psia = {vogel_q(qmax, 3100, PB):,.0f} STB/d ' f'(target 1,200 -> missed by {1200 - vogel_q(qmax, 3100, PB):,.0f})') print(f' Fetkovich no-test check = {fetkovich_qmax_no_test(0.60, PB):,.0f} STB/d (J·Pb/1.8)') print(f' Mid-point check Q(0.5·P_bar) = {vogel_factor(0.5 * PB, PB):.3f} · Qmax (true Vogel value 0.70)') print(' Eight-row IPR table (recomputed; see README for source-page variances):') for pwf in [3650, 3000, 2500, 2000, 1500, 1000, 500, 0]: print(f' Pwf {pwf:>5,} psia -> Q = {vogel_q(qmax, pwf, PB):>7,.0f} STB/d') lin = linear_aofp(0.60, PB) print(f' Linear AOFP (J = 0.60) = {lin:,.0f} STB/d -> over-prediction ' f'{(lin / qmax - 1) * 100:+.1f}% (page: +18.7%)') pwf_req = vogel_pwf_for_target(1200.0, qmax, PB) print(f' Required Pwf for 1,200 STB/d = {pwf_req:,.0f} psia (exact); page prints 1,979 psia') print(f' ESP lift differential (page) = 3,100 - 1,979 = 1,121 psi') # --- WE2: composite IPR at P_bar = 4,200 psia, measured J = 0.60 ----- J, PBAR = 0.60, 4200.0 qb, qv, qmax_c = composite_anchors(J, PBAR, PB) print('\n[WE2] Composite IPR, P_bar = 4,200 psia, Pb = 3,650, J(measured) = 0.60:') print(f' Qb = J·(P_bar - Pb) = {qb:,.0f} STB/d (page: 330)') print(f' Vogel add-on J·Pb/1.8 = {qv:,.0f} STB/d (page: 1,217)') print(f' Qmax (composite AOFP) = {qmax_c:,.0f} STB/d (page: 1,547)') q3100 = composite_q(J, PBAR, PB, 3100) print(f' Q at Pwf = 3,100 psia = {q3100:,.0f} STB/d (page: 638 — TARGET MISSED)') print(' Curve points (recomputed):') for pwf, q, seg in composite_ipr(J, PBAR, PB, [4200, 3900, 3650, 3100, 2500, 2000, 1500, 1000, 500, 0]): print(f' Pwf {pwf:>5,} psia -> Q = {q:>7,.0f} STB/d [{seg}]') lin_c = linear_aofp(J, PBAR) print(f' Linear Q at 3,100 = {J * (PBAR - 3100):,.0f} STB/d ({(J * (PBAR - 3100) / q3100 - 1) * 100:+.1f}% vs composite)') print(f' Linear AOFP = {lin_c:,.0f} STB/d -> over-prediction ' f'{(lin_c / qmax_c - 1) * 100:+.1f}% (page: +62.9%)') print(f' Required Pwf for 1,200 STB/d = {composite_pwf_for_target(1200, J, PBAR, PB):,.0f} psia ' f'(exact); page prints 1,865 psia — see README') print(' NOTE: Intro-hub variant builds this composite with Jideal = 0.926 ->') qb_i, _, qmax_i = composite_anchors(0.926, PBAR, PB) print(f' Qb = {qb_i:,.0f}, Qmax = {qmax_i:,.0f} STB/d (hub: 509 / 2,387). Decks use measured J.') # --- WE3: chalk gas well back-pressure test -------------------------- PG = 3800.0 pts = [(2800.0, 3400.0), (5100.0, 2900.0), (7600.0, 2100.0)] print('\n[WE3] Chalk gas well, P_bar = 3,800 psia, three-point back-pressure test:') for i, (q, pwf) in enumerate(pts, 1): print(f' Point {i}: Qg = {q:>5,.0f} Mscf/d @ Pwf = {pwf:,.0f} psia -> ' f'dP2 = {PG**2 - pwf**2:>12,.0f} psia2') # two-point endpoint slope, as on the page: dp2_1, dp2_3 = PG**2 - 3400**2, PG**2 - 2100**2 n_page = math.log10(7600 / 2800) / math.log10(dp2_3 / dp2_1) c_page = 2800 / dp2_1 ** 0.800 print(f' n (endpoint slope, as on page) = {n_page:.4f} -> 0.800') print(f' C (with n = 0.800) = {c_page:.5f} Mscf/d/psia^1.6 ' f'(page prints 0.01906 via rounded 146,900 denominator)') c_fit, n_fit = rawlins_schellhardt_fit(PG, pts) print(f' Least-squares fit check : n = {n_fit:.3f}, C = {c_fit:.5f}') print(f' Validation at point 2 = {rawlins_schellhardt_q(0.01906, 0.800, PG, 2900):,.0f} ' f'vs 5,100 Mscf/d measured (<1%)') aofp = backpressure_aofp(0.01906, 0.800, PG) print(f' AOFP (C = 0.01906, n = 0.800) = {aofp:,.0f} Mscf/d = {aofp / 1000:.1f} MMscf/d (page: 10.2)') # --- LIT / Jones analysis on the same test (derived, SYNTHETIC) ------ a, b = lit_fit(PG, pts) print('\n[LIT] Jones linearisation of the same three points (derived — SYNTHETIC):') print(f' A (laminar) = {a:,.1f} psia2/(Mscf/d)') print(f' B (turbulent) = {b:.4f} psia2/(Mscf/d)2') print(f' D = B/A = {b / a:.2e} (Mscf/d)^-1') print(f' LIT AOFP = {lit_aofp(a, b, PG):,.0f} Mscf/d ' f'({lit_aofp(a, b, PG) / 1000:.1f} MMscf/d)') print(f' Turbulent share of dP2 at 5,100 = {lit_turbulent_share(a, b, 5100) * 100:.0f}%') print(f' Turbulent share at AOFP = {lit_turbulent_share(a, b, lit_aofp(a, b, PG)) * 100:.0f}%') # --- assertions: lock the published anchor numbers ------------------- assert round(qmax) == 1845 assert round(vogel_q(qmax, 3100, PB)) == 467 assert round(qb) == 330 and round(qmax_c) == 1547 and round(q3100) == 638 assert round(n_page, 3) == 0.800 assert abs(c_page - 0.01906) < 2e-4 assert round(aofp / 1000, 1) == 10.2 print('\nAll anchor values verified against the Topic 2.4 brief. OK.')