다공성 매질의 세 얼굴 — Darcy·Forchheimer·Brinkman

1856년, Henry Darcy가 디종 시청에 낸 보고서#

19세기 중반, 프랑스 디종은 시민에게 깨끗한 수돗물을 공급할 새 시스템이 필요했다. 토목 기술자 Henry Darcy는 도시 외곽의 모래층으로 물을 거르는 안을 제시했다. 안을 뒷받침하려고 그는 단순한 실험을 끈질기게 반복했다. 모래를 채운 수직 관에 압력차를 걸고 통과 유량을 잰다. 둘 사이가 정확히 선형이었다.

당시까지의 유체역학은 Bernoulli, Euler, Navier–Stokes 같은 자유 흐름의 언어로만 적혀 있었다. Darcy의 한 줄은 그 옆에 빠져 있던 또 하나의 세계를 열었다.

속도가 커지면 Forchheimer의 이차항이 끼어든다. 벽이 가까워지면 Brinkman의 점성이 되살아난다. log₁₀(K) 슬라이더와 β 슬라이더가, 점성항과 관성항이 만나는 빨간 점선 u*를 좌우로 민다.

$\beta$를 0으로 두면 곡선은 한 줄의 직선이다. $\beta$를 키울수록 위로 휘어 올라간다. 자갈층 영역($\beta \sim 10^3$)에서 동작 속도가 $u^\star$를 넘는 순간, Darcy만으로 설계한 펌프는 실제로는 거의 두 배의 압력을 요구한다.

펌프 곡선이 휘어지는 순간 — Forchheimer 보정#

Darcy의 법칙은 이렇게 쓰인다.

여기서 $\mathbf{u}_D$는 Darcy 속도(공극과 고체를 함께 포함한 단면 평균 유속, m/s), $K$는 투과도(permeability, m²), $\mu$는 동점성도, $\nabla p$는 압력 구배다. 실제 공극 안 유체가 보는 속도는 $\mathbf{u} = \mathbf{u}_D / \phi$로 더 빠르다. $\phi$는 공극률(porosity, 단위 셀에서 유체가 차지하는 부피 비율).

속도가 조금만 올라가면 그림이 바뀐다. 입자 직경 $d$, 공극 속도 $u$로 정의되는 입자 Reynolds 수

가 1을 넘기 시작하면, 압력 강하 vs 유량 곡선이 선형에서 휘어진다. Forchheimer(1901)는 이 휘어짐을 이차항으로 적었다.

$\beta$는 Forchheimer 계수(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 때문에 속도가 0으로 떨어져야 하고, 그 얇은 층 안에서는 점성이 직접 작용하는 라플라시안 항이 살아 있어야 한다.

Brinkman(1947)이 다시 적어 넣은 것이 그 항이다.

$\mu_{\text{eff}}$는 유효 점성으로, 단순한 경우 $\mu/\phi$를 쓰거나 실험으로 보정한다. 비교 잣대는 Darcy 수다.

$L$은 채널 폭 같은 거시 길이 스케일이다. $\mathrm{Da} \to 0$이면 채널 거의 전 영역이 Darcy 플러그가 된다. $\mathrm{Da} \to \infty$이면 점성이 다시 살아나 Poiseuille 포물선에 가까워진다. 그 사이에 두께 $\sqrt{K}$의 얇은 Brinkman 경계층이 형성된다. 이 길이가 공극 한두 개를 가로지를 만큼 작아야 Darcy를 그대로 써도 안전하다.

공극률 가중평균 — 한 셀에 유체와 고체가 같이 살 때#

CFD 격자 한 셀 안에 모래 알(고체)과 그 사이 물(유체)이 같이 있다고 하자. 보존 방정식은 둘의 가중평균으로 적힌다. 질량 보존만 보면

면 플럭스에 $\mathbf{u}_D$를 그대로 쓰면 자연스럽게 보존형이 된다. 한쪽 셀의 $\phi$가 0이면(= 고체 벽이면) face flux도 0이 되어야 한다. 그러지 않으면 0인 공극으로 가짜 질량이 새 들어간다. 수치적으로 $\phi=0$인 셀은 대각항을 1, 비대각을 0으로 강제하거나 $1\mathrm{e}{-8}$ 같은 작은 값으로 클리핑해 솔버를 안정화한다.

에너지 방정식은 더 까다롭다. 유체와 고체가 같은 온도라는 국소 열평형 가정 아래, 유효 열전도율 $k_{\text{eff}}$를 만들어 써야 한다. 자주 쓰이는 세 혼합 모델은 다음과 같다.

입자가 평행으로 정렬돼 있으면 선형, 직렬이면 조화, 무작위면 기하평균(GM)이 실험에 가깝다. 토양·암석·다공성 금속 같은 실측에서 GM이 자주 쓰인다.

Python — 세 모델을 한 채널 위에서 비교#

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

으로 쓰인다. $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 직선만 읽으면 이 지점의 차압은 전체의 1/3이다. 절반이 되는 교점은 0.1 m/s다.

벽 쪽 log₁₀(Da)를 내리면 플러그만 남는 자리#

log₁₀(Da) 슬라이더를 −4 쪽으로 밀면 채널 중앙 속도가 한 값으로 눕는다.

$\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$와 관의 천이가 같은 문턱인지 헷갈리면 물줄기가 흐트러지는 레이놀즈 임계수를 열어 길이 척도를 대조하면 된다.