[juliom]

Dear Community;
I have read many post about the topics but I still have a question. I am working with LES and I want to verify my LES code. Basically I am doing validation using the channel flow data base from JHUTB. So, I have a 3D flow field: u(x,y,z,t),v(x,y,z,t),w(x,y,z,t). I am interesting in compute the power spectrum at each plane. I have read about defining a point a collect the temporal velocity at a single point; which make sense, but I am more interested over the space rather than time..!!!
I have read about compute the FFT from the velocity, but from Wilcox book, I read that the Spectrum Energy is the FFT of the velocity correlation tensor.
I was planning to use FORTRAN to write the script, but I have read that the majority is using Matlab, so it seems that I will lean towards Matlab. However, I want to have the picture clear before starting!!!

[FMDenaro]

I used to fix a value of y+ and then perform the FFT of the streamwise (x) and spanwise (z) velocity components along Kx and Kz at each x and z stations of the computatonal grid.
Then, I compute the Euu, Evv and Eww by the squared modulus of the Fourier coefficents performing an average along the plane.

This procedure is used for one field at a characteristic time step and is then repeated for a sufficient number of samples to perform also a time averaging

[kaya]

If you're going to compare your spectrum with an already given one, you should try to find the way they made their spectrum otherwise I think its very likely that you'll end up with something shifted on y axis.

I had the exact same question months ago. I was trying to validate my LES solver by using this AGARD database. It was an isotropic turbulence decay spectrum at a specific time and whatever I tried, I could not match my power spectrum.

I had a normalization problem and my solution to that was to compute correlations manually (instead of using ffts) and integrate over bins as defined in Pope's book (turbulent flows) section 6.5, but haven't had time to look into this again. Maybe this wouldn't be the best idea, I am not sure tbh.

ps. do yourself a favor and do this in matlab, fortran is an overkill for such thing.

[juliom]

Dear Professor,
I really appreciate your time replying to my post.!!! I have just a couple of question in regards your variables.
Are you referring to Kinetic Energy when you tiped Kx and Kz or are they the FFT of U (streamwise velocity) and W respectively?
When the average is done over the plane, is the average only for the Fourier coefficients? or from the FFT of the velocity vectors.

Thanks!

[FMDenaro]

Kx = n *2*pi/Lx and Kz = m*2*pi/Lz are the spatial frequencies...

The average is performed on the Fourier coefficients

[juliom]

Thank you very much professor!

[FMDenaro]

here a part of my matlab code:

________________
j = 41
%
%
for k=1:l
% 1) ____________ Calcolo di E_uu (stream-wise)______________
% ff è un vettore di appoggio su cui fare la FFT
for i = 1:np+1, ff(i) = u(i,j,k); end
% Trasformata discreta calcolo dei coefficienti di Fourier
c=fft(ff,np); for q = 1:np, c(q)=c(q)/np; end
% Calcolo modulo quadro. Ciclo su ogni armonica k
for q=1:np, cu(q,k)=c(q)*conj(c(q)); end
%
end
%
for q = 1:np
sum_u = 0;
for k=1:l-1
sum_u = sum_u + 0.5*(cu(q,k) + cu(q,k+1))*dz;
end
ek1_uu(q) = sum_u/lz;
end
%
%
fmax_x = pi/dx; % frequenza massima di campionamento (Nyquist) lungo x
Ncx = np/2-1; % Numero di campioni in k
kx = 2*pi*(0:Ncx)/lx; % Posizione lungo l'asse delle frequenze NB: sono Ncx + 1 componenti
yrif = (kx(1:Ncx+1)).^(-5./3.); % Pendenza retta inerziale
xc = kx(1:Ncx+1); % Posizioni in frequenza

[juliom]

Thank you very much professor, this contribution is very important to me!!!! I really appreciate it!

[juliom]

Professor, form this code I learnt that the Spectrum is done only in one direction (the homogeneous direction). I assumed that this is due to the restriction of the Velocity correlation tensor. AM I right?
I would like to know what does Np stands for? I assumed that this are the number of velocity nodes. Am I right.
Finally professor, I am trying to follow this procedure from Popes textbook.. without luck.. Is there any other texbook you recommend me to follow this explanation?

Thanks!

[juliom]

Dear professor, I have been following the procedure and I have also trying to match it with Popes; chapter 6. However, the procedure I read there was based on the two-point correlation tensor. Thus, in the homogeneous direction i compute 1/2(Phi_ii) = Fourier transform of the R_ii (two point correlation tensor) in the Homogeneous direction this is not more than the kinetic energy for i = j.
Yet, your script computes the FFT from the velocity component (streamwise) and then it computes magnitude from the Fourier Coefficients (line 5).
Then an average between k and k + 1 is done (line 9).
Is it right?
Why are the Fourier coefficient being divided by 'np' (c(q)=c(q)/np)
Finally, what does np stand for?
Thank you very much professor!!

[FMDenaro]

np is the number of cells

[sadsid]

Originally Posted by juliom

Dear professor, I have been following the procedure and I have also trying to match it with Popes; chapter 6. However, the procedure I read there was based on the two-point correlation tensor. Thus, in the homogeneous direction i compute 1/2(Phi_ii) = Fourier transform of the R_ii (two point correlation tensor) in the Homogeneous direction this is not more than the kinetic energy for i = j.
Yet, your script computes the FFT from the velocity component (streamwise) and then it computes magnitude from the Fourier Coefficients (line 5).
Then an average between k and k + 1 is done (line 9).
Is it right?
Why are the Fourier coefficient being divided by 'np' (c(q)=c(q)/np)
Finally, what does np stand for?
Thank you very much professor!!

Have you implemented the code successfully?