Flujo de Poiseuille: Derivación de la solución analítica y validación por diferencias finitas

El flujo de Poiseuille es un flujo laminar totalmente desarrollado, impulsado por un gradiente de presión constante entre dos placas paralelas o dentro de un tubo circular. Debido a que existe una solución analítica exacta, aparece en los libros de texto de todo el mundo como el primer caso de prueba para la validación de códigos CFD.

En este artículo se derivará la solución analítica para un flujo en canal (entre placas 2D) paso a paso y se verificará el orden de convergencia mediante la escritura de un código simple de diferencias finitas (FD).

Configuración del problema#

Geometría y condiciones de contorno#

Se considera un flujo laminar totalmente desarrollado entre dos placas paralelas infinitas.

• Semiancho del canal: $h$ (las placas están situadas en $y = \pm h$)
• Dirección del flujo: $x$ (se asume infinitamente largo)
• Condiciones de contorno: no-slip (sin deslizamiento) en ambas paredes $\Rightarrow$ $u(\pm h) = 0$
• Fuerza impulsora: gradiente de presión constante $\dfrac{dp}{dx} = -G < 0$ (donde $G > 0$ es una constante)

Ecuación de gobierno#

Dado que el flujo está totalmente desarrollado, $\partial/\partial x = 0$ y $v = 0$. La ecuación de Navier-Stokes en la dirección $x$ es:

$\rho \left( u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} \right) = -\frac{dp}{dx} + \mu \frac{\partial^2 u}{\partial y^2}$

Al aplicar las condiciones de flujo totalmente desarrollado ($\partial u/\partial x = 0$, $v = 0$):

$\mu \frac{d^2 u}{dy^2} = \frac{dp}{dx} = -G$

Derivación de la solución analítica#

Paso 1: Simplificación a una EDO#

$\frac{d^2 u}{dy^2} = -\frac{G}{\mu} \equiv -C \quad (C > 0)$

Paso 2: Primera integración#

$\frac{du}{dy} = -C y + A_1$

Paso 3: Segunda integración#

$u(y) = -\frac{C}{2} y^2 + A_1 y + A_2$

Paso 4: Aplicación de condiciones de contorno#

$u(+h) = 0$:

0 = -\frac{C}{2} h^2 + A_1 h + A_2 \tag{1}

$u(-h) = 0$:

0 = -\frac{C}{2} h^2 - A_1 h + A_2 \tag{2}

(1) $-$ (2): $0 = 2 A_1 h \Rightarrow A_1 = 0$ (Simetría)

(1) $+$ (2): $0 = -C h^2 + 2A_2 \Rightarrow A_2 = \dfrac{C h^2}{2}$

Solución analítica final#

$\boxed{u(y) = \frac{G}{2\mu}\bigl(h^2 - y^2\bigr)}$

Velocidad máxima en el centro del canal ($y=0$):

$u_{\max} = \frac{G h^2}{2\mu}$

Interpretación física#

Perfil de velocidad#

La solución analítica presenta una distribución parabólica. La fricción de la pared se transmite uniformemente a través de la sección transversal del flujo, resultando en un máximo en el centro y cero en las paredes.

La velocidad adimensional $u/u_{\max} = 1 - (y/h)^2$ es un perfil universal, independiente del semiancho del canal $h$ o de la viscosidad $\mu$.

Esfuerzo cortante en la pared#

$\tau_w = \mu \left.\frac{du}{dy}\right|_{y=h} = \mu \cdot \frac{G}{\mu}(-h) = -G h$

(Magnitud: $\tau_w = G h$) Cuanto mayor sea el gradiente de presión y más ancho sea el canal, mayor será la fricción en la pared.

Caudal volumétrico#

$Q = \int_{-h}^{h} u \, dy = \frac{G}{2\mu} \int_{-h}^{h} (h^2 - y^2) \, dy = \frac{G}{2\mu} \cdot \frac{4h^3}{3} = \frac{2 G h^3}{3\mu}$

Esto se conoce como la versión 2D de la ley de Hagen-Poiseuille. Dado que el caudal es proporcional a $h^3$, duplicar el ancho del canal aumenta el caudal ocho veces.

