[Tutorial de OpenFOAM] Análisis de flujo en escalón orientado hacia atrás — Capturando zonas de recirculación con simpleFoam

El tramo de 1.5 m aguas abajo son 15 alturas de escalón en la malla y 150 si cuenta el H = 0.01 m del comentario de nu.

Justo detrás de un escalón invertido la velocidad en x se vuelve negativa. La x donde esa zona termina es el punto de readhesión. Un caso de simpleFoam solo puede colocar ese punto en alturas de escalón cuando la altura de la malla y el comentario de nu nombran el mismo H.

El deslizador de Re va de 50 a 1500 y alarga la bolsa de flujo inverso detrás del escalón. El deslizador del número de partículas, de 50 a 500, solo cambia cuántas trayectorias llenan esa bolsa.

Re 400입자 220

스텝 직후 재순환 영역(Lr)이 형성된다. Re가 커지면 Lr가 길어지다가 Re ≳ 1000부근부터 비정상이 된다 — Armaly et al. 1983 실험과 일치.

Abajo, Lr queda escrito en alturas de escalón. Con Re bajo la burbuja sigue pegada al escalón. Hacia el extremo del deslizador se estira aguas abajo.

Requisitos previos#

OpenFOAM v2312 o posterior. El armado del caso es el mismo en .com y en .org.
Una terminal y las tres carpetas del caso: 0, constant y system.
Las figuras se abren en ParaView 5.x.

Paso 1: Creación del directorio del caso#

Un directorio vacío recibe 0, constant y system.

# Cargar el entorno de OpenFOAM (varía según la ruta de instalación)
source /opt/openfoam2312/etc/bashrc
 
# Copiar el directorio del tutorial o crear uno directamente
mkdir -p backwardFacingStep
cd backwardFacingStep

El árbol queda partido así.

tree backwardFacingStep/

backwardFacingStep/
├── 0/
│   ├── U          # Condiciones iniciales/de contorno de velocidad
│   ├── p          # Condiciones iniciales/de contorno de presión
│   ├── k          # Energía cinética turbulenta
│   └── epsilon    # Tasa de disipación turbulenta
├── constant/
│   ├── transportProperties
│   └── turbulenceProperties
└── system/
    ├── blockMeshDict
    ├── fvSchemes
    └── fvSolution

Paso 2: Generación de malla — blockMeshDict#

system/blockMeshDict parte el escalón en dos bloques. El hueco en y de los vértices es 1 y scale es 0,1, así que el H = 0,1 m del comentario concuerda con esta malla.

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2312                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
scale   0.1;   // Escala global (se multiplica por 0.1 para obtener el tamaño final en metros)
 
vertices
(
    // Bloque del canal de entrada (sobre el escalón)
    (0  1  0)  // 0
    (5  1  0)  // 1
    (5  2  0)  // 2
    (0  2  0)  // 3
    (0  1  1)  // 4
    (5  1  1)  // 5
    (5  2  1)  // 6
    (0  2  1)  // 7
 
    // Bloque aguas abajo del escalón (altura completa)
    (5  0  0)  // 8
    (20 0  0)  // 9
    (20 2  0)  // 10
    (5  2  0)  // 11 (= vértice 2)
    (5  0  1)  // 12
    (20 0  1)  // 13
    (20 2  1)  // 14
    (5  2  1)  // 15 (= vértice 6)
);
 
blocks
(
    // Canal de entrada
    hex (0 1 2 3 4 5 6 7)   (50 10 1) simpleGrading (1 1 1)
    // Aguas abajo del escalón (malla más densa)
    hex (8 9 10 11 12 13 14 15) (150 20 1) simpleGrading (1 1 1)
);
 
boundary
(
    inlet
    {
        type patch;
        faces ((0 3 7 4));
    }
    outlet
    {
        type patch;
        faces ((9 10 14 13));
    }
    wall
    {
        type wall;
        faces
        (
            (0 1 5 4)   // Pared superior de entrada
            (1 8 12 5)  // Superficie vertical del escalón
            (8 9 13 12) // Suelo inferior
            (3 2 6 7)   // Pared superior (entrada)
            (11 10 14 15) // Pared superior (aguas abajo)
        );
    }
    frontAndBack
    {
        type empty;  // Análisis 2D
        faces
        (
            (0 1 2 3)
            (4 5 6 7)
            (8 9 10 11)
            (12 13 14 15)
        );
    }
);

# Ejecutar generación de malla
blockMesh

Salida prevista:

Create time
 
Creating block mesh from "system/blockMeshDict"
...
Writing polyMesh with 3600 cells
End

Paso 3: Configuración de condiciones de contorno#

0/U — Condiciones de contorno de velocidad#

La entrada es 1 m/s en x, las paredes son noSlip y la salida es zeroGradient.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
dimensions      [0 1 -1 0 0 0 0];  // m/s
 
internalField   uniform (0 0 0);   // Valor inicial: estado de reposo
 
