[Hi123]

Originally Posted by juliom

Dear Community;
I have read many post about the topics but I still have a question. I am working with LES and I want to verify my LES code. Basically I am doing validation using the channel flow data base from JHUTB. So, I have a 3D flow field: u(x,y,z,t),v(x,y,z,t),w(x,y,z,t). I am interesting in compute the power spectrum at each plane. I have read about defining a point a collect the temporal velocity at a single point; which make sense, but I am more interested over the space rather than time..!!!
I have read about compute the FFT from the velocity, but from Wilcox book, I read that the Spectrum Energy is the FFT of the velocity correlation tensor.
I was planning to use FORTRAN to write the script, but I have read that the majority is using Matlab, so it seems that I will lean towards Matlab. However, I want to have the picture clear before starting!!!

Thanks community!

Hi Juliom,
I am also currently looking for writing a FORTRAN program to do the spectral analysis. AS you have mentioned, I have LES simulated data of the flow field : u(x,y,z,t) , v(x,y,z,t) , w(x,y,z,t). My ultimate aim is to apply a spatial filter to view only the large eddies. I was thinking of first doing a spectral analysis over the complete 3D domain which is a cuboid to see how the distribution of energy is.
Kindly advise me on the same.

It would also be helpful if anyone could suggest a good literature related to this.

Thanks in Advance.

[FMDenaro]

If you have a LES solution, you have data that are already filtered by the LES filter.

[Hi123]

Hi Denaro,
Thank you for your reply.
It is my understanding that "The field is filtered. The lower scale turbulence are modeled Dynamic Smagorisky and the NS equation is numerically solved. Then the resulting velocity field data has both large and small scale turbulence in it. "
Please correct me if I am wrong.

[FMDenaro]

Originally Posted by Hi123

Hi Denaro,
Thank you for your reply.
It is my understanding that "The field is filtered. The lower scale turbulence are modeled Dynamic Smagorisky and the NS equation is numerically solved. Then the resulting velocity field data has both large and small scale turbulence in it. "
Please correct me if I am wrong.

No, the resulting velocity field contains only the large eddy flow scales wherein the modelled action of the SGS model acts on. The energy spectra will not contain the content of the small scale but they finish at the wavenumber of the LES filter.

[Hi123]

Okay. I am getting it now. Thank you.
But, What should I do if I need to separate the flow fields within the large scale structures themselves ? E.g to set a lower cut off frequency than that of used in LES ?
Should I go on with a spatial filter once again ?

[FMDenaro]

Originally Posted by Hi123

Okay. I am getting it now. Thank you.
But, What should I do if I need to separate the flow fields within the large scale structures themselves ? E.g to set a lower cut off frequency than that of used in LES ?
Should I go on with a spatial filter once again ?

Yes, you can post-filter your LES solution. But if you are using the dynamic model, that is already part of your code (the test-filtering routine) so that you can directly use that field.

[kk415]

Dear All,

I would like to recall this topic with the original question: I have an Nx, Ny, Nz grid (say all are equally spaced N=32). So I have 3 component velocity vectors u, v, w, each of size 32768. I am using the FFT command of Matlab to get the FFT of the velocity component. I am not sure what is my kx, ky, kz ? The FFT gives u(k) vector of same size as 32768. Can my kx is 0 to 32768/2 ?? I think not!! Can anyone please help me on this. I would really appreciate it. Thanks!!

regards
krishna

[LuckyTran]

fft in matlab is 1-D. If your grid is 32x32x32 then you need to do fft per dimension (and then do spatial statistics). Do the n-point fft along the non-trivial dimension and average over the homogeneous directions. That means you should have 32x32=64 vectors of length 32. You probably will take the maximum number of points as a naive estimate in which case fft returns a 32 element vector (64 times).

And then I pray that given this you know how to rescale the output of the fft to get wavenumbers. The matlab documentation / help file has several good examples.

If you're at the level that you're doing 2D/3D discrete fourier transforms, then all I have to tell you is that fft is 1D and it is the wrong function in matlab to use. Still, you won't have a 32768 size vector as an input.

Also, test your fft script against a sine function of known periodicity/wavenumber!

[FMDenaro]

Originally Posted by kk415

Dear All,

