Diagnosing Delayed Instability in Steady-State RANS Liquid-Metal CFD

 

A particularly frustrating phenomenon in steady-state Reynolds-Averaged Navier-Stokes (RANS) simulations—especially in low-viscosity, high-momentum fluids like liquid metals—is delayed numerical instability. The iterative solver (such as SIMPLE or SIMPLEC) initially exhibits smooth convergence: residuals drop, and turbulent kinetic energy (k) as well as dissipation rate (ε) evolve toward physical-looking values. However, after dozens or hundreds of well-behaved iterations, k begins to grow exponentially in isolated cells, blowing up the turbulent viscosity (μ_t) and eventually crashing the simulation with floating-point errors (NaN).

Here is an in-depth breakdown of the physical and numerical mechanisms causing this delayed feedback loop and how to systematically diagnose it.

---

1. The Primary Culprit: Standard Wall Functions at y⁺ ≈ 1

If you are using a standard high-Reynolds number k-ε model (which relies on logarithmic wall functions) on a mesh refined to y⁺ ≈ 1, you are introducing a mathematical inconsistency into the boundary conditions.

• Log-Law Assumption Breakdown: Standard wall functions assume that the first cell centroid resides in the fully turbulent logarithmic layer (typically 30 < y⁺ < 300). They enforce shear stress and production of k based on the logarithmic velocity profile:

u⁺ = (1 / κ) · ln(E · y⁺)

• Over-Prediction of Turbulence Production: In the viscous sublayer (y⁺ < 5), the true velocity profile is linear (u⁺ = y⁺). Evaluating log-law formulation at y⁺ ≈ 1 causes the wall shear stress and the velocity gradient (∂u/∂y) to be drastically overestimated. This injects unphysical turbulent kinetic energy production P_k into the wall-adjacent cells:

P_k = τ_w · (∂u / ∂y) ∼ μ_t · S²

2. The Positive Feedback Loop & Stagnation Anomaly

Standard k-ε models compute eddy viscosity as:

μ_t = C_μ · ρ · (k² / ε)

Where C_μ = 0.09 is treated as a constant. This creates a non-linear positive feedback mechanism:

1. An unphysical spike in P_k increases k in near-wall or high-shear/stagnation cells.
2. Higher k raises μ_t quadratically relative to ε.
3. Increased μ_t causes excessive momentum diffusion into surrounding cells, altering local velocity gradients and local strain rate S.
4. Because standard k-ε lacks a realizability constraint, P_k = μ_t S² grows without bound, driving k and μ_t even higher in subsequent iterations.

3. Why Is the Instability "Delayed"?

Why does the calculation run smoothly for 100+ iterations before exploding?

• Under-Relaxation Suppression: Low under-relaxation factors (URFs) for momentum and turbulence equations damp the local growth rate per iteration. The numerical error accumulates slowly, effectively hiding the instability behind small iteration updates.
• Convective Transport Lag: The local numerical error generated at the wall or stagnation region needs time to advect into the main core flow. Until the local peak in μ_t alters the global flow field, global residuals may continue to fall. Once the perturbed flow feeds back into the shear layer, the system reaches a tipping point, overcoming numerical damping and causing an exponential explosion.

4. Recommended Diagnostics & Remediation Workflow

To isolate and resolve this issue, apply the following diagnostic steps:

• Monitor Spatial Maxima: Do not rely solely on volume-averaged residuals. Create monitor points for max(k), max(μt / μ), and max(Pk) per iteration. The iteration where max(μt / μ) begins its upward bend marks the exact onset of instability.
• Locate the Source Cell: Save field data right before divergence occurs. Contour plots of k, μ_t, and P_k / ε will reveal the exact boundary or recirculation zone where energy accumulation originates.
• Match Wall Treatment to Mesh:
  • If y⁺ ≈ 1: Switch from standard k-ε to a Low-Reynolds model or k-ω SST, which handles viscous sublayer physics integration directly to the wall.
  • If standard k-ε must be used: Coarsen the near-wall mesh inflation layer so that y⁺ > 30 across all solid walls.
• Enforce Production Limiters / Realizability: Enable turbulence production limiters (e.g., Kato-Launder formulation or Realizable k-ε constraint) to cap P_k in high strain rate / impinging flow regions, preventing runaway C_μ S² growth.

Have you encountered delayed turbulence runaway in low-Prandtl number CFD applications? Share your solver parameters and mesh strategies in the comments below!