[Samoza]

Hello everyone,
I'm doing my first "numerical analysis" in turbulence and I'm trying to plot all the terms of the TKE (Turbulent Kinetic Energy) Budget in a single graph to compare the curves with the data obtained by Moser, Kim, and Mansour (1999). The formula for the TKE Budget is shown in the attached image (TKE_eq (2).JPG).
The flow under consideration is a turbulent channel flow. I have successfully plotted almost all the terms (Production, Turbulence Diffusion, Viscous Diffusion, Dissipation rate), and the results align well with the data from Moser, Kim, and Mansour as you can see in the attached image (TKE Budget.png). However, I am having trouble correctly plotting the Pressure Diffusion term.
The available data includes:

• Three 3D matrices (256x128x128) named u_prime, v_prime, and w_prime, representing the fluctuations of the velocity components along the x, y, and z directions, respectively.
• A 3D matrix p_prime (256x128x128) representing the pressure field.
• Three 3D matrices x, y, and z (256x128x128) representing the spatial coordinates.

The data from Moser, Kim, and Mansour can be accessed at the following website: https://turbulence.oden.utexas.edu/MKM_1999.html. The specific file containing the data for comparison is located at chan180/balances/chan180.kbal.
In the file chan180.kbal, it states “Normalization: U_tau, nu/U_tau,” and the numerical method adopted is also mentioned: “Kim, Moin & Moser, 1987, J. Fluid Mech. vol 177, 133-166,” which can be viewed at the link http://www.oroussel.free.fr/articles/KMM87.pdf.
Given that the formula for Pressure Diffusion is shown in the image, labeled as uppercase pi, how would you write the code to correctly plot the Pressure Diffusion term?
Here is my attempt:

clc;
 
clear;

close all;




% Space Statistics

% Importing cells centers coordinates

ccx = importdata('ccx');

ccy = importdata('ccy');

ccz = importdata('ccz');




% Importing velocity field

filePath = 'U.gz';

tempDir = fullfile(tempdir, 'tempExtraction');

gunzip(filePath, tempDir);

U = importdata(fullfile(tempDir, 'U'));




% Importing pressure field

pField = importdata('p');




% Parameters

nx = 256;

ny = 128;

nz = 128;

delta = 1;

u_tau = 1;

nu = 1/180;

Re = (delta * u_tau) / nu;

rho = 1;




% Allocating values in matrices

u = zeros(nx, ny, nz);

v = zeros(nx, ny, nz);

w = zeros(nx, ny, nz);

x = zeros(nx, ny, nz);

y = zeros(nx, ny, nz);

z = zeros(nx, ny, nz);

p = zeros(nx, ny, nz);

index = 1;

for k = 1:nz

    for j = 1:ny

        for i = 1:nx

            u(i, j, k) = U(index, 1);

            v(i, j, k) = U(index, 2);

            w(i, j, k) = U(index, 3);

            x(i, j, k) = ccx(index);

            y(i, j, k) = ccy(index);

            z(i, j, k) = ccz(index);

            p(i, j, k) = pField(index);

            index = index + 1;

        end

    end

end




% Calculating mean profiles for every y

u_mean = mean(u, [1, 3]);

v_mean = mean(v, [1, 3]);

w_mean = mean(w, [1, 3]);

p_mean = mean(p, [1, 3]);




% Calculating fluctuations fields

u_prime = u - u_mean;

v_prime = v - v_mean;

w_prime = w - w_mean;

p_prime = p - p_mean;




% Extracting position vectors

x_direc = x(:, 1, 1);

y_direc = y(1, :, 1);

z_direc = z(1, 1, :);




% Calculating y_plus

y_plus = (y_direc * u_tau) / nu;




% Calculating differences along directions x, y, and z

dx = diff(x_direc);

dy = diff(y_direc);

dz = diff(z_direc);




% Importing Moser, Kim and Mansour data

fileID = fopen('chan180.kbal', 'r');

chan180kbal = textscan(fileID, '%f%f%f%f%f%f%f%f%f', 'CommentStyle', '#');

fclose(fileID);

y_plusMKM = chan180kbal{2};

pressDiffMKM = chan180kbal{6};




% Pressure diffusion term

mean_pu_prime = mean(p_prime .* u_prime, [1, 3]);

mean_pv_prime = mean(p_prime .* v_prime, [1, 3]);

mean_pw_prime = mean(p_prime .* w_prime, [1, 3]);

dpu_dx = diff(mean_pu_prime) ./ dx;

dpv_dy = diff(mean_pv_prime) ./ dy;

