# Energy spectrum-wrong spectrum shape

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

 April 21, 2014, 07:23 Energy spectrum-wrong spectrum shape #1 Senior Member   Huang Xianbei Join Date: Sep 2013 Location: CAU,China Posts: 277 Rep Power: 6 Hi,all Recently I got stuck with energy spectrum. I sampled the data across the spanwise direction(z) in channel flow, Uz, then I want to plot the energy spectrum, I adopt the method here Energy spectrum calculation as Robert said,here is the code Code: clear all; [z Ux Uy Uz]=textread('Zleft_U.xy','%f%f%f%f'); N=128; Lz=6.28; kz=(2*pi*(1:128)/Lz)'; kzz=round(kz); Uzk=fft(Uz,N)/N; %Fs=2*pi; %f = kz/2*linspace(0,1,N/2+1); E=0.5*4*((2*abs(Uzk(1:N/2+1))).^2).*pi.*(kzz(1:N/2+1)).^2 plot(kzz(1:N/2+1),E) While the result is quite far from usual,what's wrong with this method?

 April 21, 2014, 08:47 #2 Senior Member   Huang Xianbei Join Date: Sep 2013 Location: CAU,China Posts: 277 Rep Power: 6 Usually,the spectrum is averaged in time. How long should the time period be? As soon as the curve becomes smooth?

 April 21, 2014, 13:25 #3 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 2,498 Rep Power: 31 first, you should plot in a double logarithmic scale, then you must sample the velocity field in a sufficient number of file, the sampling time depending on the turnover characteristic time. At each time compute the Fourier coefficients and then do the time average

April 21, 2014, 20:35
#4
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro first, you should plot in a double logarithmic scale, then you must sample the velocity field in a sufficient number of file, the sampling time depending on the turnover characteristic time. At each time compute the Fourier coefficients and then do the time average
Thank you for your suggestion, yes, I forgot to use the loglog scale
Here, I have another question: I have read some papers about this and found that the usual method is to do the spectrum average in time and plane parallel to the wall. I'm quite confused about the average in plane(say the plane y+=10 for example). If I want to plot the spanwise spectrum(z direction in channel flow), is that means I should do average of each z position, with all the x positions(z fixed, and average the values at different x)?

April 21, 2014, 22:51
#5
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro first, you should plot in a double logarithmic scale, then you must sample the velocity field in a sufficient number of file, the sampling time depending on the turnover characteristic time. At each time compute the Fourier coefficients and then do the time average
I find a big mistake in my process. The velocity used in the energy spectrum should be the fluctuating velocity other than the temporal velocity, I think this should lead to the wrong shape

 April 22, 2014, 03:35 #6 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 2,498 Rep Power: 31 the averaging along a plane is used only for homogeneous flows (for example channel with periodic conditions) in order to make the statistics more "robust".

April 22, 2014, 04:12
#7
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro the averaging along a plane is used only for homogeneous flows (for example channel with periodic conditions) in order to make the statistics more "robust".
so If I want to plot the spanwise spectrum(z direction in channel flow), is that means I should do average of each z position, with all the x positions(z fixed, and average the values at different x)? Just as I have asked?

 April 22, 2014, 04:31 #8 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 2,498 Rep Power: 31 well, assume the channel has: x-streamwise (periodic), y-normal to walls, z-spanwise (periodic) directions. - Fix some y value (or y+ value) - Chose if you want the spectrum along Kx or Kz. - For example, along the streamwise direction the frequencies are Kx=nx*2*pi/Lx (nx=1,...namx), then for each z value of the y=constant plane perform the FFT and store the Fourier coefficients. - Perform the average of the FOurier coefficients over z, this way you have the coefficients depending only on Kx. - Repeat the procedure for all the temporal fields mgg likes this.

April 22, 2014, 06:07
#9
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro well, assume the channel has: x-streamwise (periodic), y-normal to walls, z-spanwise (periodic) directions. - Fix some y value (or y+ value) - Chose if you want the spectrum along Kx or Kz. - For example, along the streamwise direction the frequencies are Kx=nx*2*pi/Lx (nx=1,...namx), then for each z value of the y=constant plane perform the FFT and store the Fourier coefficients. - Perform the average of the FOurier coefficients over z, this way you have the coefficients depending only on Kx. - Repeat the procedure for all the temporal fields
Thank you very much for your details. I'll try it