Visualización de la solución analítica con Python#

import numpy as np
import matplotlib.pyplot as plt
 
# Parámetros
h   = 1.0    # Semiancho del canal [m]
G   = 1.0    # Gradiente de presión [Pa/m]
mu  = 1.0    # Viscosidad dinámica [Pa·s]
 
# Malla
y = np.linspace(-h, h, 400)
u_exact = G / (2 * mu) * (h**2 - y**2)
u_max   = G * h**2 / (2 * mu)
 
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
 
# Perfil de velocidad (Horizontal)
axes[0].plot(u_exact / u_max, y / h, "k-", lw=2)
axes[0].fill_betweenx(y / h, u_exact / u_max, alpha=0.15)
axes[0].axvline(0, color="gray", lw=0.8, ls="--")
axes[0].set_xlabel(r"$u / u_{\max}$")
axes[0].set_ylabel(r"$y / h$")
axes[0].set_title("Perfil de velocidad de Poiseuille")
axes[0].set_xlim(-0.05, 1.1)
 
# Distribución del esfuerzo cortante
tau = mu * np.gradient(u_exact, y)
axes[1].plot(tau / G / h, y / h, "r-", lw=2, label=r"$\tau / \tau_w$")
axes[1].axvline(0, color="gray", lw=0.8, ls="--")
axes[1].set_xlabel(r"$\tau / \tau_w$")
axes[1].set_ylabel(r"$y / h$")
axes[1].set_title("Distribución del esfuerzo cortante (Lineal)")
axes[1].legend()
 
plt.tight_layout()
plt.savefig("poiseuille_profile.png", dpi=150)
plt.show()

Validación numérica con el método de diferencias finitas#

Se discretiza $d^2u/dy^2 = -G/\mu$ utilizando un esquema de diferencia central de segundo orden. En una malla uniforme $\Delta y = 2h/(N-1)$:

$\frac{u_{i+1} - 2u_i + u_{i-1}}{(\Delta y)^2} = -\frac{G}{\mu}$

Condiciones de contorno: $u_0 = u_{N-1} = 0$.

Expresando esto como un sistema tridiagonal de ecuaciones lineales:

$\begin{pmatrix} -2 & 1 \\ 1 & -2 & 1 \\ & \ddots & \ddots & \ddots \\ & & 1 & -2 \end{pmatrix} \mathbf{u} = -\frac{G (\Delta y)^2}{\mu} \mathbf{1}$

Implementación en Python#

import numpy as np
 
def solve_poiseuille_fd(N, h=1.0, G=1.0, mu=1.0):
    """
    N : Número de puntos de malla internos (excluyendo límites)
    retorna (y, u_num, u_exact, L2_error)
    """
    dy = 2 * h / (N + 1)
    y_inner = np.linspace(-h + dy, h - dy, N)
 
    # Matriz de coeficientes tridiagonal
    diag  = -2 * np.ones(N)
    upper =  np.ones(N - 1)
    lower =  np.ones(N - 1)
    A = (np.diag(diag) + np.diag(upper, k=1) + np.diag(lower, k=-1))
 
    rhs = -G * dy**2 / mu * np.ones(N)
 
    u_num = np.linalg.solve(A, rhs)
 
    # Malla completa incluyendo límites
    y_full  = np.concatenate([[-h], y_inner, [h]])
    u_full  = np.concatenate([[0.0], u_num, [0.0]])
    u_exact = G / (2 * mu) * (h**2 - y_full**2)
 
    L2 = np.sqrt(np.mean((u_full - u_exact)**2))
    return y_full, u_full, u_exact, L2
 
# Prueba de convergencia de malla
Ns = [4, 8, 16, 32, 64, 128, 256]
errors = []
for N in Ns:
    _, _, _, err = solve_poiseuille_fd(N)
    errors.append(err)
 
# Orden de convergencia
for i in range(1, len(Ns)):
    ratio = errors[i-1] / errors[i]
    order = np.log2(ratio)
    print(f"N={Ns[i]:4d}  L2={errors[i]:.3e}  Orden de convergencia={order:.2f}")

Resultados de convergencia#

