matlab code for tremine pressure diffuseon in tke budget

[Samoza]

Bye to all,

I am trying to plot all the terms of the tke budget in a single chart to compare the curves with the data obtained from moser, kim and mansour (1999). the formula of the tke budget is shown in the attached image:

https://imgur.com/d9mbt1

the considered flow is a turbulent flow in the channel. I managed to plot almost all terms (production, turbulence diffuseon, viscous diffuseon, dissipation rate), obtaining consistent results with the data of moser, kim and mansour as seen in the image:

https://imgur.com/9w9gfz5

However, I am having difficulty plotting correctly the term pressure diffuseon.

available data include:

- three matrices 3d (256x128x128) called `u_prime`, `v_prime` and `w_prime`, which represent the fluctuations of the speed components along the directions x, y and z, respectively.
- a 3d `p_prime` matrix (256x128x128) which represents the pressure field.
- three matrices 3d `x`, `y` and `z` (256x128x128) representing the spatial coordinates.

moser, kim and mansour data can be found on the following website: dns data for turbulent channel flow up to friction reynolds number of 590 . the specific file containing the data to compare is located at the `chan180/balances/chan180.kbal` path.

here is my attempt:

clc;
clear;
close all;


% space statistics
% importing cells centers coordinates
ccx = importdata(fullfile('velfieldschannel_corsoturb', 'ccx'));
ccy = importdata(fullfile('velfieldschannel_corsoturb', 'ccy'));
ccz = importdata(fullfile('velfieldschannel_corsoturb', 'ccz'));


% importing velocity field
choosesnapshot = 3; % select a snapshot (from 1 to 22)
filepath = fullfile('velfieldschannel_corsoturb', num2str(choosesnapshot), 'u.gz');
tempdir = fullfile(tempdir, 'tempextraction');
gunzip(filepath, tempdir);
u = importdata(fullfile(tempdir, 'u'));


% importing pressure field
pfield = importdata(fullfile('velfieldschannel_corsoturb', '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 differencies along directions x, y and z
dx = diff(x_direc);
dy = diff(y_direc);
dz = diff(z_direc);


% importing moser, kim and mansour data
filepath = fullfile('chan180', 'balances', 'chan180.kbal');
fileid = fopen(filepath, 'r');
chan180kbal = textscan(fileid, '%f%f%f%f%f%f%f%f%f', 'commentstyle', '#');
fclose(fileid);
y_plusmkm = chan180kbal{2};
dissipationmkm = chan180kbal{3};
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]);
(p _ prime. * w _ prime, [1 3]);


dpu_dx = diff(mean_pu_prime) ./ dx;
dpv_dy = diff(mean_pu_prime) ./ dy;
dpw_dz = diff(mean_pu_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', '[imath]\pi[/imath]');
plot(y_plusmkm, pressdiffmkm, '--', 'displayname', '[imath]\ pi _ {mkm}[/imath]');
xlabel('[imath]y^+[/imath]',' 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 such code is the following image:

https://imgur.com/twbq9cs

and as you see the two curves don't match.

since the formula for pressure diffuseon is shown in the image, how would you write the code to properly plot the term pressure diffuseon?

for any clarification, do not hesitate to ask and I will try to answer.