Analytical Validation - Simulator vs the Ei Line-Source

Level 3

Chapter 11: Reservoir Simulation

descriptionProblem

The line-source solution for pressure drawdown in an infinite reservoir is $P(r,t) = P_i + \frac{q\mu}{4\pi kh}Ei\left(-\frac{\phi\mu c_t r^2}{4kt}\right)$ . Compare the simulator's output at early times (before the pressure front reaches the boundary) to this analytical solution. How close is the numerical result? Hint: Use scipy.special.expi() for the exponential integral function.

---

A simulator you cannot validate is just an expensive random-number generator. The gold-standard check is to reproduce a problem with a known closed-form answer. For a single constant-rate well in an effectively infinite reservoir at early time, the pressure drawdown follows the line-source (exponential-integral) solution

\Delta P(r,t) = P(r,t) - P_i = -\frac{70.6\,q\,B\,\mu}{k\,h}\;\operatorname{Ei}\!\left(-\frac{948\,\phi\,\mu\,c_t\,r^2}{k\,t}\right)

where $\operatorname{Ei}$ is the exponential integral (scipy.special.expi). This is a radial solution, so the embedded simulator here is a radial version of the implicit single-phase model, same verified numerics (harmonic-mean $kh$ transmissibility, implicit tridiagonal solve, $C = V_b\phi c_t/(5.615\,\Delta t)$ accumulation, $\beta = 1.127\times10^{-3}$ , a rate well adds $q\,B$ to the RHS) but on a logarithmic radial grid so it can be compared to $\operatorname{Ei}$ .

Reservoir OD-021 on OML 58: $P_i = 4000$ psi, $k = 80$ md, $\phi = 0.22$ , $c_t = 1.0\times10^{-5}$ psi $^{-1}$ , $\mu = 1.0$ cp, $h = 40$ ft, $B = 1.0$ . A single producer at the wellbore pulls $q = 300$ STB/d. The grid runs from $r_w = 0.25$ ft to $r_e = 20000$ ft (far enough that the pressure front never reaches it in $t = 3$ days. The reservoir is infinite-acting).

Two sign / unit conventions to respect (this is where everyone slips):

Producer rate is negative. The simulator adds rate_stbd*B to the RHS, so a

producer is add_well(0, rate_stbd=-q_stbd). Pass that same negative q to ei_drawdown so both come out negative (a drawdown, not a buildup).

The 948 constant is the field-units form for t in HOURS. The simulator

marches in days, so convert: t_hours = t_days * 24 before calling ei_drawdown.

Your tasks (the radial Reservoir1D is embedded; do not edit it):

Write ei_drawdown(r, t, q, k, phi, mu, ct, h, B) to return the line-source

$\Delta P$ above using expi. With a producer (q < 0) the result must be negative.

Write simulate_drawdown(nr=NR, r_obs_ft=R_OBS_FT):

Build Reservoir1D(nr, RW_FT, RE_FT, H_FT, K_MD, PHI, CT_PSI, MU_CP, PI_PSI, BO).
Add a producer at the innermost cell: add_well(0, rate_stbd=-Q_STBD).
Run for T_DAYS with step DT_DAYS.
Find the cell whose center radius res.r is closest to r_obs_ft

(int(np.argmin(np.abs(res.r - r_obs_ft)))).

Return a tuple (r_mid_ft, dP_sim_psi) = that cell's center radius and its

drawdown res.P[cell] - PI_PSI.

Call simulate_drawdown() and store the tuple in sim_result. Then compute the

analytical drawdown at that same radius and time (remember: hours!) into dP_ei_psi, and the relative error rel_err = abs(dP_sim_psi - dP_ei_psi) / abs(dP_ei_psi) into rel_err.

> Think about it: at $r \approx 75$ ft and $t = 3$ days the simulator and the > Ei prediction should agree to better than 1 %. Drawdown is linear in rate > (double $q$ , double $\Delta P$ ) and shrinks with radius. Why does the match > degrade if you (a) forget to convert days to hours, or (b) shrink $r_e$ so the > front reaches the boundary before $t = 3$ days?

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
from scipy.special import expi

