CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Energy Spectrum !!!

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

Like Tree14Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   April 20, 2014, 22:58
Default
  #41
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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?
shuige likes this.
huangxianbei is offline   Reply With Quote

Old   April 21, 2014, 03:17
Default
  #42
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,511
Rep Power: 31
FMDenaro will become famous soon enoughFMDenaro will become famous soon enough
Quote:
Originally Posted by huangxianbei View Post
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
FMDenaro is offline   Reply With Quote

Old   April 21, 2014, 03:33
Default
  #43
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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))
huangxianbei is offline   Reply With Quote

Old   April 21, 2014, 04:13
Default
  #44
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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
huangxianbei is offline   Reply With Quote

Old   April 21, 2014, 04:18
Default
  #45
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,511
Rep Power: 31
FMDenaro will become famous soon enoughFMDenaro will become famous soon enough
Quote:
Originally Posted by huangxianbei View Post
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???
FMDenaro is offline   Reply With Quote

Old   April 21, 2014, 04:31
Default
  #46
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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?
huangxianbei is offline   Reply With Quote

Old   April 21, 2014, 04:40
Default
  #47
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,511
Rep Power: 31
FMDenaro will become famous soon enoughFMDenaro will become famous soon enough
Quote:
Originally Posted by huangxianbei View Post
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
FMDenaro is offline   Reply With Quote

Old   April 21, 2014, 04:56
Default
  #48
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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
huangxianbei is offline   Reply With Quote

Old   April 26, 2014, 09:00
Default
  #49
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Quote:
Originally Posted by ehsankf View Post
modes should change from -N/2:N/2.
This should be -N/2:N/2-1 for an odd N.
huangxianbei is offline   Reply With Quote

Old   April 29, 2014, 15:40
Default
  #50
New Member
 
Amjad
Join Date: May 2012
Posts: 21
Rep Power: 6
Amjad Asad is on a distinguished road
hello,
could i get ur code pls, my Email is amjadju@yahoo.com
thanks
Amjad Asad is offline   Reply With Quote

Old   April 29, 2014, 21:42
Default
  #51
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Quote:
Originally Posted by Amjad Asad View Post
hello,
could i get ur code pls, my Email is amjadju@yahoo.com
thanks
you can find it in Energy spectrum-wrong spectrum shape
huangxianbei is offline   Reply With Quote

Old   July 8, 2014, 01:10
Default
  #52
Member
 
SHEIKH NASIRUDDIN
Join Date: Aug 2013
Location: DELHI
Posts: 49
Rep Power: 5
nasir.iitd is on a distinguished road
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.
nasir.iitd is offline   Reply With Quote

Old   October 27, 2014, 20:42
Default
  #53
New Member
 
Purvi
Join Date: Oct 2014
Posts: 4
Rep Power: 3
p3vdim is on a distinguished road
Can I get the code for energy spectrum?
p3vdim is offline   Reply With Quote

Old   October 27, 2014, 21:08
Default
  #54
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Quote:
Originally Posted by Amjad Asad View Post
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?
huangxianbei is offline   Reply With Quote

Old   October 27, 2014, 21:08
Default
  #55
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Quote:
Originally Posted by p3vdim View Post
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
huangxianbei is offline   Reply With Quote

Old   October 27, 2014, 21:30
Default
  #56
New Member
 
Purvi
Join Date: Oct 2014
Posts: 4
Rep Power: 3
p3vdim is on a distinguished road
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
p3vdim is offline   Reply With Quote

Old   October 28, 2014, 00:53
Default
  #57
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Quote:
Originally Posted by p3vdim View Post
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
huangxianbei is offline   Reply With Quote

Old   October 28, 2014, 05:34
Default
  #58
New Member
 
chandra shekhar pant
Join Date: Oct 2010
Posts: 6
Rep Power: 7
chandra shekhar pant is on a distinguished road
Please send it to me also ..
chandukec@gmail.com
chandra shekhar pant is offline   Reply With Quote

Old   October 29, 2014, 15:33
Default
  #59
New Member
 
Join Date: Nov 2013
Posts: 3
Rep Power: 4
Leon is on a distinguished road
Would you mind sending me the code as well? bmeyer907@gmail.com
Thx in advance!
Leon is offline   Reply With Quote

Old   October 29, 2014, 18:55
Default Energy spectra
  #60
New Member
 
Purvi
Join Date: Oct 2014
Posts: 4
Rep Power: 3
p3vdim is on a distinguished road
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!!
p3vdim is offline   Reply With Quote

Reply

Tags
energy spectrum, fft

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


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


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