Ergun packed-bed pressure drop (Darcy -> Forchheimer)¶

Validation of the full volumetric porous-medium momentum sink - the linear Darcy coefficient porous_alpha and the quadratic Forchheimer coefficient porous_beta together (models.volumetric_regions, issues #358 / #743) - against the Ergun (1952) packed-bed correlation.

On flagged region nodes the Guo source term gains

\[F_\alpha = -\alpha\,u_\alpha - \beta\,|u|\,u_\alpha.\]

In a fully-developed, gradient-free porous-filled box driven by a constant streamwise body force \(G\), the steady momentum balance is purely local:

\[G = \alpha\,u + \beta\,u^2,\]

which is the Ergun pressure-gradient law \(dP/dx = \rho G\) (\(\rho = 1\)) with

\[\alpha = 150\,\frac{(1-\varepsilon)^2}{\varepsilon^3}\,\frac{\nu}{d_p^2}, \qquad \beta = 1.75\,\frac{(1-\varepsilon)}{\varepsilon^3}\,\frac{1}{d_p}.\]

The case (validation/porous_media/03_ergun_pressure_drop/) has two parts:

``ergunGate`` (03_ergun_pressure_drop.nassu.yaml) - a triply-periodic 32^3 porous box driven at five body forces (!unroll), sweeping the superficial velocity across the linear -> quadratic regime transition. The pair \((u, dP/dx)\) must land on the analytic Ergun curve within ~3%.
``ergunRefinedSlab`` (03.1_ergun_pressure_drop_multiblock.nassu.yaml) - the same box driven at the transition point, with an inner slab refined to level 1 so the porous region straddles two refinement interfaces. Guards the issue #704 fix (porous force restored to the block-border reconstruction): the steady velocity must be spatially uniform with no sawtooth at the interfaces.

Results pending GPU execution. This notebook is a skeleton: the analytic Ergun curve and all comparison logic (and every plot) are wired against the parsed config, but the simulation outputs do not exist until the case is run on GPU. Run uv run nassu run validation/porous_media/03_ergun_pressure_drop/03_ergun_pressure_drop.nassu.yaml (and the 03.1 multiblock file) from the repo root, then execute this notebook. The data-loading cells degrade gracefully (skip + warn) while results are absent; the analytic curves render regardless so the figures are populated the moment the run exists.

Setup: load the configs and recover the bed parameters¶

[ ]:

import numpy as np
import matplotlib.pyplot as plt
from nassu.cfg.model import ConfigScheme
import nassu.viz as common

common.use_style()

CS2 = 1.0 / 3.0
PROJECT_ROOT = common.find_project_root()

CASE_DIR = "validation/porous_media/03_ergun_pressure_drop"
GATE_FILE = f"{CASE_DIR}/03_ergun_pressure_drop.nassu.yaml"
MB_FILE = f"{CASE_DIR}/03.1_ergun_pressure_drop_multiblock.nassu.yaml"

gate_cfgs = ConfigScheme.sim_cfgs_from_file_dct(GATE_FILE)
mb_cfgs = ConfigScheme.sim_cfgs_from_file_dct(MB_FILE)

# The sweep is unrolled into sim_id 0..4 (Darcy-dominated -> Forchheimer-dominated).
N_SWEEP = len([k for k in gate_cfgs if k[0] == "ergunGate"])
gate_list = [gate_cfgs["ergunGate", i] for i in range(N_SWEEP)]
mb_cfg = mb_cfgs["ergunRefinedSlab", 0]

# Recover the bed parameters from the first parsed config (all share alpha/beta).
c0 = gate_list[0]
nu = CS2 * (c0.models.LBM.tau - 0.5)
vr = c0.models.volumetric_regions[0]
ALPHA = vr.porous_alpha
BETA = vr.porous_beta
U_STAR = ALPHA / BETA  # transition velocity (Fo = 1)

# Drive forces and their target velocities (solve G = alpha u + beta u^2).
G_sweep = np.array([float(c.models.LBM.F.x) for c in gate_list])


def ergun_velocity(G, alpha=ALPHA, beta=BETA):
    ''' Invert G = alpha u + beta u^2 for the positive root u(G). '''
    return (-alpha + np.sqrt(alpha**2 + 4.0 * beta * G)) / (2.0 * beta)


U_target = ergun_velocity(G_sweep)
Fo_target = BETA * U_target / ALPHA

print(f"nu={nu:.4g}  alpha={ALPHA:.4g}  beta={BETA:.4g}  u*={U_STAR:.4g}")
for i, (g, u, fo) in enumerate(zip(G_sweep, U_target, Fo_target)):
    print(f"  sim_id={i}: G={g:.5g}  u_target={u:.5f}  Fo={fo:.3f}  Ma={u*np.sqrt(3):.4f}")

Analytic Ergun curve¶

The pressure gradient as a function of superficial velocity splits into a linear Darcy term and a quadratic Forchheimer term:

\[\frac{dP}{dx}(u) = \underbrace{\alpha\,u}_{\text{Darcy}} + \underbrace{\beta\,u^2}_{\text{Forchheimer}}.\]

The total is a straight line in the Darcy limit (\(Fo \ll 1\)) that bends upward as the quadratic term takes over (\(Fo \gg 1\)); the two components cross at the transition velocity \(u^* = \alpha/\beta\).

[ ]:

def ergun_dpdx(u, alpha=ALPHA, beta=BETA):
    ''' Analytic Ergun pressure gradient dP/dx(u) = alpha u + beta u^2 (total). '''
    return alpha * u + beta * u**2


def ergun_darcy(u, alpha=ALPHA):
    ''' Linear Darcy component alpha u (viscous term). '''
    return alpha * u


def ergun_forchheimer(u, beta=BETA):
    ''' Quadratic Forchheimer component beta u^2 (inertial term). '''
    return beta * u**2


u_fine = np.linspace(0.0, 1.05 * U_target.max(), 200)

Component decomposition of the Ergun curve¶

Plot the total Ergun curve together with its linear (Darcy) and quadratic (Forchheimer) components, and the five sweep points marked on the velocity axis. This figure renders from the analytic coefficients alone, so it is populated even before the run exists.

[ ]:

fig, ax = common.fig_single()
ax.plot(u_fine, ergun_dpdx(u_fine), **common.markers.exp_line(),
        label=r"total $\alpha u + \beta u^2$")
ax.plot(u_fine, ergun_darcy(u_fine), color=common.colors.blue, linestyle=":",
        label=r"Darcy (linear) $\alpha u$")
ax.plot(u_fine, ergun_forchheimer(u_fine), color=common.colors.green, linestyle="-.",
        label=r"Forchheimer (quadratic) $\beta u^2$")
ax.axvline(U_STAR, color=common.colors.refline, linestyle="--", linewidth=0.8)
ax.text(U_STAR, ergun_dpdx(U_STAR), r"  $u^* = \alpha/\beta$  ($Fo=1$)",
        va="bottom", ha="left", color=common.colors.refline)
for u in U_target:
    ax.axvline(u, color="grey", linestyle="-", linewidth=0.4, alpha=0.5)
ax.set_xlabel(r"superficial velocity $u$ [lattice]")
ax.set_ylabel(r"$dP/dx$ [lattice]")
ax.set_title("Ergun curve: Darcy + Forchheimer components")
ax.legend()
fig.tight_layout()
plt.show()

Measured pressure drop: steady superficial velocity per drive¶

For each sim_id, the steady superficial velocity is the domain mean of the streamwise velocity ux at the last exported snapshot. By construction the applied drive G equals the steady pressure gradient dP/dx, so the measured pair is (u_measured, G); it must lie on the analytic Ergun curve.

[ ]:

def steady_mean_ux(sim_cfg):
    ''' Domain-mean streamwise velocity at the last exported volume snapshot. '''
    src = common.FieldSource.from_cfg(sim_cfg, project_root=PROJECT_ROOT)
    step_end = src.steps[-1]
    arrays, _ = src.read_arrays(step_end, ["ux"])
    return float(np.nanmean(arrays["ux"].astype(np.float64)))


u_measured = np.full(N_SWEEP, np.nan)
for i, cfg in enumerate(gate_list):
    try:
        u_measured[i] = steady_mean_ux(cfg)
    except (FileNotFoundError, IndexError) as exc:
        print(f"sim_id={i}: results missing ({exc}) - run on GPU first.")

have_results = np.isfinite(u_measured).any()
if have_results:
    dpdx_pred = ergun_dpdx(u_measured)
    rel_err = np.abs(G_sweep - dpdx_pred) / G_sweep
    print("Ergun-curve agreement (measured u -> predicted dP/dx vs applied G):")
    for i in range(N_SWEEP):
        if np.isfinite(u_measured[i]):
            status = "PASS" if rel_err[i] < 0.03 else "FAIL"
            print(f"  sim_id={i} Fo={Fo_target[i]:.2f}: u={u_measured[i]:.5f}  "
                  f"rel_err={rel_err[i]*100:.2f}%  [{status}]")
else:
    rel_err = np.full(N_SWEEP, np.nan)
    print("No results yet; the plots below show the analytic Ergun curve only.")

Primary validation: dP/dx vs u against the Ergun curve¶

[ ]:

fig, ax = common.fig_single()
ax.plot(u_fine, ergun_dpdx(u_fine), **common.markers.exp_line(),
        label="Ergun (Darcy + Forchheimer)")
ax.plot(u_fine, ergun_darcy(u_fine), color=common.colors.blue, linestyle=":",
        label=r"Darcy-only limit $\alpha u$")
ax.plot(u_fine, ergun_forchheimer(u_fine), color=common.colors.green, linestyle="-.",
        label=r"Forchheimer-only $\beta u^2$")
ax.axvline(U_STAR, color=common.colors.refline, linestyle="--", linewidth=0.8)
ax.text(U_STAR, ergun_dpdx(U_STAR), r"  $u^* = \alpha/\beta$  ($Fo=1$)",
        va="bottom", ha="left", color=common.colors.refline)

if have_results:
    m = np.isfinite(u_measured)
    ax.plot(u_measured[m], G_sweep[m], **common.markers.sim("o"), label="Nassu")

ax.set_xlabel(r"superficial velocity $u$ [lattice]")
ax.set_ylabel(r"$dP/dx = \alpha u + \beta u^2$ [lattice]")
ax.set_title("Ergun packed-bed pressure drop: Darcy -> Forchheimer")
ax.legend()
fig.tight_layout()
plt.show()

Relative error of the measured pressure drop vs Ergun¶

The companion plot for the agreement table above: the per-sim_id relative error \(|G - dP/dx(u_{\text{meas}})| / G\) against the Forchheimer number, with the 3% pass band shaded. Renders once the run produces u_measured.

[ ]:

fig, ax = common.fig_single()
ax.axhspan(0.0, 3.0, color=common.colors.refline, alpha=0.12, label="3% pass band")
ax.axhline(3.0, color=common.colors.refline, linestyle="--", linewidth=0.8)
if have_results:
    m = np.isfinite(u_measured)
    ax.plot(Fo_target[m], rel_err[m] * 100.0, **common.markers.sim("o"),
            label="Nassu rel. error")
else:
    ax.text(0.5, 0.5, "results pending GPU execution", transform=ax.transAxes,
            ha="center", va="center", color="grey")
ax.axvline(1.0, color="grey", linestyle=":", linewidth=0.8)
ax.text(1.0, ax.get_ylim()[1], r" $Fo=1$", va="top", ha="left", color="grey")
ax.set_xlabel(r"Forchheimer number $Fo = \beta u/\alpha$")
ax.set_ylabel("relative error of $dP/dx$ [%]")
ax.set_title("Ergun-curve agreement across the sweep")
ax.legend()
fig.tight_layout()
plt.show()

Regime decomposition¶

The fraction of the pressure gradient carried by the quadratic Forchheimer term, \(\beta u^2 / (\alpha u + \beta u^2) = Fo/(1+Fo)\), crosses 0.5 at the transition velocity. This confirms the sweep spans both regimes.

[ ]:

quad_frac = Fo_target / (1.0 + Fo_target)
print("Forchheimer fraction of dP/dx across the sweep:")
for i in range(N_SWEEP):
    print(f"  sim_id={i}: Fo={Fo_target[i]:.3f}  quadratic fraction={quad_frac[i]*100:.1f}%")

fig, ax = common.fig_single()
fo_fine = np.linspace(0.05, 1.05 * Fo_target.max(), 200)
ax.plot(fo_fine, fo_fine / (1.0 + fo_fine), **common.markers.exp_line(),
        label="quadratic fraction $Fo/(1+Fo)$")
ax.plot(Fo_target, quad_frac, **common.markers.sim("s"), label="sweep points")
ax.axvline(1.0, color=common.colors.refline, linestyle="--", linewidth=0.8)
ax.set_xlabel(r"Forchheimer number $Fo = \beta u/\alpha$")
ax.set_ylabel("quadratic fraction of $dP/dx$")
ax.set_title("Regime coverage of the sweep")
ax.legend()
fig.tight_layout()
plt.show()

Multiblock interface guard (`ergunRefinedSlab`)¶

The same Ergun box driven at the transition point, with an inner slab in x refined to level 1. The porous region straddles both C2F / F2C interfaces. Because the Ergun balance is purely local and level-independent, the converged streamwise velocity must be spatially uniform - no sawtooth at the interfaces (x = 8 and x = 24), and matching the sim_id 2 gate value within ~1%. A residual jump flags the issue #704 border-force bug.

[ ]:

def streamwise_profile(sim_cfg):
    ''' Return (x, u_of_x): the y-z-averaged streamwise velocity vs x, last step. '''
    src = common.FieldSource.from_cfg(sim_cfg, project_root=PROJECT_ROOT)
    step_end = src.steps[-1]
    arrays, _ = src.read_arrays(step_end, ["ux"])
    ux = arrays["ux"].astype(np.float64)  # indexed (x, y, z)
    u_of_x = np.nanmean(ux, axis=(1, 2))
    x = np.arange(u_of_x.size)
    return x, u_of_x


# Reference: the transition-point gate value (sim_id 2, u_target ~ 0.02).
u_ref = float(U_target[2]) if N_SWEEP > 2 else U_STAR

try:
    x_mb, u_mb = streamwise_profile(mb_cfg)
    u_mean = float(np.nanmean(u_mb))
    sawtooth = float(np.nanmax(u_mb) - np.nanmin(u_mb)) / u_mean
    rel_to_gate = abs(u_mean - u_ref) / u_ref
    print(f"multiblock mean u = {u_mean:.5f}  (gate sim_id 2 target ~ {u_ref:.5f}, "
          f"rel {rel_to_gate*100:.2f}%)")
    print(f"streamwise non-uniformity (max-min)/mean = {sawtooth*100:.3f}%  "
          f"[{'PASS' if sawtooth < 0.01 else 'FAIL (#704 sawtooth?)'}]")
except (FileNotFoundError, IndexError) as exc:
    x_mb = u_mb = None
    print(f"Multiblock results missing ({exc}) - run 03.1 on GPU first.")

Summary¶

Success criteria (checked once the case is run on GPU):

``ergunGate`` (primary): the measured (u, dP/dx) pairs lie on the analytic Ergun curve \(dP/dx = \alpha u + \beta u^2\) within ~3% across the whole sweep (relative-error plot inside the 3% band), with the Forchheimer quadratic term resolved at the high-Fo end (quadratic fraction Fo/(1+Fo) crossing 0.5 at u*).
``ergunRefinedSlab`` (multiblock guard): the converged streamwise velocity is spatially uniform (flat profile, no sawtooth at x = 8, x = 24) and matches the gate transition-point value within ~1%, confirming the issue #704 border-force fix.

Scope: validates the combined linear Darcy + quadratic Forchheimer momentum sink and its survival across a refinement interface. The porosity-modified equilibrium (Guo-Zhao REV model) is out of scope.

Reference: Ergun, S. (1952), Chem. Eng. Prog. 48(2):89-94.

Version¶

[ ]:

sim_info = gate_list[0].output.read_info()
print("Version:", sim_info["version"])
print("Commit hash:", sim_info["commit"])

Configuration¶

[ ]:

from IPython.display import Code

Code(filename=GATE_FILE)

[ ]:

Code(filename=MB_FILE)