Injection-Production Pair - Pressure Support & Mass Balance

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.sparse import diags
from scipy.sparse.linalg import spsolve


# ── Verified 1D single-phase implicit simulator (do not edit) ────────────
class Reservoir1D:
    """1D single-phase implicit simulator (harmonic-mean T, real area, Peaceman wells)."""
    def __init__(self, nx, length_ft, width_ft, height_ft, k_md, phi, ct_psi, mu_cp, Pi_psi, Bo=1.0):
        self.nx = nx; self.dx = length_ft / nx
        self.area = width_ft * height_ft; self.vb = self.area * self.dx; self.height = height_ft
        self.k = np.full(nx, k_md, float) if np.isscalar(k_md) else np.asarray(k_md, float)
        self.phi = phi; self.ct = ct_psi; self.mu = mu_cp; self.Bo = Bo
        self.P = np.full(nx, Pi_psi, float); self.Pi = Pi_psi; self.beta = 1.127e-3; self.wells = []
    def _T(self, i, j):
        kh = 2*self.k[i]*self.k[j]/(self.k[i]+self.k[j])
        return self.beta*kh*self.area/(self.mu*self.dx)
    def _wi(self, w):
        re = 0.2*self.dx
        return 2*np.pi*self.beta*self.k[w["cell"]]*self.height/(self.mu*(np.log(re/w["rw"])+w["skin"]))
    def add_well(self, cell, rate_stbd=None, bhp_psi=None, rw_ft=0.354, skin=0.0):
        self.wells.append(dict(cell=cell, rate=rate_stbd, bhp=bhp_psi, rw=rw_ft, skin=skin))
    def step(self, dt):
        n = self.nx; C = self.vb*self.phi*self.ct/(5.615*dt)
        lo = np.zeros(n); di = np.zeros(n); up = np.zeros(n); b = np.zeros(n)
        for i in range(n):
            tl = self._T(i, i-1) if i > 0 else 0.0; tr = self._T(i, i+1) if i < n-1 else 0.0
            di[i] = C + tl + tr
            if i > 0: lo[i] = -tl
            if i < n-1: up[i] = -tr
            b[i] = C*self.P[i]
        for w in self.wells:
            i = w["cell"]
            if w["rate"] is not None: b[i] += w["rate"]*self.Bo
            elif w["bhp"] is not None: wi = self._wi(w); di[i] += wi; b[i] += wi*w["bhp"]
        self.P = spsolve(diags([lo[1:], di, up[:-1]], [-1, 0, 1], format="csc"), b)
    def run(self, days, dt):
        for _ in range(int(days/dt)): self.step(dt)


# ── OML-58 reservoir constants (do not edit) ─────────────────────────────
NX = 50              # grid cells
LENGTH_FT = 2500.0   # reservoir length
WIDTH_FT = 500.0     # cross-section width
HEIGHT_FT = 40.0     # net pay thickness
K_MD = 150.0         # permeability
PHI = 0.22           # porosity
CT_PSI = 1.2e-5      # total compressibility, 1/psi
MU_CP = 1.1          # oil viscosity, cp
PI_PSI = 4000.0      # initial reservoir pressure
BO = 1.2             # formation volume factor, rb/stb
DAYS = 730.0         # ~2 years
DT_DAYS = 10.0       # time step


def run_pair(q_stbd, days=DAYS, dt_days=DT_DAYS):
    """Balanced injector/producer pair on OML 58.

    Injector adds +q_stbd at cell 0, producer pulls -q_stbd at cell NX-1.
    Returns (Pavg, gradient, P_profile):
      Pavg     = mean final reservoir pressure (psi); stays ~= PI_PSI by mass balance
      gradient = P_profile[0] - P_profile[-1] (psi); a steady few-hundred-psi drop
      P_profile = final pressure array of length NX
    """
    sim = Reservoir1D(NX, LENGTH_FT, WIDTH_FT, HEIGHT_FT, K_MD, PHI,
                      CT_PSI, MU_CP, PI_PSI, Bo=BO)
    sim.add_well(0, rate_stbd=+q_stbd)        # injector at the inlet
    sim.add_well(NX - 1, rate_stbd=-q_stbd)   # producer at the outlet
    sim.run(days, dt_days)
    P_profile = sim.P
    Pavg = float(np.mean(P_profile))
    gradient = float(P_profile[0] - P_profile[-1])
    return Pavg, gradient, P_profile


pavg_psi, gradient_psi, p_profile = run_pair(350.0)

print("Pavg (psi):", pavg_psi, "  deviation from Pi:", pavg_psi - PI_PSI)
print("gradient P[0]-P[-1] (psi):", gradient_psi)

lockCopying code is a Full Access feature.