Three Faces of Porous Flow — Darcy, Forchheimer, and Brinkman

1856, Henry Darcy's report to the city of Dijon#

In the middle of the nineteenth century, Dijon needed a new system to deliver clean water to its citizens. The civil engineer Henry Darcy proposed filtering river water through a sandbank outside the city, and he ran a stubbornly simple experiment to back up the plan. Drive a pressure drop across a vertical column packed with sand, and measure the flow that comes out the bottom. The relationship turned out to be strictly linear.

Fluid mechanics up to that point spoke only the language of free flow — Bernoulli, Euler, Navier–Stokes. Darcy's single line opened up another world that had been sitting quietly next to it.

Today's post starts at that line. As the velocity rises, Forchheimer's quadratic term sneaks in. When a wall comes close, Brinkman's viscous Laplacian comes back to life. We will walk through the three faces a fluid puts on as it pushes through pores, using a pump curve and a channel profile.

When the pump curve starts to bend — the Forchheimer correction#

Darcy's law reads

where $\mathbf{u}_D$ is the Darcy velocity (the cross-section-averaged flow including both pores and solid, m/s), $K$ is the permeability (m²), $\mu$ is the dynamic viscosity, and $\nabla p$ is the pressure gradient. The fluid inside the actual pore moves faster: $\mathbf{u} = \mathbf{u}_D / \phi$, with $\phi$ the porosity (the volume fraction of pores per unit cell).

A small bump in velocity changes the picture. Define a pore-scale Reynolds number from the grain diameter $d$ and pore velocity $u$,

Once $\mathrm{Re}_p$ climbs past one, the pressure-drop curve begins to bend away from a straight line. Forchheimer (1901) wrote down that bend as a quadratic term.

The coefficient $\beta$ (1/m) sets the speed at which inertial losses — fluid having to swerve around each grain — overtake viscous losses. Gravel beds, heat-sink fins, and catalyst reactors all live on this curve in practice.

When the wall is still alive — Brinkman's return of viscosity#

Darcy and Forchheimer both absorb viscosity into a grain-scale loss term. That assumption breaks the moment the porous region touches a wall or an open fluid layer. No-slip forces the velocity to zero at the wall, and inside that thin film the bare viscous Laplacian has to be alive again.

Brinkman (1947) put that term back in.

$\mu_{\text{eff}}$ is an effective viscosity, often taken as $\mu/\phi$ or calibrated from experiment. The ruler of choice is the Darcy number

with $L$ a macroscopic length such as the channel half-width. When $\mathrm{Da} \to 0$ the channel is filled by a Darcy plug. When $\mathrm{Da} \to \infty$ viscosity comes back and the profile relaxes to Poiseuille parabolas. In between, a thin Brinkman boundary layer of thickness $\sqrt{K}$ shows up — and it has to be small enough to span only one or two pores before plain Darcy is safe.

Pore-fraction averaging — when fluid and solid share a single cell#

Imagine a CFD cell that contains both grains of sand (solid) and water in between (fluid). The conservation laws are written as pore-fraction-weighted averages. Mass conservation reads

and that form is naturally conservative if face fluxes use $\mathbf{u}_D$. If a neighbor cell has $\phi = 0$ (a solid wall), the face flux must vanish too — otherwise spurious mass leaks into a pore that does not exist. Numerically, $\phi = 0$ cells are stabilized by forcing the diagonal to one with zero off-diagonals, or by clipping $\phi$ to a tiny value like $10^{-8}$.

The energy equation is trickier. Under the local thermal equilibrium assumption (fluid and solid share the same temperature), we have to combine the conductivities into an effective $k_{\text{eff}}$. Three mixing rules show up most often.

Linear matches grains lined up in parallel, harmonic matches grains stacked in series, and the geometric mean (GM) matches random orientations best. Soils, rocks, and porous metals usually fit the GM model.

Python — three models on a single channel#

The 1-D non-dimensional momentum balance for a walled porous channel is

