Turbulence closure (RANS, Spalart–Allmaras)

One extra transported quantity, $\tilde\nu$ , produces an eddy viscosity $\mu_t$ that is added to the molecular viscosity seen by the mean flow; everything else in the solver is unchanged.

yaml

flow:
  mach: 0.2
  reynolds: 5.0e6
  turbulence: sa          # 'none' (default) or 'sa'; needs reynolds
  extra: {chi_inf: 3.0, prt: 0.9}   # optional: free-stream nu_tilde/nu, turbulent Prandtl number

With turbulence: sa the state has one more variable per node, $\rho\tilde\nu$ , stored after the energy. No-slip walls impose $\tilde\nu = 0$ and define the wall distance; the far field and inlets supply $\tilde\nu = \chi_\infty \nu$ . Outputs, gradients, error estimates and the discipline interface work as for laminar flow; chi_inf and prt are ordinary differentiable flow parameters.

The model is integrated to the wall (no wall functions), so the first element must be thin enough for the polynomial inside it to resolve the viscous sublayer: see the flat-plate study for the resolution used there.

Included: SA-neg with the 2012 modifications of Allmaras, Johnson and Spalart (the $\tilde S$ limiter that keeps the modified vorticity positive, the negative branch, the $f_{t2}$ laminar-suppression term), a constant turbulent Prandtl number, an exact curved wall distance that takes part in shape derivatives.
Not included: trip terms and transition modelling (the flow is turbulent wherever the model produces eddy viscosity), rotation/curvature and compressibility corrections, the quadratic constitutive relation, the $\tfrac23 \rho k$ term of the Boussinesq stress (neglected as is usual for this model), wall functions, other models.
Verified on a manufactured solution, a linearisation study, the zero-pressure-gradient flat plate of the NASA Turbulence Modeling Resource and adjoint-versus-finite-difference derivatives (verification); in three dimensions a spanwise extrusion of that plate reproduces its section. Separated flows, aerofoils at flight Reynolds number and turbulent flows that are three-dimensional in themselves have not been run.

The working variable is $\tilde\nu$ , with $\chi = \tilde\nu/\nu$ . Gyre solves the conservative form for $\rho\tilde\nu$ ,

\pd{(\rho\tilde\nu)}{t} + \nabla\cdot(\rho\mathbf u\,\tilde\nu) = \rho\,(P - D) + \frac{1}{\sigma}\nabla\cdot\big(\rho(\nu + \tilde\nu f_n)\nabla\tilde\nu\big) - \frac{1}{\sigma}(\nu + \tilde\nu f_n)\,\nabla\rho\cdot\nabla\tilde\nu + \frac{c_{b2}}{\sigma}\,\rho\,\abs{\nabla\tilde\nu}^2 , \tag{1}

with production $P$ and destruction $D$ that depend on the sign of $\tilde\nu$ :

	$\tilde\nu \ge 0$	$\tilde\nu < 0$
$P$	$c_{b1}(1 - f_{t2})\,\tilde S\,\tilde\nu$	$c_{b1}(1 - c_{t3})\,S\,\tilde\nu$
$D$	$\big(c_{w1} f_w - \dfrac{c_{b1}}{\kappa^2} f_{t2}\big)\dfrac{\tilde\nu^2}{d^2}$	$-c_{w1}\dfrac{\tilde\nu^2}{d^2}$
$f_n$	$1$	$\dfrac{c_{n1} + \chi^3}{c_{n1} - \chi^3}$
$\mu_t$	$\rho\tilde\nu f_{v1}$ , $\; f_{v1} = \dfrac{\chi^3}{\chi^3 + c_{v1}^3}$	$0$

$S$ is the magnitude of the vorticity and $d$ the distance to the nearest no-slip wall. The mean flow sees the viscosity $\mu + \mu_t$ and the conductivity $c_p(\mu/\mathit{Pr} + \mu_t/\mathit{Pr}_t)$ .

For $\tilde\nu \ge 0$ the modified vorticity is

\tilde S = \begin{cases} S + \bar S & \bar S \ge -c_{v2} S \\[4pt] S + \dfrac{S\,(c_{v2}^2 S + c_{v3}\bar S)}{(c_{v3} - 2c_{v2})S - \bar S} & \text{otherwise} \end{cases} \qquad \bar S = \frac{\tilde\nu f_{v2}}{\kappa^2 d^2},\quad f_{v2} = 1 - \frac{\chi}{1 + \chi f_{v1}} ,

