# Energy Spectrum !!!

 Register Blogs Members List Search Today's Posts Mark Forums Read

 April 6, 2009, 04:16 Energy Spectrum !!! #1 New Member   Yilei Shi Join Date: Apr 2009 Posts: 5 Rep Power: 15 I have a question about turbulent energy spectrum !!! Now I have the data of position and velocity at grid. the data is like this: X Y Vx Vy I have already used FFT to transform Vx and Vy to V_kx and V_ky. Then I calculate the kx = 2*pi*N / L and kx = 2*pi*N / L. (N = 1,2,....) Here is a question. N is from 0 or from 1??? Now the data is like this: kx ky V_kx V_ky How can I calculate the energy spectrum?? 1. How can I calculate the k? k = sqrt(kx^2+ky^2), is right? 2. V(k) = sqrt(V_kx^2 + V_ky^2), is right?? 3. plot(k,V(K),k^(-5/3)) please tell me the detail of the calculation procedure!!! thanks!!! Last edited by stanking; April 6, 2009 at 06:21.

 April 6, 2009, 05:48 #2 Senior Member     Daniel WEI (老魏) Join Date: Mar 2009 Location: Beijing, China Posts: 689 Blog Entries: 9 Rep Power: 20 I remember I have asked the same question, and there are indeed some good detailed answer to them, so please search the forum and see if you can get them. Daniel __________________ ~ Daniel WEI ------------- Boeing Research & Technology - China Beijing, China Email

 April 6, 2009, 06:04 #3 New Member   Yilei Shi Join Date: Apr 2009 Posts: 5 Rep Power: 15 Sorry, I have searched all of the forum, but didn't find the article of your question about turbulent energy spectrum calculation. So can you give me a link? Or please explain the last 3 steps! Thanks!

 April 6, 2009, 09:45 #4 Member   Newton KF Join Date: Mar 2009 Posts: 36 Rep Power: 16 The book from Pope (turbulent flows) have this kind of information... It is very complex, but has all energy spectrum detais... including the equations for the two point correlation... Good luck...

 April 6, 2009, 21:56 #5 New Member   Yilei Shi Join Date: Apr 2009 Posts: 5 Rep Power: 15 Someone can answer the last 3 questions???? we assume the following are right. 1. How can I calculate the k? k = sqrt(kx^2+ky^2), is right? 2. V(k) = sqrt(V_kx^2 + V_ky^2), is right?? Then we got two 32*32 matrix for k and V(k) (for example each direction has 32 points) * + ! ... The k matrix is like this. It's symmetric along the diagonal. This means k(i,j) = k(j,i). + * @ ... Now I calculate the E(k). It should be E(k) = V(k(i,j)) ^2 + V(k(j,i))^2 ( Because of k(i,j) = k(j,i) ), is right?? ! @ * ... Last edited by stanking; April 7, 2009 at 04:37.

 April 7, 2009, 16:02 #6 New Member   Nathan Grube Join Date: Mar 2009 Posts: 14 Rep Power: 16 Let's call Vx u and call Vy v. let us call the Fourier transforms U(kx,ky)=fft(u(x,y)) and V(kx,ky)=fft(v(x,y)). First I would use a 2D FFT (from FFTW of matlab or whatever you want) to get U(kx,ky) and V(kx,ky). Then you have your spectral info in terms of wavenumber vectors (kx,ky). You probably want a scalar k=sqrt(kx**2+ky**2). To get your info in terms of a scalar k, you loop through kx and ky and look at the energy in ecah mode and add it to the appropriate "k" bin. So, define an array E(k)=[0,0,0,....], loop through kx and ky, compute k from kx and ky, round to an integer k, then do E(k) = E(k) + U(kx,ky)**2+V(kx,ky)**2. This turns your 2D data into the 1D E(k) plots that you usually see. To answer your other question, usually N goes from -M/2 to (M/2)-1 or else from 0 to M where M is the umber of grid points in one direction. This varies depending on FFT package. You need to read the documentation with your particular FFT library or Matlab, etc. Then I would recommend doing an FFT of some simple sinusoids to make sure you get the results you expect and that your understanding of N is correct. I hope that is correct. I am in a hurry and doing this from memory without proofreading. Good luck. PS. (kx,ky)=(0,0) corresponds to mean flow. You may wish to ignore that mode. arvindpj, mgg, callmetao and 6 others like this.

