One-equation RANS closure
ERF offers a Reynolds-averaged (RANS) turbulence closure alongside its LES
and PBL schemes: the one-equation \(k\) model of Axell and Liungman
(2001, Environ. Fluid Mech. 1, 71-106; AL01 below), selected with
erf.rans_type = kEqn. The turbulent kinetic energy \(k\) is
prognostic and the turbulent length scale \(l\) is algebraic, built from
the distance to the ground and the local stratification. The model is
meant for coarse grids on which the boundary layer is not resolved, i.e.
horizontal spacings of tens of metres and up, where the closure carries
the whole turbulent transport. The regression cases in
Exec/CanonicalTests/Canonical_RANS exercise it on flat ground (neutral,
stable and convective) and over 2D and 3D hills.
Eddy viscosity and diffusivity
The Kolmogorov-Prandtl relations (AL01 Eqs. 11-12) give the eddy viscosity and the eddy diffusivity of heat
with the stability functions \(c_\mu\) and \(c_\mu'\) of the next
section. ERF stores \(\rho \nu_t\) in the vertical and horizontal
momentum diffusivities (plot variables Kmv, Kmh) and
\(\rho \nu_t'\) in the vertical heat diffusivity (Khv). By
default the horizontal heat diffusivity and the scalar and moisture
diffusivities are the eddy viscosity divided by \(Pr_t\) and
\(Sc_t\), as for the LES closures; with
erf.rans_consistent_diffusivities = true they all follow
\(\rho \nu_t'\), so that heat, scalars and moisture share the AL01
stability function in every direction. The diffusivity of \(k\)
itself is \(\nu_t / \sigma_k\) with \(\sigma_k = 1\) (AL01
Table I), which the code sets unless erf.sigma_k is given.
Realizability cap
With erf.rans_realizability = C_r (off by default) the eddy
viscosity is capped by the shear,
the form of the limiter in the SST model. \(\nu_t'\) keeps the ratio \(c_\mu' / c_\mu\) of the uncapped closure, and every other eddy diffusivity (heat, scalars, moisture, \(k\)) is built from the capped \(\nu_t\), so all of them are capped by the same factor; so is the shear production of \(k\), which uses the capped viscosity. \(S\) is the cell-centred average of the face and edge strains of the velocity at the start of the step, the same average the Smagorinsky closure uses. The production of \(k\) recomputes the strain at every Runge-Kutta stage, so within a step the capped stress \(\nu_t S = C_r k\) holds for the start-of-step strain only and the cap lags a fast transient by one step; applying it at each stage is left for later. In the first cell above a surface layer the wall-face strain enters that average, so \(S\) there is below the log-law value and the cap acts less than the closed form below suggests.
Why \(C_r = (c_\mu^0)^2 = 0.309\). With Eqs. 11 and 19 the closure gives
so the cap with \(C_r = (c_\mu^0)^2\) is active where
\(P / \varepsilon > c_\mu^0 / c_\mu\). In a neutral constant-stress
layer (\(c_\mu = c_\mu^0\), \(P = \varepsilon\)) the ratio is
exactly \((c_\mu^0)^2 = 0.30936\) for the default
\(c_\mu^0 = 0.5562\), so this value leaves the layer unchanged and
acts where the production exceeds the dissipation, as in the shear layer
that leaves a crest. There the uncapped closure returns a viscosity large
enough to keep the flow attached behind steep terrain, and the cap lets the
lee separate. In a stably stratified equilibrium layer \(c_\mu\) falls
faster than \(P / \varepsilon\) rises and the cap stays inactive; in a
weakly unstable one (\(-0.95 < R_t < 0\)) it trims \(\nu_t\) by up
to 0.3 %. The rounded \(C_r = 0.309\) trims a neutral equilibrium cell
by 0.1 %. GUST-CFD uses the same cap and value (hgopalan/gust-abl#51); its
authors report comparisons with the Bolund, Askervein, Big Southern Butte,
Ishihara-Zhou and RUSHIL data sets, and an earlier, unmerged
realizable-\(c_\mu\) study in kynema/kynema-sgf#2036. Values of
\(C_r\) are accepted in \([0.01, 10]\); a smaller value trims the
viscosity in the approach flow too. The key needs erf.rans_type = kEqn
and no PBL scheme (a PBL scheme overwrites the vertical diffusivities the
cap acts on); one value for the whole hierarchy (a single entry, or a
per-level list whose entries are all equal) applies to the k-eqn levels and
is reset to 0 on the others, while a per-level list with different entries
that switches the cap on at a level without the closure stops the run. The
check against a PBL scheme runs after that reset, so an LES level with a
PBL scheme beside a k-eqn level with the cap is accepted. EB terrain
cannot be combined with the k-eqn closure, so the cap never meets the
cut-cell strain. The
plotfile variable cap_ratio holds \(\nu_t / (C_r k / S)\) of the
uncapped closure: above 1 the cap is active in that cell. The
RANS_Neutral_Hill_2D_Realizability test checks the capped Kmv and
Khv of the ridge case against
\(\rho\, c_\mu k^{1/2} l / \max(\text{cap\_ratio}, 1)\) and
\(\rho\, c_\mu' k^{1/2} l / \max(\text{cap\_ratio}, 1)\) cell by cell,
with \(k^{1/2}\) recovered from diss, which the same kernel writes,
and the unit test RANSClosure.StrainMagnitudeForTheCap pins \(S\)
on simple shear and plane strain.
Transport equation for \(k\)
The prognostic equation (AL01 Eq. 13) is the Deardorff form ERF already integrates for LES,
with the shear production \(P_s = 2 \rho \nu_t S_{ij} S_{ij}\), the
buoyancy production \(P_b = (g / \theta_0) \, \overline{w' \theta'}\)
from the modelled heat flux \(-\rho \nu_t' \, \partial \theta / \partial z\)
(AL01 Eq. 15, dry air; \(\theta_0\) is erf.theta_ref or the local
\(\theta\)), and the dissipation (AL01 Eq. 19)
The proportionality constant \((c_\mu^0)^3\) is what makes the
logarithmic law hold in a constant-stress layer (AL01 Sect. 4.1). The
dissipation is explicit in the update by default;
erf.implicit_tke_dissipation = true treats it as
\((c_\mu^0)^3 k_{old}^{1/2} / l\) times the new \(k\), folded into
the time update, which removes the dissipation time-step restriction. In
the plotfile diss holds \(\rho \varepsilon\) from the start of the
last step. A floor erf.tke_floor (default: machine epsilon on
\(\rho k\)) bounds \(k\) from below.
The heat flux in \(P_b\) at the cell centre is the average of the
fluxes of the \(\theta\) diffusion at the two vertical faces of the
cell, in every cell; at the bottom of the domain the face flux is the
surface-layer flux when a surface layer is used. These are the full face
fluxes, including the stretched-mesh spacing, whether the vertical
diffusion of \(\theta\) is explicit or implicit
(erf.vert_implicit), so the source does not depend on that choice.
On a terrain-fitted mesh the stored face flux is the component along the
\(\zeta\) coordinate only: the slope corrections and the Jacobian
are applied later, when the fluxes enter the right-hand side, so over
sloping terrain both \(P_b\) and the hfx3 output omit those
terms. The Deardorff closure shares this source term. The k-eqn closure
cannot be combined with a PBL scheme that also transports the turbulent
kinetic energy (MYJ, MYNN25, MYNNEDMF, the SHOC schemes);
the Deardorff closure cannot be combined with any PBL scheme.
Length scale
Near a wall the length scale is the geometric one (AL01 Eq. 22),
with \(d\) the distance to the ground and \(z_0\) the roughness
length of the surface layer. ERF caps it harmonically at
erf.max_geom_lscale (30 m by default, about \(\kappa \times 0.1 z_i\)),
or, with erf.rans_lscale_from_pblh = true, at
\(\kappa \times 0.1 z_i\) from the surface layer’s boundary-layer height
diagnostic (erf.most.pblh_calc = MYNN25), clamped between
erf.rans_lscale_min and erf.max_geom_lscale. On a refined region
that does not run from the ground to the top of the domain, \(z_i\) is
the height diagnosed on the next coarser level, so the cap does not depend on
the extent of the refined region; each fine column takes the height of the
coarse column that holds it, so the cap is constant over each block of
refinement-ratio fine columns, as it is over each column of the coarse level.
Stratification shortens or lengthens \(l\) through the buoyancy frequency \(N^2 = (g/\theta_0) \, \partial \theta / \partial z\). In stable air (AL01 Eq. 26) the geometric and buoyancy length scales combine harmonically,
so that \(l \to c_b k^{1/2} / N\) far from the wall under strong stratification. In unstable air AL01 rewrite the same relation with the dissipation of Eq. 19 on the right-hand side (their Eq. 28),
with \(R_t < 0\) the turbulent Richardson number. ERF evaluates Eq. 28 once, with \(R_t\) from the geometric length and passed through the smoothing of the next section, which bounds the unstable length at \(l_g \, (1 + (c_\mu^0)^6 c_b^{-2} |R_t^{min}|)^{1/2}\), about \(1.31 \, l_g\) for the default constants. Iterating Eq. 28 instead (dissipation from the new length, length from the new \(R_t\)) is the fixed-point map of Eq. 26 with \(N^2 < 0\), which has no fixed point once \(l_g^2 |N^2|\) exceeds \(c_b^2 k\) (AL01 p. 78), so the length would grow without bound in strongly convective, weakly turbulent air.
Stability functions
With the length in hand the turbulent Richardson number is recomputed as \(R_t = l^2 N^2 / ((c_\mu^0)^6 k)\) and the stability functions of Launder, as recalibrated by AL01 against the Högström data set, are
(AL01 Eqs. 31-32; the turbulent Prandtl number is \(c_\mu / c_\mu'\)).
\(c_\mu\) has its pole at the larger root of its denominator,
\(R_t^{pole} = (-0.308 + \sqrt{0.308^2 - 4 \cdot 0.00837}) / (2 \cdot 0.00837)
\approx -3.59869\), and \(c_\mu'\) at \(-1/0.277 \approx -3.61011\).
Just below \(R_t^{pole}\) the denominator of \(c_\mu\) is negative, and
\(c_\mu\) (with the eddy viscosity) stays negative down to the zero of its
numerator, \(-c_\mu^0/0.108\) (\(-5.15\) for the default
\(c_\mu^0\)). Below a critical value
\(R_t^c\) (erf.Rt_crit, default \(-1\)) the smoothing of
Burchard and Petersen maps \(R_t\) onto \((R_t^{min}, R_t^c)\)
with \(R_t^{min}\) = erf.Rt_min (default \(-3\)):
ERF uses the second, algebraically identical form: the first cancels two
terms of order \(|R_t|\) and returned meaningless values once
\(|R_t|\) exceeded about \(10^{15}\), which happens wherever
\(k\) sits at its floor under an unstable \(N^2\). The inputs are
validated at start-up: \(R_t^{min} < R_t^c \le 0\) and
\(R_t^{min} > -3.5986889\), just above \(R_t^{pole}\), so the
smoothed \(R_t\) never reaches the pole, and
\(R_t^{min} > -c_\mu^0/0.108\), so the numerator of \(c_\mu\) stays
positive (this is the tighter limit when \(c_\mu^0 < 0.3887\)).
Values close to the pole are accepted but are not physically meaningful:
\(c_\mu\) is 1.53 at \(-3\), 7.26 at \(-3.5\), 78 at
\(-3.59\) and about 980 at \(-3.598\), against
\(c_\mu^0 = 0.5562\), and the turbulent Prandtl number
\(c_\mu / c_\mu'\) rises from 0.40 at \(-3.5\) to 5.9 at
\(-3.598\). Keep \(R_t^{min}\) at about \(-3.5\) or above.
The smoothed \(R_t\) and the two functions can be written to the plotfile as Rt, cmu and cmu_prime.
Wall condition on \(k\)
Under a surface layer (zlo.type = surface_layer) and with
erf.dirichlet_k = true, the first cell above the ground takes the
equilibrium value of AL01 Eq. 16,
with \(u_*\) and the buoyancy flux \(B = -(g/\theta_0) u_* \theta_*\)
from the surface-layer model and \(d_1\) the distance of the cell
centre to the ground; only a destabilising flux enters. The value is set
at the start of each step and held through every Runge-Kutta stage,
including the implicit vertical diffusion solve, and the ghost cell
below the wall carries the same value so no diffusive flux of \(k\)
crosses the wall face. Without the flag the wall flux of \(k\) is zero
and the first cell settles at about half the equilibrium value, because
the coarse first cell cannot resolve the shear production there; the
mean wind still follows the logarithmic law since the surface layer
supplies the stress, but everything fed by near-wall \(k\) is off. ERF
warns at start-up in that configuration. With
erf.init_tke_from_ustar = true the initial \(k\) profile tapers
linearly from \(u_*^2\) at the surface.
Wall distance
On a flat or vertically stretched mesh the distance \(d\) is the
height of the cell centre above the ground, \(z - z_{lo}\), where
\(z_{lo}\) is geometry.prob_lo[2]. On a terrain-fitted mesh two
choices are offered through erf.wall_dist_type:
poisson(default): the differential-equation distance of Tucker (2003, J. Comput. Phys. 190, 229-248). A Poisson equation \(\nabla^2 \phi = -1\) is solved in terrain-following coordinates with \(\phi = 0\) on the ground and zero normal gradient elsewhere, and \(d = -|\nabla \phi| + (|\nabla \phi|^2 + 2 \phi)^{1/2}\). The gradient is formed at cell centres from centred differences with the cell-centre metrics. On the canonical hills the mean error against the exact distance is about 1 %, with the largest errors, 5 to 10 %, in the first cells over a convex crest.terrain_height: the height of the cell centre above the local surface, projected on the surface normal, \(d = (z - z_s) / (1 + h_\xi^2 + h_\eta^2)^{1/2}\) with \(h_\xi, h_\eta\) the terrain slopes. It is exact on a plane, within a few percent of the true distance for slopes below about 0.3 (mean errors of 0.02 % and 0.01 % on the canonical ridge and hill), and needs no linear solve. It cannot see side walls or overhangs, which the Poisson distance can.
The distance is available as the plot variable walldist. It is
computed at initialisation, also on a restart, and again on every level
that a regrid adds or remakes.
Limitations
The buoyancy production and \(N^2\) are dry: no virtual or liquid-water potential temperature.
All levels must use the closure; a hybrid RANS-LES set-up is refused.
Embedded and thin-body boundaries are not supported by the Poisson wall distance.
Under the anelastic integrator all vertical diffusion is explicit, so the time step is bounded by \(\Delta z^2 / (2 K)\); with the eddy viscosities a convective boundary layer produces (tens of m2/s) this is the binding constraint, not the closure.
The closure is local: a convective layer keeps a superadiabatic lapse of order \(-F/K_h\) through its depth where a countergradient scheme or LES would mix it out.
Inputs
Parameter |
Definition |
Default |
|---|---|---|
erf.rans_type |
|
|
erf.Cmu0, erf.Cb |
\(c_\mu^0\), \(c_b\) |
0.5562, 0.35 |
erf.Rt_crit, erf.Rt_min |
smoothing bounds on \(R_t\); see the limits above |
-1, -3 |
erf.max_geom_lscale |
cap on \(l_g\) [m] |
30 |
erf.sigma_k |
Schmidt number of \(k\) |
1 (RANS) |
erf.theta_ref |
\(\theta_0\); 0 uses the local \(\theta\) |
0 |
erf.dirichlet_k |
wall value of \(k\) from \(u_*\), AL01 Eq. 16 |
false |
erf.init_tke_from_ustar |
initial \(k\) profile from \(u_*\) |
false |
erf.tke_floor |
floor on \(k\) [m2/s2] |
0 (epsilon) |
erf.implicit_tke_dissipation |
linearised implicit dissipation |
false |
erf.rans_consistent_diffusivities |
heat, scalar and moisture diffusivities from \(c_\mu'\) |
false |
erf.rans_lscale_from_pblh |
cap \(l_g\) at \(\kappa 0.1 z_i\) |
false |
erf.rans_lscale_min |
floor of that cap [m] |
1 |
erf.wall_dist_type |
|
|
See Inputs for the full descriptions.