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

El número que se pide primero en un generador de vapor helicoidal es el factor de transferencia de calor frente a un tubo recto. En un tubo curvo, la fuerza centrífuga pone un par de movimientos secundarios en la sección. Ese par es un vórtice de Dean. La perturbación que Dean escribió en 1928 es pequeña solo para $De \ll 1$. El punto de operación es mucho mayor.

Al mover el deslizador de De entre 5 y 400, el par de vórtices de la sección se fortalece. El deslizador del número de partículas solo cambia cuántas trayectorias se dibujan.

Con De pequeño los dos remolinos casi no se ven. Al llevar De a varios cientos, un par desplazado hacia el exterior llena la sección. La función de corriente más adelante arma ese par con $\sin\theta$.

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$

La solución de orden cero es la parábola del tubo recto. El flujo de Poiseuille que ajusta esa misma parábola con diferencias finitas es lo que queda en esta sección cuando se apaga la curvatura.

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.

Duplicar la malla divide el error casi por 4. $3.21\times 10^{-3}/8.14\times 10^{-4}\approx 3.94$, un orden de 1.98. Ese orden es el error de discretizar el forzamiento lineal en $De=10$. Con $De=10$, $De^{2}=100$, así que el término omitido en $\psi = De^{2} F + O(De^{4})$ no es una corrección pequeña.

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

$\sqrt{0.012/(2\times 2.3)}\approx 0.0511$, y por $10{,}000$ da 511. En el extremo bajo del diámetro de bobina citado, 3.36 m, $R=1.68$ m y $\sqrt{0.012/3.36}\approx 0.060$, así que $De\approx 600$. Al meter el mismo $De\approx 510$ en el ajuste, $\log_{10}510\approx 2.708$, $2.708^{4}\approx 53.8$, $0.033\times 53.8\approx 1.77$, una razón cercana a 2.77. No es el 30–60 por ciento de la frase siguiente.

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.

En De 510 la perturbación ya no es una corrección pequeña#

Para $De \ll 1$, $\psi$ se corta en $De^{2} F$. En el punto de operación $De \approx 510$, $De^{2}=260100$. Dividir por 576 deja 451. La $F$ impresa en $\tilde r=0.5$, $\theta=\pi/2$ tiene una parte polinómica cercana a 2.37, así que $\psi$ queda cerca de 1070. No es una corrección pequeña sobre la parábola de orden cero. El orden 2.00 de las diferencias finitas en $De=10$ solo dice que ambos lados resolvieron el mismo miembro derecho lineal. Subir la malla a 80 sigue dejando $O(De^{4})$ fuera de la ecuación.

El ajuste $1+0.033(\log_{10} De)^{4}$ vale cerca de 2.77 en $De=510$. No puede escribirse junto al aumento del 30–60 por ciento como si fueran el mismo número. El valor que entra en un diseño queda solo después de elegir uno de los dos. Un orden de 2.00 no convierte esa elección en una verificación. El caso en el que un mismo código dio a la vez órdenes 0.79 y 2.00 pregunta de qué ecuación es el error que la malla realmente redujo.