April 10, 2009, 03:21
#7
New Member

Yilei Shi
Join Date: Apr 2009
Posts: 5
Rep Power: 15
thx ngrube, i have tried it. But the result is not good!!
The curve has big oszillation.
Maybe can you check my code?

clc;
clear all;
close all;

E = zeros(1024,1);

L = 1;
N = 1:32;
kx = (2 * pi .* N / L)';
ky = (2 * pi .* N / L)';
kz = (2 * pi .* N / L)';

v = reshape(Data(:,4),32,32,32);
u = reshape(Data(:,5),32,32,32);
w = reshape(Data(:,6),32,32,32);

V = fftn(v);
U = fftn(u);
W = fftn(w);

for i = 1:32,
for j = 1:32,
for m = 1:32,
k = round(sqrt(kx(i)^2 + ky(j)^2 + kz(m)^2));
E(k) = E(k) + abs(V(i,j,m))^2 + abs(U(i,j,m))^2 + abs(W(i,j,m))^2;
end
end
end

index = find(E~=0);
E_index = E(index);
loglog(index,E_index,index,index.^(-5/3));
Attached Images
 result.jpg (16.9 KB, 1307 views)

 April 10, 2009, 20:54 check your k's #8 New Member   Nathan Grube Join Date: Mar 2009 Posts: 14 Rep Power: 16 I don't have Matlab (or its manual) in front of me right now, but I think your kx, ky, kz look wrong. You need to check the documentation to see what the wavenumbers are and how the FFT'd data is arranged because these vary from program to program. If you have 32 points in one direction, then the largest magnitude of kx should be only 16 (if you don't see why, google "Nyquist limit") So, your kx ky kz look like they go from 0 to 32 rather than -15 to 16 (the exact upper and lower limit can be handled different ways since kx=16 is redundant with kx=-16... ). I am not sure what MatLab uses exactly. If you can't find it in the documentation, just play around and try an FFT of sin(5*x) or something to see where k=(+/-)5 is in the resulting output. To see what I mean, suppose N = 8. Some common output orders would be arranged with k's as follows: (-4,-3,-2,-1,0,1,2,3) or (0,1,2,3,-4,-3,-2,-1) or (-3,-2,-1,0,1,2,3,4) etc. Note that in place of -4 I could write +4 and it would be equivalent. Some programs (but probably not MatLab) might even give you (-4,-3,-2,-1,0,1,2,3,4) but then -4 and +4 will each be half of the usual kx=4 contribution. I think that MatLab will use one of my first 3 suggestions. You can just use trial and error if you like.

 April 11, 2009, 15:39 #9 Senior Member   N/A Join Date: Mar 2009 Posts: 189 Rep Power: 16 You need to do fftshift to center the mean flow frequency. Also fft in matlab gives you frequencies from -N/2:N/2-1, scaled accordingly. If you can give your email address, I will send the matlab code that I had used. Amjad Asad, mgg and mmmn036 like this.

 April 12, 2009, 15:32 #10 New Member   Yilei Shi Join Date: Apr 2009 Posts: 5 Rep Power: 15 thx harishg, this is my email : stanking1983@hotmail.com I'm waiting for your code!!!

 April 24, 2009, 23:08 #11 New Member   Nathan Join Date: Apr 2009 Posts: 3 Rep Power: 15 Yilei, I was looking at the same problem actually. Could you please send me the code Harish sent you? I would like to try it out. Would help a lot. Thanks! nathanpetrelli77@yahoo.com Nathan

April 28, 2009, 11:02
#12
Member

vishwanath somashekar
Join Date: Apr 2009
Posts: 41
Rep Power: 15
Hi,
i am interested in calculating the spectrum as well. My geometry however, in not regular..

