| huangxianbei |
April 20, 2014 22:58 |
Quote:
Originally Posted by FMDenaro
(Post 406843)
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? |
