[Samoza]

Hello everyone,
I am having trouble with my very first numerical analysis. The problem is about calculating the dissipation term in the TKE budget.

Consider a rectangular channel like in the image below. Let delta be = 1:

The meshgrid for this channel consists of 256x128x128 points along the x, y, and z directions, respectively. At each point, 6 values are stored: 3 velocity components (u along x, v along y, and w along z) and 3 spatial coordinates (x, y, and z). Therefore, there are 6 matrices of size 256x128x128: 3 velocity matrices (u(:, :, :), v(:, :, :) , and w(:, :, :)) and 3 position matrices (x(:, :, :), y(:, :, :), and z(:, :, :)). Next, we calculate the velocity fluctuation matrices u', v', and w' (denoted as `u_prime`, `v_prime`, and `w_prime`, respectively), which are simply the difference between the velocity matrices and their mean values: u' = u - u_mean, v' = v - v_mean, and w' = w - w_mean.
Assume the flow is steady and homogeneous along the x and z directions, so the partial derivatives with respect to x and z are zero.
The dissipation term in the TKE budget is given by the formula:

and I guess the formula calculated along the y direction becomes like this:

Assume \nu = 1/180. Here is my attempt to implement this equation:

u_mean = mean(u, [1, 3]); %calculating u mean along the y-direction
u_prime = u - u_mean; %calculating u'
[~, du_dy_prime, ~] = gradient(u_prime); %calculating the partial derivative of u' with respect to y
epsilon = - nu * mean(du_dy_prime .* du_dy_prime, [1, 3]); %calculating the dissipation term

To validate my results, I need to compare them with data from a study by Moser, Kim, and Mansour. Here’s the comparison: the curves are plotted along the y-axis. The blue curve represents my results, while the orange curve corresponds to the data from the study by Moser, Kim, and Mansour.

All other terms in the TKE budget match well, except for the dissipation term, as you can see. I’m confident that the input data for the dissipation calculation are accurate, but I suspect the issue lies in the implementation of the formula. This assumption is based on the fact that all the other terms align closely with the data from Moser, Kim, and Mansour.

So, finally, how would you implement the dissipation function ε in light of the information provided?
Please feel free to ask for any further clarifications if needed.

[FMDenaro]

A part of my very old Matlab code (non uniform grid along y):


% ------ Calcolo dissipazione:
% SGS =- tau_ij*S_ij=2*ni_sgs_ij*S_ij*s_ij
% molecolare=(2/Re_tau)*S_ij*s_ij
%


cx1 = 0.5*dx; cz1 = 0.5*dz;
for j = 2:ny-1,
for k = 2:nz-1, for i = 2:nx-1,
%
ux = cx1*(u(i+1,j,k)-u(i-1,j,k)); vx = cx1*(v(i+1,j,k)-v(i-1,j,k)); wx = cx1*(w(i+1,j,k)-w(i-1,j,k));
%
tj= -(yg(j-1)^2-yg(j+1)^2-2*yg(j)*yg(j-1)+2*yg(j)*yg(j+1))/(yg(j-1)*yg(j)^2- ...
yg(j-1)*yg(j+1)^2+yg(j)*yg(j+1)^2+yg(j+1)*yg(j-1)^2-yg(j+1)*yg(j)^2-yg(j)*yg(j-1)^2);
tjm1=-(-yg(j)^2-yg(j+1)^2+2*yg(j)*yg(j+1))/(-yg(j)*yg(j+1)^2-yg(j-1)*yg(j)^2+ ...
yg(j-1)*yg(j+1)^2+yg(j)*yg(j-1)^2-yg(j+1)*yg(j-1)^2+yg(j+1)*yg(j)^2);
tjp1=-(yg(j)^2+yg(j-1)^2-2*yg(j)*yg(j-1))/(yg(j+1)*yg(j)^2-yg(j+1)*yg(j-1)^2+ ...
yg(j)*yg(j-1)^2+yg(j-1)*yg(j+1)^2-yg(j-1)*yg(j)^2-yg(j)*yg(j+1)^2);
%
uy = tj*u(i,j,k)+tjm1*u(i,j-1,k)+tjp1*u(i,j+1,k);
vy = tj*v(i,j,k)+tjm1*v(i,j-1,k)+tjp1*v(i,j+1,k);
wy = tj*w(i,j,k)+tjm1*w(i,j-1,k)+tjp1*w(i,j+1,k);
%
uz = cz1*(u(i,j,k+1)-u(i,j,k-1)); vz = cz1*(v(i,j,k+1)-v(i,j,k-1)); wz = cz1*(w(i,j,k+1)-w(i,j,k-1));
%
traccias3 = (ux + vy + wz)/3.d0;
s11 = ux-traccias3; s22 = vy-traccias3 ; s33 = wz-traccias3;
s12 = 0.5d0*(uy+vx);s13 = 0.5d0*(uz+wx); s23 = 0.5d0*(vz+wy);
%
epsilon(i,j,k)= 2*sgsqm(i,j,k)*(s11*s11+s22*s22+s33*s33+2.d0*(s12* s12+s13*s13+s23*s23));
epsilon_m(i,j,k)= 2*(s11*s11+s22*s22+s33*s33+2.d0*(s12*s12+s13*s13+s 23*s23))/reynolds;
end, end, end

[Samoza]

Thanks! (Grazie mille, adesso provo ad adattarlo al mio codice e vedo se ne vengo fuori)