Turbulent kinetic energy equation

The TurbulentKineticEnergyEquation module provides diagnostics for the turbulent kinetic energy (TKE) budget. While the KineticEnergyEquation module deals with the total (resolved) kinetic energy, this module focuses on the kinetic energy of velocity perturbations relative to a specified mean flow. TKE is defined as

\[e = \tfrac{1}{2} u_i' u_i'\]

where $u_i' = u_i - U_i$ is the velocity perturbation from the mean $U_i$. Mean velocities default to zero (in which case TKE reduces to KE) and can be set via the U, V, W keyword arguments.

A key term in the TKE budget is the shear production rate, which quantifies the transfer of kinetic energy from the mean flow to the turbulence through Reynolds stresses acting on the mean shear:

\[P = -u_i' u_j' \partial_j U_i\]

The module provides both directional components ($P_x$, $P_y$, $P_z$) and the total shear production. It also provides an isotropic dissipation rate diagnostic that computes $\varepsilon = 2\nu S'_{ij} S'_{ij}$ using the perturbation strain rate tensor.

All diagnostics are computed at (Center, Center, Center).

Example

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(grid; closure=ScalarDiffusivity(ν=1e-4));

julia> tke = TurbulentKineticEnergyEquation.TurbulentKineticEnergy(model)
TurbulentKineticEnergy KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: turbulent_kinetic_energy_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField")
└── computes: turbulent kinetic energy  ½uᵢ′uᵢ′

julia> SP = TurbulentKineticEnergyEquation.ShearProductionRate(model)
TurbulentKineticEnergyShearProductionRate KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: shear_production_rate_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField")
└── computes: total TKE shear production  -uᵢ′uⱼ′ ∂ⱼUᵢ

Turbulent kinetic energy

Oceanostics.TurbulentKineticEnergyEquation.TurbulentKineticEnergyType
TurbulentKineticEnergy(model, u, v, w; U, V, W, location)

Calculate the turbulent kinetic energy of model.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(grid);

julia> TKE = TurbulentKineticEnergyEquation.TurbulentKineticEnergy(model)
TurbulentKineticEnergy KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: turbulent_kinetic_energy_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField")
└── computes: turbulent kinetic energy  ½uᵢ′uᵢ′
source

Isotropic dissipation rate

Oceanostics.TurbulentKineticEnergyEquation.TurbulentKineticEnergyIsotropicDissipationRateFunction
TurbulentKineticEnergyIsotropicDissipationRate(
    u,
    v,
    w,
    args;
    U,
    V,
    W,
    location
)

Calculate the Turbulent Kinetic Energy Isotropic Dissipation Rate, defined as

ε = 2 ν S'ᵢⱼS'ᵢⱼ,

where S'ᵢⱼ is the strain rate tensor, for a fluid with an isotropic turbulence closure (i.e., a turbulence closure where ν (eddy or not) is the same for all directions.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(grid; closure=ScalarDiffusivity(ν=1e-4));

julia> TurbulentKineticEnergyEquation.IsotropicDissipationRate(model)
KineticEnergyIsotropicDissipationRate KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: isotropic_viscous_dissipation_rate_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "NamedTuple")
└── computes: isotropic kinetic energy dissipation rate  ε = 2νSᵢⱼSᵢⱼ
source

Shear production rates

Oceanostics.TurbulentKineticEnergyEquation.TurbulentKineticEnergyXShearProductionRateType
TurbulentKineticEnergyXShearProductionRate(
    u′,
    v′,
    w′,
    U,
    V,
    W;
    grid,
    location
)

Calculate the shear production rate in the model's x direction:

XSHEAR = uᵢ′u′∂x(Uᵢ)

where uᵢ′ is the velocity perturbation in the i direction, u′ is the velocity perturbation in the x direction, Uᵢ is the background velocity in the i direction, and ∂x is the horizontal derivative.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(grid);

julia> XSHEAR = TurbulentKineticEnergyEquation.XShearProductionRate(model)
TurbulentKineticEnergyXShearProductionRate KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: shear_production_rate_x_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField")
└── computes: TKE shear production (x)  -uᵢ′u′ ∂ₓUᵢ
source
Oceanostics.TurbulentKineticEnergyEquation.TurbulentKineticEnergyYShearProductionRateType
TurbulentKineticEnergyYShearProductionRate(
    u′,
    v′,
    w′,
    U,
    V,
    W;
    grid,
    location
)

Calculate the shear production rate in the model's y direction:

YSHEAR = uᵢ′v′∂y(Uᵢ)

where uᵢ′ is the velocity perturbation in the i direction, v′ is the velocity perturbation in the y direction, Uᵢ is the background velocity in the i direction, and ∂y is the vertical derivative.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(grid);

julia> YSHEAR = TurbulentKineticEnergyEquation.YShearProductionRate(model)
TurbulentKineticEnergyYShearProductionRate KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: shear_production_rate_y_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField")
└── computes: TKE shear production (y)  -uᵢ′v′ ∂_yUᵢ
source
Oceanostics.TurbulentKineticEnergyEquation.TurbulentKineticEnergyZShearProductionRateType
TurbulentKineticEnergyZShearProductionRate(
    u′,
    v′,
    w′,
    U,
    V,
    W;
    grid,
    location
)

Calculate the shear production rate in the model's z direction:

ZSHEAR = uᵢ′w′∂z(Uᵢ)

where uᵢ′ is the velocity perturbation in the i direction, w′ is the vertical velocity perturbation, Uᵢ is the background velocity in the i direction, and ∂z is the vertical derivative.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(grid);

julia> ZSHEAR = TurbulentKineticEnergyEquation.ZShearProductionRate(model)
TurbulentKineticEnergyZShearProductionRate KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: shear_production_rate_z_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField")
└── computes: TKE shear production (z)  -uᵢ′w′ ∂_zUᵢ
source
Oceanostics.TurbulentKineticEnergyEquation.TurbulentKineticEnergyShearProductionRateType
TurbulentKineticEnergyShearProductionRate(
    u′,
    v′,
    w′,
    U,
    V,
    W;
    grid,
    location
)

Calculate the total shear production rate (sum of the shear production rates in the model's x, y and z directions):

SHEAR = XSHEAR + YSHEAR + ZSHEAR = uᵢ′uⱼ′∂ⱼ(Uᵢ)

where XSHEAR, YSHEAR and ZSHEAR are the shear production rates in the x, y and z directions, respectively, uᵢ′ and uⱼ′ are the velocity perturbations in the i and j directions, respectively, and ∂ⱼ is the derivative in the j direction.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(grid);

julia> SHEAR = TurbulentKineticEnergyEquation.ShearProductionRate(model)
TurbulentKineticEnergyShearProductionRate KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: shear_production_rate_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField", "Oceananigans.Fields.ZeroField")
└── computes: total TKE shear production  -uᵢ′uⱼ′ ∂ⱼUᵢ
source