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

Energy spectrum-wrong spectrum shape

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

Like Tree2Likes
  • 1 Post By FMDenaro
  • 1 Post By huangxianbei

Reply
 
LinkBack Thread Tools Display Modes
Old   April 21, 2014, 07:23
Default Energy spectrum-wrong spectrum shape
  #1
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 211
Rep Power: 4
huangxianbei is on a distinguished road
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?
huangxianbei is offline   Reply With Quote

Old   April 21, 2014, 08:47
Default
  #2
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 211
Rep Power: 4
huangxianbei is on a distinguished road
Usually,the spectrum is averaged in time. How long should the time period be? As soon as the curve becomes smooth?
huangxianbei is offline   Reply With Quote

Old   April 21, 2014, 13:25
Default
  #3
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,247
Rep Power: 17
FMDenaro will become famous soon enough
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
FMDenaro is offline   Reply With Quote

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

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

Old   April 22, 2014, 03:35
Default
  #6
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,247
Rep Power: 17
FMDenaro will become famous soon enough
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".
FMDenaro is offline   Reply With Quote

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

Old   April 22, 2014, 04:31
Default
  #8
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,247
Rep Power: 17
FMDenaro will become famous soon enough
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.
FMDenaro is offline   Reply With Quote

Old   April 22, 2014, 06:07
Default
  #9
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 211
Rep Power: 4
huangxianbei is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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
huangxianbei is offline   Reply With Quote

Old   April 24, 2014, 06:43
Unhappy
  #10
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 211
Rep Power: 4
huangxianbei is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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');
[xorder]=textread('x.txt','%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)
huangxianbei is offline   Reply With Quote

Old   April 24, 2014, 09:11
Default
  #11
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,247
Rep Power: 17
FMDenaro will become famous soon enough
have you normalized correctly? why don't you try the spectrum?
FMDenaro is offline   Reply With Quote

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

Old   April 26, 2014, 09:08
Default
  #13
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 211
Rep Power: 4
huangxianbei is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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');
%[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?
mgg likes this.
huangxianbei is offline   Reply With Quote

Old   April 26, 2014, 13:15
Default
  #14
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,247
Rep Power: 17
FMDenaro will become famous soon enough
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...
FMDenaro is offline   Reply With Quote

Old   April 27, 2014, 05:30
Default
  #15
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 211
Rep Power: 4
huangxianbei is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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
huangxianbei is offline   Reply With Quote

Old   April 30, 2014, 00:30
Default
  #16
New Member
 
amjad
Join Date: May 2012
Posts: 12
Rep Power: 4
Amjad Asad is on a distinguished road
Hello,
could anyone explain the steps of this code please
thanks
Amjad Asad is offline   Reply With Quote

Old   June 8, 2014, 20:37
Default
  #17
New Member
 
ehk
Join Date: Sep 2012
Posts: 21
Rep Power: 4
ehsankf is on a distinguished road
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 View Post
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?
ehsankf is offline   Reply With Quote

Old   June 8, 2014, 23:29
Default
  #18
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 211
Rep Power: 4
huangxianbei is on a distinguished road
Quote:
Originally Posted by ehsankf View Post
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
huangxianbei is offline   Reply With Quote

Old   June 8, 2014, 23:54
Default
  #19
New Member
 
ehk
Join Date: Sep 2012
Posts: 21
Rep Power: 4
ehsankf is on a distinguished road
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 View Post
velocity is made dimensionless before FFT, so this is the nondimensional spectrum
ehsankf is offline   Reply With Quote

Old   June 9, 2014, 02:29
Default
  #20
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 211
Rep Power: 4
huangxianbei is on a distinguished road
Quote:
Originally Posted by ehsankf View Post
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]
huangxianbei is offline   Reply With Quote

Reply

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
ATTENTION! Reliability problems in CFX 5.7 Joseph CFX 14 April 20, 2010 15:45
energy in sonicFoam joern OpenFOAM Running, Solving & CFD 0 June 16, 2009 15:53
udf error srihari FLUENT 0 February 9, 2009 09:00
Energy Spectrum in LES Simulation Anindya FLUENT 0 March 17, 2006 14:20
Why FVM for high-Re flows? Zhong Lei Main CFD Forum 23 May 14, 1999 13:22


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