[OpenFOAM Tutorial] Backward Facing Step Flow Analysis — Capturing Recirculation Zones with simpleFoam

The 1.5 m downstream block is 15 step heights on the mesh and 150 if the viscosity comment's H = 0.01 m is the one that counts.

Right behind a backward-facing step the streamwise speed goes negative. The x where that region ends is the reattachment point. A simpleFoam run can place that point in step heights only when the mesh height and the viscosity comment name the same H.

The Re slider runs from 50 to 1500 and stretches the reverse-flow pocket behind the step. The particle-count slider, from 50 to 500, only changes how many paths fill that pocket.

Re 400입자 220

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

What moves on screen is that pocket, and the Lr readout under the picture, counted in step heights. Low on the slider the bubble stays against the step. Toward the top it strings out downstream.

Prerequisites#

OpenFOAM v2312 or later. The case layout is the same on .com and on .org.
A terminal, plus the three case directories 0, constant, and system.
Pictures open in ParaView 5.x.

Step 1: Create Case Directory#

An empty directory receives 0, constant, and system.

# Load OpenFOAM environment (varies by installation path)
source /opt/openfoam2312/etc/bashrc
 
# Copy tutorial directory or create directly
mkdir -p backwardFacingStep
cd backwardFacingStep

The tree has to split like this.

tree backwardFacingStep/

backwardFacingStep/
├── 0/
│   ├── U          # Velocity initial/boundary conditions
│   ├── p          # Pressure initial/boundary conditions
│   ├── k          # Turbulent kinetic energy
│   └── epsilon    # Turbulent dissipation rate
├── constant/
│   ├── transportProperties
│   └── turbulenceProperties
└── system/
    ├── blockMeshDict
    ├── fvSchemes
    └── fvSolution

Step 2: Mesh Generation — blockMeshDict#

system/blockMeshDict cuts the step into two blocks. The y gap on the vertices is 1 and scale is 0.1, so the comment's H = 0.1 m agrees with this mesh.

