Grid-Refinement Study - When Does the Answer Stop Moving?

Level 3

Chapter 11: Reservoir Simulation

descriptionProblem

Run the 1D simulator from this chapter with 10, 25, 50, 100, and 200 grid cells. Plot the pressure profile at t=180 days for each case. At what grid resolution does the solution stop changing significantly? Why does this matter for practical simulation work?

---

A reservoir model has no "correct" number of grid cells. Coarse grids run fast but smear the answer; fine grids are accurate but slow. The only honest way to pick a grid is a refinement study: run the same problem at increasing resolution and watch the answer settle. When refining the grid further stops moving the result by much, the solution has converged and the coarsest grid that is "close enough" is the one you ship.

You will run a depletion case on OML-58: a 2000-ft × 1000-ft × 100-ft sand (k = 100 md, phi = 0.2, ct = 1.5e-5 /psi, mu = 2 cp, Pi = 4000 psi, Bo = 1.2) with a single producer at cell 0 held on constant BHP of 2500 psi. Each run advances 180 days with dt = 1 day. The Reservoir1D simulator is embedded; do not re-derive it.

Your tasks:

Write metric_for_nx(nx).

Build `Reservoir1D(nx, length_ft=LENGTH_FT, width_ft=WIDTH_FT,

height_ft=HEIGHT_FT, k_md=K_MD, phi=PHI, ct_psi=CT_PSI, mu_cp=MU_CP, Pi_psi=PI_PSI, Bo=BO)`.

Add the producer: res.add_well(0, bhp_psi=BHP_PSI).
Run the depletion: res.run(DAYS, DT_DAYS).
Return the average reservoir pressure at the end of the run,

float(res.P.mean()) (in psi). This single scalar is our convergence metric: the number we watch as the grid gets finer.

Write refine(nx_list) returning a dict mapping each nx to its metric:

{nx: metric_for_nx(nx) for nx in nx_list}.

Call refine([20, 40, 80]) and store it in metric_by_nx.

> Think about it: every run depletes the same reservoir, so the average > pressure they report should agree more and more as the cells shrink. The test > of convergence is that the successive differences shrink: the gap between > the nx=80 and nx=40 answers is smaller than the gap between nx=40 and nx=20. > A coarse grid over-predicts how fast pressure falls near the well; refining > resolves that boundary layer and the answer creeps toward the true value by a > diminishing amount each time you double the cells. That diminishing change is > exactly what "the grid is fine enough" looks like, and why nobody runs a > million-cell model when ten thousand cells already stopped moving the number.

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) ────────────
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 depletion case: one BHP producer at cell 0 ────────────────────
LENGTH_FT = 2000.0
WIDTH_FT = 1000.0
HEIGHT_FT = 100.0
K_MD = 100.0
PHI = 0.2
CT_PSI = 1.5e-5
MU_CP = 2.0
PI_PSI = 4000.0
BO = 1.2
BHP_PSI = 2500.0     # producer at cell 0 held on constant BHP
DAYS = 180.0
DT_DAYS = 1.0
NX_LIST = [20, 40, 80]

def metric_for_nx(nx):
    """Run the OML-58 depletion at resolution `nx` and return the convergence
    metric: the average reservoir pressure (psi) at the end of the run."""
    # TODO: res = Reservoir1D(nx, LENGTH_FT, WIDTH_FT, HEIGHT_FT, K_MD, PHI,
    #                         CT_PSI, MU_CP, PI_PSI, Bo=BO)
    # TODO: res.add_well(0, bhp_psi=BHP_PSI)
    # TODO: res.run(DAYS, DT_DAYS)
    # TODO: return float(res.P.mean())
    pass

def refine(nx_list):
    """Dict mapping each grid resolution nx -> its convergence metric (psi)."""
    # TODO: return {nx: metric_for_nx(nx) for nx in nx_list}
    pass

metric_by_nx = None  # TODO: refine(NX_LIST)

print("metric (avg pressure, psi) by nx:", metric_by_nx)

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) ────────────
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 depletion case: one BHP producer at cell 0 ────────────────────
LENGTH_FT = 2000.0
WIDTH_FT = 1000.0
HEIGHT_FT = 100.0
K_MD = 100.0
PHI = 0.2
CT_PSI = 1.5e-5
MU_CP = 2.0
PI_PSI = 4000.0
BO = 1.2
BHP_PSI = 2500.0     # producer at cell 0 held on constant BHP
DAYS = 180.0
DT_DAYS = 1.0
NX_LIST = [20, 40, 80]


def metric_for_nx(nx):
    """Run the OML-58 depletion at resolution `nx` and return the convergence
    metric: the average reservoir pressure (psi) at the end of the run."""
    res = Reservoir1D(nx, LENGTH_FT, WIDTH_FT, HEIGHT_FT, K_MD, PHI,
                      CT_PSI, MU_CP, PI_PSI, Bo=BO)
    res.add_well(0, bhp_psi=BHP_PSI)
    res.run(DAYS, DT_DAYS)
    return float(res.P.mean())


def refine(nx_list):
    """Dict mapping each grid resolution nx -> its convergence metric (psi)."""
    return {nx: metric_for_nx(nx) for nx in nx_list}


metric_by_nx = refine(NX_LIST)

print("metric (avg pressure, psi) by nx:", metric_by_nx)

lockCopying code is a Full Access feature.

arrow_back

11.10 Capstone - A 2D Quarter-Five-Spot Pressure Solver

11.2 Permeability Heterogeneity - A Low-k Streak Throttles the Flow

arrow_forward