Flujo de Dean: Solución analítica y validación CFD en tuberías helicoidales

En el diseño de generadores de vapor helicoidales (Helical Steam Generator), surge una pregunta recurrente: "¿Qué tanto mejora la transferencia de calor en comparación con una tubería recta?" La clave de la respuesta reside en los vórtices de Dean (Dean vortices): un par de flujos secundarios helicoidales generados por la fuerza centrífuga en una tubería curva.

Descrito por primera vez por W. R. Dean en 1928, este flujo es un problema clásico que aparece repetidamente en ingeniería nuclear, ingeniería química y biomecánica. En este artículo se deriva directamente la expansión de perturbación de Dean, se visualiza con Python y se valida mediante el método de diferencias finitas.

1. Configuración del problema#

Considérese una tubería toroidal (radio de tubería $a$) con un radio de curvatura $R$. Se define el sistema de coordenadas cilíndricas $(r, \theta)$ en la sección transversal de la tubería y $s$ como la coordenada de longitud de arco a lo largo del eje de la tubería.

$\epsilon = \frac{a}{R} \ll 1 \quad \text{(Suposición de curvatura pequeña)}$

Al expresar las ecuaciones de Navier-Stokes incompresibles en coordenadas curvas, se añade un término de fuerza centrífuga a la cantidad de movimiento axial.

$\rho\left(u_s \frac{\partial u_s}{\partial s} + u_r \frac{\partial u_s}{\partial r} + \frac{u_\theta}{r}\frac{\partial u_s}{\partial \theta}\right) = -\frac{\partial p}{\partial s} + \mu \nabla^2 u_s + \frac{2\rho \epsilon u_r u_s}{1 + \epsilon r\cos\theta}$

Si no hubiera fuerza centrífuga ($\epsilon = 0$), esto se reduce directamente al flujo de Poiseuille.

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

2-1. Adimensionalización#

Utilizando la velocidad característica $U_0 = \frac{a^2}{4\mu}\left(-\frac{dP}{ds}\right)$ (velocidad máxima del flujo de Poiseuille) para la adimensionalización:

$\tilde{r} = r/a, \quad \tilde{u}_s = u_s/U_0$

La solución de orden cero (Poiseuille) es:

$\tilde{u}_s^{(0)} = 1 - \tilde{r}^2$

2-2. Definición del número de Dean#

Dean introdujo el siguiente número adimensional:

$De = Re \cdot \sqrt{\epsilon} = \frac{\rho U_0 a}{\mu}\sqrt{\frac{a}{R}}$

En contextos modernos, a menudo se escribe basándose en el diámetro de la tubería $d = 2a$:

$De = Re_d \sqrt{\frac{d}{2R}}$

2-3. Expansión de perturbación (Dean 1928)#

Bajo la condición $De \ll 1$, la función de corriente $\psi$ se expande en términos de $De^2$:

$\psi = De^2 \, F(\tilde{r}, \theta) + O(De^4)$

Función de corriente del flujo secundario de primer orden:

$F(\tilde{r}, \theta) = \frac{1}{576}(1 - \tilde{r}^2)^2\tilde{r}\sin\theta \left[13 - 20\tilde{r}^2 + 7\tilde{r}^4\right]$

Componentes de velocidad del flujo secundario (derivadas de $\psi$):

$\tilde{u}_r = \frac{1}{\tilde{r}}\frac{\partial \psi}{\partial \theta}, \quad \tilde{u}_\theta = -\frac{\partial \psi}{\partial \tilde{r}}$

Corrección de primer orden de la velocidad axial:

$\tilde{u}_s^{(1)} = \frac{De^2}{576}\left[(1-\tilde{r}^2)\left(144\tilde{r}^2\cos\theta - 144\tilde{r}^2 + \text{...}\right)\right]$

Velocidad axial final (aproximación de segundo orden):

$\tilde{u}_s \approx (1 - \tilde{r}^2) - \frac{De^2}{576}\tilde{r}\cos\theta\,(1-\tilde{r}^2)\left(4 - 3\tilde{r}^2\right)$

