RasmahRasmah

Turbulence Closures and Hydrodynamic Stability

using Rasmah
using LinearAlgebra

Overview

Turbulence is the chaotic, fluctuating motion of fluids at high Reynolds number. Engineering models avoid resolving every eddy and instead transport statistical quantities: the Reynolds-stress model (RSM) transports the six stress components directly, and the γ-Reθ transition model predicts where a boundary layer turns from laminar to turbulent. Complementary linear stability theory predicts when a smooth flow first becomes unstable — the onset of Rayleigh–Bénard convection, the Rayleigh–Taylor and Kelvin–Helmholtz instabilities, and the Taylor–Green decaying vortex. Each is verified below.

The Reynolds-stress model

The RSM transports R_ij = ⟨u_i′ u_j′⟩. Its production, the LRR pressure–strain redistribution (trace-free), and the isotropic dissipation are:

R = [2.0 0.3 0.0; 0.3 1.0 0.0; 0.0 0.0 0.8]
G = [0.0 5.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
P = rsm_production(R, G)
Π = lrr_pressure_strain(R, tr(R)/2, 1.0, P)
println("P_11 = ", P[1, 1], "  P_12 = ", P[1, 2], "  (homogeneous shear)")
println("tr(Π) = ", tr(Π), "  (trace-free redistribution)")
P_11 = -3.0  P_12 = -5.0  (homogeneous shear)
tr(Π) = -2.220446049250313e-16  (trace-free redistribution)

In homogeneous shear the RSM reaches a self-similar equilibrium with the classic anisotropy b_ij (the LRR-IP benchmark):

res = rsm_homogeneous_shear(5.0)
println("b_11 = ", round(res.b[1, 1], digits=3), "  b_12 = ", round(res.b[1, 2], digits=3), "  Sk/ε = ", round(res.Skε, digits=2))
b_11 = 0.193  b_12 = -0.185  Sk/ε = 5.68

Laminar–turbulent transition

Transition onset is predicted by the Abu-Ghannam–Shaw correlation for the momentum-thickness Reynolds number, combined with the Blasius growth:

for Tu in (0.5, 1.0, 3.0, 5.0)
    println("Tu = ", Tu, "%   Re_x,t = ", round(flat_plate_transition_rex(Tu), digits=0))
end
# higher free-stream turbulence triggers earlier transition
Tu = 0.5%   Re_x,t = 1.347479e6
Tu = 1.0%   Re_x,t = 641027.0
Tu = 3.0%   Re_x,t = 102773.0
Tu = 5.0%   Re_x,t = 65338.0

Hydrodynamic stability

The Rayleigh–Bénard critical Rayleigh number (convection onset in a heated layer) and the Rayleigh–Taylor / Kelvin–Helmholtz growth rates:

c = rayleigh_benard_critical()
println("Ra_c = ", round(c.Ra_c, digits=2), "   α_c = ", round(c.alpha_c, digits=3))

A = atwood_number(2.0, 1.0)
println("RT growth σ = ", round(rayleigh_taylor_growth_rate(2.0, 1.0, 9.81, 10.0), digits=3), "  (√(Agk))")
println("KH growth σ = ", round(kelvin_helmholtz_growth_rate(2.0, 1.0, 5.0, 3.0), digits=3))
Ra_c = 657.51   α_c = 2.221
RT growth σ = 5.718  (√(Agk))
KH growth σ = 7.071

The Taylor–Green decaying vortex

The 2D Taylor–Green vortex is an exact solution of the Navier–Stokes equations, with kinetic energy E(t) = ¼ e^{−4νt}:

ν = 0.01
println("E(0) = ", taylor_green_energy(0.0, ν), "   ε(0) = ", taylor_green_dissipation(0.0, ν), "  (= ν)")
println("E(1) = ", round(taylor_green_energy(1.0, ν), digits=4), "  = ¼ e^{-4ν}")
E(0) = 0.25   ε(0) = 0.01  (= ν)
E(1) = 0.2402  = ¼ e^{-4ν}

Next steps