[stanking]

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]

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

[stanking]

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!

[NewtonKF]

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]

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]

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]

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;

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;
end
end
end

index = find(E~=0);
E_index = E(index);
loglog(index,E_index,index,index.^(-5/3));

[ngrube]

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]

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]

thx harishg, this is my email : stanking1983@hotmail.com

I'm waiting for your code!!!

[nathanpetrelli]

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

[vishwa]

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

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

[germangsilva]

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

[NicoG]

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

[uttom]

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

[jet]

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!

[lakeat]

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]

could you please also send me the Matlab code?

jitahz@yahoo.com


tnx!

[risku9]

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!

[ghoshi1983]

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