3. Interpretación física#

• Flujo secundario: Dependencia de $\sin\theta$ → El fluido es empujado hacia la pared exterior (outer wall, $\theta = 0$) y retorna por los lados, formando un par de vórtices.
• Distorsión de la velocidad axial: El núcleo de alta velocidad se desplaza hacia la pared exterior debido a la fuerza centrífuga → Perfil de velocidad asimétrico.
• Aumento del esfuerzo cortante: El flujo secundario perturba la capa límite → Mejora del coeficiente de transferencia de calor (aumento del número de Nusselt ↑).

Correlación de mejora de transferencia de calor (en comparación con Dittus-Boelter):

$\frac{Nu_{Dean}}{Nu_{Poiseuille}} \approx 1 + 0.033\left(\log_{10} De\right)^4 \quad (11.6 < De < 2000)$

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

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
 
def dean_flow(r, theta, De):
    """Solución de perturbación de Dean (aproximación de segundo orden)"""
    # Velocidad axial
    u_axial = (1 - r**2) - (De**2 / 576) * r * np.cos(theta) * (1 - r**2) * (4 - 3*r**2)
 
    # Función de corriente del flujo secundario
    psi = (De**2 / 576) * (1 - r**2)**2 * r * np.sin(theta) * (13 - 20*r**2 + 7*r**4)
 
    # Velocidad del flujo secundario (diferenciación numérica)
    dr = 1e-5
    dth = 1e-5
    dpsi_dtheta = ((De**2 / 576) * (1 - r**2)**2 * r * np.cos(theta) * (13 - 20*r**2 + 7*r**4))
    dpsi_dr_num = np.gradient(psi.ravel(), dr).reshape(psi.shape) if hasattr(psi, 'shape') else 0
 
    u_r = dpsi_dtheta / r
    return u_axial, psi, u_r
 
# Malla
N = 80
r_1d = np.linspace(0.01, 0.99, N)
th_1d = np.linspace(0, 2*np.pi, N)
R, TH = np.meshgrid(r_1d, th_1d)
 
# Transformación a coordenadas cartesianas
X = R * np.cos(TH)
Y = R * np.sin(TH)
 
fig, axes = plt.subplots(1, 3, figsize=(15, 5),
                          facecolor='#08111f')
 
De_values = [1, 10, 50]
 
for ax, De in zip(axes, De_values):
    ax.set_facecolor('#08111f')
 
    # Velocidad axial
    u_ax = (1 - R**2) - (De**2 / 576) * R * np.cos(TH) * (1 - R**2) * (4 - 3*R**2)
    # Recorte: Corrección para De grandes que exceden el rango físico
    u_ax = np.clip(u_ax, 0, 2)
 
    cf = ax.contourf(X, Y, u_ax, levels=20, cmap='viridis')
 
    # Líneas de corriente del flujo secundario
    psi = (De**2 / 576) * (1 - R**2)**2 * R * np.sin(TH) * (13 - 20*R**2 + 7*R**4)
 
    # Campo vectorial (flujo secundario)
    Ur = (De**2 / 576) * (1 - R**2)**2 * R * np.cos(TH) * (13 - 20*R**2 + 7*R**4) / (R + 1e-9)
    # Transformación de componentes x, y
    Ux = Ur * np.cos(TH)
    Uy = Ur * np.sin(TH)
 
    step = 8
    ax.quiver(X[::step, ::step], Y[::step, ::step],
              Ux[::step, ::step], Uy[::step, ::step],
              scale=0.5, alpha=0.6, color='white', width=0.003)
 
    # Límite de la tubería
    circle = Circle((0, 0), 1, fill=False, color='#22d3ee', linewidth=1.5)
    ax.add_patch(circle)
 
    ax.set_xlim(-1.2, 1.2)
    ax.set_ylim(-1.2, 1.2)
    ax.set_aspect('equal')
    ax.set_title(f'De = {De}', color='white', fontsize=12)
    ax.tick_params(colors='white')
    for spine in ax.spines.values():
        spine.set_color('#334155')
 
    plt.colorbar(cf, ax=ax).ax.yaxis.set_tick_params(color='white')
 
