""" nondarcy.py — Non-Darcy (rate-dependent) skin toolkit for Topic 3.2, Well Productivity Programme · Course 01 · Module 03 (Abalt Solutions). Reference implementation behind Lectures 3.2a ("Total Skin — Breaking the Well Test Number Apart") and 3.2b ("Calculating D — Theory, Correlations, and Field Practice"). Physics ------- Forchheimer (1901): -dP/dr = (mu/k)*u + beta*rho*u**2 Total skin: S' = S + D*q (D in d/Mscf gas, d/stb oil) Two-rate separation: D = (S'2 - S'1)/(q2 - q1); S = S'1 - D*q1 Theoretical D (BP CDM Eq 4a.25-26): D = 2.22e-15 * beta * gamma_g * k * h / (hp**2 * rw * mu_wf) [d/Mscf] (oil analogue: gamma_o, ko, mu_o -> d/stb) Beta correlations: Katz (Eq 4a.23): beta = 4.11e10 / k**1.33 [1/ft] Dake (Eq 4a.24): beta = 2.73e10 / k**1.10 [1/ft] Effective open interval from the D ratio: hp_eff = hp_nominal / sqrt(D_measured / D_calculated) Canonical cases (TECH-BRIEF-m03 §1b / §2 T3.2) ---------------------------------------------- FK-7 (Fulmar gas well): kg 12 mD, h 65 ft, hp 52 ft, rw 0.35 ft, gamma_g 0.62, mu_wf 0.022 cP, T 620 R. Two-rate DST: S' = +11.2 at 8 MMscf/d, +18.4 at 18 MMscf/d -> D = 0.720 d/MMscf, S = +5.44, D_theory = 0.086 d/MMscf, ratio 8.4x -> hp_eff ~= 18 ft (~35% open), Dq at 25 MMscf/d design rate = 18.0. GK-22 (Gashaka oil well): ko 85 mD, h = hp = 42 ft, rw 0.35 ft, gamma_o 0.85, mu_o 1.8 cP, q 782 stb/d -> D ~= 9.4e-7 d/stb, Dq ~= 0.001 (negligible: S = S' = +14). NOTE on published beta values: the course pages print beta_Katz(12 mD) = 4.11e10/24.8 = 1.66e9 1/ft and beta_Katz(85 mD) = 4.11e10/265.3 = 1.55e8 1/ft. Exact evaluation gives 12**1.33 = 27.25 (-> 1.51e9) and 85**1.33 = 368.3 (-> 1.12e8). The decks and this demo carry the PUBLISHED values so the canonical chain (0.086 / 8.4x / 18 ft) reproduces exactly; beta_katz() itself evaluates the formula exactly. See README.md. """ import math # ---------------------------------------------------------------------- # Forchheimer pressure gradient (consistent units) # ---------------------------------------------------------------------- def forchheimer(u, mu, k, beta, rho): """Forchheimer pressure gradient -dP/dr = (mu/k)*u + beta*rho*u**2. Consistent-unit form (e.g. SI). Returns (total, viscous, inertial) so callers can inspect the share of each term. """ viscous = (mu / k) * u inertial = beta * rho * u ** 2 return viscous + inertial, viscous, inertial # ---------------------------------------------------------------------- # Beta correlations (k in mD, beta in 1/ft) # ---------------------------------------------------------------------- def beta_katz(k_md): """Katz correlation (BP CDM Eq 4a.23): beta = 4.11e10 / k^1.33 [1/ft].""" return 4.11e10 / k_md ** 1.33 def beta_dake(k_md): """Dake correlation (BP CDM Eq 4a.24): beta = 2.73e10 / k^1.10 [1/ft].""" return 2.73e10 / k_md ** 1.10 # ---------------------------------------------------------------------- # Theoretical D from beta # ---------------------------------------------------------------------- def d_from_beta(beta, gamma, k_md, h_ft, hp_ft, rw_ft, mu_cp): """Theoretical non-Darcy coefficient D (BP CDM Eq 4a.25-26). D = 2.22e-15 * beta * gamma * k * h / (hp^2 * rw * mu) Gas: gamma_g, kg, mu at flowing wellbore conditions -> D in d/Mscf. Oil analogue: gamma_o, ko, mu_o -> D in d/stb. Note D ~ 1/hp^2: halving the open interval quadruples D. """ return 2.22e-15 * beta * gamma * k_md * h_ft / (hp_ft ** 2 * rw_ft * mu_cp) # ---------------------------------------------------------------------- # Multi-rate separation of S and D # ---------------------------------------------------------------------- def multirate_fit(rates, total_skins): """Least-squares straight line through (q, S') points: S' = S + D*q. Returns (D, S). With exactly two points this reduces to the two-rate formulas D = (S'2-S'1)/(q2-q1), S = S'1 - D*q1. Units of D follow the rate units supplied (MMscf/d -> d/MMscf). """ n = len(rates) if n < 2: raise ValueError("need at least two (rate, S') points") sx = sum(rates) sy = sum(total_skins) sxx = sum(q * q for q in rates) sxy = sum(q * s for q, s in zip(rates, total_skins)) d = (n * sxy - sx * sy) / (n * sxx - sx * sx) s = (sy - d * sx) / n return d, s def hp_eff_from_d_ratio(hp_nominal_ft, d_measured, d_calculated): """Effective flowing interval: hp_eff = hp_nom / sqrt(D_meas/D_calc). Because D ~ 1/hp^2, a measured D exceeding the clean-perforation theoretical D by a factor R implies only hp_nom/sqrt(R) of the interval is actually taking flow (plugged perforations). """ return hp_nominal_ft / math.sqrt(d_measured / d_calculated) # ---------------------------------------------------------------------- # Demo: reproduces the Module 03 brief's worked numbers # ---------------------------------------------------------------------- if __name__ == "__main__": print("=" * 68) print("FK-7 (Fulmar gas well) — full Topic 3.2 workup") print("=" * 68) # -- two-rate separation ------------------------------------------- rates = [8.0, 18.0] # MMscf/d skins = [11.2, 18.4] # S' (total skin) D_meas, S_true = multirate_fit(rates, skins) # d/MMscf print(f"Two-rate test: S' = {skins[0]} @ {rates[0]:.0f} MMscf/d, " f"S' = {skins[1]} @ {rates[1]:.0f} MMscf/d") print(f" D measured = {D_meas:.3f} d/MMscf (target 0.720)") print(f" S true skin = +{S_true:.2f} (target +5.44)") assert abs(D_meas - 0.720) < 1e-9 assert abs(S_true - 5.44) < 1e-9 for q in (8, 18, 25): dq = D_meas * q sp = S_true + dq tag = "design forecast" if q == 25 else "measured" print(f" q = {q:>2} MMscf/d: Dq = {dq:5.2f}, S' = {sp:5.2f} " f"({dq / sp * 100:.0f}% turbulence, {tag})") assert abs(D_meas * 25 - 18.0) < 1e-9 # Dq @ 25 MMscf/d # -- theoretical D (clean perforations) ---------------------------- BETA_PUBLISHED = 1.66e9 # 1/ft, as printed (4.11e10/"24.8") beta_exact = beta_katz(12) print(f"\nKatz beta at kg = 12 mD: published 1.66e+09 1/ft " f"(exact formula {beta_exact:.2e}; see README quirk note)") print(f"Dake beta at kg = 12 mD: {beta_dake(12):.2e} 1/ft " f"(correlations bracket ~2-5x; run both)") D_calc_Mscf = d_from_beta(BETA_PUBLISHED, gamma=0.62, k_md=12, h_ft=65, hp_ft=52, rw_ft=0.35, mu_cp=0.022) D_calc = D_calc_Mscf * 1000.0 # d/Mscf -> d/MMscf print(f"Theoretical D = {D_calc_Mscf:.2e} d/Mscf = " f"{D_calc:.3f} d/MMscf (target 0.086)") assert abs(D_calc - 0.086) < 0.001 # -- plugged-perforation diagnosis --------------------------------- ratio = D_meas / D_calc hp_eff = hp_eff_from_d_ratio(52.0, D_meas, D_calc) print(f"D_measured / D_theoretical = {ratio:.1f}x (target 8.4x)") print(f"hp,eff = 52/sqrt({ratio:.1f}) = {hp_eff:.1f} ft " f"= {hp_eff / 52 * 100:.0f}% of perforations open " f"(target ~18 ft, ~35%)") assert abs(ratio - 8.4) < 0.05 assert abs(hp_eff - 18.0) < 0.2 # -- post-workover forecast ---------------------------------------- dq18_post = D_calc * 18 print(f"Post-TCP-reperforation (D -> {D_calc:.3f}): Dq @ 18 MMscf/d " f"= {dq18_post:.2f} vs {D_meas * 18:.2f} now " f"(-{D_meas * 18 - dq18_post:.1f} skin units, no acid)") print() print("=" * 68) print("GK-22 (Gashaka oil well) — why Dq is negligible") print("=" * 68) BETA_GK22_PUBLISHED = 1.55e8 # 1/ft, as printed (4.11e10/"265.3") print(f"Katz beta at ko = 85 mD: published 1.55e+08 1/ft " f"(exact formula {beta_katz(85):.2e})") D_oil = d_from_beta(BETA_GK22_PUBLISHED, gamma=0.85, k_md=85, h_ft=42, hp_ft=42, rw_ft=0.35, mu_cp=1.8) dq = D_oil * 782 print(f"D_oil = {D_oil:.2e} d/stb (target 9.43e-07)") print(f"Dq at 782 stb/d = {dq:.2e} ~= 0.001 " f"({dq / 14 * 100:.3f}% of S' = +14) -> S = S' = +14") assert abs(D_oil - 9.4e-7) < 0.1e-7 assert dq < 0.001 print("\nAll canonical checks passed: 0.72 / +5.44 / 8.4x / 18 ft / 18.0")