MATLAb tremine code Pressão Difusão em Orçamento TKE

[Samoza]

Adeus a todos.

Estou tentando traçar todos os termos do orçamento do tke em um único gráfico para comparar as curvas com os dados obtidos de muser, kim e mansour (1999). a fórmula do orçamento do tke é mostrada na imagem anexada:

https://imgur.com/D9mbt1

o fluxo considerado é um fluxo turbulento no canal. Consegui traçar quase todos os termos (produção, difusão de turbulência, difusão viscoso, taxa de dissipação), obtendo resultados consistentes com os dados de muser, kim e mansour como visto na imagem:

https://imgur.com/9w9gf

No entanto, estou tendo dificuldade em traçar corretamente o termo difusão de pressão.

dados disponíveis incluem:

- três matrizes 3d (256x128x128) chamado `u_prime`, `v_prime` e `w_prime`, que representam as flutuações dos componentes de velocidade ao longo das direções x, y e z, respectivamente.
- uma matriz 3d `p_prime` (256x128x128) que representa o campo de pressão.
- três matrizes 3d `x`, `y` e `z` (256x128x128) representando as coordenadas espaciais.

moser, kim e dados do mansour podem ser encontrados no seguinte site: dns dados para o fluxo de canal turbulento até reynolds de atrito número de 590 . o arquivo específico contendo os dados para comparar está localizado no caminho `chan180/balances/chan180.kbal`.

Aqui está a minha tentativa:

clc;
claro;
fechar todos;


% estatísticas espaciais
% importando coordenadas de centros de células
ccx = import (fullfile('velfieldschannel_corsoturb', 'ccx'));
ccy = importdata(fullfile('velfieldschannel_corsoturb', 'ccy'));
ccz = import (fullfile('velfieldschannel_corsoturb', 'ccz');


% campo de velocidade de importação
escolhesnapshot = 3; % seleciona um instantâneo (de 1 a 22)
filepath = fullfile('velfieldschannel_corsoturb', num2str(choosesnapshot), 'u.gz');
tempdir = fullfile(tempdir, 'tempextraction');
gunzip (filepath, tempdir);
u = importdata(fullfile(tempdir, 'u'));


% campo de pressão de importação
pfield = importdata(fullfile('velfieldschannel_corsoturb', 'p'));


% dos parâmetros
nx = 256;
ny = 128;
nz = 128;
delta = 1;
= 1;
n = 1/180;
re = (delta * u_tau) / nu;
rho = 1;


% alocando valores em matrizes
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);
índice = 1;
para k = 1:nz
para j = 1:ny
para 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);
índice = índice + 1;
fim
fim
fim


% calculando perfis médios para cada y
u_mean = média (u, [1, 3]);
v_mean = média (v, [1, 3]);
w_mean = média (w, [1, 3]);
p_mean = média (p, [1, 3]);


% campos de cálculo das flutuações
u_prime = u - u_mean;
v_prime = v - v_mean;
w_prime = 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, : );


% calculando y_plus
y_plus = (y_direc * u_tau) / nu;


% calculando diferenças ao longo de direções x, y e z
dx = diff(x_direc);
dy = diff(y_direc);
dz = diff(z_direc);


% importação de dados de muser, kim e mansour
filepath = fullfile('chan180', 'balances', 'chan180.kbal');
fileid = fopen(filepath, 'r');
chan180kbal = textcan(fileid, %f%f%f%f%f%f%f%f%f%f, 'commentstyle', '#');
fclose(fileid);
y_plusmkm = chan180kbal{2};
dissipaçãomkm = chan180kbal{3};
pressdiffmkm = chan180kbal{6};


% termo de difusão da pressão
significa_pu_prime = média (p_prime .* u_prime, [1 3]);
significa_pv_prime = média (p_prime.* v_prime, [1 3]);
Primeiro. [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;




% plotando termos de difusão de pressão
figura;
Espera.
plot(y_plus, pressdiff/re, 'displayname', '[imath]\pi[/imath]');
plot(y_plusmkm, pressdiffmkm, '-', 'displayname', '[imath]\ pi _ {mkm}[/imath]');
xlabel( '[imath]Sim.[/imath]',' intérprete ',' látex ');
xlim ([0 max(y_plusMKM)/2]);
plottitle = termo de difusão de pressão;
title(plottitle, 'interpreter', 'latex');
legend('show', 'location', 'northeast', 'interpreter', 'latex');
grade sobre;
caixa em;
Espera.

a saída de tal código é a seguinte imagem:

https://imgur.com/Twbq9cs

e como você vê as duas curvas não combinam.

uma vez que a fórmula para difusão de pressão é mostrada na imagem, como você escreveria o código para traçar corretamente o difusão de pressão do termo?

para qualquer esclarecimento, não hesite em perguntar e vou tentar responder.