#!/usr/bin/env python3 """ sand_control.py — Sand-control & gravel-pack skin reference implementation for WPP Topic 3.6 (Gravel Pack & Sand Control Skin, Sg). Well Productivity Programme · Course 01 · Module 03 · Topic 3.6 Running PBL case: Gashaka GK-22 (Agbada Fm, Niger Delta — Module 03 skin-audit PBL). Implements: * mohr_coulomb_cdp — simplified Niger Delta critical-drawdown-pressure rule * saucier_sizing — Saucier (1974) gravel-median rule + commercial-mesh check * gravelpack_sg — gravel-pack skin (Furui simplified + detailed Darcy form) * productivity_ratio— PR = ln(re/rw) / (ln(re/rw) + Sd + Sg) * flow_efficiency — FE = 7/(7+S') (Standing, ln≈7 approximation) * radial_q — Darcy radial oil inflow at a given S' (canon constants) * selection_score — qualitative sand-control method screening Canon (TECH-BRIEF-m03 §1a/§1e, QUICK CANON CARD): GK-22 k 85 mD · h 42 ft · μo 1.8 cp · Bo 1.32 rb/stb · p̄R 4,200 psi · pwf 2,500 psi · rw 0.35 ft · re 1,650 ft · q 782 stb/d · UCS 0.8 MPa · Sw 0.22 · D50 sand 215 μm. Worked outputs the __main__ demo reproduces: CDP 62 psi · 12/20 mesh · Sg +0.12. Where a quantity is illustrative (degraded-pack scenarios, ageing rates) it is labelled SYNTHETIC — TRAINING USE ONLY in the CSVs and README. © 2026 Abalt Solutions Ltd. Licensed for personal educational use only. """ import math # ---------------------------------------------------------------------------- # GK-22 canonical constants # ---------------------------------------------------------------------------- K = 85.0 # mD, formation permeability H = 42.0 # ft, net pay (fully perforated, hp = h) MU = 1.8 # cp, oil viscosity BO = 1.32 # rb/stb, oil FVF PR = 4200.0 # psi, average reservoir pressure PWF = 2500.0 # psi, flowing bottomhole pressure RW = 0.35 # ft, wellbore radius RE = 1650.0 # ft, drainage radius Q_AUDIT = 782.0 # stb/d, DST/audit rate at S' = +14 UCS = 0.8 # MPa, unconfined compressive strength (Agbada) SW = 0.22 # current water saturation D50_SAND = 215.0 # micron, formation median grain size LN_APPROX = 7.0 # ln(0.472 re/rw) ≈ 7 (FE / J-ratio work) LN_FULL = math.log(0.472 * RE / RW) # ≈ 7.709 (full log, FE/skin work) LN_OPENHOLE = math.log(RE / RW) # ≈ 8.46 (open-hole PR basis) # ---------------------------------------------------------------------------- # Geomechanics — critical drawdown pressure # ---------------------------------------------------------------------------- def mohr_coulomb_cdp(ucs_mpa=UCS, sw=SW, const=100.0): """Simplified Niger Delta critical-drawdown-pressure rule (TECH-BRIEF §1a). CDP (psi) ≈ const × UCS (MPa) × (1 − Sw) Calibrated regional shortcut for the full Coates & Denoo (1981) form. Returns the largest safe drawdown (psi) before disaggregation-driven sand failure. """ return const * ucs_mpa * (1.0 - sw) def drawdown_ratio(drawdown_psi, ucs_mpa=UCS, sw=SW): """Ratio of an applied drawdown to the critical drawdown pressure.""" return drawdown_psi / mohr_coulomb_cdp(ucs_mpa, sw) # ---------------------------------------------------------------------------- # Gravel sizing — Saucier (1974) # ---------------------------------------------------------------------------- # Commercial gravel meshes: (label, low_micron, high_micron, d50_micron). # D50 values are the published representative medians (TECH-BRIEF §1e); these # are the supplier-quoted D50s, which for 16/30 (840) sit below the geometric # midpoint (890) — we carry the published value so deck/CSV/code agree. GRAVEL_MESHES = ( ('20/40', 420.0, 850.0, 635.0), ('16/30', 600.0, 1180.0, 840.0), ('12/20', 850.0, 1680.0, 1265.0), ) def saucier_window(d50_sand=D50_SAND, lo=5.0, hi=6.0): """Saucier (1974) target gravel-median window: D50,gravel = 5–6 × D50,sand.""" return lo * d50_sand, hi * d50_sand def saucier_sizing(d50_sand=D50_SAND, meshes=GRAVEL_MESHES, lo=5.0, hi=6.0): """Evaluate each commercial mesh against the Saucier window. Returns a list of dicts: mesh, d50_gravel, ratio, verdict (str). 'IDEAL' when 5 ≤ ratio ≤ 6; 'UNDERSIZED' below; 'OVERSIZED' above. """ out = [] for label, glo, ghi, d50_g in meshes: ratio = d50_g / d50_sand if ratio < lo: verdict = 'UNDERSIZED' elif ratio > hi: verdict = 'OVERSIZED' else: verdict = 'IDEAL' out.append({'mesh': label, 'd50_gravel': d50_g, 'ratio': ratio, 'verdict': verdict}) return out # ---------------------------------------------------------------------------- # Gravel-pack skin # ---------------------------------------------------------------------------- def furui_sg(k=K, kg=300000.0, h=H, hp=H, lp=0.75, np_count=168.0, rw=RW): """Furui et al. (2003) simplified Sg (lower bound — ignores the screen): Sg = (k/kg) × (h/hp) × (Lp/rw) × (1/Np) """ return (k / kg) * (h / hp) * (lp / rw) * (1.0 / np_count) def gravelpack_sg(q=Q_AUDIT, mu=MU, bo=BO, kg=800000.0, h=H, lp=0.75, rp=0.021, spf=8.0, hp=H, sscreen=0.12): """Detailed Darcy gravel-pack skin (TECH-BRIEF §1e working formula): Sg = 241.1·q·μ·B/(kg·h) × Lp/(2·rp·Np) + Sscreen, Np = spf × hp Returns (sg_total, tunnel_term, sscreen). """ np_count = spf * hp lead = 241.1 * q * mu * bo / (kg * h) geom = lp / (2.0 * rp * np_count) tunnel = lead * geom return tunnel + sscreen, tunnel, sscreen # ---------------------------------------------------------------------------- # Productivity, flow efficiency, inflow # ---------------------------------------------------------------------------- def productivity_ratio(sd, sg, ln_term=LN_OPENHOLE): """PR = ln(re/rw) / (ln(re/rw) + Sd + Sg) — fraction of open-hole ideal.""" return ln_term / (ln_term + sd + sg) def flow_efficiency(s_total, ln_term=LN_APPROX): """FE = ln_term / (ln_term + S') (Standing; ln≈7 by default).""" return ln_term / (ln_term + s_total) def radial_q(s_total, pr=PR, pwf=PWF, k=K, h=H, mu=MU, bo=BO, rw=RW, re=RE): """Darcy radial oil rate (stb/d) at total skin S', full-log ln term.""" ln_term = math.log(0.472 * re / rw) return (0.00708 * k * h * (pr - pwf)) / (mu * bo * (ln_term + s_total)) def radial_q_approx(s_total, pr=PR, pwf=PWF): """Inflow on the 7-approx canon: J·Δp with J = Jideal·7/(7+S'). Jideal (7-approx) = 1.380 stb/d/psi (TECH-BRIEF §1a). Gives the audit and design points exactly: 782 @ S'=14, 2,022 @ S'=+1.12, 2,346 @ S=0. Note the acid-only state uses the published J = 1.208 (q 2,056), not 1.380·7/8 — see TECH-BRIEF Inconsistency #2; for that state use canon_q(). """ j_ideal = 1.380 j = j_ideal * LN_APPROX / (LN_APPROX + s_total) return j * (pr - pwf) # Locked inflow rates at fixed pwf 2,500 (QUICK CANON CARD / TECH-BRIEF §1a). # The published J carries to ~4 sig figs; these are the canon stb/d the decks, # CSVs and narration all quote (Inconsistency #2 standardises the 7-approx set). CANON_Q = {14.0: 782.0, 1.12: 2022.0, 1.00: 2056.0, 0.0: 2346.0} def canon_q(s_total): """Locked inflow (stb/d) at pwf 2,500 for the standard skin states. Reproduces the QUICK CANON CARD exactly: 782 / 2,022 / 2,056 / 2,346 stb/d at S' = 14 / 1.12 / 1.00 / 0, so the 34-stb/d gravel-pack cost lands. """ return CANON_Q[s_total] # ---------------------------------------------------------------------------- # Method selection (qualitative screen) # ---------------------------------------------------------------------------- # (method, sg_low, sg_high, exclusion, note) — TECH-BRIEF §1e selection matrix. SELECTION_MATRIX = ( ('Rate-limit / none', 0.0, 0.0, 'none', 'Consolidated rock only; not suitable for UCS < 1 MPa.'), ('SAS (standalone screen)', 0.5, 2.0, 'moderate', 'Gap/erosion risk in fines-laden sand; possible but risky at GK-22.'), ('ICHGP (cased-hole GP)', 1.0, 3.0, 'good', 'RECOMMENDED for GK-22 (designed Sg 0.12–0.4, below generic range).'), ('OHGP (open-hole GP)', 0.5, 1.5, 'excellent', 'Excellent in high-k; needs open-hole completion.'), ('Frac-pack', -1.0, 1.0, 'excellent', 'Best for k < 20 mD; over-engineered/risky at 85 mD & UCS 0.8 MPa.'), ) def selection_score(k=K, ucs_mpa=UCS, cased=True): """Return the selection matrix annotated with a simple GK-22 fit flag. Heuristic only — encodes the qualitative reasoning on the source page: frac-pack penalised above 20 mD; SAS penalised in weak, fines-prone sand; ICHGP preferred for a cased, weak, high-k completion. """ rows = [] for method, lo, hi, excl, note in SELECTION_MATRIX: fit = 'consider' if method.startswith('Frac-pack') and k > 20: fit = 'over-engineered' elif method.startswith('SAS') and ucs_mpa < 1.0: fit = 'risky' elif method.startswith('Rate-limit'): fit = 'unsuitable' elif method.startswith('ICHGP') and cased and k > 50 and ucs_mpa < 1.0: fit = 'RECOMMENDED' elif method.startswith('OHGP') and cased: fit = 'needs open hole' rows.append({'method': method, 'sg_range': (lo, hi), 'exclusion': excl, 'fit': fit, 'note': note}) return rows # ---------------------------------------------------------------------------- # Demo — reproduces the Topic 3.6 worked / locked numbers # ---------------------------------------------------------------------------- if __name__ == '__main__': print('WPP Topic 3.6 — sand-control reference demo (Gashaka GK-22)') print('=' * 66) # 1 · Critical drawdown pressure cdp = mohr_coulomb_cdp() cdp_bt = mohr_coulomb_cdp(sw=0.80) dd = PR - PWF print(f'CDP (Sw 0.22) = {cdp:5.1f} psi (expect 62.4 ≈ 62)') print(f'CDP at breakthrough (0.80) = {cdp_bt:5.1f} psi (expect 16)') print(f'Current drawdown = {dd:5.0f} psi') print(f'Drawdown / CDP = {drawdown_ratio(dd):5.1f}x (expect 27.4x)') print(f'Post-acid pwf 1,500 = {drawdown_ratio(PR-1500):5.1f}x (expect ~44x)') # 2 · Saucier gravel sizing lo, hi = saucier_window() print(f'\nSaucier window (5-6x 215um)= {lo:.0f}-{hi:.0f} um (expect 1,075-1,290)') for r in saucier_sizing(): print(f" {r['mesh']:>6} D50 {r['d50_gravel']:6.0f} um " f"ratio {r['ratio']:.2f}x -> {r['verdict']}") print(' expect: 20/40 UNDERSIZED (2.95x), 16/30 UNDERSIZED (3.9x), ' '12/20 IDEAL (5.9x)') # 3 · Gravel-pack skin chain sg8, tun8, scr = gravelpack_sg(spf=8.0, kg=800000.0) sg4, tun4, _ = gravelpack_sg(spf=4.0, kg=300000.0, sscreen=0.20) print(f'\nSg 8 spf, kg 800k = {sg8:.4f} (tunnel {tun8:.5f} + ' f'screen {scr:.2f}) -> +0.12') print(f'Sg 4 spf, kg 300k = {sg4:.4f} (expect ~0.20)') print(f'Furui simplified (4 spf) = {furui_sg():.3e} (expect ~3.61e-06)') # 4 · Productivity ratio pr_now = productivity_ratio(sd=14.0, sg=0.0) pr_post = productivity_ratio(sd=1.0, sg=0.12) print(f'\nPR today (Sd 14, no pack) = {pr_now:.3f} (expect 0.377)') print(f'PR post (Sd 1 + Sg 0.12) = {pr_post:.3f} (expect 0.883)') print(f'PR gain = {pr_post/pr_now:.2f}x (expect 2.34x)') # 5 · Flow efficiency and inflow on the canon print(f'\nFE S\' +14 = {flow_efficiency(14.0):.3f} (expect 0.333)') print(f'FE S\' +1.12 = {flow_efficiency(1.12):.3f} (expect 0.862)') print(f'FE S +1.00 = {flow_efficiency(1.00):.3f} (expect 0.875)') print(f'q S\' +14 (locked canon) = {canon_q(14.0):6.0f} stb/d (expect 782)') print(f'q S\' +1.12 (locked canon) = {canon_q(1.12):6.0f} stb/d (expect 2,022)') print(f'q S +1.00 (locked canon) = {canon_q(1.00):6.0f} stb/d (expect 2,056)') print(f'q S 0 (locked canon) = {canon_q(0.0):6.0f} stb/d (expect 2,346)') print(f'GP cost vs acid-only = {canon_q(1.00)-canon_q(1.12):4.0f} ' f'stb/d (expect 34)') # 6 · Selection screen print('\nSand-control selection screen (GK-22: 85 mD, UCS 0.8 MPa, cased):') for r in selection_score(): print(f" {r['method']:<26} Sg {r['sg_range'][0]:+.1f}..{r['sg_range'][1]:+.1f} " f"{r['exclusion']:<10} -> {r['fit']}")