# Energy Spectrum !!!

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 LinkBack Thread Tools Display Modes
April 20, 2014, 22:58
#41
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro just consider a 3d Cartesian mesh and (j,k) some fixed positions. This is some lines of my Matlab code for the energy spectrum % ____________ E_uu (stream-wise)______________ % ff is a dump array for using the FFT of the field u(i,j,k) for i = 1:np+1, ff(i) = u(i,j,k); end % Discrete transform: computation of the Fourier coefficient c=fft(ff,np); for q = 1:np, c(q)=c(q)/np; end % Computation of the modulus squared. Cycle for all the wavenumbers for q=1:np, cu(q,k)=c(q)*conj(c(q)); end Then, you must simply plot your power spectrum. Generally, the computed spectrum are also computed along the whole homogeneous directions (e.g. for all the k) and then averaged to have a statistically meaning.
Hi,Filippo
Sorry for the second reply.I followed your advice, here is my code in matlab
Code:
```clear all;
[x y z Ux Uy Uz]=textread('U_constantPlane.raw','%f%f%f%f%f%f','headerlines',2);
N=32768;
Lz=6.28;

kz=(2*pi*(1:256)/Lz);
kzz=round(kz);
w=reshape(Uz,256,128);
wk=fft(w);
for i = 1:256,
for j=1:128,
wk(i,j)=wk(i,j)/N;
end
end
for i = 1:256,
for j=1:128,
R(i,j,kzz)=wk(i,j)*conj(wk(i,j));
end
end
for j=1:128,
E(kzz)=0.5.*R(j,j,kzz(j))*4*pi.*kzz(j)^2;
end
plot(kzz,E(kzz))```
I have the velocity on a plane and the points are 256*128, so I turned the U into a 256*128 matrix. By doing the whole program, the result is zero! I don't know why, could you please help me?

April 21, 2014, 03:17
#42
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,511
Rep Power: 31
Quote:
 Originally Posted by huangxianbei Hi,Filippo I think the third line should be: cu(q)=c(q)*conj(c(q)) As known, the velocity spectrum tensor is phi[i,j](k)=cu(k),the energy spectrum is integrated to E(k)=1/2*cu(k)*dS(k) where S(k) is the sphere's surface area with radius mag(k) is this means to calculate E(k)=1/2*cu(k)*4*pi*k^2？
yes, the second integer k is just becasue I stored in a matrix for all the points along the z-axis

April 21, 2014, 03:33
#43
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro yes, the second integer k is just becasue I stored in a matrix for all the points along the z-axis
Hi,Filippo,
Thank you for your reply, could you please have a look at my previous reply for the code?
Also, I tried another case which only contains 1 line of velocity in z direction.The result is still incorrect
Code:
```clear all;
[z Ux Uy Uz]=textread('Zleft_U.xy','%f%f%f%f');
N=128;
Lz=6.28;
kz=(1:N)*2*pi/Lz;
k=round(kz);
c=fft(Uz,128);
for q = 1:N,
c(q)=c(q)/N;
end
for q=1:N,
cu(q)=c(q)*conj(c(q));
end
%up to now, the velocity-spectrum is obtained,however, the q~=k,
%how to match them?
for q=1:N,
E(k)=0.5*cu(q)*4*pi*k.^2;
end
plot(k,E(k))```

April 21, 2014, 04:13
#44
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro yes, the second integer k is just becasue I stored in a matrix for all the points along the z-axis
I managed to match them with this:
E=0.5.*cu.*4*pi.*k.^2;
Again,the result is incorrect

April 21, 2014, 04:18
#45
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,511
Rep Power: 31
Quote:
 Originally Posted by huangxianbei Hi,Filippo Sorry for the second reply.I followed your advice, here is my code in matlab Code: ```clear all; [x y z Ux Uy Uz]=textread('U_constantPlane.raw','%f%f%f%f%f%f','headerlines',2); N=32768; Lz=6.28; kz=(2*pi*(1:256)/Lz); kzz=round(kz); w=reshape(Uz,256,128); wk=fft(w); for i = 1:256, for j=1:128, wk(i,j)=wk(i,j)/N; end end for i = 1:256, for j=1:128, R(i,j,kzz)=wk(i,j)*conj(wk(i,j)); end end for j=1:128, E(kzz)=0.5.*R(j,j,kzz(j))*4*pi.*kzz(j)^2; end plot(kzz,E(kzz))``` I have the velocity on a plane and the points are 256*128, so I turned the U into a 256*128 matrix. By doing the whole program, the result is zero! I don't know why, could you please help me?

w=reshape(Uz,256,128);
wk=fft(w)

this way, w is a matrix???

