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.

As the velocity rises, Forchheimer's quadratic term sneaks in. When a wall comes close, Brinkman's viscous Laplacian comes back to life. The log10(K) slider and the β slider move the red dashed line, the speed where the viscous and inertial terms meet.

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.

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.

Set $(\mu/K)\,u$ equal to $\beta\rho u^{2}$ and the crossing speed is one division, $u^{\star}=\mu/(K\beta\rho)$. The gravel-bed numbers in the snippet are $\mu=10^{-3}$, $K=10^{-8}$, $\beta=10^{3}$, $\rho=1000$, so the denominator is $10^{-2}$ and $u^{\star}=0.1$ m/s. Double only $\beta$ and $u^{\star}$ drops to $0.05$ m/s, which puts the same $0.1$ m/s duty point on the inertial side. Multiply $K$ by ten and $u^{\star}$ moves with it.

At $u=0.05$ m/s the inertia-to-viscosity ratio is $u/u^{\star}=1/2$, and the bend is still mild. At $u=0.2$ m/s the ratio is 2, so the total drop is three times the viscous term and a Darcy-only reading keeps one third. The half-drop speed is $0.1$ m/s, not $0.2$ m/s.

A grain-diameter Reynolds number calls the quadratic term in near 1, the same split as how a Reynolds number divides the flow. Only the length changes, from a pipe diameter to a grain diameter $d$.

With $G=1$, mid-channel speed is $\mathrm{Da}$ times $(1-1/\cosh(0.5/\sqrt{\mathrm{Da}}))$. At $\mathrm{Da}=10^{-4}$, $0.5/\sqrt{\mathrm{Da}}=50$, so $\cosh$ runs away and the center speed sticks to $\mathrm{Da}$ itself, $10^{-4}$. At $\mathrm{Da}=1$, $0.5/\sqrt{\mathrm{Da}}=0.5$, the denominator only steps a little off 1, and the plug no longer spans the channel.

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 Darcy-only reading at that speed keeps about one third of the true drop. The half-drop crossing sits at 0.1 m/s.

Where lowering log10(Da) leaves only a plug#

Move the log10(Da) slider toward −4 and the mid-channel speed flattens to one value.

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.

Where a straight pump curve halves the gravel-bed drop#

On the side $\mathrm{Re}_p < 1$ the straight line is the right face, and once $\mathrm{Re}_p \gtrsim 1$ the quadratic term joins. Leave $u^{\star}=\mu/(K\beta\rho)$ off the duty sheet and a straight Darcy curve still gets picked after the point has already moved to the right. On the gravel-bed numbers in the snippet that boundary is 0.1 m/s. At 0.2 m/s the straight line points at one third of the true drop.

A wall that puts the channel in $\sqrt{\mathrm{Da}} \lesssim L$ splits the choice again. Drop the viscous Laplacian before measuring the layer of thickness $\sqrt{K}$, and the wall velocity never falls to zero. If the pore $\mathrm{Re}_p$ and the pipe transition look like the same threshold, open the Reynolds number where a water jet breaks up and compare the length that each one uses.