多孔質媒体の三つの顔 — Darcy・Forchheimer・Brinkman

1856年、Henry Darcy がディジョン市に提出した報告書#

19世紀半ば、フランスのディジョン市は市民に清浄な水道を届ける新しい仕組みを必要としていました。土木技術者の 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$ は空隙率(porosity、単位セル内で空隙が占める体積比)です。

速度が少し上がるだけで構図は変わります。粒径 $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 格子の 1 セルに砂粒(固体)とその間の水(流体)が同居していると考えてみましょう。保存方程式は両者の加重平均で書かれます。質量保存だけを見ると

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

面フラックスに $\mathbf{u}_D$ を使えば自然に保存形になります。隣セルの $\phi$ が 0(= 固体壁)なら、面フラックスも 0 になっていなければなりません。さもなければ存在しない空隙に偽の質量が流れ込みます。数値的には $\phi = 0$ のセルは対角項を 1、非対角を 0 にして強制するか、 $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 数の二つがどの領域に落ちるかで三つの顔のどれか一つに切り替わります。シミュレーションを始める前にこの二つだけ押さえれば、どのモデルを使うかの答えは九割方決まります。