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

運転点を決める前に u* = μ/(Kβρ) を置くと、直線モデルが差圧を削る側にいるかが先に分かります。

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

19世紀半ば、フランスのディジョン市は市民に清浄な水道を届ける新しい仕組みを必要としていました。土木技術者の 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$ は空隙率(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)は、粒の周りを迂回することで生まれる慣性損失が、粘性損失に追いつき始める速度を決めます。砂利層、ヒートシンクのフィン、触媒反応器 — 現場で出会う多孔質はほぼこの曲線の上で生きています。

粘性項 $(\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 格子の 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 の直線だけを読むと、この速度での差圧は全体の 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$ と管内の遷移が同じ閾値に見えるときは、水筋が乱れるレイノルズ臨界を開いて、長さの物差しの違いを並べると判断が戻ります。