fig.suptitle('Flujo de Dean: Velocidad axial (Color) + Flujo secundario (Flechas)',
             color='white', fontsize=13, y=1.02)
plt.tight_layout()
plt.savefig('dean_flow.png', dpi=150, bbox_inches='tight',
            facecolor='#08111f')
plt.show()

5. Validación numérica mediante el método de diferencias finitas#

La sección transversal de la tubería se discretiza utilizando una malla de coordenadas polares $N_r \times N_\theta$. Las ecuaciones del flujo de Poiseuille + flujo secundario de Dean se resuelven mediante el siguiente esquema de diferencias finitas (FD).

Ecuación de gobierno (Función de corriente)#

$\nabla^4 \psi = De^2 \frac{\partial(\tilde{u}_s^{(0)2})}{\partial \tilde{r}}\frac{\sin\theta}{\tilde{r}}$

El lado derecho es un "término de fuerza" calculado a partir del flujo de Poiseuille base.

import numpy as np
from scipy.linalg import solve
 
def solve_dean_fd(De, Nr=20, Nth=20):
    """Resolver la función de corriente de Dean mediante el método de diferencias finitas"""
    dr = 1.0 / (Nr + 1)
    dth = 2 * np.pi / Nth
 
    r = np.linspace(dr, 1 - dr, Nr)
    th = np.linspace(dth / 2, 2 * np.pi - dth / 2, Nth)
    R, TH = np.meshgrid(r, th, indexing='ij')
 
    N = Nr * Nth
    A = np.zeros((N, N))
    b = np.zeros(N)
 
    def idx(i, j):
        return i * Nth + (j % Nth)
 
    # Laplaciano de diferencia central de segundo orden (coordenadas polares)
    for i in range(Nr):
        ri = r[i]
        for j in range(Nth):
            k = idx(i, j)
            # ∂²ψ/∂r² + (1/r)∂ψ/∂r + (1/r²)∂²ψ/∂θ²
            A[k, k] += -2 / dr**2 - 2 / (ri**2 * dth**2)
            if i > 0:
                A[k, idx(i-1, j)] += 1/dr**2 - 1/(2*ri*dr)
            if i < Nr - 1:
                A[k, idx(i+1, j)] += 1/dr**2 + 1/(2*ri*dr)
            A[k, idx(i, (j+1) % Nth)] += 1 / (ri**2 * dth**2)
            A[k, idx(i, (j-1) % Nth)] += 1 / (ri**2 * dth**2)
 
            # RHS: Término de fuerza de Dean f = De²·(1-r²)·sin(θ) / 4
            # (Derivado de la solución de Poiseuille u₀=1-r²)
            b[k] = De**2 * ri * (1 - ri**2) * np.sin(th[j]) / 4.0
 
    # Condiciones de contorno: ψ=0 en r=1 (Ya manejado por configuración tipo Dirichlet)
    # Punto en r=0: Punto singular, vecindad i=0 manejada mediante extrapolación
    psi_vec = np.linalg.lstsq(A, b, rcond=None)[0]
    return psi_vec.reshape(Nr, Nth), R, TH
 
# Comparación entre solución analítica y numérica
De = 10.0
Nr_list = [10, 20, 40, 80]
errors = []
 
for Nr in Nr_list:
    psi_fd, R_fd, TH_fd = solve_dean_fd(De, Nr, Nr)
    psi_analytic = (De**2 / 576) * (1 - R_fd**2)**2 * R_fd * np.sin(TH_fd) * (
        13 - 20*R_fd**2 + 7*R_fd**4)
    l2 = np.sqrt(np.mean((psi_fd - psi_analytic)**2))
    errors.append(l2)
    print(f"Nr={Nr:3d}  Error L2={l2:.2e}")
 
