Turbulence: Spalart–Allmaras RANS

println("Julia ", VERSION)

Julia 1.12.6

Overview

Directly resolving every turbulent eddy is far too expensive for engineering flow, so industry uses Reynolds-averaged Navier–Stokes (RANS): the velocity is split into a mean and a fluctuating part, and the effect of the fluctuations is folded into a turbulent (eddy) viscosity μ_t. Rasmah implements the Spalart–Allmaras one-equation model (Spalart & Allmaras 1992): a transport equation for a modified eddy viscosity ν̃,

\[(u\cdot\nabla)\tilde\nu = c_{b1}\tilde S\,\tilde\nu + \frac{1}{\sigma}\left[\nabla\cdot((\nu+\tilde\nu)\nabla\tilde\nu) + c_{b2}(\nabla\tilde\nu)^2\right] - c_{w1}f_w\left(\frac{\tilde\nu}{d}\right)^2,\]

coupled to the momentum equation through μ_eff = μ + ρ ν̃ f_v1. Here d is the distance to the nearest wall, which drives the destruction term.

The flat-plate test

The canonical validation is a turbulent flat-plate boundary layer: uniform freestream flow over a plate, where the skin-friction coefficient follows the classical 1/5-power-law correlation

\[C_f \approx 0.0592\,Re_x^{-1/5} \qquad \text{(White 2006)}.\]

We use a mesh graded toward the plate (y=0) so the thin boundary layer is resolved:

using Rasmah
using LinearAlgebra

function flat_plate_mesh(nx, ny, L, H; grad = 1.15)
    nv = (nx + 1) * (ny + 1)
    verts = Matrix{Float64}(undef, 3, nv)
    ys = zeros(ny + 1)
    h = H * (1 - 1 / grad) / (1 - 1 / grad^ny)
    for j in 2:(ny + 1)
        ys[j] = ys[j - 1] + h * grad^(j - 2)
    end
    ys .= ys ./ ys[end] .* H
    for j in 0:ny, i in 0:nx
        idx = i + j * (nx + 1) + 1
        verts[1, idx] = (i / nx) * L; verts[2, idx] = ys[j + 1]
    end
    faces = Matrix{Int}(undef, 3, 2nx * ny)
    k = 0
    for j in 0:(ny - 1), i in 0:(nx - 1)
        a = i + j * (nx + 1) + 1; b = a + 1; c = a + nx + 1; d = c + 1
        k += 1; faces[:, k] = [a, b, c]
        k += 1; faces[:, k] = [b, d, c]
    end
    TriangleMesh(verts, faces)
end

flat_plate_mesh (generic function with 1 method)

The operator is three-field [U, P, N] (velocity, pressure, eddy viscosity). wall_distance builds the wall-distance field from the plate nodes. Boundary conditions: freestream on the inlet/top, no-slip on the plate, ν̃ = 0 on the plate and ν̃ = 3ν (freestream) on the inlet/top. We solve with Newton from a freestream initial guess:

function flat_plate_cf(nx, ny, Re_L)
    L, H = 1.0, 0.1
    m = flat_plate_mesh(nx, ny, L, H)
    U = FESpace(m, ReferenceFE(2, 2); valuetype = VectorValue{2})
    P = FESpace(m, ReferenceFE(2, 1))
    N = FESpace(m, ReferenceFE(2, 1))
    μ = 1.0 / Re_L; ρ = 1.0; Uinf = 1.0; ν = μ / ρ
    wall_nodes = [i for i in 1:size(m.vertices, 2) if m.vertices[2, i] < 1e-12]
    op = spalart_allmaras_operator(U, P, N, SpalartAllmaras(), μ, ρ, wall_distance(m, wall_nodes))

    nU = ndofs(U); nP = ndofs(P); nN = ndofs(N); n = nU + nP + nN
    ufixed = Int[]; ug = Float64[]; nufixed = Int[]; nug = Float64[]
    for i in 1:U.nnodes
        x = U.nodes[1, i]; y = U.nodes[2, i]
        onwall = y < 1e-12; onin = x < 1e-12; ontop = y > H - 1e-12
        if onwall || onin || ontop
            append!(ufixed, 2i - 1, 2i)
            onin || ontop ? (append!(ug, Uinf, 0.0)) : (append!(ug, 0.0, 0.0))
        end
    end
    for i in 1:N.nnodes
        y = N.nodes[2, i]; x = N.nodes[1, i]
        if y < 1e-12
            push!(nufixed, i); push!(nug, 0.0)
        elseif x < 1e-12 || y > H - 1e-12
            push!(nufixed, i); push!(nug, 3 * ν)
        end
    end
    fixed = vcat(ufixed, nU .+ [1], nU + nP .+ nufixed)
    g = vcat(ug, [0.0], nug)
    x0 = zeros(n)
    for i in 1:U.nnodes; x0[2i - 1] = Uinf; end
    x0[(nU + nP + 1):end] .= 3 * ν
    free = setdiff(1:n, fixed)
    f = wf -> begin
        w = zeros(eltype(wf), n); w[fixed] .= g; w[free] .= wf
        assemble_residual(op, w)[free]
    end
    jac = wf -> begin
        w = zeros(Float64, n); w[fixed] .= g; w[free] .= wf
        sparse_jacobian(op, w)[free, free]
    end
    res = newton_solve(f, x0[free]; jac = jac, maxiter = 40)
    w = zeros(n); w[fixed] .= g; w[free] .= res.solution

    # skin friction τ_w = μ du/dy|_wall at x ≈ 0.8L
    xstation = 0.8 * L
    wn = [i for i in 1:U.nnodes if abs(U.nodes[1, i] - xstation) < 0.02 && U.nodes[2, i] < 1e-9]
    τw = 0.0; cnt = 0
    for i in wn
        best = nothing; besty = Inf
        for j in 1:U.nnodes
            if abs(U.nodes[1, j] - U.nodes[1, i]) < 1e-9 && 1e-9 < U.nodes[2, j] < besty
                besty = U.nodes[2, j]; best = j
            end
        end
        best === nothing && continue
        τw += μ * (w[2best - 1] / besty); cnt += 1
    end
    cf = (τw / cnt) / (0.5 * ρ * Uinf^2)
    return cf, maximum(w[(nU + nP + 1):end]), res.converged
end

flat_plate_cf (generic function with 1 method)

Verifying against the 1/5-power law

Run the plate at Re_L = 10^4 and compare the skin friction at x = 0.8L against the correlation (a coarse boundary-layer mesh lands within ~15%):

cf, numax, conv = flat_plate_cf(30, 15, 1e4)
Re_x = 1e4 * 0.8
cf_corr = 0.0592 * Re_x^(-0.2)

println("converged      = ", conv)
println("turbulence on  = ", numax > 0.0)
println("C_f computed   = ", round(cf, digits = 5))
println("C_f ≈ 0.0592 Reₓ^(-1/5) = ", round(cf_corr, digits = 5))

converged      = true
turbulence on  = true
C_f computed   = 0.01099
C_f ≈ 0.0592 Reₓ^(-1/5) = 0.00981

Next steps

Navier–Stokes — the laminar baseline the turbulence model augments.
Two-phase flow — another closure (surface tension) on the same incompressible momentum equation.