CFD Training Roadmap -- Learning Sequence and Progression Guide

Introduction

This note provides a structured learning pathway for CFD using OpenFOAM, progressing from fundamental theory to practical solver selection, targeting PhD-level researchers who must understand what the solver does internally.

The sequence assumes familiarity with continuum mechanics and numerical methods. Skip Governing Equations if starting from zero.

Learning Sequence

The recommended progression is Governing Equations, FVM, Discretization Schemes, Case Setup, Solver Selection, and Post-Processing. Each stage builds on the previous. Skipping steps causes problems—you might write a C++ program that compiles without understanding pointers, but you will not know why it works.

Stage 1: Governing Equations

Understand the physics before touching OpenFOAM. The Navier-Stokes equations explain why your simulation crashes at t = 0.03 seconds.

The Reynolds-averaged Navier-Stokes (RANS) equations form the backbone of most engineering simulations:

$$\frac{\partial \bar{u}_i}{\partial x_i} = 0$$

$$\frac{\partial (\rho \bar{u}_i)}{\partial t} + \frac{\partial (\rho \bar{u}_i \bar{u}_j)}{\partial x_j} = -\frac{\partial \bar{p}}{\partial x_i} + \frac{\partial}{\partial x_j}\left[\mu_{eff}\left(\frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i}\right)\right]$$

Review Governing Equations for full derivations including the energy equation and turbulent closure models.

Stage 2: Finite Volume Method

FVM conserves fluxes at each cell face. Use FVM with a mass source term and check the outlet flux. The integral form over a control volume ensures exact conservation. The divergence theorem converts surface integrals of fluxes into volume integrals of sources. OpenFOAM uses the fvSchemes file to control discretization of each term.

See FVM for discrete operator treatment covering divergence, gradient, and Laplacian schemes.

Stage 3: Discretization Schemes

Spatial discretization determines accuracy and stability. Second-order upwind is the industrial CFD workhorse. First-order upwind is the default, making everyone produce wrong results in predictable ways.

Higher-order schemes capture sharper gradients, benefiting shock fronts and boundary layers. They can also produce oscillations. First-order schemes are diffusive and smear everything. See Numerical Schemes for formal stability analysis.

Stage 4: Case Setup

Case setup in OpenFOAM follows the directory structure paradigm. Each file has a purpose: boundaryField in 0/U or 0/p defines conditions at each patch, constant/fvSolution contains solver configurations. Understanding this map prevents getting lost. See Case Setup for a complete walkthrough.

Stage 5: Solver Selection

OpenFOAM solcers fall into categories by physics: incompressible (simpleFoam, pimpleFoam), compressible steady (rhoSimpleFoam), compressible PIMPLE (rhoPimpleFoam), compressible central-upwind (rhoCentralFoam for high-speed shocks), and two-phase (interFoam, multiphaseInterFoam). See Solver Selection for decision trees.

Stage 6: Post-Processing

After a solver finishes, analyze results using the postProcess utility and paraFoam for visualization. Common tasks include surface force and pressure integration, profile extraction with sample, field visualization via .vtk export for ParaView, and Q-criterion vortex identification. See Post-Processing for a detailed workflow.

Mesh Effects on Gradient Computation and Convergence

Mesh quality determines whether your simulation converges or produces wrong results that can look plausible. Structured hexahedral meshes use 30-50% fewer cells than unstructured tetrahedral for the same accuracy, with better gradient accuracy, skewness tolerance, boundary-layer resolution, and parallel scaling. Hexahedral inflation creates prism layers for boundary layers. Tetrahedral elements near walls produce poor gradient resolution, yielding unreliable results. Hexahedral mesh generation is slower (manual), while tetrahedral auto-mesher is faster.

Non-Orthogonality and Pressure Convergence

Non-orthogonal meshes require correction iterations. The pressure equation correction depends on the non-orthogonalCorrectors parameter in fvSolution:

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-07;
        relTol          0.1;
        nNonOrthogonalCorrectors 2;
    }
}

Set nNonOrthogonalCorrectors to 3 or 4 for highly non-orthogonal meshes. Use an orthogonality value below 0.3 as the threshold. Set the value to 0 or 1 for near-orthogonal meshes. Use 2 as the default for engineering-quality meshes.

Each additional non-orthogonal corrector adds 10-20% to the wall-clock time per iteration. The question is whether that cost improves accuracy enough. Most industrial cases do not justify the cost. You will discover this when your mesh has 0.15 orthogonality and you leave nNonOrthogonalCorrectors at 0.

Solver Accuracy Comparison: Sod Shock Tube

The Sod shock tube problem is the canonical test for compressible solvers. A diaphragm separates high-pressure gas from low-pressure gas. State 1 has p1 = 10^5 Pa and rho1 = 1.0 kg/m^3. State 2 has p2 = 10^4 Pa and rho2 = 0.125 kg/m^3. The diaphragm ruptures at t = 0. The flow reaches a Mach 2.5 shock. The shock propagates into the low-pressure region.