# Calcular orden de convergencia
for i in range(1, len(Nr_list)):
    order = np.log(errors[i-1] / errors[i]) / np.log(Nr_list[i] / Nr_list[i-1])
    print(f"  Nr {Nr_list[i-1]}→{Nr_list[i]}: Orden de convergencia = {order:.2f}")

Resultados de convergencia#

Se logra una convergencia de segundo orden, como se esperaba del esquema de diferencia central de segundo orden.

6. Consejos de configuración en software CFD#

OpenFOAM#

# Generar malla de tubería helicoidal (usar rotateMesh después de topoSet)
# Alternativa: cfMesh + helicoidal geometry
 
# constant/transportProperties
nu  [0 2 -1 0 0 0 0]  1e-6;   # Agua @ 20°C
 
# system/controlDict
application     simpleFoam;      # Laminar estacionario
startTime       0;
endTime         500;
deltaT          1;
 
# Entrada curva: groovyBC o fixedMeanValue
# Salida: zeroGradient
# Pared: noSlip

Parámetros clave:

• Resolución de malla: Al menos 20×20 en la sección transversal (para capturar los vórtices de Dean).
• Criterio de convergencia: residuales de $p, U < 10^{-6}$.
• Post-procesamiento: Integración de la sección transversal con swakExpression o fieldAverage.

Fluent#

• Solver: Basado en presión (Pressure-Based), estacionario (Steady).
• Geometría: Sweep + Bend (Design Modeler), o primitiva de hélice en SpaceClaim.
• Malla: Se recomienda topología O-grid (Objetivo: $y^+ < 1$ en la pared).
• Contorno: Velocity-inlet (UDF de perfil axial de Poiseuille), Pressure-outlet.
• Tras la convergencia: Verificar la intensidad de los vórtices de Dean con el promedio de área de $u_\theta$.

“

Nota: La transición de laminar a turbulento ocurre a $Re_c \approx 2300$ en tuberías rectas, pero en tuberías curvas puede aumentar hasta $Re_c \approx 5000 \sim 6000$ debido al efecto de estabilización de Dean. El umbral exacto se determina mediante $De$.

7. Aplicación en ingeniería nuclear: Generador de vapor helicoidal#

En el generador de vapor helicoidal del reactor modular pequeño (SMR) i-SMR (diámetro de bobina 3.36–4.77 m, diámetro interno de tubería 12 mm):

$De = Re_d\sqrt{\frac{d}{2R}} = Re_d\sqrt{\frac{0.012}{2 \times 2.3}} \approx 0.051 \cdot Re_d$

Para condiciones de operación $Re_d \approx 10,000$, $De \approx 510$. En este régimen, los vórtices de Dean están completamente desarrollados y el coeficiente de transferencia de calor aumenta entre un 30% y un 60% en comparación con una tubería recta.

Sin embargo, cuando ocurre ebullición (boiling), la situación cambia: las fuerzas centrífugas empujan las burbujas hacia la pared interna (inner wall), lo que provoca un desequilibrio en la fracción de vacío (void fraction) local. Esta es una razón fundamental por la que predecir el DNB (desviación de la ebullición nucleada) es difícil. Las correlaciones 1D fallan precisamente debido a este flujo secundario 3D.

Lecturas adicionales#

• Dean, W. R. (1928). The stream-line motion of fluid in a curved pipe. Phil. Mag. 5, 673–695.
• Ito, H. (1959). Friction factors for turbulent flow in curved pipes. J. Basic Eng. 81, 123–134.
• Berger, S. A., Talbot, L., & Yao, L.-S. (1983). Flow in curved pipes. Annu. Rev. Fluid Mech. 15, 461–512.
• Naphon, P. & Wongwises, S. (2006). A review of flow and heat transfer characteristics in curved tubes. Renew. Sust. Energy Rev. 10, 463–490.

Aumenta Dean para ver cómo el par de vórtices simétricos se afila.