Las tres caras del flujo en medios porosos — Darcy, Forchheimer y Brinkman

1856, el informe de Henry Darcy al ayuntamiento de Dijon#

A mediados del siglo XIX, la ciudad francesa de Dijon necesitaba un sistema nuevo para llevar agua limpia a sus habitantes. El ingeniero civil Henry Darcy propuso filtrar el agua del río a través de un banco de arena a las afueras y respaldó la idea con un experimento obstinadamente simple. Aplicar una caída de presión sobre una columna vertical llena de arena y medir el caudal que sale por debajo. La relación entre ambos resultó estrictamente lineal.

Hasta entonces la mecánica de fluidos hablaba solo el idioma del flujo libre: Bernoulli, Euler, Navier–Stokes. La línea de Darcy abrió otro mundo que dormía justo al lado.

Hoy partimos de esa línea. Cuando la velocidad sube, el término cuadrático de Forchheimer se cuela. Cuando se acerca una pared, el laplaciano viscoso de Brinkman vuelve a la vida. Recorreremos las tres caras que adopta un fluido al atravesar los poros, con una curva de bomba y un perfil de canal.

Cuando la curva de bomba empieza a doblarse — la corrección de Forchheimer#

La ley de Darcy se escribe

\mathbf{u}_D = -\frac{K}{\mu}\,\nabla p

donde $\mathbf{u}_D$ es la velocidad de Darcy (la velocidad media en la sección, que incluye tanto poros como sólido, m/s), $K$ es la permeabilidad (m²), $\mu$ es la viscosidad dinámica y $\nabla p$ es el gradiente de presión. El fluido dentro del poro real va más rápido: $\mathbf{u} = \mathbf{u}_D / \phi$ , con $\phi$ la porosidad (la fracción de volumen ocupada por los poros).

Un pequeño aumento de velocidad cambia la escena. Definiendo un Reynolds de poro con el diámetro de grano $d$ y la velocidad de poro $u$ ,

\mathrm{Re}_p = \frac{\rho\,u\,d}{\mu}

en cuanto $\mathrm{Re}_p$ supera la unidad, la curva de caída de presión frente al caudal se aparta de la recta y comienza a doblarse. Forchheimer (1901) escribió esa curvatura como un término cuadrático.

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

El coeficiente $\beta$ (1/m) marca la velocidad a la que las pérdidas inerciales — el fluido teniendo que esquivar cada grano — alcanzan a las pérdidas viscosas. Lechos de grava, aletas de disipador, reactores catalíticos: los medios porosos reales viven casi siempre sobre esta curva.

Cuando la pared sigue viva — el regreso de Brinkman#

Tanto Darcy como Forchheimer absorben la viscosidad en una pérdida a escala de grano. Esa hipótesis se rompe en el momento en que el medio poroso toca una pared o una capa de fluido abierto. La condición no-slip obliga a que la velocidad caiga a cero en la pared, y dentro de esa película fina el laplaciano viscoso desnudo tiene que volver a estar activo.

Brinkman (1947) lo escribió de nuevo.

-\nabla p = \frac{\mu}{K}\,\mathbf{u}_D - \mu_{\text{eff}}\,\nabla^2\mathbf{u}_D

$\mu_{\text{eff}}$ es una viscosidad efectiva, que en el caso simple se toma como $\mu/\phi$ o se calibra con experimentos. La regla que sirve para comparar es el número de Darcy

\mathrm{Da} = \frac{K}{L^2}

con $L$ una longitud macroscópica como la semi-anchura del canal. Cuando $\mathrm{Da} \to 0$ el canal queda casi todo lleno por un tapón de Darcy. Cuando $\mathrm{Da} \to \infty$ la viscosidad recupera el control y el perfil se acerca a una parábola de Poiseuille. Entre ambos límites aparece una capa límite delgada de Brinkman de espesor $\sqrt{K}$ , y tiene que ser tan pequeña como para abarcar uno o dos poros antes de que usar Darcy puro sea seguro.

Promediado por fracción de poro — cuando fluido y sólido viven en la misma celda#

Imaginemos una celda CFD que contiene a la vez granos de arena (sólido) y agua entre ellos (fluido). Las ecuaciones de conservación se escriben como promedios ponderados por la fracción de poro. La conservación de masa, por sí sola, queda

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