“

Nota: Dado que el flujo de Poiseuille es un polinomio cuadrático, el esquema de diferencia central de segundo orden arroja resultados al nivel de la precisión de máquina (machine precision). El error real no escala con $\Delta y^2$, sino que se satura al nivel del error de redondeo de punto flotante.

Esta es precisamente la razón por la que el flujo de Poiseuille es un caso de validación CFD ideal: proporciona respuestas casi perfectas incluso con resoluciones de malla bajas, lo que permite una verificación rápida de la precisión fundamental del código.

Gráfico comparativo de resultados#

import matplotlib.pyplot as plt
 
fig, ax = plt.subplots(figsize=(6, 5))
 
colors = plt.cm.viridis(np.linspace(0.2, 0.9, 4))
Ns_plot = [4, 8, 16, 64]
 
for N, c in zip(Ns_plot, colors):
    y, u_num, _, _ = solve_poiseuille_fd(N)[:3]
    _, _, u_exact, _ = solve_poiseuille_fd(N)
    ax.plot(u_num, y, "o--", color=c, ms=4, label=f"FD N={N}")
 
# Solución analítica
y_fine = np.linspace(-1, 1, 400)
ax.plot(0.5 * (1 - y_fine**2), y_fine, "k-", lw=2, label="Solución Analítica")
 
ax.set_xlabel(r"$u$ (Adimensional)")
ax.set_ylabel(r"$y/h$")
ax.set_title("Comparación: FD vs Solución Analítica")
ax.legend(fontsize=9)
plt.tight_layout()
plt.savefig("poiseuille_fd_vs_exact.png", dpi=150)
plt.show()

Consejos de configuración en OpenFOAM#

Al configurar el flujo de Poiseuille con icoFoam de OpenFOAM, los parámetros clave son:

0/U — Condiciones de contorno de velocidad#

boundaryField
{
    walls
    {
        type    noSlip;         // Sin deslizamiento en la pared
    }
    inlet
    {
        type    fixedValue;
        value   uniform (1 0 0); // O parabolicVelocity
    }
    outlet
    {
        type    zeroGradient;
    }
}

0/p — Condiciones de contorno de presión#

boundaryField
{
    inlet   { type fixedValue; value uniform 1; }
    outlet  { type fixedValue; value uniform 0; }
    walls   { type zeroGradient; }
}

constant/transportProperties#

nu              nu [ 0 2 -1 0 0 0 0 ] 0.01;  // Ajustar viscosidad cinemática

Lista de verificación para la validación#

1. Condición de flujo totalmente desarrollado: La longitud del canal debe ser suficientemente mayor que la longitud de desarrollo de la entrada ($L_{dev} \approx 0.05 \, Re \cdot h$).
2. Independencia de malla: Un conteo de malla de 16 o más en la dirección $y$ asegura una precisión de segundo orden.
3. Métricas de validación: Comparar la velocidad en la línea central $u_c = G h^2 / (2\mu)$ y el esfuerzo cortante en la pared $\tau_w = G h$.

Conclusión#

El flujo de Poiseuille puede parecer simple, pero este benchmark permite verificar varios aspectos:

• Si el esquema de diferencia central de segundo orden está correctamente implementado.
• Si las condiciones de contorno de presión y de no-slip en la pared se manejan adecuadamente.
• Si el tensor de esfuerzo viscoso ($\mu \nabla^2 \mathbf{u}$) se calcula con precisión en el código.

El principio de validar comenzando por problemas con soluciones analíticas (V&V) es una cultura compartida en todos los campos de CFD, incluyendo la aeroespacial, nuclear y meteorología. En el próximo artículo, se abordará el desafío de la validación de soluciones analíticas no estacionarias mediante el primer problema de Stokes (placa que arranca impulsivamente).

“

Referencias

• Batchelor, G. K. An Introduction to Fluid Dynamics. Cambridge, 1967.
• Ferziger, J. H., Perić, M. Computational Methods for Fluid Dynamics. Springer, 2002.
• Documentación de OpenFOAM — icoFoam tutorial: cavity

Modifica gradiente de presión y viscosidad para ver el perfil parabólico.