Thanks for the code!
However, I'm still not clear on how I should calculate the spectrum for a 3d, periodic box with homogeneous isotropic turbulence.
Is it really just averaging over all the directions? However, the result is wrong.
Hence, I did it like the following


%Number of points in each direction
np=64;

%

%Box Length
L=1.;

%Minimal wave number
k0=2.*pi/L;

%Wave numbers
kx=k0*(-np/2:1:np/2-1);
ky=k0*(-np/2:1:np/2-1);
kz=k0*(-np/2:1:np/2-1);


ux=reshape(u(:,1),np,np,np);
uy=reshape(u(:,2),np,np,np);
uz=reshape(u(:,3),np,np,np);

Ux=fftn(ux)/np;
Ux=fftshift(Ux);
Uy=fftn(uy)/np;
Uy=fftshift(Uy);
Uz=fftn(uz)/np;
Uz=fftshift(Uz);


%% 3D example


%Wave number for 1d plot
nK=np; %This should be np/2
k_1d=k0.*(1:nK);



% Computing the Energy for each wave-number
% and sorting everything into Bins centered around n*k0 for n=1,..,N/2
E=zeros(nK,1);
E_t=zeros(nK,1);
for i=1:np %k
for j=1:np %ky
for m=1:np %kz


k_abs =sqrt(kx(i)^2+ky(j)^2+kz(m)^2 );
for n=1:numel(k_1d)

%Sort into appropriate Wavenumber bin and sum it up
if( k_abs >= (k0*(n-0.5)) && k_abs < (k0*(n+0.5) ) )
E(n)=E(n)+0.5*(abs(Ux(i,j,m))^2+abs(Uy(i,j,m))^2+a bs(Uz(i,j,m))^2);
end

end


end %kz
end %ky
end %kx
loglog(k_1d,E_t), hold on

However, I only get the 'correct' shape if I include all the wavenumber in the box and not just the one's in the sphere with radius pi*np/L.
Any idea what I might be doing wrong?


Any help would greatly be appreciated!

fft will give you the symmetry spectrum about the y axis after fftshift, half that (positive) is ok. The isotropic turbulent flow would have to average over the whole domain, I think a volume average is more appropriate than an algebra average, <u>=add(u*dx*dy*dz)/volume. In fact, in my case, I realize that a plane average is more appropriate than algebra average too.

Dear all,

I am also willing to plot Energy spectra for my supersonic turbulent boundary layer (STBL) DNS data. I have some questions:
1- if this Energy spectra should be in time or space?
2- if it is in time what would be the characteristic time for gathering the data? if it is in space what would you suggest for space dimensions in a 3D Flat plate simulation?
3- If any of you can send me the code for calculating the Energy spectra? (my email address is: msshadloo@google.com)

Thank you in advance,
Best,
Mostafa

Data and Energy Spectrum
 

I have a code written to calculate a 2D energy spectrum. However, I would like to confirm my code is working. Does anyone have 2D velocity data and an energy spectrum in which i can confirm my code is working properly.

my email address is
rehager.edu@gmail.com

Thank you

I have 3D velocity fields

Well I don't have a code that works for 3D yet, but I was planning to extend my code eventually so any data will help. That would be great if you have some data I can test my code on. My email is rehager.edu@gmail.com
Thank you!

Hi, I'm having some trouble too, I'd appreciate if someone could send me the matlab code

david.t.fitzgerald@gmail.com

Thanks!

Hi Mr. Harish, could you please share the Matlab script to compute the PSD from a 3D flow field.?? I would appreciate it!!!
THanks!

Dear Mr. Harish,
I am doing the same task you discussed in the forum. I read that you shared the script with other members. Could you please share it with me too. I will be very thankful with you all!!...

Energy spectrum - matlab code
 

Yilei or anyone who has the code,
I was looking at the same problem and the matlab code would be really helpful! Could you post it here or send it to cfdstudent07@outlook.com ?
Thanks for your collaboration!

someone should do a wikipage for this.

Spectral Density Calculation
 

Hi everyone,

I followed the explanation in Roland B. Stull book (An Introduction for Boundary Layer Meteorology), there you can find a really good example with which you can check that you calculate correctly the Spectral Density whose integral should the variance of your signal.

In my case I calculated from the experiments, here it is my code:

U = WS-mean(WS); % WS: Wind Speed
n = length(U)/2;
Fs = 1/10; % sampling frequency
Dt = 1/Fs;
tt = n*Dt;
t = 0:1/Fs:n/Fs;
nfft=2^nextpow2(length(U));
Fu = fft(U,nfft);
Fhu = Fu(1:end/2);
f=Fs/2*linspace(0,1,nfft/2);
Eu = abs(Fhu).^2;

%%%%%%%%%% Plot loglog
hfig1=figure(1);
loglog(f,Eu)
xlabel('Frequency [Hz]','fontsize',fs)
ylabel('Power Spectral Density','fontsize',fs)
legend('Sprectra','-2/3 slope line')
set(gca,'fontsize',fs-2)
title('Power Spectral Density','fontsize',fs)

