다공성 매질의 세 얼굴 — 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을 넘기 시작하면, 압력 강하 vs 유량 곡선이 선형에서 휘어진다. Forchheimer(1901)는 이 휘어짐을 이차항으로 적었다.

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

$\beta$ 는 Forchheimer 계수(1/m). 알갱이를 돌아 나가는 관성 손실이 점성 손실을 따라잡기 시작하는 속도를 결정한다. 자갈층, 발열소자의 핀, 촉매 반응기 — 실무에서 만나는 다공질은 거의 이 곡선 위에 산다.

벽이 살아 있는 다공질 — Brinkman 항의 부활#

Darcy도 Forchheimer도 점성을 입자 단위로 흡수해 버린다. 하지만 다공질이 벽이나 자유 유동과 맞닿는 순간 그 가정이 깨진다. 벽에서는 no-slip 때문에 속도가 0으로 떨어져야 하고, 그 얇은 층 안에서는 점성이 직접 작용하는 라플라시안 항이 살아 있어야 한다.

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이면(= 고체 벽이면) face flux도 0이 되어야 한다. 그러지 않으면 0인 공극으로 가짜 질량이 새 들어간다. 수치적으로 $\phi=0$ 인 셀은 대각항을 1, 비대각을 0으로 강제하거나 $1\mathrm{e}{-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 — 세 모델을 한 채널 위에서 비교#

벽이 있는 채널의 1D 무차원 운동량 방정식은

\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  →  3중대각 시스템
    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의 두 값이 어느 영역에 떨어지는지에 따라 셋 중 하나의 얼굴로 바뀐다. 시뮬레이션을 시작하기 전에 이 두 숫자만 먼저 잡아도 — 어떤 모델을 써야 할지 9할은 정해진다.