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

先把 u* = μ/(Kβρ) 放在工作点旁边,再决定直线泵曲线能不能用。

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

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

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

速度一抬高,Forchheimer 的二次项就挤进来;壁面一靠近,Brinkman 的粘性拉普拉斯算子又活过来。log₁₀(K) 滑块和 β 滑块会左右推动红虚线 u*,也就是粘性项与惯性项相抵的速度。

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 设计的泵在现实中需要将近两倍的压力。

泵曲线开始弯曲的那一刻 — 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)决定了"绕颗粒而行"的惯性损失何时追上粘性损失。砾石层、散热器翅片、催化反应器 — 工程上遇到的多孔介质几乎都活在这条曲线上。

粘性项 $(\mu/K)\,u$ 与惯性项 $\beta\rho u^{2}$ 相等的速度,一次除法就落出来。 $u^{\star}=\mu/(K\beta\rho)$ 。代码里的砾石层数值 $\mu=10^{-3}$ 、 $K=10^{-8}$ 、 $\beta=10^{3}$ 、 $\rho=1000$ 代进去,分母是 $10^{-2}$ ,于是 $u^{\star}=0.1$ m/s。只把 $\beta$ 加倍, $u^{\star}$ 被拉到 $0.05$ m/s。同一个 $0.1$ m/s 的工作点从此站到惯性一侧。把 $K$ 放大十倍, $u^{\star}$ 也放大十倍,直线又盖住泵曲线上更长的一段。

$u=0.05$ m/s 时,惯性与粘性之比是 $u/u^{\star}=1/2$ ,曲线还缓。 $u=0.2$ m/s 时这个比是 2。总压降是粘性项的 3 倍,只读 Darcy 直线就只看到压降的 1/3。两项相等、直线正好指向真实值一半的位置,不是 $0.2$ m/s,而是 $0.1$ m/s。

用粒径定义的 Reynolds 数在 1 附近叫出二次项,和Reynolds 数把流动世界切开的方式是同一件事。变的只是长度,从管径换成粒径 $d$ 。管内在转捩之前可以长时间留在层流,粒与粒之间这个数一过 1,泵曲线就弯。

$G=1$ 时,通道中央速度是 $\mathrm{Da}$ 乘上 $(1-1/\cosh(0.5/\sqrt{\mathrm{Da}}))$ 。 $\mathrm{Da}=10^{-4}$ 时 $0.5/\sqrt{\mathrm{Da}}=50$ , $\cosh$ 实际上发散,中央速度贴在 $\mathrm{Da}$ 自己身上,就是 $10^{-4}$ 。 $\mathrm{Da}=1$ 时 $0.5/\sqrt{\mathrm{Da}}=0.5$ ,分母只离开 1 一点点,平推盖不住整个通道。

壁面尚存的多孔介质 — 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 直线来读,这个速度上的压降只剩全部的 1/3。刚好一半的交点在 0.1 m/s。

把 log₁₀(Da) 拉低后通道里只剩平推的地方#

把 log₁₀(Da) 滑块推向 −4,通道中央的速度会摊成同一个值。

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 边界层)仍然存活。它正是多孔介质与自由流动交界面的指纹。

直线泵曲线把砾石层压降砍成一半的位置

$\mathrm{Re}_p < 1$ 时直线是对的, $\mathrm{Re}_p \gtrsim 1$ 时二次项加进来。不把 $u^{\star}=\mu/(K\beta\rho)$ 放在工作点旁边,点已经到了右侧,仍会选中 Darcy 直线。按代码里的砾石层数字,这条边界是 0.1 m/s。到了 0.2 m/s,直线指出的压降只有全部的 1/3。

壁面贴上来、通道落在 $\sqrt{\mathrm{Da}} \lesssim L$ 时,选择再分一次。还没量厚度 $\sqrt{K}$ 的那一层,就把粘性拉普拉斯拿掉,壁面速度不会落到 0。孔隙 $\mathrm{Re}_p$ 和圆管转捩看起来像同一道门槛时,打开水流开始散开的雷诺临界,把两套长度对照一下。