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

Dos placas, separación $2h$. La presión cae solo en $x$. En régimen desarrollado, $v=0$ y $u$ depende solo de $y$. Un tubo circular lleva el mismo nombre, flujo de Poiseuille, pero las fórmulas de abajo son el canal plano. Si el código se aparta de esta parábola, el término viscoso o la pared sin deslizamiento está mal.

Al mover el control −dp/dx o μ cambia solo el $u_{\max}$ impreso. La curva dibujada no se mueve.

La curva del lienzo es $u/u_{\max} = 1 - (y/h)^2$. El control Re solo desplaza el número de la esquina. Con $h=1$, el pico mostrado es el valor de −dp/dx dividido por $2\mu$.

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

$u(-h) = 0$:

(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.

Con $h=G=\mu=1$, las mismas fórmulas dan $u_{\max}=1/2$ y $Q=2/3$. Dividir por la altura $2h$ deja la media en $1/3$, dos tercios del pico. Integrar $1-(y/h)^2$ de $-h$ a $h$ y dividir por $2h\,u_{\max}$ devuelve el mismo $2/3$. $h$ y $\mu$ se cancelan.

Al duplicar solo $h$, $u_{\max}$ crece por 4 ($h^2$), $Q$ por 8 ($h^3$) y el esfuerzo $Gh$ por 2. Leer el esfuerzo con el factor del caudal deja la pared cuatro veces alta. El tubo laminar que fija el coeficiente de energía cinética en 2.0000 está en el coeficiente delante de la altura de velocidad.

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.

Con $N$ bajo la parábola ya se pega a la fórmula. Ahí se ve si el término viscoso sigue vivo. La columna de orden no sirve para ese juicio.

La discretización escribe $\Delta y = 2h/(N-1)$. La tabla no. Para $N=4,8,16,32,64$ los espaciados son 0.4000, 0.2222, 0.1176, 0.0606 y 0.0308. Con $h=1$ coinciden con $2/(N+1)$: $2/5=0.4$, $2/9 \approx 0.2222$, $2/17 \approx 0.1176$, $2/33 \approx 0.0606$, $2/65 \approx 0.0308$. Las dos paredes salen de las incógnitas, así que el paso es $2h/(N+1)$.

Los cinco $L_2$ son $2.78\times 10^{-16}$, $1.67\times 10^{-16}$, $2.22\times 10^{-16}$, $2.78\times 10^{-16}$ y $2.22\times 10^{-16}$. Los logaritmos en base 2 de los cocientes siguen cerca de 0.7, negativo, negativo y 0.3. No es el ~2.0 impreso. Un polinomio cuadrático en $y$ tiene derivada tercera nula, así que el error de truncamiento de segundo orden es cero. Esta columna no mide el orden.

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

Longitud de desarrollo y velocidad central#

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

¿Mirar la columna de orden cuando L2 ya es 1e-16?#

En $N=4$, $L_2$ ya vale $2.78\times 10^{-16}$. Llevar la malla hasta 64 solo hace oscilar esa columna dentro de la banda de $10^{-16}$. El ~2.0 impreso no salió de esos cocientes. Sobre un polinomio cuadrático el esquema centrado tiene error de truncamiento nulo, así que dar por bueno el orden con esta tabla es la decisión equivocada. Que $\mu \nabla^2 \mathbf{u}$ siga vivo se ve en la velocidad central y en el esfuerzo de pared.

Los números a contrastar son $Gh^2/(2\mu)$ y $Gh$. El orden que de verdad se parte vive en una solución con discontinuidad. El caso en que el mismo código da a la vez 0.79 y 2.00 mide aparte que esa columna puede no ser el orden del esquema. Con el residuo pegado al redondeo, una malla más fina tampoco produce un orden.

“

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