CFD Online Discussion Forums

CFD Online Discussion Forums (
-   Main CFD Forum (
-   -   Energy Spectrum !!! (

stanking April 6, 2009 04:16

Energy Spectrum !!!
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!!!

lakeat April 6, 2009 05:48

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.


stanking April 6, 2009 06:04

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!


NewtonKF April 6, 2009 09:45

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...

stanking April 6, 2009 21:56

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??
! @ * ...

ngrube April 7, 2009 16:02

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.

stanking April 10, 2009 03:21

1 Attachment(s)
thx ngrube, i have tried it. But the result is not good!!
The curve has big oszillation.
Maybe can you check my code?

clear all;
close all;

Data = load('tg3dn');

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;

index = find(E~=0);
E_index = E(index);

ngrube April 10, 2009 20:54

check your k's
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.

harishg April 11, 2009 15:39

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.

stanking April 12, 2009 15:32

thx harishg, this is my email :

I'm waiting for your code!!!

nathanpetrelli April 24, 2009 23:08


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.


vishwa April 28, 2009 11:02

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..


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

germangsilva April 8, 2010 14:21


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.

NicoG April 16, 2010 16:31

Info for every one ...
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,


uttom July 4, 2010 07:50


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.

jet August 26, 2010 05:14

calculating energy spectra
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?


lakeat August 26, 2010 10:02

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.

jitahz2 January 16, 2011 12:57

code request
could you please also send me the Matlab code?


risku9 March 13, 2012 14:46

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?

Thank you! :)

ghoshi1983 May 4, 2012 12:55

Calculating Energy Spectra
Hi Guys,
I am also having problems with plotting energy spectrum. Could somebody plz send me the code ( also so that i can compare it with mine and figure out what i am doing wrong. Thanks in advance!!

All times are GMT -4. The time now is 00:22.