which is positive whenever $S$ is and has a continuous first derivative at the switch. The wall function of the destruction term is

f_w = g\left(\frac{1 + c_{w3}^6}{g^6 + c_{w3}^6}\right)^{1/6},\qquad g = r + c_{w2}(r^6 - r),\qquad r = \min\!\Big(\frac{\tilde\nu}{\tilde S\kappa^2 d^2},\, r_\text{lim}\Big),

and $f_{t2} = c_{t3}\exp(-c_{t4}\chi^2)$ . Constants: $c_{b1} = 0.1355$ , $c_{b2} = 0.622$ , $\sigma = 2/3$ , $\kappa = 0.41$ , $c_{w1} = c_{b1}/\kappa^2 + (1 + c_{b2})/\sigma$ , $c_{w2} = 0.3$ , $c_{w3} = 2$ , $c_{v1} = 7.1$ , $c_{v2} = 0.7$ , $c_{v3} = 0.9$ , $c_{t3} = 1.2$ , $c_{t4} = 0.5$ , $c_{n1} = 16$ , $r_\text{lim} = 10$ ; $\mathit{Pr}_t = 0.9$ and $\chi_\infty = 3$ by default (prt, chi_inf).

A polynomial that represents the sharp edge of a boundary layer undershoots, so $\tilde\nu$ becomes slightly negative there. The original model then loses its footing: the eddy viscosity turns negative and the diffusion coefficient can change sign. SA-neg continues the equation to $\tilde\nu < 0$ so that it stays well posed: no eddy viscosity, a diffusion coefficient that remains positive through $f_n$ , and source terms that drive $\tilde\nu$ back towards zero. The positive and the negative branch join with continuous first derivatives, which is what Newton’s method and the adjoint need.

Equation (1) goes through the machinery of the mean-flow equations (DGSEM, BR2):

Convection. In the volume, the two-point mass flux of the mean flow times the arithmetic mean of $\tilde\nu$ ; on faces, upwinding of $\tilde\nu$ with the interface mass flux $f_\rho^*$ of the Riemann solver, $f^*_{\rho\tilde\nu} = \tfrac12 f_\rho^*(\tilde\nu_L + \tilde\nu_R) - \tfrac12\sqrt{(f_\rho^*)^2 + \delta^2}\,(\tilde\nu_R - \tilde\nu_L)$ .
Diffusion. The conservative diffusion term is one more row of the viscous flux and takes the BR2 path unchanged. $\tilde\nu$ and $\rho$ join the velocity and the temperature as gradient variables.
Sources. Production, destruction and the two non-conservative gradient terms are evaluated with the BR2 gradient $q = G + \sum_f r_f$ , the same lifted gradient the viscous volume flux uses, so that they see the jumps between elements.

The Gauss–Lobatto solution nodes include points on the wall, where $d = 0$ and $\tilde\nu/d$ is undefined for a discrete $\tilde\nu$ that satisfies its wall condition only weakly. The source is therefore integrated with Gauss quadrature: state and BR2 gradient are interpolated to the $(p+1)^d$ Gauss points of each element, where $d > 0$ , the source is evaluated there, and the result is tested against the nodal basis,

(J R)_{i} \mathrel{-}= \frac{1}{w_i}\sum_g \ell_i(\xi_g)\, w_g\, J_g\, s_t(U_g, q_g, d_g).

A source that depends on the gradient of the solution is adjoint consistent only if it is evaluated with a gradient that includes the interface contributions (Oliver and Darmofal). Using the element-local gradient $G$ alone would leave the primal scheme consistent but make the discrete adjoint inconsistent with the continuous one, and functionals would lose their superconvergence. The BR2 gradient $q$ already contains the face liftings, so the same array serves the viscous flux and the source.

Newton’s method and the adjoint need a residual with continuous first derivatives. The model’s own switches (the $\tilde S$ limiter and the sign branch) have them by construction. The remaining kinks are smoothed:

