[CFDrookie]

Hi all 1st post so be gentle!

I have a 3D (turbulent) velocity field. For the sake of example lets say this is contained in 3 n^3 arrays, ux, uy, uz.

I want the 3d energy spectrum, my method at the moment is to calculate fft of these arrays (everything below is in matlab syntax) i.e.

>fux=fftn(ux)/(n^3);
>fuy=fftn(uy)/(n^3);
>fuz=fftn(uz)/(n^3);

I then get real and imaginary parts:

>energyr=real(fux).^2+real(fuy).^2+real(fuz).^2;
>energyi=imag(fux).^2+imag(fuy).^2+imag(fuz).^2;

Define the maximum wavenumber

>midpt=n/2+1;

Now I compute energy spectra in shells,

spect(1:1.5*n)=0.;
for i=1:n
for j=1:n
for k=1:n
ii=i; jj=j; kk=k;
if ii>midpt ; ii=n-ii+1; ; end ;
if jj>midpt ; jj=n-jj+1; ; end ;
if kk>midpt ; kk=n-kk+1; ; end ;
r=round(sqrt(ii^2+jj^2+kk^2));
spect(r)=spect(r)+energyr(i,j,k)+energyi(i,j,k);
end
end
end

My question is, whilst this works well, I do not get the energy at the box scale (k=1). Where am I going wrong?

Many thanks in advance

[CFDrookie]

Anyone got any ideas?

[saeedi]

Hey,

You know you should first calculate the Autocorrelation or autocovariance, and then take the fft of that one not the vlocity itslf.

[CFDstudent07]

Hi
Did you figure it out? I may need the same code too

Originally Posted by CFDrookie

Hi all 1st post so be gentle!

I have a 3D (turbulent) velocity field. For the sake of example lets say this is contained in 3 n^3 arrays, ux, uy, uz.

I want the 3d energy spectrum, my method at the moment is to calculate fft of these arrays (everything below is in matlab syntax) i.e.

>fux=fftn(ux)/(n^3);
>fuy=fftn(uy)/(n^3);
>fuz=fftn(uz)/(n^3);

I then get real and imaginary parts:

>energyr=real(fux).^2+real(fuy).^2+real(fuz).^2;
>energyi=imag(fux).^2+imag(fuy).^2+imag(fuz).^2;

Define the maximum wavenumber

>midpt=n/2+1;

Now I compute energy spectra in shells,

spect(1:1.5*n)=0.;
for i=1:n
for j=1:n
for k=1:n
ii=i; jj=j; kk=k;
if ii>midpt ; ii=n-ii+1; ; end ;
if jj>midpt ; jj=n-jj+1; ; end ;
if kk>midpt ; kk=n-kk+1; ; end ;
r=round(sqrt(ii^2+jj^2+kk^2));
spect(r)=spect(r)+energyr(i,j,k)+energyi(i,j,k);
end
end
end

My question is, whilst this works well, I do not get the energy at the box scale (k=1). Where am I going wrong?

Many thanks in advance