# ── Embedded RADIAL implicit simulator (do not edit) ─────────────────────
# Same verified scheme as the Cartesian Reservoir1D - harmonic-mean kh
# transmissibility, implicit tridiagonal spsolve, accumulation
# C = Vb*phi*ct/(5.615*dt), beta = 1.127e-3, a rate well adds q_stbd*B to the
# RHS - but on a logarithmic RADIAL grid so it can be checked against the
# radial line-source (Ei) analytical solution. res.r holds the block-center
# radii (ft); res.P holds the cell pressures (psi).
class Reservoir1D:
    def __init__(self, nr, rw_ft, re_ft, h_ft, k_md, phi, ct_psi, mu_cp, Pi_psi, Bo=1.0):
        self.nr = nr; self.h = h_ft
        self.k = np.full(nr, k_md, float)
        self.phi = phi; self.ct = ct_psi; self.mu = mu_cp; self.Bo = Bo
        self.P = np.full(nr, Pi_psi, float); self.Pi = Pi_psi; self.beta = 1.127e-3
        u = np.log(re_ft / rw_ft) / nr
        self.rb = rw_ft * np.exp(u * np.arange(nr + 1))          # block boundaries
        self.r = np.sqrt(self.rb[:-1] * self.rb[1:])             # block-center radii
        self.vb = np.pi * (self.rb[1:] ** 2 - self.rb[:-1] ** 2) * self.h  # annular Vb
        self.wells = []

def _T(self, i, j):
        kh = 2 * self.k[i] * self.k[j] / (self.k[i] + self.k[j])
        ro, ri = self.r[max(i, j)], self.r[min(i, j)]
        return 2 * np.pi * self.beta * kh * self.h / (self.mu * np.log(ro / ri))

def add_well(self, cell, rate_stbd):
        self.wells.append(dict(cell=cell, rate=rate_stbd))

def step(self, dt):
        n = self.nr; 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[i] + tl + tr
            if i > 0: lo[i] = -tl
            if i < n - 1: up[i] = -tr
            b[i] = C[i] * self.P[i]
        for w in self.wells:
            b[w["cell"]] += w["rate"] * self.Bo
        self.P = spsolve(diags([lo[1:], di, up[:-1]], [-1, 0, 1], format="csc"), b)

def run(self, days, dt):
        for _ in range(int(round(days / dt))):
            self.step(dt)

# ── OML 58 reservoir OD-021: single constant-rate producer, infinite-acting ─
PI_PSI = 4000.0     # initial reservoir pressure
K_MD = 80.0         # permeability
PHI = 0.22          # porosity
CT_PSI = 1.0e-5     # total compressibility, 1/psi
MU_CP = 1.0         # oil viscosity
H_FT = 40.0         # net pay
BO = 1.0            # formation volume factor
RW_FT = 0.25        # wellbore (innermost) radius
RE_FT = 20000.0     # outer radius - far enough to stay infinite-acting at T_DAYS
NR = 300            # radial cells
Q_STBD = 300.0      # producer rate (added as a NEGATIVE rate at the wellbore)
T_DAYS = 3.0        # run time
DT_DAYS = 0.02      # time step
R_OBS_FT = 75.0     # mid-field observation radius

def ei_drawdown(r, t, q, k, phi, mu, ct, h, B):
    """Line-source (exponential-integral) drawdown dP = P(r,t) - Pi.

t is in HOURS (the 948 constant is the field-units hours form). Pass q < 0
    for a producer so the returned drawdown is NEGATIVE.
    """
    # TODO: return -70.6*q*B*mu/(k*h) * expi(-948.0*phi*mu*ct*r**2/(k*t))
    pass

def simulate_drawdown(nr=NR, r_obs_ft=R_OBS_FT):
    """Run the radial simulator and return (r_mid_ft, dP_sim_psi) at the cell
    whose center radius is closest to r_obs_ft."""
    # TODO: res = Reservoir1D(nr, RW_FT, RE_FT, H_FT, K_MD, PHI, CT_PSI, MU_CP, PI_PSI, BO)
    # TODO: res.add_well(0, rate_stbd=-Q_STBD)     # producer = NEGATIVE rate
    # TODO: res.run(T_DAYS, DT_DAYS)
    # TODO: cell = int(np.argmin(np.abs(res.r - r_obs_ft)))
    # TODO: return float(res.r[cell]), float(res.P[cell] - PI_PSI)
    pass

sim_result = None        # TODO: simulate_drawdown()
r_mid_ft = None          # TODO: sim_result[0]
dP_sim_psi = None        # TODO: sim_result[1]
dP_ei_psi = None         # TODO: ei_drawdown(r_mid_ft, T_DAYS*24.0, -Q_STBD, K_MD, PHI, MU_CP, CT_PSI, H_FT, BO)
rel_err = None           # TODO: abs(dP_sim_psi - dP_ei_psi) / abs(dP_ei_psi)

print("r_mid (ft):", r_mid_ft)
print("dP_sim (psi):", dP_sim_psi, " dP_ei (psi):", dP_ei_psi)
print("relative error:", rel_err)

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
from scipy.special import expi


