[huangxianbei]

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 http://www.cfd-online.com/Forums/mai...lculation.html as Robert said,here is the 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]

Usually,the spectrum is averaged in time. How long should the time period be? As soon as the curve becomes smooth?

[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

[huangxianbei]

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)?

[huangxianbei]

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

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

[huangxianbei]

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?

[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

[huangxianbei]

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

[huangxianbei]

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

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)

[FMDenaro]

have you normalized correctly? why don't you try the spectrum?

[huangxianbei]

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?

[huangxianbei]

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.

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?

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

[huangxianbei]

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

[Amjad Asad]

Hello,
could anyone explain the steps of this code please
thanks

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

Originally Posted by huangxianbei

After numerous trying, I think I get the right(seems) spectrum.

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?

[huangxianbei]

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

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

Originally Posted by huangxianbei

velocity is made dimensionless before FFT, so this is the nondimensional spectrum

[huangxianbei]

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]