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.

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. Los controles log₁₀(K) y β desplazan la línea roja discontinua, la velocidad en la que se igualan el término viscoso y el inercial.

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.

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

La ley de Darcy se escribe

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$,

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.

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.

Igualar $(\mu/K)\,u$ con $\beta\rho u^{2}$ deja la velocidad de cruce en una división, $u^{\star}=\mu/(K\beta\rho)$. Los números del lecho de grava en el fragmento son $\mu=10^{-3}$, $K=10^{-8}$, $\beta=10^{3}$, $\rho=1000$, así que el denominador vale $10^{-2}$ y $u^{\star}=0.1$ m/s. Doblar solo $\beta$ lleva $u^{\star}$ a $0.05$ m/s y el mismo punto de $0.1$ m/s queda del lado inercial. Multiplicar $K$ por diez arrastra $u^{\star}$ con él.

En $u=0.05$ m/s la razón inercia/viscosidad es $u/u^{\star}=1/2$ y la curva aún es suave. En $u=0.2$ m/s la razón vale 2 y la caída total triplica el término viscoso, así que Darcy solo conserva un tercio. La mitad exacta está en $0.1$ m/s, no en $0.2$ m/s.

Un Reynolds con diámetro de grano llama al término cuadrático cerca de 1, como cuenta cómo un número de Reynolds parte el flujo. Solo cambia la longitud, de la tubería al grano $d$.

Con $G=1$, la velocidad central es $\mathrm{Da}$ por $(1-1/\cosh(0.5/\sqrt{\mathrm{Da}}))$. En $\mathrm{Da}=10^{-4}$, $0.5/\sqrt{\mathrm{Da}}=50$, así que $\cosh$ se dispara y la velocidad central se pega al propio $\mathrm{Da}$, $10^{-4}$. En $\mathrm{Da}=1$, $0.5/\sqrt{\mathrm{Da}}=0.5$, el denominador solo se aparta un poco de 1 y el tapón ya no cubre el canal.

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.

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

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

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

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

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 lectura solo con Darcy se queda ahí con cerca de un tercio de la caída real. El cruce a la mitad está en 0.1 m/s.

Donde bajar log₁₀(Da) deja solo el tapón#

Al llevar el control log₁₀(Da) hacia −4, la velocidad en el centro del canal se aplana a un solo valor.

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.

Donde una curva de bomba recta deja en la mitad la caída del lecho de grava#

Con $\mathrm{Re}_p < 1$ la recta es la cara correcta, y con $\mathrm{Re}_p \gtrsim 1$ entra el término cuadrático. Sin colocar $u^{\star}=\mu/(K\beta\rho)$ junto al punto de operación, la curva recta de Darcy sigue elegida aunque el punto ya esté a la derecha. Con los números del lecho de grava en el fragmento, ese límite es 0.1 m/s. En 0.2 m/s la recta señala un tercio de la caída real.

Una pared que deja el canal en $\sqrt{\mathrm{Da}} \lesssim L$ parte la decisión otra vez. Quitar el laplaciano viscoso antes de medir la capa de espesor $\sqrt{K}$ hace que la velocidad en la pared no caiga a cero. Si el $\mathrm{Re}_p$ del poro y la transición en una tubería parecen el mismo umbral, conviene abrir el Reynolds en el que un chorro de agua se desarma y contrastar la longitud que usa cada uno.