/*--------------------------------*- 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;   // Global scale (multiplied by 0.1 for final size in meters)
 
vertices
(
    // Inlet channel block (above the step)
    (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
 
    // Downstream block (full height)
    (5  0  0)  // 8
    (20 0  0)  // 9
    (20 2  0)  // 10
    (5  2  0)  // 11 (= vertex 2)
    (5  0  1)  // 12
    (20 0  1)  // 13
    (20 2  1)  // 14
    (5  2  1)  // 15 (= vertex 6)
);
 
blocks
(
    // Inlet channel
    hex (0 1 2 3 4 5 6 7)   (50 10 1) simpleGrading (1 1 1)
    // Downstream (finer mesh)
    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)   // Top inlet wall
            (1 8 12 5)  // Step vertical face
            (8 9 13 12) // Bottom floor
            (3 2 6 7)   // Top wall (inlet)
            (11 10 14 15) // Top wall (downstream)
        );
    }
    frontAndBack
    {
        type empty;  // 2D Analysis
        faces
        (
            (0 1 2 3)
            (4 5 6 7)
            (8 9 10 11)
            (12 13 14 15)
        );
    }
);

# Execute mesh generation
blockMesh

Expected output:

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

Step 3: Boundary Condition Setup#

0/U — Velocity Boundary Conditions#

The inlet is 1 m/s along x, the walls are noSlip, and the outlet is 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);   // Initial value: rest state
 
boundaryField
{
    inlet
    {
        type        fixedValue;
        value       uniform (1 0 0);  // Uniform inflow of 1 m/s in x-direction
    }
    outlet
    {
        type        zeroGradient;  // Outlet: free outflow
    }
    wall
    {
        type        noSlip;        // Wall: no-slip condition
    }
    frontAndBack
    {
        type        empty;
    }
}

0/p — Pressure Boundary Conditions#

Pressure here is kinematic. Only the outlet is fixed at 0.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
 
dimensions      [0 2 -2 0 0 0 0];  // Kinematic pressure (p/rho), m²/s²
 
internalField   uniform 0;
 
boundaryField
{
    inlet
    {
        type        zeroGradient;
    }
    outlet
    {
        type        fixedValue;
        value       uniform 0;  // Reference pressure of 0 at outlet
    }
    wall
    {
        type        zeroGradient;
    }
    frontAndBack
    {
        type        empty;
    }
}

0/k — Turbulent Kinetic Energy#

A 5% inlet intensity gives $k = 1.5 \times (0.05)^2 = 0.00375$ , the value written in the file.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      k;
}
 
dimensions      [0 2 -2 0 0 0 0];  // m²/s²
 
// Inlet turbulence intensity 5%, velocity 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;
    }
}

Step 4: Transport Properties and Turbulence Model Setup#

constant/transportProperties#

ν lives in this file. Whether the H that produces Re 1000 matches the mesh height is checked after the results table.

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

constant/turbulenceProperties#

RAS stays on, with kEpsilon.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      turbulenceProperties;
}
 
simulationType  RAS;
 
RAS
{
    RASModel        kEpsilon;   // k-epsilon two-equation model
    turbulence      on;
    printCoeffs     on;
}

Step 5: Numerical Schemes and Solvers Setup#

system/fvSchemes#

Time is steadyState, momentum convection is linearUpwind, and k and epsilon stay upwind.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
 
ddtSchemes      { default steadyState; }  // Steady-state analysis
 
gradSchemes
{
    default         Gauss linear;
    grad(U)         cellLimited Gauss linear 1;  // Gradient limiter
}
 
divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);  // 2nd-order accuracy
    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; }           // Pressure relaxation factor (to prevent divergence)
    equations
    {
        U       0.7;
        k       0.7;
        epsilon 0.7;
    }
}

Step 6: Run Analysis#

Pressure relaxation 0.3, and 0.7 on U, k, and epsilon, is the coupling that carries the log to iteration 500 with residuals under 1e-4. How that pressure interpolation hides a checkerboard is written up separately in SIMPLE and Rhie–Chow.

# Check mesh quality (optional but recommended)
checkMesh
 
# Run simpleFoam
simpleFoam 2>&1 | tee log.simpleFoam

Convergence log:

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

The stop is every residual under 1e-4. This log crosses that line at iteration 500.

Step 7: Results Visualization (ParaView)#

paraFoam opens the case.

# Create dummy file for ParaView and execute
touch backwardFacingStep.foam
paraFoam

In ParaView the order is this.

Under Coloring, pick U and paint Magnitude.
Glyph sprays the vectors so the reverse pocket shows up.
Stream Tracer draws the streamlines.

Once the fields are VTK, this command is enough for the numbers.

foamToVTK -latestTime

The x where near-wall ux changes sign is taken as reattachment. The script hard-wires H at 0.01 m.

import numpy as np
 
# Read bottom wall (y=0) data from foamToVTK output file
# (Adjust path according to your environment)
data = np.loadtxt("VTK/backwardFacingStep_500/internal.vtk",
                  skiprows=10)  # Exclude VTK header
 
x = data[:, 0]
ux = data[:, 3]  # x-direction velocity
 
# Filter cells near the bottom wall (y ≈ 0)
wall_mask = data[:, 1] < 0.005
x_wall = x[wall_mask]
ux_wall = ux[wall_mask]
 
# Position of sign change = reattachment point
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  # Step height (m)
    print(f"Reattachment Length: {x_reattach/H:.1f}H")

Common Errors and Solutions#

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

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

0/epsilon is missing. Build one on the same skeleton as 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: Analysis Divergence (Floating point exception)#

Floating point exception (core dumped)

Too much relaxation, or a coarse mesh, and the run blows up. Lower the factors in fvSolution.

relaxationFactors
{
    fields      { p 0.2; }   // Decrease 0.3 → 0.2
    equations   { U 0.5; }   // Decrease 0.7 → 0.5
}

Error 3: Residuals Oscillating Instead of Decreasing#

nNonOrthogonalCorrectors is 0 while the mesh is far from orthogonal, and the residual wobbles. Raise the correction count.

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

Expected Results Summary#

Item	Value
Reynolds Number ( $Re = UH/\nu$ )	1,000
Reattachment Length ( $x_r / H$ )	6–8
Max Reverse Velocity	$\approx -0.1$ m/s
Convergence Iterations	300–600

Armaly et al. 1983, at $Re = 800$ , sits near 6H, next to the table's 6–8.

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

The step runs from y = 0 to y = 1. Times scale 0.1, H is 0.1 m. The nu comment has $Re = 1 \times 0.01 / 10^{-5} = 1000$ , so its H is 0.01 m. The table's 1000 is that comment. Mesh height gives $1 \times 0.1 / 10^{-5} = 10^{4}$ . Then 6–8H is 0.6–0.8 m, or 0.06–0.08 m.

Downstream x is 5 to 20. Face at 0.5 m, outlet at 2.0 m, gap 1.5 m, which is 15 mesh heights, so 0.6–0.8 m fits. By the comment's H the gap is 150H, and 6H is 0.06 m past the face, at x = 0.56 m. The script uses H = 0.01 m too. A bottom cell is $(2/20) \times 0.1 = 0.01$ m. y under 0.005 keeps half of it.

k matches $1.5 \times (0.05)^2 = 0.00375$ . With $C_\mu = 0.09$ , $C_\mu^{3/4} k^{3/2} \approx 3.77 \times 10^{-5}$ , and over 0.0009 the length is 0.042 m. $L \approx 0.1H$ is 0.01 m or 0.001 m, so epsilon would be 0.00377 or 0.0377, not the file's 0.0009. Reverse flow near $-0.1$ m/s is 10% of 1 m/s, twice the 0.05 m/s of 5% intensity.

The Step Height That Stretches 6H by Ten#

Mesh H makes the table's Re 1000 false. nu must become $10^{-4}$ . The comment's H forces scale down to 0.01, and 6H shrinks to 0.06 m on the 1.5 m run.

Which of those is H? If inlet 1000 is a pipe-transition Reynolds number, is k-epsilon already too heavy, or does the step separate first? Still open in where a pipe stream breaks up.