how to calculate turbulent kinetic energy of resolved and sgs components.

[dshen1993]

I want to calculate the total turbulent kinetic energy from resolved and sgs components.

Here is code i used to calculate the TKE.

u_center=@at (Center, Center, Center) model.velocities.u
v_center=@at (Center, Center, Center) model.velocities.v
w_center=@at (Center, Center, Center) model.velocities.w

U = Field(Average(u_center, dims=(1, 2)))
V = Field(Average(v_center, dims=(1, 2)))
W = Field(Average(w_center, dims=(1, 2)))

u_prime=@at (Center, Center, Center) Field(u_center-U)
v_prime=@at (Center, Center, Center) Field(v_center-V)
w_prime=@at (Center, Center, Center) Field(w_center-W)

TKE_pointwise=Field(0.5*(u_prime^2 + v_prime^2 + w_prime^2))
TKE=Average(TKE_pointwise, dims=(1, 2))
####################################################
The TKE from this calculation is too small, less than 10^-32.

So, if i want to get the internal tke that is used to calculated the eddy viscosity in les turbulence closure. How to get this internal tke.

By the way, how do i calculate the tke from resolved component?

Thanks.

[glwagner]

@dshen1993, we suggest that you use https://github.com/tomchor/Oceanostics.jl for this computation.

Also, you likely want to interpolate to cell centers after computing u_prime^2, etc. The reason is because this corresponds to a discrete TKE that is feasible to derive from the discrete momentum equations. The "other" TKE computed by averaging first is more difficult to compute.

Also note that the average of w should be 0 if it is valid to average over the entire horizontal dimension (eg with dims=(1, 2)).

Oceanostics will do this for you. It will look something like

U = Field(Average(u, dims=(1, 2)))
V = Field(Average(v, dims=(1, 2)))

u_prime = u - U
v_prime = v - V

e = @at (Center, Center, Center) sqrt(u_prime^2 + v_prime^2 + w^2) / 2
E = Field(Average(e, dims=(1, 2)))
compute!(E)

As for why your TKE is too small, I am not sure. We need all of your code to be able to reproduce the issue.

PS use triple backticks "```" to enclose your code so that it formatted correctly (you can also opt for syntax highlighting this way).

[glwagner]

Here's a sample calculation:

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))

u = XFaceField(grid)
v = YFaceField(grid)
w = ZFaceField(grid)

set!(u, (x, y, z) -> 1 + randn())
set!(v, (x, y, z) -> 2 + randn())
set!(w, (x, y, z) -> randn())

U = Field(Average(u, dims=(1, 2)))
V = Field(Average(v, dims=(1, 2)))

u_prime = u - U
v_prime = v - V

e = @at (Center, Center, Center) sqrt(u_prime^2 + v_prime^2 + w^2) / 2
E = Field(Average(e, dims=(1, 2)))
compute!(E)

executing this gives me

julia> compute!(E)
1×1×4 Field{Nothing, Nothing, Center} reduced over dims = (1, 2) on RectilinearGrid on CPU
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (1, 1, 4)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── operand: Average of BinaryOperation at (Center, Center, Center) over dims (1, 2)
├── status: time=0.0
└── data: 1×1×10 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:7) with eltype Float64 with indices 1:1×1:1×-2:7
    └── max=1.00921, min=0.883394, mean=0.937948