Code MATLAb tremine Diffusion de pression dans le budget TKE

[Samoza]

Au revoir à tous,

J'essaie de tracer tous les termes du budget en un seul graphique pour comparer les courbes avec les données obtenues de moser, kim et mansour (1999). la formule du budget est indiquée dans l'image ci-jointe:

https://imgur.com/d9mbt1

le débit considéré est un débit turbulent dans le canal. J'ai réussi à tracer presque tous les termes (production, diffusion de turbulences, diffusion visqueuse, débit de dissipation), en obtenant des résultats cohérents avec les données de moser, kim et mansour telles que vues dans l'image:

https://imgur.com/9w9gfz5

Cependant, j'ai de la difficulté à tracer correctement le terme diffuseur de pression.

Les données disponibles comprennent:

- trois matrices 3d (256x128x128) appelées `u_prime`, `v_prime` et `w_prime`, qui représentent les fluctuations des composantes de vitesse le long des directions x, y et z, respectivement.
- une matrice 3d `p_prime` (256x128x128) qui représente le champ de pression.
- trois matrices 3d `x`, `y` et `z` (256x128x128) représentant les coordonnées spatiales.

Les données moser, kim et mansour se trouvent sur le site web suivant: dns pour un débit turbulent jusqu'au nombre de réynolds de frottement de 590 . le fichier spécifique contenant les données à comparer est situé sur le chemin `chan180/balance/chan180.kbal`.

Voici ma tentative :

ccc;
clair;
fermer tout;


% de statistiques spatiales
% des coordonnées des centres de cellules importatrices
ccx = importdata(fullfile('velfieldschannel_corsoturb', 'ccx'));
ccy = importdata(fullfile('velfieldschannel_corsoturb', 'ccy'));
ccz = importdata(fullfile('velfieldschannel_corsoturb', 'ccz'));


% du champ de vitesse d'importation
sélectionne un instantané = 3; % sélectionne un instantané (de 1 à 22)
fichierpath = fullfile('velfieldschannel_corsoturb', num2str(choosesnapshot), 'u.gz');
tempdir = fullfile(tempdir, 'tempextraction');
gunzip(filepath, tempdir);
u = importdata(fullfile(tempdir, 'u'));


% du champ de pression d'importation
pfield = importdata(fullfile('velfieldschannel_corsoturb', 'p'));


Paramètres %
nx = 256;
ny = 128;
nz = 128;
delta = 1;
u_tau = 1;
nu = 1/180;
re = (delta * u_tau) / nu;
rho = 1;


% des valeurs attribuées dans les matrices
u = zéros (nx, ny, nz);
v = zéros (nx, ny, nz);
w = zéros (nx, ny, nz);
x = zéros (nx, ny, nz);
y = zéros (nx, ny, nz);
z = zéros (nx, ny, nz);
p = zéros (nx, ny, nz);
indice = 1;
pour k = 1:nz
pour j = 1:ny
pour 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);
indice = indice + 1;
fin
fin
fin


% calcul des profils moyens pour chaque y
u_moyen = moyenne(u, [1, 3]);
v_moyen = moyenne(v, [1, 3]);
w_medium = moyenne(w, [1, 3]);
p_medium = moyenne(p, [1, 3]);


% calcul des champs de fluctuations
u_prime = u - u_moyen;
v_prime = v - v_moyen;
w_prime = w - w_moyen;
p_prime = p - p_moyen;


Pourcentage de vecteurs de position extractibles
x_direc = x, 1, 1);
y_direc = y(1, :, 1);
z_direc = z(1, 1, : );


% de calcul y_plus
y_plus = (y_direc * u_tau) / nu;


% calcul des différences selon les directions x, y et z
dx = diff(x_direc);
dy = diff(y_direc);
dz = diff(z_direc);


% d'importation de données de moser, kim et mansour
fichierpath = fullfile('chan180', 'balance', 'chan180.kbal');
fichier = 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};
pressiondiffmkm = chan180kbal{6};


Durée de diffusion de la pression en %
moyenne_pu_prime = moyenne(p_prime .* u_prime, [1 3]);
moyenne_pv_prime = moyenne(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]«);
xétiquette( '[imath]Y^+[/imath]',' interprète ',' latex ');
xlim ([0 max(y_plusMKM)/2]);
timing = 'terme de diffusion de pression';
titre (plottitle, 'interprète', 'latex');
légende('show', 'location', 'nord-est', 'interprète', 'latex');
sur le réseau;
boîte allumée;
Accroche-toi;

la sortie de ce code est l'image suivante:

https://imgur.com/Twbq9cs

Et comme vous voyez, les deux courbes ne correspondent pas.

Puisque la formule pour le diffuseur de pression est affichée dans l'image, comment écrire le code pour tracer correctement le terme diffuseur de pression?

pour toute clarification, n'hésitez pas à demander et je vais essayer de répondre.