# ── Embedded RADIAL implicit simulator (do not edit) ─────────────────────
# Same verified scheme as the Cartesian Reservoir1D - harmonic-mean kh
# transmissibility, implicit tridiagonal spsolve, accumulation
# C = Vb*phi*ct/(5.615*dt), beta = 1.127e-3, a rate well adds q_stbd*B to the
# RHS - but on a logarithmic RADIAL grid so it can be checked against the
# radial line-source (Ei) analytical solution. res.r holds the block-center
# radii (ft); res.P holds the cell pressures (psi).
class Reservoir1D:
    def __init__(self, nr, rw_ft, re_ft, h_ft, k_md, phi, ct_psi, mu_cp, Pi_psi, Bo=1.0):
        self.nr = nr; self.h = h_ft
        self.k = np.full(nr, k_md, float)
        self.phi = phi; self.ct = ct_psi; self.mu = mu_cp; self.Bo = Bo
        self.P = np.full(nr, Pi_psi, float); self.Pi = Pi_psi; self.beta = 1.127e-3
        u = np.log(re_ft / rw_ft) / nr
        self.rb = rw_ft * np.exp(u * np.arange(nr + 1))          # block boundaries
        self.r = np.sqrt(self.rb[:-1] * self.rb[1:])             # block-center radii
        self.vb = np.pi * (self.rb[1:] ** 2 - self.rb[:-1] ** 2) * self.h  # annular Vb
        self.wells = []

    def _T(self, i, j):
        kh = 2 * self.k[i] * self.k[j] / (self.k[i] + self.k[j])
        ro, ri = self.r[max(i, j)], self.r[min(i, j)]
        return 2 * np.pi * self.beta * kh * self.h / (self.mu * np.log(ro / ri))

    def add_well(self, cell, rate_stbd):
        self.wells.append(dict(cell=cell, rate=rate_stbd))

    def step(self, dt):
        n = self.nr; 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[i] + tl + tr
            if i > 0: lo[i] = -tl
            if i < n - 1: up[i] = -tr
            b[i] = C[i] * self.P[i]
        for w in self.wells:
            b[w["cell"]] += w["rate"] * self.Bo
        self.P = spsolve(diags([lo[1:], di, up[:-1]], [-1, 0, 1], format="csc"), b)

    def run(self, days, dt):
        for _ in range(int(round(days / dt))):
            self.step(dt)


# ── OML 58 reservoir OD-021: single constant-rate producer, infinite-acting ─
PI_PSI = 4000.0     # initial reservoir pressure
K_MD = 80.0         # permeability
PHI = 0.22          # porosity
CT_PSI = 1.0e-5     # total compressibility, 1/psi
MU_CP = 1.0         # oil viscosity
H_FT = 40.0         # net pay
BO = 1.0            # formation volume factor
RW_FT = 0.25        # wellbore (innermost) radius
RE_FT = 20000.0     # outer radius - far enough to stay infinite-acting at T_DAYS
NR = 300            # radial cells
Q_STBD = 300.0      # producer rate (added as a NEGATIVE rate at the wellbore)
T_DAYS = 3.0        # run time
DT_DAYS = 0.02      # time step
R_OBS_FT = 75.0     # mid-field observation radius


def ei_drawdown(r, t, q, k, phi, mu, ct, h, B):
    """Line-source (exponential-integral) drawdown dP = P(r,t) - Pi.

    t is in HOURS (the 948 constant is the field-units hours form). Pass q < 0
    for a producer so the returned drawdown is NEGATIVE.
    """
    return -70.6 * q * B * mu / (k * h) * expi(-948.0 * phi * mu * ct * r ** 2 / (k * t))


def simulate_drawdown(nr=NR, r_obs_ft=R_OBS_FT):
    """Run the radial simulator and return (r_mid_ft, dP_sim_psi) at the cell
    whose center radius is closest to r_obs_ft."""
    res = Reservoir1D(nr, RW_FT, RE_FT, H_FT, K_MD, PHI, CT_PSI, MU_CP, PI_PSI, BO)
    res.add_well(0, rate_stbd=-Q_STBD)            # producer = NEGATIVE rate
    res.run(T_DAYS, DT_DAYS)
    cell = int(np.argmin(np.abs(res.r - r_obs_ft)))
    return float(res.r[cell]), float(res.P[cell] - PI_PSI)


sim_result = simulate_drawdown()
r_mid_ft = sim_result[0]
dP_sim_psi = sim_result[1]
dP_ei_psi = ei_drawdown(r_mid_ft, T_DAYS * 24.0, -Q_STBD, K_MD, PHI, MU_CP, CT_PSI, H_FT, BO)
rel_err = abs(dP_sim_psi - dP_ei_psi) / abs(dP_ei_psi)

print("r_mid (ft):", r_mid_ft)
print("dP_sim (psi):", dP_sim_psi, " dP_ei (psi):", dP_ei_psi)
print("relative error:", rel_err)

lockCopying code is a Full Access feature.

arrow_back

11.6 Fractional Flow vs Mobility Ratio

11.8 Boundary Conditions - No-Flow vs Aquifer Support

arrow_forward