April 21, 2014, 04:31
#46
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro w=reshape(Uz,256,128); wk=fft(w) this way, w is a matrix???
yes,because z values are ranged in a periodic way(256 values one loop), so 32768 values corresponding to 128 loops, representing the whole value on a plane. Is this resonable?
I'm now confused about the wavanumber k, is this write as:
(1:N)*2*pi/Lz
each one is connected to a cu value in order?

April 21, 2014, 04:40
#47
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,511
Rep Power: 31
Quote:
 Originally Posted by huangxianbei yes,because z values are ranged in a periodic way(256 values one loop), so 32768 values corresponding to 128 loops, representing the whole value on a plane. Is this resonable? I'm now confused about the wavanumber k, is this write as: (1:N)*2*pi/Lz each one is connected to a cu value in order?
therefore, fft(w) returns the Fourier transform of each column of the matrix and store in wc.
Anyway, I suggest checking for a simple test: generate a sin(k*x) function (try for some k value) and store in a vector, after doing the fft plot the spectrum and check if the energy density is located at the right frequency

April 21, 2014, 04:56
#48
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro therefore, fft(w) returns the Fourier transform of each column of the matrix and store in wc. Anyway, I suggest checking for a simple test: generate a sin(k*x) function (try for some k value) and store in a vector, after doing the fft plot the spectrum and check if the energy density is located at the right frequency
All right, I'll try it. Another question is still needed to be solved as I have asked:
what's definition formula of k?
2*pi*n/Lz?

Best reguards

April 26, 2014, 09:00
#49
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by ehsankf modes should change from -N/2:N/2.
This should be -N/2:N/2-1 for an odd N.

 April 29, 2014, 15:40 #50 New Member   Amjad Join Date: May 2012 Posts: 21 Rep Power: 6 hello, could i get ur code pls, my Email is amjadju@yahoo.com thanks

April 29, 2014, 21:42
#51
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by Amjad Asad hello, could i get ur code pls, my Email is amjadju@yahoo.com thanks
you can find it in Energy spectrum-wrong spectrum shape

 July 8, 2014, 01:10 #52 Member   SHEIKH NASIRUDDIN Join Date: Aug 2013 Location: DELHI Posts: 49 Rep Power: 5 Dear Harishg, I 'm also facing the same problem and want to plot energy spectra using matlab. Can you send me the code you have used. My e-mail is:nasir.iitd@gmail.com Your help is highly appreciated. nasiruddin.

 October 27, 2014, 20:42 #53 New Member   Purvi Join Date: Oct 2014 Posts: 4 Rep Power: 3 Can I get the code for energy spectrum?

October 27, 2014, 21:08
#54
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by Amjad Asad hello, could i get ur code pls, my Email is amjadju@yahoo.com thanks
Sorry that I just check my thread. Do u still need it?

October 27, 2014, 21:08
#55
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by p3vdim Can I get the code for energy spectrum?
Yes, most of my code is in this thread, if you need, I'll send you a .m file with email

 October 27, 2014, 21:30 #56 New Member   Purvi Join Date: Oct 2014 Posts: 4 Rep Power: 3 Yes, that will be great and if you can send me understanding of it too!!! Thanks Purvi my email address is : p3vdim@yahoo.co.in

October 28, 2014, 00:53
#57
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by p3vdim Yes, that will be great and if you can send me understanding of it too!!! Thanks Purvi my email address is : p3vdim@yahoo.co.in
Hi,
I have sent you the code. The explaination is not so much for the limited time, hoping this can help you.

Xianbei

 October 28, 2014, 05:34 #58 New Member   chandra shekhar pant Join Date: Oct 2010 Posts: 6 Rep Power: 7 Please send it to me also .. chandukec@gmail.com

 October 29, 2014, 15:33 #59 New Member   Join Date: Nov 2013 Posts: 3 Rep Power: 4 Would you mind sending me the code as well? bmeyer907@gmail.com Thx in advance!

 October 29, 2014, 18:55 Energy spectra #60 New Member   Purvi Join Date: Oct 2014 Posts: 4 Rep Power: 3 Hi Guys, I am also having problems with plotting energy spectrum. Could somebody plz send me the code. My email address is : p3vdim@yahoo.co.in Thanks in advance!!

 Tags energy spectrum, fft

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post cfd.newbie Main CFD Forum 6 September 24, 2015 16:13 cfd.newbie CD-adapco 1 June 19, 2008 23:48 cfd.newbie FLUENT 0 June 18, 2008 18:34 Fabian Main CFD Forum 4 October 18, 2005 02:04 Emad Khalifa Main CFD Forum 3 June 30, 2003 16:03

All times are GMT -4. The time now is 04:21.

 Contact Us - CFD Online - Top