boundaryField
{
    inlet
    {
        type        fixedValue;
        value       uniform (1 0 0);  // Entrada uniforme de 1 m/s en dirección x
    }
    outlet
    {
        type        zeroGradient;  // Salida: flujo libre
    }
    wall
    {
        type        noSlip;        // Pared: condición de no deslizamiento
    }
    frontAndBack
    {
        type        empty;
    }
}

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

La presión es cinemática. Solo la salida queda fija en 0.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
 
dimensions      [0 2 -2 0 0 0 0];  // Presión cinemática (p/rho), m²/s²
 
internalField   uniform 0;
 
boundaryField
{
    inlet
    {
        type        zeroGradient;
    }
    outlet
    {
        type        fixedValue;
        value       uniform 0;  // Presión de referencia 0 en la salida
    }
    wall
    {
        type        zeroGradient;
    }
    frontAndBack
    {
        type        empty;
    }
}

0/k — Energía cinética turbulenta#

Una intensidad de entrada del 5% da $k = 1.5 \times (0.05)^2 = 0.00375$ , el valor escrito en el archivo.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      k;
}
 
dimensions      [0 2 -2 0 0 0 0];  // m²/s²
 
// Intensidad de turbulencia de entrada 5%, velocidad 1 m/s → k = 1.5*(I*U)² = 1.5*(0.05*1)² ≈ 0.00375
internalField   uniform 0.00375;
 
boundaryField
{
    inlet
    {
        type        fixedValue;
        value       uniform 0.00375;
    }
    outlet
    {
        type        zeroGradient;
    }
    wall
    {
        type        kqRWallFunction;
        value       uniform 0.00375;
    }
    frontAndBack
    {
        type        empty;
    }
}

Paso 4: Configuración de propiedades de transporte y modelo de turbulencia#

constant/transportProperties#

ν está en este archivo. Si el H que produce Re 1000 coincide con la altura de la malla se contrasta después de la tabla de resultados.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      transportProperties;
}
 
transportModel  Newtonian;
nu              1e-5;   // Coeficiente de viscosidad cinemática (m²/s), Re = U*H/nu = 1*0.01/1e-5 = 1000

constant/turbulenceProperties#

RAS queda encendido, con kEpsilon.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      turbulenceProperties;
}
 
simulationType  RAS;
 
RAS
{
    RASModel        kEpsilon;   // Modelo de dos ecuaciones k-epsilon
    turbulence      on;
    printCoeffs     on;
}

Paso 5: Configuración de esquemas numéricos#

system/fvSchemes#

El término temporal es steadyState, la convección de la velocidad es linearUpwind, y k y epsilon quedan en upwind.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
 
ddtSchemes      { default steadyState; }  // Análisis de estado estacionario
 
gradSchemes
{
    default         Gauss linear;
    grad(U)         cellLimited Gauss linear 1;  // Limitador de gradiente
}
 
divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);  // Precisión de segundo orden
    div(phi,k)      bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}
 
laplacianSchemes { default Gauss linear corrected; }
 
interpolationSchemes { default linear; }
 
snGradSchemes { default corrected; }

system/fvSolution#

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
 
solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-7;
        relTol          0.01;
    }
    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
    }
    "(k|epsilon)"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
    }
}
 
SIMPLE
{
    nNonOrthogonalCorrectors    0;
    pRefCell                    0;
    pRefValue                   0;
 
    residualControl
    {
        p           1e-4;
        U           1e-4;
        "(k|epsilon)" 1e-4;
    }
}
 
relaxationFactors
{
    fields      { p 0.3; }           // Factor de relajación de presión (para evitar divergencia)
    equations
    {
        U       0.7;
        k       0.7;
        epsilon 0.7;
    }
}

Paso 6: Ejecución del análisis#

La relajación 0,3 de la presión y el 0,7 de U, k y epsilon son el acople que lleva el registro a la iteración 500 con residuos bajo 1e-4. Cómo esa interpolación de presión esconde un tablero de ajedrez queda aparte, en SIMPLE y Rhie–Chow.

# Verificación de la calidad de la malla (opcional, pero recomendado)
checkMesh
 
# Ejecutar simpleFoam
simpleFoam 2>&1 | tee log.simpleFoam

Registro de convergencia:

Time = 1
 
smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 0.0234, No Iterations 3
smoothSolver:  Solving for Uy, Initial residual = 0.412, Final residual = 0.00891, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.987, Final residual = 0.00821, No Iterations 15
...
 
Time = 500
 
smoothSolver:  Solving for Ux, Initial residual = 1.23e-05, ...
GAMG:  Solving for p, Initial residual = 8.91e-05, ...
 
SIMPLE solution converged in 500 iterations

La parada es todo residuo bajo 1e-4. Este registro cruza esa línea en la iteración 500.

Paso 7: Visualización de resultados (ParaView)#

paraFoam abre el caso.

# Crear un archivo ficticio para ParaView y ejecutar
touch backwardFacingStep.foam
paraFoam

En ParaView el orden es este.