I am running LES using fluent and have saved the velocity fields in the midplane. what would be the procedure to calculate the spectrum...

I have read the many posts here and in books like pope. It really isn't very clear to me. I would appreciate if you could guide me in calculating the spectrum..

The .dat file that i saved for velocity, contains the following header..

Quote:
 TITLE = "title" VARIABLES = X, Y, Z, "z-velocity", "y-velocity", "x-velocity" ZONE T="midplane", N=63808, E=60912, ET=QUADRILATERAL, F=FEBLOCK
Regards,
vishwa

 April 8, 2010, 14:21 #13 New Member   Germán González Silva Join Date: Mar 2010 Location: Campinas-SP Posts: 20 Rep Power: 15 Yilei, I was looking at the same problem actually. Could you please send me the code Harish sent you? I would like to try it out. Would help a lot. Thanks! germangsilva@gmail.com

 April 16, 2010, 16:31 Info for every one ... #14 New Member   Join Date: Apr 2010 Posts: 1 Rep Power: 0 Yilei or anyone who has the code, it would be very nice for every newbie to have the matlab code posted here (if it is not too long). In fact, Pope explains things in details but example on howto implement such as thing in matlab would be helpfull. Thanks in advance for the help, Nico

 July 4, 2010, 07:50 #15 New Member   Sumon Saha Join Date: Jul 2010 Location: Melbourne Posts: 1 Rep Power: 0 Yilei, I was looking at the same problem actually. Could you please send me the code Harish sent you? I would like to try it out. Would help a lot. Thanks! sumonmebuet@gmail.com

 August 26, 2010, 05:14 calculating energy spectra #16 New Member   Join Date: Nov 2009 Posts: 12 Rep Power: 15 hi, i want to calculate an energy spectrum out of the velocity (or better out of the velocity fluctuations in one direction) measured at one point in a turbulent pipe flow calculated with OpenFoam (LES). i get a velocity value each timestep (0.00004 s). and here is my way to calculate the energy spectra: 1) calculating the autocorrelation of the velocity vector 2) fftshift of the autocorrelation 3) calculating abs(fft(fftshift(signal))) --> the energy spectrum values for plotting the spectra i need a frequency-axis and i don't know how to calculate it! my sampling freq should be 1/0.00004 = 25000 [Hz] ? what should be my first value on the freq axis? what is my delta freq? and what is the last axis value? Thanks!

 August 26, 2010, 10:02 #17 Senior Member     Daniel WEI (老魏) Join Date: Mar 2009 Location: Beijing, China Posts: 689 Blog Entries: 9 Rep Power: 20 I guess, when we based on Taylor frozen hypothesis, the fluid particle will move a certain distance in the time step = 0.00004s, and this distance is equivalent to spacial flow structure length in its homogeneous direction (w/o Taylor frozen hypothesis), and hence it is related to the wave number. Just some thoughts. randolph likes this. __________________ ~ Daniel WEI ------------- Boeing Research & Technology - China Beijing, China Email

 January 16, 2011, 11:57 code request #18 New Member   Join Date: Jan 2011 Posts: 1 Rep Power: 0 could you please also send me the Matlab code? jitahz@yahoo.com tnx!

 March 13, 2012, 13:46 Doubts!! #19 New Member   Clara Join Date: Sep 2010 Location: Bruxelles Posts: 16 Rep Power: 14 Hi everyone!! I've been visiting this web for a long time trying to understand how to plot the turbulence spectra without good results!! Can anybody send me the code to see what am I doing wrong? gsclara@gmail.com Thank you!

 May 4, 2012, 12:55 Calculating Energy Spectra #20 New Member   Join Date: Jun 2009 Posts: 21 Rep Power: 15 Hi Guys, I am also having problems with plotting energy spectrum. Could somebody plz send me the code (zietapk@hotmail.com) also so that i can compare it with mine and figure out what i am doing wrong. Thanks in advance!!

 Tags energy spectrum, fft