April 24, 2014, 06:43
#10
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro well, assume the channel has: x-streamwise (periodic), y-normal to walls, z-spanwise (periodic) directions. - Fix some y value (or y+ value) - Chose if you want the spectrum along Kx or Kz. - For example, along the streamwise direction the frequencies are Kx=nx*2*pi/Lx (nx=1,...namx), then for each z value of the y=constant plane perform the FFT and store the Fourier coefficients. - Perform the average of the FOurier coefficients over z, this way you have the coefficients depending only on Kx. - Repeat the procedure for all the temporal fields
Hi,Filippo:
Sorry for my late. I tried similar as you say to calculate the two-point correlation(just to simply see if it's right), however, the correlation exceeds 1
Code:
clear all;
[x y z Ux Uy Uz]=textread('U_constantPlane.txt','%f%f%f%f%f%f');
[xx yy zz Uxm Uym Uzm]=textread('UMean_constantPlane.txt','%f%f%f%f%f%f');
Uxmatrix=reshape(Ux,128,256);
Uxmmatrix=reshape(Uxm,128,256);
ux=Uxmatrix-Uxmmatrix;
uxplane=mean(ux);
for i=1:256,
corrx(i)=uxplane(1)*uxplane(i)/(uxplane(1)*uxplane(1));
end
plot(xorder,corrx)

 April 24, 2014, 09:11 #11 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 2,498 Rep Power: 31 have you normalized correctly? why don't you try the spectrum?

April 24, 2014, 12:21
#12
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro have you normalized correctly? why don't you try the spectrum?
I think It's due to the data structure. My plane is 128*128 grid. While the output contains twice the points value(32768 values). I'll try to figure out a better way to output. I have just tried the spectrum(one temporal data), the result seems incorrect too.
Q1:Is it need to do like this:
for i=1:N,
fftcoefficient/N
as someone suggested ?

Q2:Also, I have another question. What's the mean velocity's definition? If we do the average on the plane, is that means the average should be the fff(U)dxdzdt?
In the present reply , I use the temporal velocity U and time-averaged UMean, is this right or just as I said before, should do further average?

April 26, 2014, 09:08
#13
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro have you normalized correctly? why don't you try the spectrum?
After numerous trying, I think I get the right(seems) spectrum.
Code:
clear all;
[x y z Ux Uy Uz]=textread('Uzyplus35.txt','%f%f%f%f%f%f');
[xx yy zz Uxm Uym Uzm]=textread('UzMeanyplus35.txt','%f%f%f%f%f%f');
N=128;
Lz=6.28;%length in the z direction
kz=2*pi*(-N/2:1:N/2-1)/Lz;%frequencies, or wavenumber
Uzre=reshape(Uz,N,N);%change the column of Ux into a matrix,each column corresponding to one z
Uzmre=reshape(Uzm,N,N);
Uzmreave=mean(Uzmre(:));%field average
Uzmfield=zeros(N,N);
for i=1:N,
for j=1:N,
Uzmfield(i,j)=Uzmreave;
end
end
w=Uzre-Uzmfield;%get the velocity fluctuation=instannous velocity - timeaveraged velocity
wk=fft(w);%each column corresponds to a z value
wkk=fftshift(wk);
E=zeros(N,N);
for i=1:N,
for j=1:N,
E(i,j)=(wkk(i,j)*conj(wkk(i,j)));%equal to abs(wkmean(i))^2
end
end
Emean=mean(E,2);
loglog(kz,Emean);

The latter one is from a paper with DNS data in channel flow,Re=2666, while my simulation has Re=2900, the deviation is somehow large, do you think it's reasonable or my method is wrong?

 April 26, 2014, 13:15 #14 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 2,498 Rep Power: 31 the shape seems correct, however I see one order of magnitude in the lowest wavenumber... check it carefully, I suspect you are doing the streamwise component of the spectrum while the DNS reports the spanwise...

April 27, 2014, 05:30
#15
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by FMDenaro the shape seems correct, however I see one order of magnitude in the lowest wavenumber... check it carefully, I suspect you are doing the streamwise component of the spectrum while the DNS reports the spanwise...
I have clarified the method by comparing with the database of Moin and Kim with Ret=180. In fact , the velocity in the calculation should be non-dimensional divided by ut. So in the fft, non-dimensional velocity is used and also, the amplitude obtained from fft should divided by N.

Now, the result is correct

 April 30, 2014, 00:30 #16 New Member   Amjad Join Date: May 2012 Posts: 21 Rep Power: 6 Hello, could anyone explain the steps of this code please thanks

June 8, 2014, 20:37
#17
New Member

ehk
Join Date: Sep 2012
Posts: 21
Rep Power: 6
What is the dimension of Energy Spectra calculated in this code? Is it (m^4/s^3)?
As FFT(m/s)=m^2/s.
Quote:
 Originally Posted by huangxianbei After numerous trying, I think I get the right(seems) spectrum. Code: clear all; [x y z Ux Uy Uz]=textread('Uzyplus35.txt','%f%f%f%f%f%f'); [xx yy zz Uxm Uym Uzm]=textread('UzMeanyplus35.txt','%f%f%f%f%f%f'); %[xorder]=textread('z.txt','%f');,just ignore this line N=128; Lz=6.28;%length in the z direction kz=2*pi*(-N/2:1:N/2-1)/Lz;%frequencies, or wavenumber Uzre=reshape(Uz,N,N);%change the column of Ux into a matrix,each column corresponding to one z Uzmre=reshape(Uzm,N,N); Uzmreave=mean(Uzmre(:));%field average Uzmfield=zeros(N,N); for i=1:N, for j=1:N, Uzmfield(i,j)=Uzmreave; end end w=Uzre-Uzmfield;%get the velocity fluctuation=instannous velocity - timeaveraged velocity wk=fft(w);%each column corresponds to a z value wkk=fftshift(wk); E=zeros(N,N); for i=1:N, for j=1:N, E(i,j)=(wkk(i,j)*conj(wkk(i,j)));%equal to abs(wkmean(i))^2 end end Emean=mean(E,2); loglog(kz,Emean); The latter one is from a paper with DNS data in channel flow,Re=2666, while my simulation has Re=2900, the deviation is somehow large, do you think it's reasonable or my method is wrong?

June 8, 2014, 23:29
#18
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by ehsankf What is the dimension of Energy Spectra calculated in this code? Is it (m^4/s^3)? As FFT(m/s)=m^2/s.
velocity is made dimensionless before FFT, so this is the nondimensional spectrum

June 8, 2014, 23:54
#19
New Member

ehk
Join Date: Sep 2012
Posts: 21
Rep Power: 6
As I mentioned FFT[m/s]=m^2/s. If you apply FFT for nondimensional variable, the result would have the dimension of [m], as there is an integration over space to derive FFT coefficients. If you also check the literature for plots of Energy Spectra, it is seen the dimension is like [m^3/s^2]=FFT(Kinetic Energy)=FFT([m^2/s^2]). The way to non-dimensionalize Energy Spectra is dividing by (\eta u_*) where \eta is the kolmogrov length scale.

Quote:
 Originally Posted by huangxianbei velocity is made dimensionless before FFT, so this is the nondimensional spectrum

June 9, 2014, 02:29
#20
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
Quote:
 Originally Posted by ehsankf As I mentioned FFT[m/s]=m^2/s. If you apply FFT for nondimensional variable, the result would have the dimension of [m], as there is an integration over space to derive FFT coefficients. If you also check the literature for plots of Energy Spectra, it is seen the dimension is like [m^3/s^2]=FFT(Kinetic Energy)=FFT([m^2/s^2]). The way to non-dimensionalize Energy Spectra is dividing by (\eta u_*) where \eta is the kolmogrov length scale.
While, u(x,t)=sum(e^ikx*u(k,t)) ,here k[/m],x[m], so the u(k,t)[m/s] which is the same as u(x,t). So if energy spectrum is calculated, that's E(k)=u(k,t)*u-(k,t)[m^2/s^2]. Considering the FFT, as you say, [m] will be introduced into the result, try think like this, we can easily consider the length parameter to be dimensionless, for example, L[m]/1[m], this makes it dimensionless, thus ,the FFT won't bring in [m]

 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 Joseph CFX 14 April 20, 2010 15:45 joern OpenFOAM Running, Solving & CFD 0 June 16, 2009 15:53 srihari FLUENT 0 February 9, 2009 10:00 Anindya FLUENT 0 March 17, 2006 15:20 Zhong Lei Main CFD Forum 23 May 14, 1999 13:22

All times are GMT -4. The time now is 17:29.