En Coloring se elige U y se pinta Magnitude.
Glyph esparce los vectores y la bolsa inversa aparece.
Stream Tracer dibuja las líneas de corriente.

Pasado el campo a VTK, este comando basta para los números.

foamToVTK -latestTime

La x donde ux cerca de la pared cambia de signo se toma como readhesión. El script deja H clavado en 0,01 m.

import numpy as np
 
# Leer datos de la pared inferior (y=0) del archivo de salida de foamToVTK
# (Ajustar la ruta según el entorno)
data = np.loadtxt("VTK/backwardFacingStep_500/internal.vtk",
                  skiprows=10)  # Excluir cabecera VTK
 
x = data[:, 0]
ux = data[:, 3]  # Velocidad en dirección x
 
# Filtrado de celdas cerca de la pared inferior (y ≈ 0)
wall_mask = data[:, 1] < 0.005
x_wall = x[wall_mask]
ux_wall = ux[wall_mask]
 
# Posición del cambio de signo = punto de reajuste
sign_changes = np.where(np.diff(np.sign(ux_wall)))[0]
if len(sign_changes) > 0:
    x_reattach = x_wall[sign_changes[0]]
    H = 0.01  # Altura del escalón (m)
    print(f"Longitud de reajuste: {x_reattach/H:.1f}H")

Errores comunes y soluciones#

Error 1: `FOAM FATAL ERROR: Cannot find file "0/epsilon"`#

--> FOAM FATAL ERROR:
    cannot find file
    "backwardFacingStep/0/epsilon"

Falta 0/epsilon. Se arma uno con el mismo esqueleto que k.

// 0/epsilon
dimensions      [0 2 -3 0 0 0 0];  // m²/s³
 
// epsilon = C_mu^(3/4) * k^(3/2) / L, L ≈ 0.1*H
internalField   uniform 0.0009;
 
boundaryField
{
    inlet       { type fixedValue; value uniform 0.0009; }
    outlet      { type zeroGradient; }
    wall        { type epsilonWallFunction; value uniform 0.0009; }
    frontAndBack { type empty; }
}

Error 2: Divergencia del análisis (Floating point exception)#

Floating point exception (core dumped)

Con demasiada relajación, o con una malla gruesa, el caso revienta. Se bajan los factores en fvSolution.

relaxationFactors
{
    fields      { p 0.2; }   // Reducción de 0.3 → 0.2
    equations   { U 0.5; }   // Reducción de 0.7 → 0.5
}

Error 3: El residuo no disminuye y oscila#

nNonOrthogonalCorrectors está en 0 y la malla se aleja de la ortogonalidad, así que el residuo oscila. Se sube el número de correcciones.

SIMPLE
{
    nNonOrthogonalCorrectors    2;  // 0 → 2
    ...
}

Resumen de resultados previstos#

Ítem	Valor
Número de Reynolds ( $Re = UH/\nu$ )	1,000
Longitud de reajuste ( $x_r / H$ )	6–8
Velocidad máxima inversa	$\approx -0.1$ m/s
Número de iteraciones de convergencia	300–600 veces

Armaly et al. 1983, con $Re = 800$ , queda cerca de 6H, al lado del 6–8 de la tabla.

x_r \approx 6\text{–}8 \cdot H \quad (Re \sim 10^3)

El escalón va de y = 0 a y = 1. Con scale 0,1, H es 0,1 m. El comentario de nu da $Re = 1 \times 0.01 / 10^{-5} = 1000$ , con H = 0,01 m. La malla da $1 \times 0.1 / 10^{-5} = 10^{4}$ . 6–8H son 0,6–0,8 m, o 0,06–0,08 m.

x va de 5 a 20. Cara en 0,5 m, salida en 2,0 m, tramo de 1,5 m: 15H, y 0,6–0,8 m cabe. Con el H del comentario el tramo es 150H y 6H queda a 0,06 m, en x = 0,56 m. El script usa H = 0,01 m. La celda de fondo mide $(2/20) \times 0.1 = 0,01$ m. Con y bajo 0,005 se toma la mitad.

k coincide con $1.5 \times (0.05)^2 = 0.00375$ . Con $C_\mu = 0.09$ sale $C_\mu^{3/4} k^{3/2} \approx 3.77 \times 10^{-5}$ , y sobre 0,0009 la longitud es 0,042 m. $L \approx 0.1H$ es 0,01 m o 0,001 m: epsilon 0,00377 o 0,0377, no el 0,0009 del archivo. El inverso de $-0,1$ m/s es el 10% de 1 m/s, el doble de 0,05 m/s.

La altura de escalón que estira 6H por diez#

Con la altura de malla el Re 1000 es falso. nu debe subir a $10^{-4}$ . Con la del comentario, scale baja a 0,01 y 6H queda en 0,06 m, al inicio de 1,5 m.

¿Cuál es el H del caso? Si el 1000 de entrada es el Reynolds de transición en tubería, el k-épsilon ya encendido puede sobrar, o el escalón separarse antes. Sigue abierto en el instante en que un chorro de tubería se desordena.