rhoCentralFoam uses a central-upwind flux scheme from Kurganov and Tadmor. This scheme captures shocks well. rhoPimpleFoam handles transonic cases. You must select deltaT carefully to avoid CFL-based timestep crashes. rhoSimpleFoam is a steady-state solver. It cannot capture transient shock phenomena.

Typical rhoCentralFoam configuration:

divSchemes
{
    default         none;
   (phi,U)         flux(U,phi) limitedLinearV 1;
   (phi,rho)       flux(rho,phi) bounded limitedLinearV 1;
   (phi,k)         flux(k,phi) bounded limitedLinearV 1;
   (phi,omega)     flux(omega,phi) bounded limitedLinearV 1;
   (phi,E)         flux(E,phi) limitedLinearV 1;
   ((nabla|U),rho) limitedLinearV 1;
   ((nabla|U),rhoE) limitedLinearV 1;
   (((rho*g)*h)Snb) Gauss linear;
   ("(nabla|U).*") Gauss linear;
   ("(nabla|T).*") Gauss linear;
   ("(nabla|k).*") Gauss linear;
   ("(nabla|omega).*") Gauss linear;
   ("(nabla|e).*") Gauss linear;
    "(nabla*p)"     Gauss linear;
    "(nabla*rho)"   Gauss linear;
    "(nabla|fv)"    Gauss linear;
    "(nabla|fvE)"   Gauss linear;
    default         none;
}

Under-Relaxation Behavior

Under-relaxation factors (URFs) control how much of the new solution the solver accepts. The solver retains the old solution otherwise. OpenFOAM uses under-relaxation to prevent the solver from handling sudden changes.

Typical URFs for simpleFoam (steady, incompressible):

The under-relaxation formula:

$$\phi^{new} = \alpha \phi_{calculated}^{new} + (1 - \alpha)\phi^{old}$$

where alpha is the under-relaxation factor. A lower alpha retains more of the old solution. Lower alpha trades convergence speed for stability.

Start with conservative URFs. Set p to 0.2 and U to 0.7. Observe residuals. Increase p toward 0.3 for faster convergence. Do not change all URFs simultaneously. That wastes 12 hours of wall clock time on debugging.

Tolerance and Convergence Criteria

Monitoring convergence requires looking at multiple indicators:

Residuals alone are insufficient. A residual of 1e-6 does not prove correctness. It only proves that the solver stopped changing the solution at that rate. Monitor physical quantities alongside residuals. Track forces, mass flow rates, and temperature differences.

High-Speed Aerodynamics: Steady vs Unsteady vs Local Time-Stepping

For high-speed external aerodynamics:

Local time-stepping uses a CFL-based local timestep. The adjustTimeStep and maxDeltaT framework controls it. The method advances toward a steady state. Use this method when the final state is steady but the path to reach it involves multiple time scales. The simulation reaches the steady solution. The journey to that solution does not matter.

Learning Progression Pathway

The recommended progression mirrors increasing physical complexity:

1. **Laminar flow** -- Solve laminarFoam cases. Understand the exact solution validity. 2. **Turbulent flow** -- Transition to simpleFoam with kOmegaSRT or kEpsilon. Compare with DNS and experimental data. 3. **Multiphase flows** -- Study interFoam for VOF and reactingTwoPhaseEulerFoam for Euler-Euler. 4. **Compressible flow** -- Begin with rhoSimpleFoam for subsonic and transonic cases. Then use rhoCentralFoam for shock-capturing. 5. **Combustion** -- Use reactive solvers reactingFoam and rhoChemFoam. Combustion runs can produce fire instead of convergence.

Validate each stage with at least one case study. Compare against analytical solutions like Poiseuille flow and Couette flow. Compare against experimental data like the Garnier DNS database.

References

• Versteeg, H.E. and Malalasekera, W. An Introduction to Computational Fluid Dynamics: The Finite Volume Method. 2nd Edition, Pearson, 2007.
• Hirsch, C. Numerical Computation of Internal and External Flows. 2nd Edition, Wiley, 2007.
• Ferziger, J.H. and Peric, M. Computational Methods for Fluid Dynamics. 3rd Edition, Springer, 2002.
• Jansson, N. and Nikishkov, Y. et al. OpenFOAM Extension Releases User Guide, 2024.
• Weller, H.G., Tabor, G.R., Jasak, H. and Fureby, C. "A Tensorial Approach to Computational Modelling Using Object-Oriented Finite Volume Methods," Computational Physics, Vol. 12, 1998.
• Anderson, J.D. Computational Fluid Dynamics: The Basics with Applications. McGraw-Hill, 1995.
• Panton, R.L. Incompressible Flow. 4th Edition, Wiley, 2013.