Turbulence Closures

In turbulence.jl we specify turbulence closures. Currently, pointwise models of the eddy viscosity/eddy diffusivity type are supported for turbulent shear and tracer diffusivity. Methods currently supported are:
ConstantViscosityWithDivergence
SmagorinskyLilly
Vreman
AnisoMinDiss

Note

Usage: This is a quick-ref guide to using turbulence models as a subcomponent of AtmosModel
$\nu$ is the kinematic viscosity, $C_smag$ is the Smagorinsky Model coefficient,

  • turbulence=ConstantViscosityWithDivergence(ν)
  • turbulence=SmagorinskyLilly(C_smag)
  • turbulence=Vreman(C_smag)
  • turbulence=AnisoMinDiss(C_poincare)
using DocStringExtensions
using CLIMAParameters.Atmos.SubgridScale: inv_Pr_turb
export ConstantViscosityWithDivergence, SmagorinskyLilly, Vreman, AnisoMinDiss
export turbulence_tensors

Abstract Type

We define a TurbulenceClosure abstract type and default functions for the generic turbulence closure which will be overloaded with model specific functions. Minimally, overloaded functions for the following stubs must be defined for a turbulence model.

abstract type TurbulenceClosure end


vars_state_gradient((::TurbulenceClosure, FT) = @vars()
vars_state_gradient_flux(::TurbulenceClosure, FT) = @vars()
vars_state_auxiliary(::TurbulenceClosure, FT) = @vars()

function atmos_init_aux!(
    ::TurbulenceClosure,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
) end
function compute_gradient_argument!(
    ::TurbulenceClosure,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function compute_gradient_flux!(
    ::TurbulenceClosure,
    ::Orientation,
    diffusive,
    ∇transform,
    state,
    aux,
    t,
) end

The following may need to be addressed if turbulence models require additional state variables or auxiliary variable updates (e.g. TKE based models)

vars_state_conservative(::TurbulenceClosure, FT) = @vars()
function atmos_nodal_update_auxiliary_state!(
    ::TurbulenceClosure,
    ::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end

Eddy-viscosity Models

The following function provides an example of a stub for an eddy-viscosity model. Currently, scalar and diagonal tensor viscosities and diffusivities are supported.

ClimateMachine.Atmos.turbulence_tensorsFunction
ν, D_t, τ = turbulence_tensors(
                ::TurbulenceClosure,
                orientation::Orientation,
                param_set::AbstractParameterSet,
                state::Vars,
                diffusive::Vars,
                aux::Vars,
                t::Real
            )

Compute the kinematic viscosity (ν), the diffusivity (D_t) and SGS momentum flux tensor (τ) for a given turbulence closure. Each closure overloads this method with the appropriate calculations for the returned quantities.

Arguments

  • ::TurbulenceClosure = Struct identifier for turbulence closure model
  • orientation = AtmosModel.orientation
  • param_set = AtmosModel.param_set
  • state = Array of prognostic (state) variables. See vars_state_conservative in AtmosModel
  • diffusive = Array of diffusive variables
  • aux = Array of auxiliary variables
  • t = time
source

Generic math functions for use within the turbulence closures such as the principal tensor invariants, symmetric tensors and tensor norms have been included.

Pricipal Invariants

\[\textit{I}_{1} = \mathrm{tr(X)} \\ \textit{I}_{2} = (\mathrm{tr(X)}^2 - \mathrm{tr(X)^2}) / 2 \\ \textit{I}_{3} = \mathrm{det(X)} \\\]

Symmetrize

\[\frac{\mathrm{X} + \mathrm{X}^{T}}{2} \\\]

2-Norm

Given a tensor X, return the tensor dot product

\[\sum_{i,j} S_{ij}^2\]

Strain-rate Magnitude

By definition, the strain-rate magnitude, as defined in standard turbulence modelling is computed such that

\[|\mathrm{S}| = \sqrt{2 \sum_{i,j} \mathrm{S}_{ij}^2}\]

where

\[\vec{S}(\vec{u}) = \frac{1}{2} \left(\nabla\vec{u} + \left( \nabla\vec{u} \right)^T \right)\]

\mathrm{S} is the rate-of-strain tensor. (Symmetric component of the velocity gradient). Note that the skew symmetric component (rate-of-rotation) is not currently computed.

"""
    strain_rate_magnitude(S)
Given the rate-of-strain tensor `S`, computes its magnitude.
"""
function strain_rate_magnitude(S::SHermitianCompact{3, FT, 6}) where {FT}
    return sqrt(2 * norm2(S))
end

Constant Viscosity Model

ConstantViscosityWithDivergence requires a user to specify the constant viscosity (kinematic) and appropriately computes the turbulent stress tensor based on this term. Diffusivity can be computed using the turbulent Prandtl number for the appropriate problem regime.

\[\tau = - 2 \nu \mathrm{S}\]

Smagorinsky-Lilly

The Smagorinsky turbulence model, with Lilly's correction to stratified atmospheric flows, is included in ClimateMachine. The input parameter to this model is the Smagorinsky coefficient. For atmospheric flows, the coefficient C_smag typically takes values between 0.15 and 0.23. Flow dependent C_smag are currently not supported (e.g. Germano's extension). The Smagorinsky-Lilly model does not contain explicit filtered terms.

Equations

\[\nu = (C_{s} \mathrm{f}_{b} \Delta)^2 \sqrt{|\mathrm{S}|}\]

with the stratification correction term

\[f_{b} = \begin{cases} 1 & \mathrm{Ri} \leq 0 ,\\ \max(0, 1 - \mathrm{Ri} / \mathrm{Pr}_{t})^{1/4} & \mathrm{Ri} > 0 . \end{cases}\]
\[\mathrm{Ri} = \frac{N^2}{{|S|}^2}\]
\[N = \left( \frac{g}{\theta_v} \frac{\partial \theta_v}{\partial z}\right)^{1/2}\]

Here, $\mathrm{Ri}$ and $\mathrm{Pr}_{t}$ are the Richardson and turbulent Prandtl numbers respectively. $\Delta$ is the mixing length in the relevant coordinate direction. We use the DG metric terms to determine the local effective resolution (see src/Mesh/Geometry.jl), and modify the vertical lengthscale by the stratification correction factor $\mathrm{f}_{b}$ so that $\Delta_{vert} = \Delta z f_b$.

ClimateMachine.Atmos.SmagorinskyLillyType
SmagorinskyLilly <: TurbulenceClosure

Fields

  • C_smag

    Smagorinsky Coefficient [dimensionless]

Smagorinsky Model Reference

article{doi:10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2,
  author = {Smagorinksy, J.},
  title = {General circulation experiments with the primitive equations},
  journal = {Monthly Weather Review},
  volume = {91},
  number = {3},
  pages = {99-164},
  year = {1963},
  doi = {10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2},
  URL = {https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2},
  eprint = {https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2}
  }

Lilly Model Reference

article{doi:10.1111/j.2153-3490.1962.tb00128.x,
  author = {LILLY, D. K.},
  title = {On the numerical simulation of buoyant convection},
  journal = {Tellus},
  volume = {14},
  number = {2},
  pages = {148-172},
  doi = {10.1111/j.2153-3490.1962.tb00128.x},
  url = {https://onlinelibrary.wiley.com/doi/abs/10.1111/j.2153-3490.1962.tb00128.x},
  eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.2153-3490.1962.tb00128.x},
  year = {1962}
  }

Brunt-Väisälä Frequency Reference

Brunt-Vaisala frequency N² defined as in equation (1b) in
  Durran, D.R. and J.B. Klemp, 1982:
  On the Effects of Moisture on the Brunt-Väisälä Frequency.
  J. Atmos. Sci., 39, 2152–2158,
  https://doi.org/10.1175/1520-0469(1982)039<2152:OTEOMO>2.0.CO;2
source

Vreman Model

Vreman's turbulence model for anisotropic flows, which provides a less dissipative solution (specifically in the near-wall and transitional regions) than the Smagorinsky-Lilly method. This model relies of first derivatives of the velocity vector (i.e., the gradient tensor). By design, the Vreman model handles transitional as well as fully turbulent flows adequately. The input parameter to this model is the Smagorinsky coefficient - the coefficient is modified within the model functions to account for differences in model construction.

Equations

\[\nu_{t} = 2.5 C_{s}^2 \sqrt{\frac{B_{\beta}}{u_{i,j}u_{i,j}}},\]

where ($i,j, m = (1,2,3)$)

\[\begin{align} B_{\beta} &= \beta_{11}\beta_{22} + \beta_{11}\beta_{33} + \beta_{22}\beta_{33} - (\beta_{13}^2 + \beta_{12}^2 + \beta_{23}^2) \\ \beta_{ij} &= \Delta_{m}^2 u_{i, m} u_{j, m} \\ u_{i,j} &= \frac{\partial u_{i}}{\partial x_{j}}. \end{align}\]
ClimateMachine.Atmos.VremanType
Vreman{FT} <: TurbulenceClosure

Filter width Δ is the local grid resolution calculated from the mesh metric tensor. A Smagorinsky coefficient is specified and used to compute the equivalent Vreman coefficient.

  1. ν_e = √(Bᵦ/(αᵢⱼαᵢⱼ)) where αᵢⱼ = ∂uⱼ∂uᵢ with uᵢ the resolved scale velocity component.
  2. βij = Δ²αₘᵢαₘⱼ
  3. Bᵦ = β₁₁β₂₂ + β₂₂β₃₃ + β₁₁β₃₃ - β₁₂² - β₁₃² - β₂₃²

βᵢⱼ is symmetric, positive-definite. If Δᵢ = Δ, then β = Δ²αᵀα

Fields

  • C_smag

    Smagorinsky Coefficient [dimensionless]

Reference

@article{Vreman2004,
  title={An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications},
  author={Vreman, AW},
  journal={Physics of fluids},
  volume={16},
  number={10},
  pages={3670--3681},
  year={2004},
  publisher={AIP}
}
source

Anisotropic Minimum Dissipation

This method is based Vreugdenhil and Taylor's minimum-dissipation eddy-viscosity model. The principles of the Rayleigh quotient minimizer are applied to the energy dissipation terms in the conservation equations, resulting in a maximum dissipation bound, and a model for eddy viscosity and eddy diffusivity.

\[\nu_e = (\mathrm{C}\delta)^2 \mathrm{max}\left[0, - \frac{\hat{\partial}_k \hat{u}_{i} \hat{\partial}_k \hat{u}_{j} \mathrm{\hat{S}}_{ij}}{\hat{\partial}_p \hat{u}_{q} \hat{\partial}_p \hat{u}_{q}} \right]\]
ClimateMachine.Atmos.AnisoMinDissType
AnisoMinDiss{FT} <: TurbulenceClosure

Filter width Δ is the local grid resolution calculated from the mesh metric tensor. A Poincare coefficient is specified and used to compute the equivalent AnisoMinDiss coefficient (computed as the solution to the eigenvalue problem for the Laplacian operator).

Fields

  • C_poincare

Reference

@article{
    doi:10.1063/1.5037039,
    author = {Vreugdenhil,Catherine A.  and Taylor,John R. },
    title = {Large-eddy simulations of stratified plane Couette flow using the anisotropic minimum-dissipation model},
    journal = {Physics of Fluids},
    volume = {30},
    number = {8},
    pages = {085104},
    year = {2018},
    doi = {10.1063/1.5037039},
    URL = {https://doi.org/10.1063/1.5037039}
}
source