Hope this helps..... =)

Clara

Could anyone send the code to fperezsegovia@gmail.com ?

Thank you in advance!!!!

Look at Professor Denaro code at:
http://www.cfd-online.com/Forums/mai...tml#post565863

You will need to add just few lines but the bulk of the code is right there.!!!

Density of state of a 3D lattice
 

Hi,

I am working with matlab and I have to calculate the density of state of a 3D lattice as the fourier transform of the velocity autocorrelation function. I havethe velocity vector inthe 3 drections (vx,vy,vz) and the time vector t (0:50000:6000000). I found the autocorrelation velocity vector V. Then I do abs(fftshiftV) but the resut is not good. How to build the vector frequency and the fourier tranfor?
I don't know also how to plot the fourier spectrum wrt the frequencies in THz. If you have any code that will be helpful please email me.
Thanks a lot for the help.
Jess1

Jess1, look at Professor Denaro script, this is everything you need.

http://www.cfd-online.com/Forums/mai...tml#post565863

But in that code the energy spectrum is computed based on the velocity field by computing the kinetic energy. Your case is different though, because you are using the two point correlation. I can share with you my code (that I wrote based on prof Denaro script, I need to acknowledge his contribution) based on Popes' book.

Turbulence energy spectrum for a 3D curved pipe
 

Hi,
I need to make turbulence energy spectrum for a 3D curved pipe. I would like to see how LES is able to capture the turbulence flow details with my current mesh.
I followed many threads and posts. However, there are still many questions in my mind. The most important one is: should I store the data that belongs to a cross section (a plane, e.g. outlet) or entire domain (a volume e.g. entire pipe) or along the curved pipe (a line e.g. the pipe centerline)? which one is preferred for my case? I think I should use a plane but since people use a domain, I am confused now!
Then, since my case is 3D, should I store all three velocity components (u,v,w) of the cells? Do I need to calculate fluctuating velocity components of the cells later in my matlab code?
Please someone help me with this issue?
Many thanks in advance

My guess is that your are using cylindrical coordiantes, I have not seen an Energy spectrum under that coordinates.
Also, remember that you need to compute the FFT to the data set which is only possible with periodic boundary conditions that I think is not your case. Off course you can compute it for non-periodic BC but the procedure is quite cumbersome and Pope's book is a very good reference.
Nonetheless, I am looking forward to read professor Denaro answer.
If you were using Cartesian coordinates, then you need to average the velocity field by planes and you will use that average to compute the spectrum. Like I said in a previous post, professor Denaro explained it very well before.

If you have inflow/outflow BC.s, then the only periodic direction is the azimuthal one. Thus, I suppose you can proceed as follow:
1) establish some longitudinal positions far from the inflow;
2) for each of these positions fix some radial values r+ (note that r+=0 is the wall);
3) perform the FFT of longitudinal, radial and azimuthal velocity components along the azimuthal direction theta [0:2pi];
4) perform the ensemble spatial average
5) repeat steps 1)-4) for several time units
6) perform the ensemble temporal average

Dear Mr. juliom and Mr. Filippo
Thank you very much for your reply.
Let me explain my case a bit more.
I am using Fluent and my inlet and outlet BCs are velocity-inlet and pressure-outlet. Hopefully, I am using Cartesian coordinates as my case is a helical pipe which has 3 turns (it is long and it has pitch and torsion). I assume it is somehow impossible to apply cylindrical coordinate for such a geometry.
On the other hand, because of existence of secondary flow effects, periodic BC is not my case. So if you monitor streamlines at a certain cross-section far from the inlet you will see the pairs of secondary vortices and Dean Vortices.
What I am trying to do is to find a relationship between turbulence energy and heat transfer from the wall at any point far from the inlet using Turbulence Energy Spectrum.
This is my understanding from what Mr Filippo said:
1. I establish a cross section somewhere far from the inlet (Can I do this with one cross-section only?).
2. As I want to test capability of LES in capturing turbulence details, Turbulence Energy Spectrum of boundary layer is as important as that of mean flow area. So am I allowed to fix some radial values r+? If I do this, I assume that I will lose high frequency part of the spectrum which is associated with the low energy containing eddies. Won’t I?
3. I will have to perform the FFT for x-, y-, z-velocity components. I should use 3 velocity components in the same energy spectrum. By doing this, I will get 1 energy spectrum using 3 velocity components or 3 separate energy spectrum diagrams? Please correct me if I am mistaken.
4. I should take a spatial average over all the cells in the cross section which gives me a single value for each velocity component. I assume that this is to calculate fluctuating velocity component for all the cells? Isn’t it?
5. Repeat steps 1-4 for several time steps, e.g. 2s, 2.5s, 3s, 3.5, 4s, 4.5s, 5s?
6. I should average over those time steps.
I really appreciate your contribution.