I would like to recall this topic with the original question: I have an Nx, Ny, Nz grid (say all are equally spaced N=32). So I have 3 component velocity vectors u, v, w, each of size 32768. I am using the FFT command of Matlab to get the FFT of the velocity component. I am not sure what is my kx, ky, kz ? The FFT gives u(k) vector of same size as 32768. Can my kx is 0 to 32768/2 ?? I think not!! Can anyone please help me on this. I would really appreciate it. Thanks!!

regards
krishna

There is no doubt about your question:


You have a 1D periodicity lenght L subdivided by N=32 steps h=L/N. Each 1D FFT will give you N/2 Fourier components. The spatial frequency is k=n*2*pi/L=n*(2/N)*(pi/h). Your Nyquist frequency is k_max=pi/h reached at n=N/2.

[kk415]

Dear Lucky and Prof. Filippo,

Thanks for your reply. I did some more research on the topic and found out that I had to do reshaping of my velocity vector.

So my u vector of size 32768 is now a matrix of 32*32*32.

I then use the FFT to get u(k) a matrix of 32*32*32.

e11=abs(u(k))/(32^3).

E=0.5*(e11^2+e22^2+e33^2).

Not sure, though, if this is correct!!

A summation step is missing. A simple sum function will give a matrix of 1*32*32. Is that what you meant Lucky?

I originally got my vector field for doing synthetic homogenous isotropic turbulence using the Tony Saad algorithm. And trying to validate my obtained velocity field. So Length of the domain is 2*pi. kmin=2*pi/L=1 and kmax=2*pi/delta=32.

so my wave number vector is kx=kmin+(kmax-kmin)*m, where m is from 1 to ??(N/2 ??)

kmag = sqrt(kx^2+ky^2+kz^2).
And I have to multiply E to 4*pi*kmag*kmag also.

So I am stuck in the summation step and maximum value of modes of FFT. Kindly help me or suggest if I am wrong somewhere.

Regards
krishna

[FMDenaro]

Originally Posted by kk415

Dear Lucky and Prof. Filippo,

Thanks for your reply. I did some more research on the topic and found out that I had to do reshaping of my velocity vector.

So my u vector of size 32768 is now a matrix of 32*32*32.

I then use the FFT to get u(k) a matrix of 32*32*32.

e11=abs(u(k))/(32^3).

E=0.5*(e11^2+e22^2+e33^2).

Not sure, though, if this is correct!!

A summation step is missing. A simple sum function will give a matrix of 1*32*32. Is that what you meant Lucky?

I originally got my vector field for doing synthetic homogenous isotropic turbulence using the Tony Saad algorithm. And trying to validate my obtained velocity field. So Length of the domain is 2*pi. kmin=2*pi/L=1 and kmax=2*pi/delta=32.

so my wave number vector is kx=kmin+(kmax-kmin)*m, where m is from 1 to ??(N/2 ??)

kmag = sqrt(kx^2+ky^2+kz^2).
And I have to multiply E to 4*pi*kmag*kmag also.

So I am stuck in the summation step and maximum value of modes of FFT. Kindly help me or suggest if I am wrong somewhere.

Regards
krishna

k_max is pi/h not 2*pi/h

[kk415]

Dear Community,

I finally found good links for creating the synthetic isotropic turbulence and calculating spectrum, respectively. It's an OpenFoam program code.

Creating Turbulence Field:
https://www.openfoam.com/documentati...8C_source.html

Creating Energy Spectrum:
https://www.openfoam.com/documentati...8C_source.html

What I am trying to do now is create a similar code in Matlab. And instead of creating the CBC spectrum, I am trying to generate the Batchelor-Towend spectrum. I am attaching my results for two grids with each side division of 32 and 64. Results are poor . But wondering if it is the best we get by above-said algorithms!!

The only difference I made is in the creation of turbulence field code they are using the second trigonometry function inside the first trigonometry function. However, I kept the trigonometry functions separately.

OpenFoam code:

    return vector
    (
        Foam::sin(thetam*Foam::cos(phim)),
        Foam::sin(thetam*Foam::sin(phim)),
        Foam::cos(thetam)
    );

Matlab code:

Khatxm=sin(thetham).*cos(phim);
Khatym=sin(thetham).*sin(phim);
Khatzm=cos(thetham);