A modern RANS k-ω turbulence solver for DOLFINx/FEniCSx 0.10+.
A finite element RANS solver implementing the Wilcox (2006) standard k-ω and Menter (1994) k-ω SST turbulence models. Use it as-is for turbulent channel flow, or as a validated foundation for extending to other geometries.
- Why This Exists
- Quick Start
- Features
- Installation
- Configuration Reference
- Governing Equations
- Numerical Method
- Stabilization: What Makes FEM RANS Hard
- Architecture
- Validation Against Nek5000
- Extending to Other Geometries
- Python API
- Output Files
- References
- License
DOLFINx is one of the most capable open-source finite element frameworks — automated variational forms, parallel assembly, flexible function spaces. Yet RANS turbulence solvers for it are scarce. Users who need turbulence modeling often start from scratch, re-discovering the same FEM-specific pitfalls (positivity violations, ν_t feedback instability, convection oscillations) that make RANS in finite elements harder than in finite volumes.
This project provides a working, validated k-ω RANS solver that exposes both the mathematics and the engineering decisions needed to make FEM RANS converge.
RANS turbulence models solve time-averaged equations with a turbulent viscosity
| Model | Solves for | Wall behavior |
|---|---|---|
| k-ε | k (TKE), ε (dissipation) | Needs wall functions or damping |
| k-ω | k (TKE), ω (specific dissipation) | Natural wall BC: |
| k-ω SST | Blends k-ω (wall) with k-ε (freestream) | Best of both worlds |
Decision: We chose k-ω because it has the simplest wall boundary conditions of any two-equation model. The asymptotic wall value
Limitation: Standard k-ω is sensitive to freestream
Turbulent channel flow is the canonical RANS validation case because it removes all geometric ambiguity: one inlet condition (periodic), one driving mechanism (body force), and an exact momentum balance (
Decision: We validate against Nek5000's RANS channel tutorial rather than DNS because the goal is cross-code RANS verification (do two codes with the same model produce the same answer?), not model-vs-physics validation (does RANS match DNS?). Nek5000 uses a well-tested finite-volume-style RANS implementation, so agreement confirms our FEM implementation is correct.
- DOLFINx users who need RANS turbulence and want a working starting point
- Researchers exploring FEM-specific turbulence numerics
- Students learning how RANS models are implemented in variational form
# 1. Activate FEniCSx environment
conda activate fenicsx
# 2. Install dolfinx-rans
cd dolfinx-rans && uv pip install -e .
# 3. Run canonical case (Re = 100,000 channel flow)
./run_channel.sh
# Results appear in channel/results/:
# final_fields.png — Final field/profile composite
# fields0000000.png — Channel snapshot at logged iterations
# convergence.png — Live residual and dt history
# history.csv — Per-iteration residual/metrics log
# profiles.csv — Wall-normal benchmark profile
#
# For BFS case outputs and one canonical file only:
# ./run_bfs.sh
# bfs/results/
# bfs_fields0000000.png — BFS field snapshots (iteration tagged)
# bfs_cf0000000.png — BFS wall-friction snapshots (iteration tagged)
# bfs_profiles0000000.png — BFS profile snapshots (iteration tagged)The runner supports MPI parallelism:
./run_channel.sh 4 # 4 MPI processes
./run_channel.sh 8 path/to/config.jsonc # 8 processes, custom configCase layout is intentionally strict:
channel/channel.jsoncandbfs/bfs.jsoncare the working canonical configschannel/results/andbfs/results/are the only case outputs
Clean outputs before a fresh rerun:
rm -rf channel/results bfs/results Turbulent Channel Flow Setup (nondimensional)
y
^
y = 2δ | top wall: u = 0, k = 0, ω = ω_wall
| ┌──────────────────────────────────────────────┐
x = 0 | │ → → → → → → → → → → → → → │ x = Lx
(periodic) | │ │ (periodic)
↔ | │ ------------- centerline: y = δ --------- │ ↔
| │ → → → → → → → → → → → → → │
y = 0 | └──────────────────────────────────────────────┘
| bottom wall: u = 0, k = 0, ω = ω_wall
+------------------------------------------------------> x
streamwise driving term: f_x = +1
Physics
- Standard k-ω model (Wilcox 2006) with cross-diffusion term
- k-ω SST model (Menter 1994) with F1/F2 blending functions
- Optional Durbin realizability limiter for stagnation regions
- Wilcox 2006 stress limiter ($\tilde{\omega} = \max(\omega, C_\text{lim} |S| / \sqrt{\beta^*})$)
Numerics
- Pseudo-transient continuation to steady state
- Adaptive time stepping with hysteresis to prevent oscillation
- Picard iteration with under-relaxation for k, ω, and ν_t
- IPCS fractional-step pressure-velocity coupling
- Wall-refined meshes with geometric or tanh stretching
- Half-channel (symmetry) or full-channel mode
Infrastructure
- Periodic boundary conditions via dolfinx_mpc
- Config-driven JSON workflow — one file controls everything
- Profile CSV export for cross-code benchmark comparison
- Config-driven regression gates (bulk velocity, wall shear, profile RMSE)
- VTX snapshots for ParaView time series
- Live convergence plots during iteration
- MPI-parallel execution
DOLFINx and dolfinx_mpc require conda-forge (not pip-installable). After setting up the FEniCSx environment, install dolfinx-rans with uv or pip.
# Create dedicated environment with DOLFINx 0.10.0+
conda create -n fenicsx python=3.11
conda activate fenicsx
conda install -c conda-forge fenics-dolfinx dolfinx_mpc petsc4py mpi4py# Activate your FEniCSx environment first
conda activate fenicsx
# Install dolfinx-rans (editable for development)
cd dolfinx-rans
uv pip install -e .
# Or install directly from GitHub
uv pip install git+https://github.com/ricardofrantz/dolfinx-rans.gitAlternative: pip install
conda activate fenicsx
pip install -e .
# or
pip install git+https://github.com/ricardofrantz/dolfinx-rans.git| Package | Source | Notes |
|---|---|---|
| DOLFINx 0.10.0+ | conda-forge | FEniCSx finite element library |
| dolfinx_mpc | conda-forge | Multi-point constraints (periodic BCs) |
| petsc4py | conda-forge | PETSc linear solvers |
| mpi4py | conda-forge | MPI parallelization |
| numpy ≥ 1.24 | pip/conda | Numerical arrays |
| matplotlib ≥ 3.7 | pip/conda | Plotting |
Click to expand
For turbulent channel flow, we need periodic boundary conditions in the streamwise direction:
The problem: Standard FEM boundary conditions can only set values (Dirichlet) or flux (Neumann) at boundaries — neither can express "right boundary equals left boundary", which is a constraint between DOFs.
How dolfinx_mpc solves it:
Standard FEM: With MPC:
[K]{u} = {f} [K']{u'} = {f'}
DOFs: [u₀, u₁, ..., uₙ] DOFs: [u₀, u₁, ..., uₘ] where m < n
Right-boundary DOFs eliminated and
replaced by left-boundary DOFs
MPC modifies the stiffness matrix to enforce
Why periodic + body force instead of inlet/outlet?
| Approach | Pros/Cons |
|---|---|
| Inlet/outlet BCs | Unknown inlet profile, entrance effects, long domain needed |
| Periodic + body force | Body force |
A single JSON file controls the entire simulation. Core sections (geom, nondim, turb, solve) are required. The benchmark section is optional and enables regression gates. Keys starting with _ are ignored, allowing JSON comments.
Full example (the canonical channel case, channel/channel.jsonc):
| Parameter | Type | Description |
|---|---|---|
Lx |
float | Channel length (streamwise) |
Ly |
float | Channel height ( |
Nx |
int | Mesh cells in x |
Ny |
int | Mesh cells in y |
mesh_type |
str |
"triangle" or "quad"
|
y_first |
float | First cell height from wall. For low-Re BCs: |
growth_rate |
float | Geometric stretching ratio (>1 for wall refinement, 1 for uniform) |
stretching |
str |
"geometric" or "tanh" (default: "geometric") |
y_first_tol_rel |
float | Relative tolerance for mesh/BC consistency check (hard fail if exceeded) |
use_symmetry |
bool |
true = half-channel with symmetry at top; false = full channel with walls on both sides |
Decision: We provide both geometric and tanh stretching. Geometric stretching (y_{i+1} = r \cdot y_i) is standard in CFD and gives explicit control via growth rate. Tanh stretching (y(\eta) = H[1 - \tanh(\beta(1-\eta))/\tanh(\beta)]) clusters more points near the wall for the same cell count. The solver computes y_first from the mesh geometry at runtime and hard-fails if it disagrees with the config — this catches silent mesh/BC mismatches that would produce wrong
Wall refinement example:
For Re_τ ≈ 1116 (Re = 100,000), need y⁺ < 2.5:
y_first < 2.5 / 1116 ≈ 2.24e-3
With y_first = 1.60e-3 and tanh stretching:
y⁺ ≈ 1.79 ✓
| Parameter | Type | Description |
|---|---|---|
Re_tau |
float | Friction Reynolds number |
use_body_force |
bool |
true = body force |
Decision: Everything is nondimensionalized by
| Parameter | Type | Description | Default rationale |
|---|---|---|---|
model |
str |
"wilcox2006" or "sst"
|
Wilcox 2006 is simpler and sufficient for channel flow |
beta_star |
float | k-ω constant |
Standard across all k-ω variants |
nu_t_max_factor |
float | Max |
Safety cap; see ν_t limiting |
omega_min |
float | ω floor to prevent |
See ω floor |
k_min |
float | k floor for positivity |
1e-10 prevents division by zero |
k_max |
float | k cap for safety | Prevents runaway during early iterations |
C_lim |
float | Durbin limiter: $\nu_t \le C_\text{lim} \cdot k / (\sqrt{6} | S |
| Parameter | Type | Description | Default rationale |
|---|---|---|---|
dt |
float | Initial pseudo-time step | Small for stability during startup transient |
dt_max |
float | Maximum dt | Prevents excessively large steps |
dt_growth |
float | dt multiplier when converging | See adaptive dt |
dt_growth_threshold |
float | Hysteresis: only grow dt if |
0.8 means residual must be decreasing steadily |
t_final |
float | Max pseudo-time (safety limit) | Usually not reached — convergence check stops iteration |
max_iter |
int | Max outer iterations | Safety limit |
steady_tol |
float | Convergence tolerance on |
1e-3 is typical for RANS |
picard_max |
int | Inner Picard iterations per time step | See Picard iteration |
picard_tol |
float | Picard convergence tolerance | Inner loop exits early if k change is below this |
under_relax_k_omega |
float | Under-relaxation for k and ω | See feedback loop |
under_relax_nu_t |
float | Under-relaxation for ν_t |
0.4 is aggressive but needed for stability |
log_interval |
int | Print diagnostics every N iterations | |
snapshot_interval |
int | Save VTX + PNG every N iterations (0 = disabled) | |
out_dir |
str | Output directory for results | |
min_iter |
int | Minimum iterations before accepting steady convergence | |
min_dt_ratio |
float | Convergence is blocked until dt >= min_dt_ratio * dt_initial
|
Decision: Under-relaxation factors of 0.6 (k, ω) and 0.4 (ν_t) were determined empirically. Higher values cause oscillation in the ν_t feedback loop; lower values converge but take more iterations. The ν_t under-relaxation is more aggressive than k/ω because ν_t amplifies small changes through the momentum equation.
| Parameter | Type | Description |
|---|---|---|
gate_u_bulk_bounds |
[min, max] |
Acceptable range for final |
gate_tau_wall_bounds |
[min, max] |
Acceptable range for final |
reference_profile_csv |
str | Path to Nek reference CSV (columns: y_over_delta,u_over_ubulk) |
u_plus_rmse_max |
float | Max RMSE for velocity profile vs reference |
k_plus_rmse_max |
float | Max RMSE for TKE profile vs reference |
All equations are nondimensionalized by the half-channel height
Dropping the
The body force
The Reynolds-averaged Navier-Stokes momentum equation with the Boussinesq eddy-viscosity hypothesis:
where
Turbulent kinetic energy (k):
Specific dissipation rate (ω):
where:
-
$P_k = \nu_t S^2$ is the TKE production, with$S^2 = 2 S_{ij} S_{ij}$ and$S_{ij} = \frac{1}{2}(\partial u_i / \partial x_j + \partial u_j / \partial x_i)$ - The last term in the ω-equation is the Wilcox 2006 cross-diffusion term, which reduces freestream sensitivity
Decision: The ω-equation production term uses
| Constant | Symbol | Value | Notes |
|---|---|---|---|
| k destruction | 0.09 | Standard across all k-ω variants | |
| ω destruction | 0.0708 | Updated from 0.075 (1998) | |
| k diffusion | 0.6 | Updated from 0.5 (1998) | |
| ω diffusion | 0.5 | ||
| von Kármán | 0.41 | ||
| Production ratio |
|
Chosen to yield correct log-layer slope | |
| Cross-diffusion | 0.125 (= 1/8) | New in Wilcox 2006 | |
| Stress limiter | 0.875 (= 7/8) | New in Wilcox 2006 |
Decision: We use the Wilcox 2006 constants (not 1998). The 2006 revision added the cross-diffusion term
The k-ω SST (Shear Stress Transport) model blends Wilcox k-ω near walls with transformed k-ε in the freestream using a blending function
where
Inner-layer (k-ω) constants:
Outer-layer (k-ε transformed) constants:
The cross-diffusion term for SST uses the blending:
Wilcox 2006 with stress limiter:
The stress limiter prevents excess
SST with shear-stress transport limiter:
where
Walls (
| Variable | BC | Value |
|---|---|---|
| Dirichlet (no-slip) | ||
| Dirichlet | ||
| Dirichlet |
where
Symmetry plane (
-
$v = 0$ (no wall-normal velocity) -
$\partial u / \partial y = 0$ (natural Neumann, not explicitly imposed) -
$\partial k / \partial y = 0$ (natural Neumann) -
$\partial \omega / \partial y = 0$ (natural Neumann)
Streamwise (
- Periodic via dolfinx_mpc multi-point constraints
Decision: The
The equations are implemented in DOLFINx's UFL (Unified Form Language). Here is how the k-equation maps from math to code.
Mathematical form (semi-discrete in time):
UFL implementation (from solver.py):
# Diffusion coefficient
D_k = nu_c + sigma_k_c * nu_t_
# Production: P_k = ν_t * 2 * S_ij * S_ij
S_tensor = sym(grad(u_n))
S_sq = 2.0 * inner(S_tensor, S_tensor)
P_k = nu_t_ * S_sq
# Weak form: find k_trial such that F_k(k_trial, phi_k) = 0
F_k = (
(k_trial - k_n) / dt_c * phi_k * dx # time derivative
+ dot(u_n, grad(k_trial)) * phi_k * dx # convection
+ D_k * inner(grad(k_trial), grad(phi_k)) * dx # diffusion (integrated by parts)
+ beta_star_c * omega_safe * k_trial * phi_k * dx # destruction
- P_k * phi_k * dx # production (RHS)
)The ω-equation follows the same pattern with an additional cross-diffusion source term:
# Cross-diffusion: σ_d/ω * max(∇k·∇ω, 0) — explicit treatment
grad_k_dot_grad_w = dot(grad(k_n), grad(omega_n))
grad_kw_positive = ufl.conditional(
ufl.gt(grad_k_dot_grad_w, 0.0), grad_k_dot_grad_w, 0.0
)
cross_diff = sigma_d_c / omega_safe * grad_kw_positiveDecision: The convection term lhs() and rhs() from UFL and solve a single linear system per equation per Picard iteration. Fully implicit treatment would require a nonlinear solver and is not worth the complexity for pseudo-transient RANS.
The Incremental Pressure Correction Scheme (IPCS) splits the Navier-Stokes equations into three sequential linear solves per time step:
Step 1 — Tentative velocity (momentum without pressure gradient correction):
where
Step 2 — Pressure correction (enforce divergence-free velocity):
Step 3 — Velocity correction (project onto divergence-free space):
Decision: IPCS was chosen over monolithic (fully coupled) Navier-Stokes because it reuses the same matrix structure across iterations (only the RHS changes when
We do not seek time accuracy — the
The outer loop structure is:
for step in 1..max_iter:
for picard_iter in 1..picard_max:
solve momentum (IPCS steps 1-3)
under-relax velocity: u_n = α·u_new + (1-α)·u_old
solve k-equation, clip, under-relax
solve ω-equation, clip, under-relax
update ν_t with under-relaxation
check Picard convergence
compute outer residual
adapt dt based on residual trend
check steady-state convergence
Within each pseudo-time step, the nonlinear coupling between momentum and turbulence is resolved by Picard (fixed-point) iteration. Each Picard iteration:
- Solves momentum with current
$\nu_t$ (IPCS steps 1-3) - Under-relaxes velocity: $\mathbf{u}^n \leftarrow \alpha_u \mathbf{u}\text{new} + (1 - \alpha_u) \mathbf{u}\text{old}$
- Updates SST blending functions (if using SST model)
- Solves k with current velocity and
$\omega$ - Clips k to
$[k_\text{min}, k_\text{max}]$ , under-relaxes - Solves ω with current velocity and k
- Clips ω to
$[\omega_\text{min}, \omega_\text{max}]$ , under-relaxes - Computes
$\nu_t$ from updated k and ω, under-relaxes
The inner loop exits early if the relative change in k drops below picard_tol.
Decision: Under-relaxation is applied at three levels: velocity (under_relax_k_omega = 0.6), and ν_t (under_relax_nu_t = 0.4). The ν_t relaxation is the most critical because ν_t couples all equations — a 10% change in ν_t causes ~10% change in effective viscosity everywhere, which changes production, which changes k and ω, which changes ν_t again. Without under-relaxation, this feedback loop diverges.
Each linear subsystem uses a tuned PETSc solver:
| System | KSP | Preconditioner | Rationale |
|---|---|---|---|
| Momentum (step 1) | BiCGStab | Jacobi | Non-symmetric system; Jacobi is cheap and MPI-friendly |
| Pressure (step 2) | CG | Hypre BoomerAMG | SPD Laplacian; AMG gives mesh-independent convergence |
| Velocity correction (step 3) | CG | Jacobi | Mass matrix; cheap and well-conditioned |
| k-equation | BiCGStab | Hypre BoomerAMG | Convection-diffusion; AMG handles anisotropic meshes |
| ω-equation | BiCGStab | Hypre BoomerAMG | Same as k |
All solvers use rtol = 1e-8.
Decision: We use algebraic multigrid (BoomerAMG) for the pressure Poisson and turbulence transport equations because wall-refined meshes create highly anisotropic cells. Direct solvers or simple preconditioners (Jacobi, ILU) struggle with the resulting ill-conditioning. AMG automatically constructs a hierarchy that handles the anisotropy. Note: ILU is incompatible with MPI parallelism; BoomerAMG works in parallel out of the box.
Most RANS solvers use finite volume methods (FVM), which naturally suit conservation laws. FEM offers advantages — higher-order accuracy, complex geometries, rigorous error analysis — but faces specific challenges with turbulence. This section documents what we tried, what failed, and what works.
The problem:
What happens without it: Negative
Solution: Hard clipping after each linear solve, before under-relaxation:
k_.x.array[:] = np.clip(k_.x.array, k_min, k_max)
omega_.x.array[:] = np.clip(omega_.x.array, omega_min, omega_max_limit)where omega_max_limit = 10 * omega_wall scales the cap with the current case.
Decision: Clipping is mathematically inelegant (it violates the variational principle) but pragmatically essential. More sophisticated approaches exist — DG with slope limiters, variational inequalities — but they add enormous complexity. Clipping works reliably and the clipped values are small perturbations on the converged solution. The key insight is that k_min = 1e-10 (not 0) prevents 0/0 = NaN, and the ω floor has a much larger impact on solution quality (see below).
The problem: At high Re, the k and ω transport is dominated by convection (
FVM handles this naturally with upwind schemes that add stabilizing numerical diffusion in the flow direction.
FEM workarounds:
- SUPG stabilization — adds artificial streamline diffusion (complex to tune, especially with turbulence)
-
Fine mesh — keep local Péclet number
$Pe_h = |u| h / (2D)$ moderate -
Pseudo-transient stepping — the
$k/\Delta t$ mass-matrix term acts as isotropic numerical diffusion
Decision: We rely on mesh refinement and pseudo-transient damping rather than SUPG. This is simpler to implement and debug, and the wall-refined mesh already provides the resolution where gradients are steepest. SUPG would help with coarser meshes but introduces its own tuning parameters and can interfere with turbulence production near walls.
The problem:
ω ↑ → ν_t ↓ → less diffusion → steeper gradients → more production → k ↑ → ν_t ↑ → ...
What failed: No under-relaxation → oscillation and divergence within 5-10 iterations. Under-relaxing only k/ω but not ν_t → still oscillates because ν_t amplifies the changes. Using the same relaxation factor for all fields → either too aggressive (oscillates) or too conservative (takes thousands of iterations).
What works: Three-level under-relaxation with different factors:
# Level 1: velocity (mild)
u_n = 0.7 * u_solved + 0.3 * u_old
# Level 2: turbulence quantities (moderate)
k_new = 0.6 * k_solved + 0.4 * k_old
w_new = 0.6 * w_solved + 0.4 * w_old
# Level 3: eddy viscosity (aggressive)
nu_t_new = 0.4 * nu_t_computed + 0.6 * nu_t_oldDecision: The ν_t under-relaxation factor (0.4) is the most critical. This means each iteration only moves ν_t 40% of the way toward its "correct" value, which slows convergence but prevents the feedback loop from overshooting. The factor was determined empirically: 0.5 oscillates on the Re = 100,000 case, 0.3 is safe but slow. 0.4 is the sweet spot for this class of problems.
The problem: Fixed
What works: Hysteresis-based adaptation. The dt only grows when the residual ratio
residual_ratio = residual / residual_prev
if residual_ratio < dt_growth_threshold: # residual clearly decreasing
dt = min(dt * dt_growth, dt_max)
elif residual_ratio > 1.0 / dt_growth_threshold: # residual clearly increasing
dt = max(dt / dt_growth, dt_initial)
# Otherwise: hold dt steady (hysteresis band)Decision: The dt_growth_threshold = 0.8 means dt only grows when the residual dropped by at least 20%. The dead band between 0.8 and dt_growth = 1.05 is conservative (5% increase per step) to avoid overshooting the stability boundary.
The problem: Where ω is small (channel core),
Solution: Two complementary limits:
-
ω floor (
omega_min): Prevents$\omega$ from dropping below a minimum value. This directly limits$\nu_t$ in the channel core. -
ν_t cap (
nu_t_max_factor): Hard cap at$\nu_t \le \text{factor} \cdot \nu$ . Safety net for extreme cases. -
Wilcox 2006 stress limiter (built-in):
$\tilde{\omega} = \max(\omega, C_\text{lim} |S| / \sqrt{\beta^*})$ . This limits ν_t in regions of high strain rate without affecting the log layer.
Decision: The omega_min = 1.0 in the canonical config is a compromise. Higher values (10-100) produce a better C_lim = 0.0 in config) because it can interact with the ω floor in unexpected ways.
Low-Re wall treatment (direct integration to the wall, no wall functions) requires the first cell center in the viscous sublayer:
For the canonical case (
This forces thin wall cells → high aspect ratios (>100:1 on a 192 × 166 mesh) → potential conditioning issues in the linear solvers. AMG preconditioning handles this well, but simpler preconditioners (ILU, Jacobi) may struggle.
Decision: Tanh stretching concentrates more cells near the wall than geometric stretching for the same Ny, achieving lower y_first_tol_rel) catches cases where the config specifies one y_first but the mesh generator produces a different one — this mismatch would silently corrupt the
dolfinx-rans/
├── src/dolfinx_rans/
│ ├── solver.py — Core solver: equations, BCs, IPCS loop, k-ω/SST models
│ ├── cli.py — CLI entry point: config parsing, case setup, post-processing
│ ├── utils.py — Config loading, diagnostics, wall shear stress, CSV logging
│ ├── plotting.py — Mesh, field, convergence, and profile plots
│ └── validation/
│ └── nek_poiseuille_profile.py — Nek5000 profile extraction and comparison
├── bfs/
│ ├── bfs.jsonc — Canonical BFS case (single config file)
│ ├── run.file — Human-readable run transcript
│ └── results/ — BFS outputs (plots, history, snapshots, etc.)
├── channel/
│ ├── channel.jsonc — Canonical channel case configuration
│ ├── run.file — Human-readable run transcript
│ └── results/ — Channel outputs (plots, history, snapshots, ...)
├── run_channel.sh — One-command runner for channel case
├── run_bfs.sh — One-command runner for BFS case
└── pyproject.toml — Package metadata (hatchling build system)
| Variable | Space | Element | Degree | Rationale |
|---|---|---|---|---|
| Velocity |
V | Lagrange (vector) | 2 | Taylor-Hood P2-P1 for inf-sup stability |
| Pressure |
Q | Lagrange (scalar) | 1 | Taylor-Hood P2-P1 |
|
|
S | Lagrange (scalar) | 1 | P1 matches pressure space; avoids interpolation issues |
Decision: P2-P1 (Taylor-Hood) velocity-pressure is the standard inf-sup stable pair. Using P1 for turbulence quantities (same as pressure) avoids the need to project between different-degree spaces when computing
┌──────────────────────────────────┐
│ Picard Iteration │
│ │
u_n, p_ ─────────────►│ 1. IPCS momentum → u_new │
│ 2. Under-relax u → u_n │
│ 3. [SST: update F1, blend coeffs] │
k_n, omega_n, nu_t_ ──►│ 4. Solve k → k_new │
│ 5. Clip + under-relax k │
│ 6. Solve ω → ω_new │
│ 7. Clip + under-relax ω │
│ 8. Compute ν_t → ν_t_new │
│ 9. Under-relax ν_t │
│ 10. Check Picard convergence │
└──────────────────────────────────┘
│
▼
Outer residual check
Adaptive dt update
Logging / snapshots
Nek5000 is a well-established spectral-element CFD code with a validated RANS implementation. Comparing against Nek5000 is cross-code verification: do two independent implementations of the same turbulence model produce the same answer? This is more rigorous than comparing against DNS, which tests the model itself (and all RANS models disagree with DNS in known ways).
The canonical validation case matches the Nek5000 RANS tutorial:
| Parameter | Value |
|---|---|
| Bulk Reynolds number |
100,000 |
| Friction Reynolds number |
~1,116 |
| Domain |
|
| Mesh | 192 × 166 quads, tanh stretching |
| ~1.79 | |
| Turbulence model | Wilcox 2006 k-ω |
| Driving | Body force |
- Run dolfinx-rans:
./run_channel.sh→ produceschannel/profiles.csv - Extract Nek5000 profile:
python -m dolfinx_rans.validation.nek_poiseuille_profile \ --nek-case-dir path/to/poiseuille_RANS \ --out-dir nek_re100k \ --dolfinx-profiles channel/profiles.csv - Configure regression gate: Set
benchmark.reference_profile_csvto the Nek reference CSV - Re-run with gate:
./run_channel.sh— the post-run check compares profiles automatically
The comparison uses outer scaling (
The run_channel.sh script includes post-run regression checks driven by the benchmark config section:
-
Bulk velocity gate:
$U_\text{bulk}$ must fall withingate_u_bulk_bounds -
Wall shear gate:
$\tau_\text{wall}$ must fall withingate_tau_wall_bounds -
Profile RMSE gate: RMSE of
$U/U_\text{bulk}$ vs reference must be belowu_plus_rmse_max -
TKE RMSE gate: RMSE of
$k$ profile vs reference must be belowk_plus_rmse_max
These gates catch regressions when modifying the solver — if a code change breaks the solution quality, the runner exits with an error.
The solver is structured so the turbulence model and numerical method are geometry-independent. To adapt to a new geometry:
What to change:
-
Mesh: Replace
create_channel_mesh()with your geometry (external mesh import viadolfinx.io.gmshioor similar) -
Boundary conditions: Update
mark_boundaries()and BC definitions for your wall/inlet/outlet/symmetry surfaces -
Wall distance: Replace
compute_wall_distance_channel()(SST only) with a general wall-distance computation - Initial conditions: Provide appropriate ICs for your geometry
-
Body force / driving: Replace
$f_x = 1$ with your driving mechanism (pressure BC, inlet velocity, etc.)
What to keep:
- The IPCS splitting and Picard iteration loop
- Under-relaxation strategy and factors (start with the channel values)
- Clipping and ν_t limiting
- Turbulence model constants
- PETSc solver configuration
Tuning caveat: The under-relaxation factors and ω floor were tuned for channel flow. More complex geometries (recirculation zones, separation, stagnation) may require:
- Lower under-relaxation factors (0.3-0.5 for k/ω instead of 0.6)
- Higher ω floor initially, reduced as the solution develops
- The Durbin stress limiter (
C_lim > 0) for stagnation regions - SST model instead of standard k-ω for freestream-sensitive flows
from dolfinx_rans.solver import (
ChannelGeom,
NondimParams,
TurbParams,
SolveParams,
create_channel_mesh,
solve_rans_kw,
)
# Define geometry (full channel, tanh stretching)
geom = ChannelGeom(
Lx=1.0, Ly=2.0, Nx=192, Ny=166,
mesh_type="quad", y_first=0.001604628,
growth_rate=1.0, stretching="tanh",
y_first_tol_rel=0.2, use_symmetry=False,
)
# Nondimensional parameters
nondim = NondimParams(Re_tau=1115.82, use_body_force=True)
# Turbulence model (Wilcox 2006)
turb = TurbParams(
model="wilcox2006", beta_star=0.09,
nu_t_max_factor=2000.0, omega_min=1.0,
k_min=1e-10, k_max=20.0, C_lim=0.0,
)
# Solver settings
solve = SolveParams(
dt=0.005, dt_max=0.01, dt_growth=1.05,
dt_growth_threshold=0.8, cfl_target=0.25,
t_final=10000.0, max_iter=3000, steady_tol=1e-6,
enable_physical_convergence=False, physical_u_bulk_rel_tol=1e-4,
physical_tau_wall_rel_tol=2.5e-3, physical_convergence_start_iter=10,
picard_max=6, picard_tol=1e-4,
under_relax_k_omega=0.6, under_relax_nu_t=0.4,
log_interval=10, snapshot_interval=50, out_dir="results",
min_iter=50,
)
# Create mesh and run
domain = create_channel_mesh(geom, Re_tau=nondim.Re_tau)
u, p, k, omega, nu_t, V, Q, S, domain, step, t = solve_rans_kw(
domain, geom, turb, solve, Path("results"), nondim,
)The solver can also be run from the command line:
dolfinx-rans channel/channel.jsonc # Serial
mpirun -np 4 dolfinx-rans channel/channel.jsonc # ParallelResults are saved to the case results/ directory specified by solve.out_dir (defaults to "results" in repo configs):
| File | Description |
|---|---|
final_fields.png |
2D contours (u, k, ν_t/ν) + 1D profiles (U/U_bulk, k, ω, ν_t/ν) |
history.csv |
Per-iteration convergence log: residual, dt, U_bulk, τ_wall, field ranges |
profiles.csv |
Wall-normal profile: y, y_over_delta, u, u_over_ubulk, k, omega, nu_t_over_nu |
convergence.png |
Residual and dt history plots |
mesh.png |
Mesh visualization with stretching |
config_used.json |
Exact config snapshot for reproducibility |
run_info.json |
Environment metadata (Python version, git SHA, timestamp) |
snps/*.bp |
VTX time series for ParaView (ADIOS2 format, if snapshot_interval > 0) |
fields0000000.png |
Channel iteration field snapshots (numbered) |
bfs_fields0000000.png |
BFS iteration field snapshots (numbered) |
bfs_cf0000000.png |
BFS wall friction snapshots (numbered) |
bfs_profiles0000000.png |
BFS profiles snapshots (numbered) |
- Wilcox, D.C. (2006). Turbulence Modeling for CFD, 3rd ed. DCW Industries. — Standard k-ω model with cross-diffusion and stress limiter updates.
- Menter, F.R. (1994). "Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications." AIAA Journal 32(8):1598-1605. — SST model derivation and constants.
- Durbin, P.A. (1996). "On the k-ε stagnation point anomaly." Int. J. Heat and Fluid Flow 17:89-90. — Realizability limiter for stagnation regions.
- Wilcox, D.C. (1998). Turbulence Modeling for CFD, 2nd ed. DCW Industries. — Earlier k-ω constants (superseded by 2006 edition for this solver).
- Carrier, A., Errera, M.-P., Gilles, F., and Mouriaux, S. (2021). "Finite element implementation of k−ω SST with automatic wall treatment and target-function approach." Int. J. Numer. Methods Fluids 93:3598-3627. — FEM-specific k-ω implementation strategies.
- Codina, R. (1998). "Comparison of some finite element methods for solving the diffusion-convection-reaction equation." Comput. Methods Appl. Mech. Eng. 156:185-210. — Stabilization for convection-dominated transport.
- Donea, J. and Huerta, A. (2003). Finite Elements Methods for Flow Problems. Wiley. — IPCS and fractional-step methods.
- Logg, A., Mardal, K.-A., and Wells, G.N. (2012). Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book. Springer. — FEniCS framework and Navier-Stokes implementation.
- Nek5000 RANS Tutorial (periodic channel) — High-Re RANS setup and cross-code reference.
- Moser, R.D., Kim, J., and Mansour, N.N. (1999). "Direct numerical simulation of turbulent channel flow up to Re_τ = 590." Phys. Fluids 11(4):943-945. — Classic DNS reference (for model-vs-physics comparison, not used for solver verification).
- DOLFINx Documentation
- dolfinx_mpc Documentation
- CFD-Wiki: SST k-omega model
- CFD-Wiki: Wilcox k-omega model
This solver did not emerge in a vacuum. The following projects, tutorials, and discussions shaped its design and helped debug implementation choices. Listed roughly in order of influence.
FEniCS ecosystem — Navier-Stokes foundations:
- Dokken, J.S. — The FEniCSx Tutorial: Navier-Stokes — The IPCS implementation in this solver follows Dokken's DOLFINx tutorial closely. The Adams-Bashforth convection extrapolation and
nabla_gradconvention come directly from here. - Dokken, J.S. — dolfinx_mpc (docs) — Periodic boundary conditions via multi-point constraints. The topological vs geometrical constraint paths and the
backsubstitutionpattern in this solver mirror the dolfinx_mpc demos. - Mortensen, M. and Valen-Sendstad, K. — Oasis (paper) — High-performance FEniCS Navier-Stokes solver. Our IPCS splitting, P2-P1 function spaces, and solver architecture are influenced by Oasis's IPCS_ABCN scheme. Oasis demonstrated that FEniCS can match dedicated CFD codes in performance.
- Langtangen, H.P. and Logg, A. — FEniCS Tutorial (channels ft07, ft08) — The classic IPCS channel flow and cylinder examples that introduced many users to FEniCS CFD. Our initial velocity and pressure stepping follow this lineage.
- Kuchta, M. — ns-fenics — Collection of FEniCS Navier-Stokes solvers comparing IPCS, mixed, and Yanenko schemes. Useful for understanding trade-offs between splitting approaches.
RANS turbulence in FEniCS — prior art:
- Mortensen, M. and Langtangen, H.P. — CBC.RANS (paper) — The original FEniCS RANS framework (2011). Demonstrated that k-ε and other RANS models could be implemented compactly in FEniCS's variational language. Our approach to decoupling turbulence equations from momentum (Picard iteration with under-relaxation) parallels CBC.RANS, though we target k-ω on modern DOLFINx rather than k-ε on legacy DOLFIN.
- Langtangen, H.P. — cbc.pdesys — Langtangen's PDE system framework for FEniCS, which supported full RANS turbulence hierarchies: from Spalart-Allmaras through k-ε, k-ω SST, v2f, up to Reynolds stress models. Demonstrated how to decompose large coupled PDE systems into segregated subsystems with flexible linearization (Picard/Newton). The architectural pattern of "subsystem per equation + outer coupling loop" influenced our solver design.
- Valen-Sendstad, K. and Mortensen, M. — Implementing a k-ε Turbulence Model in FEniCS — Early paper on FEM RANS in FEniCS. Documents the same positivity and stability challenges we address here (clipping, under-relaxation, convection dominance).
- joove123 — k-epsilon — Lam-Bremhorst k-ε implementation in FEniCS (master's thesis, Univ. Southern Denmark). Validates on channel flow and backward-facing step — the same canonical geometries. Shows how wall-damping functions can be handled in FEniCS's variational framework.
- FEniCS 2024 Conference — Implementation of k-ε in FEniCS (Marcibál, Univ. Southern Denmark) — Recent work confirming continued community interest in FEniCS RANS. Uses DOLFINx with dolfinx_mpc for periodic hills — similar stack to this solver.
- FEniCS Discourse — FEniCSx/DOLFINx implementation for turbulent — Thread where users ask whether FEniCSx is suitable for RANS. Dokken points to the k-ε repo above as a starting point and notes it "should not be hard to convert to DOLFINx." This thread captures the community gap this solver aims to fill.
-
FEniCS Discourse — Implementing k-ε: variational terms depending on a solution — Community discussion on the core challenge of RANS in FEM: how to handle terms like
$P_k = \nu_t |S|^2$ where both$\nu_t$ and$|S|$ depend on the solution. - FEniCS Discourse — BC application with dolfinx_mpc nonlinear solver — RANS on a periodic hill with dolfinx_mpc; boundary condition handling challenges that parallel our periodic channel setup.
Turbulence model specifications and verification:
- NASA Turbulence Modeling Resource — Wilcox k-ω 2006 — Definitive reference for model constants, equations, and verification test cases. We cross-checked our constants against this page.
-
OpenFOAM — kOmega2006 source (header) — OpenFOAM's C++ implementation of the same Wilcox 2006 model. Useful for verifying constant values (
$\beta_0 = 0.0708$ ,$\sigma_k = 0.6$ ,$\sigma_d = 1/8$ ,$C_\text{lim} = 7/8$ ) and cross-diffusion term structure. - CFD-Wiki — Wilcox's modified k-ω model — Documents the 2006 updates (cross-diffusion, stress limiter) vs the 1998 baseline.
- CFD Notes — Enhancements to the k-ω model — Clear explanation of cross-diffusion and stress-limiter terms with derivations.
Cross-code validation target:
- Nek5000 RANS Channel Tutorial — Our primary validation reference. Nek5000 uses a regularized k-ω in spectral elements; comparing our FEM k-ω output against Nek's RANS output confirms implementation correctness independent of spatial discretization.
- Nek5000/NekExamples — Source for the
poiseuille_RANScase files.
MIT License. See LICENSE.
{ "_meta": { "purpose": "Canonical channel benchmark", // single production case "reference": "Nek poiseuille_RANS benchmark parity for Re_τ=1115.8187", // reference target "notes": "Keep only this file in channel/ and clear channel/results to rerun." // workflow note }, "geom": { "Lx": 1.0, // domain length "Ly": 2.0, // full-channel height "Nx": 192, // streamwise cells "Ny": 166, // wall-normal cells "mesh_type": "quad", // element family "y_first": 0.001604628, // first wall-normal spacing "growth_rate": 1.0, // wall-growth ratio "stretching": "tanh", // spacing profile "y_first_tol_rel": 0.2, // y_first consistency tolerance "use_symmetry": false // full-height channel }, "nondim": { "Re_tau": 1115.818661288065, // target friction Reynolds "use_body_force": true // body-force-driven channel }, "turb": { "model": "wilcox2006", // two-equation model "beta_star": 0.09, // Wilcox model constant "nu_t_max_factor": 2000.0, // ν_t limiter "omega_min": 1.0, // ω floor "k_min": 1e-10, // k floor "k_max": 20.0, // k cap "C_lim": 0.0 // Durbin limiter disabled }, "solve": { "dt": 0.005, // initial pseudo-time step "dt_max": 0.01, // pseudo-time step cap "dt_growth": 1.05, // dt growth factor on good convergence "dt_growth_threshold": 0.8, // residual threshold to grow dt "cfl_target": 0.25, // mapped from legacy stability target "t_final": 10000.0, // outer-loop safety time cap "max_iter": 3000, // outer-iteration cap "steady_tol": 1e-6, // residual tolerance "enable_physical_convergence": false, // keep residual-based convergence "physical_u_bulk_rel_tol": 1e-4, // physical gate (unused unless enabled) "physical_tau_wall_rel_tol": 2.5e-3, // physical gate (unused unless enabled) "physical_convergence_start_iter": 10, // gate start iteration "picard_max": 6, // inner nonlinear sweep count "picard_tol": 1e-4, // inner nonlinear stopping "under_relax_k_omega": 0.6, // k/ω under-relaxation "under_relax_nu_t": 0.4, // ν_t under-relaxation "log_interval": 10, // print + log every N iters "snapshot_interval": 50, // numbered PNG + VTX intervals "out_dir": "results", // writes into case/results "min_iter": 50 // block early convergence }, "benchmark": { "reference_profile_csv": "" // optional benchmark overlay CSV } }