Quantity	Model	Gyre	Size of the change
vorticity magnitude	$\abs{\boldsymbol\omega}$	$\sqrt{\abs{\boldsymbol\omega}^2 + \epsilon^2} - \epsilon$	at most $\epsilon = 10^{-6}$ (units $u_\infty/L_\text{ref}$ ); exactly zero in irrotational flow
limiter of $r$	$\min(r, r_\text{lim})$	$C^1$ blend on $r_\text{lim} \pm 1$	none for $r \le 9$
$f_{v1}$ and $\mu_t$ at $\tilde\nu = 0$	$\chi^3$	$\max(\chi, 0)^3$	none ( $C^2$ already)
interface upwinding	$\abs{f_\rho^*}$	$\sqrt{(f_\rho^*)^2 + \delta^2}$	a dissipation floor $\delta = 10^{-3}\bar\rho$

Condition	$\tilde\nu$
no-slip walls	$\tilde\nu = 0$ , in the convective flux and in the BR2 lifting (Dirichlet)
`farfield`	exterior state $\chi_\infty\mu_\infty$ ; the mass-flux upwinding takes it on inflow and the interior value on outflow
`total_inlet`	$\chi_\infty\,\mu(T_b)$ at the boundary state
slip walls, symmetry, `pressure_outlet`	the interior value; no diffusive flux

The wall distance is computed once per discretization for every Gauss point: a k-d tree over samples of the wall faces proposes candidate faces, and a projected Gauss–Newton iteration on each candidate’s polynomial (of the geometry degree) finds the closest point. The result is stored as interpolation weights of the foot point on the geometry nodes of its wall face, not as a number.

Because the foot point is stored as weights on mesh nodes, the distance is re-evaluated inside the differentiable geometry,

d(X) = \abs{x_g(X) - x_\text{foot}(X)} ,

whenever the node coordinates $X$ change. By the envelope theorem the derivative of a minimum with respect to a parameter is the partial derivative at the minimiser, so holding the foot-point parameters fixed gives the exact first derivative of the distance with respect to $X$ , as long as the closest face does not change. Shape derivatives of turbulent cases therefore include the wall-distance sensitivity; the adjoint study checks them against finite differences that recompute the distance.

In the units of the solver the mean-flow equations are of order one and the turbulence equation of order $\tilde\nu \sim \chi/\mathit{Re}$ , six or seven orders smaller at flight Reynolds number. Unweighted, an inexact Krylov solve and the convergence monitor would not see it. The rows of the turbulence equation are therefore weighted by $\mathit{Re}/\chi_\text{ref}$ with $\chi_\text{ref} = 100$ : the Newton residual norm is $\norm{s\,R_w}$ , the forward Krylov solver works on $(SA)x = Sb$ and the adjoint one on $(S^{-1}A\T)x = S^{-1}b$ . The solutions are those of the unscaled systems; only the norm the tolerances refer to changes. The implicit solver, the block-Jacobi preconditioner and p-multigrid are otherwise the same as for laminar flow, with element blocks one variable larger.

The eddy viscosity enters the stable time step and the block preconditioner, but there is no positivity enforcement for $\rho$ , $p$ or $\mu + \mu_t$ beyond what the model provides.
Start-up from a uniform free stream relies on pseudo-transient continuation; there is no dedicated turbulence start-up strategy and no measured robustness on separated flows.
Wall-resolved meshes have element aspect ratios of the order of a thousand at the wall. Element-block Jacobi does not see that anisotropy: used alone it leaves linear solves on such meshes slow or short of their tolerance. Two remedies exist. The piecewise-constant coarse level of the p-multigrid hierarchy (solver.levels at its default, which includes $p = 0$ ) makes them converge with block Jacobi (solver statistics of the adjoint study). The line smoother (solver.smoother: line, how it works) inverts the coupling along the wall-normal lines of elements exactly and needs a fraction of the Krylov iterations (measured). There is no incomplete-factorisation preconditioner, and nothing has been measured on a turbulent mesh larger than the flat plate of the preset.
The entropy analysis of the two-point fluxes covers the mean-flow convective terms only; nothing is proven for the model equation.

References: Allmaras, Johnson and Spalart (2012); Oliver and Darmofal (2009); see the bibliography.

Turbulence closure (RANS, Spalart–Allmaras)

Model equations#

Discretization#

Wall distance#

Solver#