dpw_dz = diff(mean_pw_prime) ./ dz;

pressDiff_pre = -(1/rho) * (dpu_dx + dpv_dy + dpw_dz);

pressDiff_pre = pressDiff_pre(1, :, 1);

pressDiff = NaN(1, 128);

pressDiff(1:127) = pressDiff_pre;




% Plotting pressure diffusion terms

figure;

hold on;

plot(y_plus, pressDiff/Re, 'DisplayName', '$\Pi$');

plot(y_plusMKM, pressDiffMKM, '--', 'DisplayName', '$\Pi_{MKM}$');

xlabel('$y^+$', 'Interpreter', 'latex');

xlim([0 max(y_plusMKM)/2]);

plotTitle = 'Pressure Diffusion Term';

title(plotTitle, 'Interpreter', 'latex');

legend('show', 'Location', 'northeast', 'Interpreter', 'latex');

grid on;

box on;

hold off;

The output of this code is attached (PressDiffTerm.png). I will provide the files via WeTransfer if anyone wants to try performing the calculation: https://we.tl/t-3UWWWrHyo0
For any clarifications, please do not hesitate to ask, and I will do my best to respond.

[sulli72]

Hi Samoza,

Just wanted to ask some clarifying questions and leave some thoughts.

First question: Are these budgets being computed from one snapshot, that is averaged in the spanwise and streamwise directions?

if so, you might see better convergence in the pressure term if multiple snapshots are considered when computing the budgets.

Second question: How has pressure been non-dimensionalized in the original DNS code? You may need to carry this non-dimensionalization factor through the derivation of the pressure diffusion term, resulting in some factor that pre-multiples the derivative terms you compute numerically. I've run into this issue before when working with compressible TKE budgets, where lack of proper non-dimensional scaling can make certain profiles be off by orders of magnitude while others look ok.

Hope this helps, and happy budgeting!

[Samoza]

Originally Posted by sulli72

Hi Samoza,

Just wanted to ask some clarifying questions and leave some thoughts.

First question: Are these budgets being computed from one snapshot, that is averaged in the spanwise and streamwise directions?

if so, you might see better convergence in the pressure term if multiple snapshots are considered when computing the budgets.

Second question: How has pressure been non-dimensionalized in the original DNS code? You may need to carry this non-dimensionalization factor through the derivation of the pressure diffusion term, resulting in some factor that pre-multiples the derivative terms you compute numerically. I've run into this issue before when working with compressible TKE budgets, where lack of proper non-dimensional scaling can make certain profiles be off by orders of magnitude while others look ok.

Hope this helps, and happy budgeting!

Hi sulli72,

Here are my responses to your questions:

First response: I was provided with 22 snapshots of the turbulent flow, and I arbitrarily chose snapshot number 3. Since I didn’t encounter any issues when plotting the other TKE terms using a single snapshot, I thought that wasn’t the problem.

Second response: If by the original DNS code you mean the numerical method used by Moser et al., then unfortunately, I don't know how the pressure was non-dimensionalized. Or rather, I’m not entirely sure. In the file `chan180.kbal`, it states “Normalization: U_tau, nu/U_tau,” which, as written, isn't very clear to me. The numerical method adopted is also mentioned: “Kim, Moin & Moser, 1987, J. Fluid Mech. vol 177, 133-166,” which can be viewed at the link [http://www.oroussel.free.fr/articles/KMM87.pdf](http://www.oroussel.free.fr/articles/KMM87.pdf).

I hope I’ve answered your questions.

[LuckyTran]

dpu_dx = diff(mean_pu_prime) ./ dx;

dpv_dy = diff(mean_pu_prime) ./ dy;

dpw_dz = diff(mean_pu_prime) ./ dz;

should be

dpu_dx = diff(mean_pu_prime) ./ dx;

dpv_dy = diff(mean_pv_prime) ./ dy;

dpw_dz = diff(mean_pw_prime) ./ dz;

[Christofer]

Oh yeah, it's a bit of a pain when you manage to track down all but one term! It's often helped me out with stuff like this.

[Samoza]

Originally Posted by LuckyTran

dpu_dx = diff(mean_pu_prime) ./ dx;

dpv_dy = diff(mean_pu_prime) ./ dy;

dpw_dz = diff(mean_pu_prime) ./ dz;

should be

dpu_dx = diff(mean_pu_prime) ./ dx;

dpv_dy = diff(mean_pv_prime) ./ dy;

dpw_dz = diff(mean_pw_prime) ./ dz;

Oops, you're right. My bad, but still, in the plot, the two curves don't match. I will edit the original post.