多孔介质的三张面孔 — Darcy、Forchheimer 与 Brinkman

1856 年,Henry Darcy 递交给第戎市政厅的报告#

十九世纪中叶,法国第戎市需要一套新的供水系统,把干净的水送到每一户居民家中。土木工程师 Henry Darcy 提出借助城郊的砂层来过滤河水,并用一个朴素得近乎枯燥的实验为方案撑腰。在装满砂子的竖直管两端施加压差,再测出穿过管子的流量。两者之间的关系,严格地呈直线。

那个时代的流体力学,只用 Bernoulli、Euler、Navier–Stokes 这些"自由流动"的语言书写。Darcy 的一行公式,打开了它身旁那片被忽略的世界。

今天的文章就从这一行开始。当速度抬高,Forchheimer 的二次项会悄悄挤进来;当壁面靠近,Brinkman 的粘性拉普拉斯算子又会复活。我们用泵曲线和通道剖面,看看穿过孔隙的流体到底有几张面孔。

泵曲线开始弯曲的那一刻 — Forchheimer 修正#

Darcy 定律写作

\mathbf{u}_D = -\frac{K}{\mu}\,\nabla p

其中 $\mathbf{u}_D$ 是 Darcy 速度(将孔隙和固体一并计入的截面平均流速,m/s), $K$ 是渗透率(permeability,m²), $\mu$ 是动力粘度, $\nabla p$ 是压力梯度。真正在孔隙内运动的流体,速度更快, $\mathbf{u} = \mathbf{u}_D / \phi$ ,其中 $\phi$ 是孔隙率(单位体积内孔隙所占的比例)。

只要速度略微提高,情况就会变化。用颗粒直径 $d$ 与孔内速度 $u$ 定义的颗粒 Reynolds 数

\mathrm{Re}_p = \frac{\rho\,u\,d}{\mu}

一旦越过 1,压力降随流量的曲线就会从直线开始弯曲。Forchheimer(1901)用二次项写下了这道弯曲。

-\nabla p = \frac{\mu}{K}\,\mathbf{u}_D + \beta\,\rho\,\lvert\mathbf{u}_D\rvert\,\mathbf{u}_D

系数 $\beta$ (1/m)决定了"绕颗粒而行"的惯性损失何时追上粘性损失。砾石层、散热器翅片、催化反应器 — 工程上遇到的多孔介质几乎都活在这条曲线上。

壁面尚存的多孔介质 — Brinkman 项的复活#

Darcy 和 Forchheimer 都把粘性吸收进了颗粒尺度的阻力。但只要多孔区域贴上壁面或者邻接自由流动,这一假设就会破裂。壁面处的 no-slip 条件强迫速度降为零,而那一薄层之内,赤裸的粘性拉普拉斯算子必须依旧活着。

Brinkman(1947)把这一项写了回来。

-\nabla p = \frac{\mu}{K}\,\mathbf{u}_D - \mu_{\text{eff}}\,\nabla^2\mathbf{u}_D

$\mu_{\text{eff}}$ 是等效粘度,简单情形取 $\mu/\phi$ ,工程上常由实验校正。比较的尺子是 Darcy 数

\mathrm{Da} = \frac{K}{L^2}

$L$ 是通道半宽这样的宏观长度。 $\mathrm{Da} \to 0$ 时通道几乎全部退化为 Darcy 平直流; $\mathrm{Da} \to \infty$ 时粘性重新主宰,剖面靠近 Poiseuille 抛物线。两者之间,出现厚度为 $\sqrt{K}$ 的薄薄的 Brinkman 边界层。只有当这个厚度小到只跨越一两个孔时,直接用 Darcy 才算稳妥。

孔隙率加权 — 一格之内流体和固体同住#

设想一个 CFD 网格单元里同时住着砂粒(固体)和它们之间的水(流体)。守恒方程要写成两者的加权平均。只看质量守恒:

\frac{\partial (\phi\,\rho_f)}{\partial t} + \nabla \!\cdot (\rho_f\,\mathbf{u}_D) = 0

面通量取 $\mathbf{u}_D$ 时,这个形式天然守恒。如果相邻单元 $\phi = 0$ (就是一堵固体墙),面通量也必须为零,否则虚假的质量会流入根本不存在的孔隙。数值上, $\phi = 0$ 的单元通常用对角项设为 1、非对角项设为 0 来强制约束,或者把 $\phi$ 截断为 $10^{-8}$ 这种小量以稳住求解器。

能量方程更敏感。在局部热平衡假设下(流体和固体共享同一温度),需要把两者的导热系数合成为等效 $k_{\text{eff}}$ 。常用的三种混合规则是:

\begin{aligned} k_{\text{lin}} &= \phi\,k_f + (1-\phi)\,k_s \\ k_{\text{harm}}^{-1} &= \frac{\phi}{k_f} + \frac{1-\phi}{k_s} \\ k_{\text{GM}} &= k_f^{\phi}\,k_s^{1-\phi} \end{aligned}

颗粒并列排布时线性混合最贴近,串联时调和混合更合适,无规取向时则用几何平均(GM)。土壤、岩石和多孔金属在实测中通常更接近 GM 模型。

Python — 在同一条通道里比较三种模型#

带壁通道的 1 维无量纲动量方程写作

\frac{d^2 U}{dY^2} - \frac{1}{\mathrm{Da}}\,U + G = 0,\qquad U(0)=U(1)=0

$G$ 是外加的无量纲压力驱动, $\mathrm{Da}$ 是 Darcy 数, $U$ 是无量纲速度。解析解一行就能写出 $U(Y) = G\,\mathrm{Da}\,\bigl(1 - \cosh((Y-0.5)/\sqrt{\mathrm{Da}}) / \cosh(0.5/\sqrt{\mathrm{Da}})\bigr)$ ,但用数值方法走一遍,结构会更清楚。

import numpy as np
 
def brinkman_channel_velocity(Da: float, G: float = 1.0, N: int = 200):
    """带壁 1D 多孔通道的 Brinkman 稳态解。Da = K/L²。"""
    Y = np.linspace(0.0, 1.0, N)
    dY = Y[1] - Y[0]
    # U'' - U/Da + G = 0  →  三对角系统
    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   # 左壁 U=0
    A[-1, -1] = 1.0; rhs[-1] = 0.0   # 右壁 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):
    """均匀截面多孔柱内的压力降: 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 剖面: 扫描 Da, 观察中心速度
    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) 砾石层的泵曲线 (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}")

运行后会看到, $\mathrm{Da}$ 越小,中心速度越平坦,通道几乎成了一整块 Darcy 平直流。砾石层的泵曲线在 $u \sim 0.05$ m/s 之前几乎是直线, $0.2$ m/s 之后二次弯曲就明显可见。只用 Darcy 设计泵,正是在这里把压力损失估错了将近一半。

用滑块看泵曲线与靠壁剖面#

下面的模拟里可以亲手调一调参数。

Pump curve · ΔP/L vs Darcy velocity u

log₁₀(K) = -8.0 (K = 1.0e-8 m²)β = 1500 1/m (Forchheimer coefficient)

Solid cyan: ΔP/L = (μ/K) u + βρ u². Dashed gray: Darcy-only (μ/K) u. Red dashed line: u* where the two terms balance.

把 $\beta$ 调到 0,曲线就是一条直线; $\beta$ 越大,曲线越往上翘。砾石层区段( $\beta \sim 10^3$ )里,工作点一旦越过 $u^{\star}$ ,只按 Darcy 设计的泵在现实中需要将近两倍的压力。

1D channel · Brinkman steady profile U(Y) for U″ − U/Da + 1 = 0, U(0)=U(1)=0

log₁₀(Da) = -2.00 (Da = K/L² = 1.00e-2)

Drop Da below ~10⁻³: the profile flattens into a Darcy plug. Increase Da past 1: it relaxes back to Poiseuille (dashed gray). The red dashed lines mark the Brinkman boundary-layer thickness δ ≈ √Da.

把 $\mathrm{Da}$ 拉到 $10^{-4}$ ,通道中央就变成一段平直的平台。从平台边缘到壁面,厚度大约 $\sqrt{\mathrm{Da}}$ 的薄粘性层(Brinkman 边界层)仍然存活。它正是多孔介质与自由流动交界面的指纹。

多孔流的三张面孔 — 一行总结#

Darcy( $\mathrm{Re}_p < 1$ ):压力梯度 = 粘性损失。直线就是真相。
Forchheimer( $\mathrm{Re}_p \gtrsim 1$ ):绕颗粒的惯性损失叠加。泵曲线开始弯。
Brinkman( $\sqrt{\mathrm{Da}} \lesssim L$ ):一旦壁面或自由流动逼近,粘性拉普拉斯算子复活,厚度 $\sqrt{K}$ 的边界层重新出现。

同样的砂或砾石,Reynolds 数与 Darcy 数落在哪个区段,决定了它呈现三张面孔中的哪一张。开始仿真前,把这两个数估清楚,模型选择的九成答案就已经写好。