Boundary Conditions - No-Flow vs Aquifer Support

Level 3

Chapter 11: Reservoir Simulation

descriptionProblem

The simulator uses no-flow boundaries by default (no flux at x=0 and x=L). Modify it to support a constant-pressure boundary at x=L (representing an aquifer). How does this change the pressure profile and the well's production rate over time?

---

The 1D simulator embedded below is sealed: it has no-flow boundaries at both ends (x=0 and x=L), so the only thing that moves fluid is a well. On block OML 58 we have producer OD-021 sitting at the heel (cell 0), pulling a steady q_stbd = -100 STB/d. With both ends sealed, the only energy available is the rock-and-fluid expansion of a finite, closed tank; the reservoir simply depletes.

Real reservoirs are rarely sealed. A water-bearing aquifer downdip acts like an effectively infinite tank held at the initial pressure $P_i$ : as the reservoir draws down, the aquifer feeds water across the boundary and supports the pressure. The cheap, standard way to bolt an aquifer onto a finite-difference model is to place a high-productivity BHP well in the last cell, pinned to $P_i$ . That cell can take or give as much fluid as needed to stay at $P_i$ , which is exactly a constant-pressure (Dirichlet) boundary at x=L.

The simulator (Reservoir1D, with add_well(cell, rate_stbd=.., bhp_psi=.., rw_ft=.., skin=..)) is embedded; do not re-derive it. The reservoir is set up for you (make_reservoir() returns a fresh, identical Reservoir1D each call). Your job is the boundary condition.

Your tasks:

Write run_bc(aquifer: bool, days=720.0, dt_days=10.0) -> (Pavg_psi, Pbnd_psi).

Start from sim = make_reservoir().
Always add the producer: sim.add_well(cell=0, rate_stbd=Q_PROD_STBD)

(Q_PROD_STBD = -100.0).

If aquifer is True, add the constant-pressure boundary: a BHP well in

the last cell (sim.nx - 1) pinned to bhp_psi = sim.Pi, using rw_ft = AQ_RW_FT and skin = AQ_SKIN (a near-wellbore radius gives this well a large Peaceman index, so the cell is clamped to $P_i$ ). If aquifer is False, add no boundary well: the model stays sealed.

Run with sim.run(days, dt_days).
Return a tuple (float(sim.P.mean()), float(sim.P[-1])): the final average

reservoir pressure and the pressure in the last (boundary) cell, both in psi.

Call run_bc(False) and store it in sealed (a (Pavg, Pbnd) tuple).

Call run_bc(True) and store it in supported.

> Think about it: with the aquifer the reservoir should deplete less. > supported's final $P_{avg}$ must be higher than sealed's, and the > boundary cell should sit right at $P_i$ (≈ 4000 psi) because the aquifer holds > it there. The sealed case has nowhere to draw energy from, so it bleeds down > hard. This single line (aquifer or no aquifer) is often the difference > between a field that needs water injection on day one and one that flows on its > own for years.

lightbulbHints (0/3)

Stuck? Reveal hints one at a time — they progress from nudge to near-solution.

codeYour solution

main.py

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve

# ── Embedded 1D single-phase implicit simulator (do not edit) ────────────
# Sealed (no-flow) on BOTH ends by default. Harmonic-mean transmissibility,
# real cross-sectional area, Peaceman well index. It conserves mass.
class Reservoir1D:
    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 setup: a finite tank with producer OD-021 at the heel ─────────
Q_PROD_STBD = -100.0   # producer at cell 0 (negative = withdrawal), STB/d
AQ_RW_FT = 0.30        # near-wellbore radius for the aquifer boundary well
AQ_SKIN = 0.05         # small skin -> large Peaceman index -> cell clamped to Pi

def make_reservoir():
    """Fresh, identical Reservoir1D for OML 58 (sealed on both ends)."""
    return Reservoir1D(nx=50, length_ft=5000.0, width_ft=500.0, height_ft=50.0,
                       k_md=100.0, phi=0.20, ct_psi=1.0e-5, mu_cp=1.2, Pi_psi=4000.0, Bo=1.2)

def run_bc(aquifer, days=720.0, dt_days=10.0):
    """Run OD-021 with a sealed (no-flow) or aquifer (constant-P) outer boundary.

Returns (Pavg_psi, Pbnd_psi): final average reservoir pressure and the
    pressure in the last (boundary) cell.
    """
    sim = make_reservoir()
    # TODO: always add the producer at cell 0:
    #   sim.add_well(cell=0, rate_stbd=Q_PROD_STBD)
    # TODO: if aquifer, add a high-PI BHP well at the LAST cell pinned to Pi:
    #   sim.add_well(cell=sim.nx - 1, bhp_psi=sim.Pi, rw_ft=AQ_RW_FT, skin=AQ_SKIN)
    # TODO: sim.run(days, dt_days)
    # TODO: return (float(sim.P.mean()), float(sim.P[-1]))
    pass

sealed = None     # TODO: run_bc(False)
supported = None  # TODO: run_bc(True)

print("sealed   (Pavg, Pbnd):", sealed)
print("aquifer  (Pavg, Pbnd):", supported)

visibilityReveal reference solutionexpand_more

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


# ── Embedded 1D single-phase implicit simulator (do not edit) ────────────
# Sealed (no-flow) on BOTH ends by default. Harmonic-mean transmissibility,
# real cross-sectional area, Peaceman well index. It conserves mass.
class Reservoir1D:
    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 setup: a finite tank with producer OD-021 at the heel ─────────
Q_PROD_STBD = -100.0   # producer at cell 0 (negative = withdrawal), STB/d
AQ_RW_FT = 0.30        # near-wellbore radius for the aquifer boundary well
AQ_SKIN = 0.05         # small skin -> large Peaceman index -> cell clamped to Pi


def make_reservoir():
    """Fresh, identical Reservoir1D for OML 58 (sealed on both ends)."""
    return Reservoir1D(nx=50, length_ft=5000.0, width_ft=500.0, height_ft=50.0,
                       k_md=100.0, phi=0.20, ct_psi=1.0e-5, mu_cp=1.2, Pi_psi=4000.0, Bo=1.2)


def run_bc(aquifer, days=720.0, dt_days=10.0):
    """Run OD-021 with a sealed (no-flow) or aquifer (constant-P) outer boundary.

    Returns (Pavg_psi, Pbnd_psi): final average reservoir pressure and the
    pressure in the last (boundary) cell.
    """
    sim = make_reservoir()
    sim.add_well(cell=0, rate_stbd=Q_PROD_STBD)
    if aquifer:
        # constant-pressure (aquifer) boundary: high-PI BHP well at the last
        # cell pinned to Pi. The near-wellbore rw gives a large Peaceman index,
        # clamping the cell to Pi.
        sim.add_well(cell=sim.nx - 1, bhp_psi=sim.Pi, rw_ft=AQ_RW_FT, skin=AQ_SKIN)
    sim.run(days, dt_days)
    return (float(sim.P.mean()), float(sim.P[-1]))


sealed = run_bc(False)
supported = run_bc(True)

print("sealed   (Pavg, Pbnd):", sealed)
print("aquifer  (Pavg, Pbnd):", supported)

lockCopying code is a Full Access feature.

arrow_back

11.7 Analytical Validation - Simulator vs the Ei Line-Source

11.9 Corey Exponent Sensitivity

arrow_forward