ML-Enhanced Turbulence Modeling in OpenFOAM

Study the Reynolds-Averaged Navier-Stokes equations and understand the closure problem: the Reynolds stress tensor τ_ij = -ρ⟨u_i'u_j'⟩ contains 6 unknowns that must be modeled. RANS models (k-ω SST, etc.) use the Boussinesq hypothesis: τ_ij ≈ 2μ_t S_ij - (2/3)kδ_ij, where S_ij is the strain rate tensor and μ_t is the eddy viscosity.

The Boussinesq hypothesis is exact only for simple shear flows. Read Pope (2000) "Turbulent Flows" Chapter 11 for the theoretical limitations, and Weatheritt & Sandberg (2016) for the data-driven correction approach you'll implement.

Start with fully developed turbulent channel flow at Re_τ = 395 — the most studied turbulent flow, with high-quality DNS reference data from Moser et al. (1999). Run OpenFOAM k-ω SST RANS and a wall-resolved LES with the Smagorinsky subgrid model on the same geometry.

Compare the mean velocity profile, Reynolds shear stress, and turbulent kinetic energy profiles between RANS, LES, and DNS reference data. Quantify the RANS errors — this is what the ML model must learn to correct. Use the built-in OpenFOAM postprocessing utilities (turbulenceFields, fieldAverage) to extract time-averaged statistics from the LES.

The target for the ML model is the discrepancy: Δτ_ij = τ_ij^LES - τ_ij^RANS, computed field by field at every grid point. Post-process both simulations in Python using the fluidfoam library to load OpenFOAM field data into NumPy arrays. Compute the full Reynolds stress tensor from LES data and compare against RANS.

Decompose the discrepancy into: eddy viscosity correction (scalar) and Reynolds stress anisotropy correction (traceless symmetric tensor). The anisotropy tensor is a better ML target because it separates the magnitude correction (easy) from the directional correction (hard and physically meaningful).

The ML model must predict the stress correction from local flow quantities that are available during RANS simulation. Standard feature choices (following Ling & Templeton 2015): the 10 invariants of the strain rate tensor S_ij, vorticity tensor Ω_ij, and their products; dimensionless pressure gradient; dimensionless wall distance; turbulence intensity k/||U||².

These inputs are invariant to rotation and reflection — an important physical constraint called Galilean invariance. Using frame-invariant features ensures the model works regardless of how the coordinate system is oriented relative to the flow.

Train a neural network mapping [local flow features] → [Reynolds stress anisotropy correction]. Use the LES-derived correction field as training labels. Architecture: 4–6 fully connected layers with 50–100 neurons, tanh activation. Train on channel flow data and validate on a different channel flow at a different Reynolds number (Re_τ = 550).

Alternatively, implement the Tensor Basis Neural Network (TBNN) of Ling et al. (2016) — a physics-constrained architecture that guarantees the output satisfies the symmetry and trace conditions of the Reynolds stress tensor by construction, significantly improving generalization.

Integrate the trained correction into OpenFOAM using the Python runtime interface or by compiling a custom turbulence model in C++ that calls the neural network at each iteration. Run the corrected RANS solver on the backward-facing step (a harder flow not in the training data) and compare against LES reference data.

Quantify the improvement: does the ML-corrected RANS capture the reattachment length better than standard k-ω SST? Does it improve velocity profiles in the separated region? Create a comparison figure showing DNS/LES reference, standard RANS, and ML-corrected RANS — this is the standard figure in data-driven turbulence modeling papers.