where $G$ is the non-dimensional pressure forcing, $\mathrm{Da}$ is the Darcy number, and $U$ is the non-dimensional velocity. A closed-form solution exists, $U(Y) = G\,\mathrm{Da}\,\bigl(1 - \cosh((Y-0.5)/\sqrt{\mathrm{Da}}) / \cosh(0.5/\sqrt{\mathrm{Da}})\bigr)$, but a numerical solve makes the structure clearer.

import numpy as np
 
def brinkman_channel_velocity(Da: float, G: float = 1.0, N: int = 200):
    """Steady 1-D Brinkman profile in a walled porous channel. Da = K/L²."""
    Y = np.linspace(0.0, 1.0, N)
    dY = Y[1] - Y[0]
    # U'' - U/Da + G = 0  →  tridiagonal system
    main = -2.0 / dY**2 - 1.0 / Da
    off  =  1.0 / dY**2
    A = np.zeros((N, N))
    rhs = np.full(N, -G)
    A[0, 0] = 1.0;   rhs[0]  = 0.0   # left wall U=0
    A[-1, -1] = 1.0; rhs[-1] = 0.0   # right wall U=0
    for i in range(1, N - 1):
        A[i, i - 1] = off
        A[i, i]     = main
        A[i, i + 1] = off
    return Y, np.linalg.solve(A, rhs)
 
 
def darcy_forchheimer_pressure_drop(u: np.ndarray, mu: float, K: float,
                                    beta: float, rho: float):
    """Pressure drop in a uniform porous column: Darcy + Forchheimer."""
    visc_loss    = (mu / K) * u
    inertia_loss = beta * rho * u * np.abs(u)
    return visc_loss + inertia_loss
 
 
if __name__ == "__main__":
    # 1) Brinkman profile: sweep Da, watch the mid-channel velocity
    for Da in [1e-4, 1e-2, 1e0]:
        Y, U = brinkman_channel_velocity(Da)
        print(f"Da = {Da:7.0e}   U_mid = {U[len(U)//2]:.4f}")
 
    # 2) Gravel-bed pump curve (K=1e-8 m², β=1e3 1/m)
    mu, K, beta, rho = 1e-3, 1e-8, 1.0e3, 1000.0
    u = np.linspace(0.0, 0.5, 6)
    dp = darcy_forchheimer_pressure_drop(u, mu, K, beta, rho)
    print("\n u [m/s]    ΔP/L [Pa/m]")
    for ui, dpi in zip(u, dp):
        print(f"  {ui:5.2f}      {dpi:.3e}")

Run it and the mid-channel velocity flattens as $\mathrm{Da}$ shrinks — most of the channel is now a Darcy plug. The pump curve stays nearly linear up to $u \approx 0.05$ m/s, then bends visibly above $0.2$ m/s. A pump designed on Darcy alone would mis-size the required pressure by almost a factor of two right there.

Sliders for the bending pump curve and the wall-bound profile#

Play with the parameters in the simulations below.

Set $\beta$ to zero and the curve is a straight line. Push $\beta$ up and the curve bends upward. Once your design point in the gravel-bed regime ($\beta \sim 10^3$) crosses $u^{\star}$, a Darcy-only design under-predicts the actual pressure loss by nearly half.

Drop $\mathrm{Da}$ to $10^{-4}$: the center of the channel becomes a flat plug. From the plug edge to the wall, a thin viscous layer (the Brinkman boundary layer) of thickness $\sqrt{\mathrm{Da}}$ survives. That layer is the fingerprint of the interface between the porous medium and the free flow.

Three faces of porous flow — a line each#

• Darcy ($\mathrm{Re}_p < 1$): pressure gradient balances viscous loss. The straight line is the truth.
• Forchheimer ($\mathrm{Re}_p \gtrsim 1$): inertial losses around grains pile on. The pump curve bends.
• Brinkman ($\sqrt{\mathrm{Da}} \lesssim L$): a wall or an open flow lives nearby, so the viscous Laplacian wakes back up. A boundary layer of thickness $\sqrt{K}$ reappears.

The same sand or gravel switches between these three faces depending on where Reynolds and Darcy numbers land. Pinning those two numbers down before any simulation already settles nine-tenths of the model-choice question.