y la forma es naturalmente conservativa si los flujos de cara usan $\mathbf{u}_D$ . Si la celda vecina tiene $\phi = 0$ (una pared sólida), el flujo de cara también debe anularse; de lo contrario se cuela una masa falsa hacia un poro que no existe. Numéricamente, las celdas con $\phi = 0$ se estabilizan forzando la diagonal a uno con extra-diagonales nulas, o recortando $\phi$ a un valor mínimo como $10^{-8}$ .

La ecuación de energía es más sensible. Bajo la hipótesis de equilibrio térmico local (fluido y sólido comparten la misma temperatura) hay que combinar las conductividades en un $k_{\text{eff}}$ efectivo. Las tres reglas de mezcla más usadas son

\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}

La lineal encaja con granos alineados en paralelo, la armónica con granos en serie y la media geométrica (GM) ajusta mejor las orientaciones aleatorias. Suelos, rocas y metales porosos suelen seguir el modelo GM en mediciones reales.

Python — tres modelos en un mismo canal#

La ecuación de cantidad de movimiento 1-D, adimensional, para un canal poroso con paredes es

\frac{d^2 U}{dY^2} - \frac{1}{\mathrm{Da}}\,U + G = 0,\qquad U(0)=U(1)=0

con $G$ el forzamiento adimensional de presión, $\mathrm{Da}$ el número de Darcy y $U$ la velocidad adimensional. Existe una solución cerrada, $U(Y) = G\,\mathrm{Da}\,\bigl(1 - \cosh((Y-0.5)/\sqrt{\mathrm{Da}}) / \cosh(0.5/\sqrt{\mathrm{Da}})\bigr)$ , pero resolverla numéricamente deja la estructura más a la vista.

import numpy as np
 
def brinkman_channel_velocity(Da: float, G: float = 1.0, N: int = 200):
    """Perfil 1D estacionario de Brinkman en un canal poroso con paredes. Da = K/L²."""
    Y = np.linspace(0.0, 1.0, N)
    dY = Y[1] - Y[0]
    # U'' - U/Da + G = 0  →  sistema tridiagonal
    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   # pared izquierda U=0
    A[-1, -1] = 1.0; rhs[-1] = 0.0   # pared derecha 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):
    """Caída de presión en una columna porosa uniforme: Darcy + Forchheimer."""
    visc_loss    = (mu / K) * u
    inertia_loss = beta * rho * u * np.abs(u)
    return visc_loss + inertia_loss
 
 
if __name__ == "__main__":
    # 1) Perfil de Brinkman: barrer Da y observar la velocidad central
    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) Curva de bomba para lecho de grava (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}")

Al ejecutarlo, la velocidad central se va aplanando conforme $\mathrm{Da}$ disminuye: casi todo el canal pasa a ser un tapón de Darcy. En la curva de bomba del lecho de grava, hasta $u \sim 0.05$ m/s la línea es casi recta, pero por encima de $0.2$ m/s la curvatura cuadrática se vuelve clara. Una bomba diseñada solo con Darcy se equivoca exactamente ahí, casi por un factor de dos.

Sliders para la curva de bomba y el perfil junto a la pared#

En las simulaciones de abajo se pueden mover los parámetros directamente.

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.

Con $\beta = 0$ la curva es una línea recta. Al subir $\beta$ se dobla hacia arriba. Una vez que el punto de operación en el régimen de grava ( $\beta \sim 10^3$ ) cruza $u^{\star}$ , una bomba diseñada solo con Darcy subestima la pérdida de presión casi a la mitad.

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.

Si se baja $\mathrm{Da}$ hasta $10^{-4}$ , el centro del canal pasa a ser un tapón plano. Desde el borde de ese tapón hasta la pared sobrevive una capa viscosa fina de espesor $\sqrt{\mathrm{Da}}$ (la capa límite de Brinkman). Esa capa es la huella dactilar de la interfaz entre el medio poroso y el flujo libre.

Las tres caras del flujo poroso — un renglón cada una#

Darcy ( $\mathrm{Re}_p < 1$ ): gradiente de presión = pérdida viscosa. La recta es la verdad.
Forchheimer ( $\mathrm{Re}_p \gtrsim 1$ ): se suman las pérdidas inerciales al rodear los granos. La curva de bomba se dobla.
Brinkman ( $\sqrt{\mathrm{Da}} \lesssim L$ ): basta con que una pared o un flujo abierto estén cerca para que el laplaciano viscoso vuelva. Reaparece una capa límite de espesor $\sqrt{K}$ .

La misma arena o la misma grava cambia entre estas tres caras dependiendo de dónde caen Reynolds y Darcy. Fijar ambos números antes de cualquier simulación ya resuelve nueve décimas partes de la elección del modelo.