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

How to plot time dependent energy spectrum

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

Like Tree3Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 8, 2021, 14:16
Smile How to plot time dependent energy spectrum
  #1
New Member
 
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4
Falu is on a distinguished road
Hey CFD-Community,

I've done a LES free stream simulation and now i want to plot my data. I have the velocity data of 24 points of my simulation environment for the x, y and z components and i want to plot my data in Matlab. The last plot i need is an energy spectrum plot.

I'm imagining 24 energy spectrum plots for my 24 observer points. But at first i just wanted to a basic energy spectrum plot for one velocity component to get an idea how to do it. I used a method i found here in the forum:

u=fft(VelocityU(:,1));
p=(u.*conj(u))/N^2;

But im not getting any viable data, my energy spectrum has a peak at point 0 and thats it. I have uploaded my result and the origin velocity data. Can you tell me what im doing wrong?

Thank you very much in advance, and have a good day.
Attached Images
File Type: jpg EnergySpectrum.jpg (12.5 KB, 29 views)
File Type: jpg VelocityU6_Zeitbereich_Normal.jpg (33.4 KB, 26 views)
Falu is offline   Reply With Quote

Old   January 8, 2021, 15:26
Default
  #2
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,252
Rep Power: 67
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
First of all, use the log scale but it is also mandatory that you plot the range of frequencies only up to Nyquist cut-off (pi/dt). In a time-frequency analsysis you should adopt a time-windowing technique, why don't you start first by doing a spatial analysis?
FMDenaro is offline   Reply With Quote

Old   January 9, 2021, 06:01
Default
  #3
New Member
 
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4
Falu is on a distinguished road
Thanks FMDenaro for your reply.

I still have a few questions.
I have implemented the Han-Window and I'm using the log scale now. I'm curious, why should i start with a spatial analysis when I'm not interested in it?(1) Anyway, i did one for my first datarow for my velocity u.
Sidenote: My simulation time is 1 Second with dt=0.0001. Because i want to view a stationary process, i cutout the first 0.2 Seconds, so my data starts with second 0.2001. To fullfil the 2^x law i added 192 zeros to my data (to get to 8192 entries) to use the FFT algorithm. But i think this is a mistake, isnt it? In my time-energy spectrum plot i have small fluctuations for high frequencies, this is due to my addition, isnt it?(2)
Additional, I have a more general question to this transformation. Am I really getting the frequencies, or am i getting the wavenumber here?(3)
My last question for now applies to the spational energy analysis. I have 24 Monitoring points with equal distance to each other. Then my nyquist frequency is just 12?(4)
Attached Images
File Type: jpg EnergySpectrum_Space_VelU_Time1.jpg (15.7 KB, 21 views)
File Type: jpg EnergySpectrumVelU6.jpg (20.4 KB, 30 views)
File Type: jpg VelocityU6_Zeitbereich_Hann-Window.jpg (29.1 KB, 15 views)
File Type: jpg VelocityU6_Zeitbereich_Normal.jpg (33.4 KB, 16 views)
Falu is offline   Reply With Quote

Old   January 9, 2021, 06:20
Default
  #4
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,252
Rep Power: 67
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Some answers were already addressed in similar posts.
However, the window technique should allow you to get some samples of "periodic" temporal signal so that you can do the spectral analysis of each one and then performing a statistical averaging of the computed spectra.


The FFT must be done only for the resolved frequencies, until the Nyquist frequency pi/dt. Use your computational dt to compute.



For a statistically steady flow you have to wait until the numerical transient is finished, only after you collect the time-samples.


The relation between frequency (f) and wavenumber (n) is f=n*2*pi/T. Therefore you need to evaluate the period T.
Falu likes this.
FMDenaro is offline   Reply With Quote

Old   January 9, 2021, 09:03
Default
  #5
New Member
 
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4
Falu is on a distinguished road
Thanks again for the respond. I've read a lot of posts in this forum, but sometimes i wasnt able to understand it properly, so i thought i post my specific problem to understand the functinality of the spectrum analysis on the basis of an example i have spend a good amount of time with.

Quote:
Originally Posted by FMDenaro View Post
Some answers were already addressed in similar posts.
However, the window technique should allow you to get some samples of "periodic" temporal signal so that you can do the spectral analysis of each one and then performing a statistical averaging of the computed spectra.
Okay, so my mistake isnt the addition of the zeros in my velocity array, but the neglection of the averaging process. I generally know the functionality of a window function, but im confused of the statistical averaging process. In my first impression i thought you mean a time averaging process by this. So I computed a variable, which is defined as the sum of my fft-transformed velocityvector multiplicated by dt and finally divided with T. T in my case is equal to 0.8192 (N(=8219)*dt(0.0001). But I think this is not what you meant?


Quote:
Originally Posted by FMDenaro View Post
The FFT must be done only for the resolved frequencies, until the Nyquist frequency pi/dt. Use your computational dt to compute.
I dont get your definition of the nyquist frequency. The nyquist frequency is f_sample/2, and f_sample i can define as 1/Delta_t to sample every datapoint. Am I wrong here? So im getting a nyquist frency of 5.000Hz(Delta_t=0.0001 -> f_sample=10.000Hz). My confusion with the "spatial nyquist frequency" stems from the fact, that if im considering the space and I "hold the time", my data is reduced to the 24 monitor points in my simulation, in contrast to the 10.000 values if I "hold the space" and analyse the time.

Quote:
Originally Posted by FMDenaro View Post
The relation between frequency (f) and wavenumber (n) is f=n*2*pi/T. Therefore you need to evaluate the period T.
Yeah i found thas definition in this forum, too. So my x-axes is wrong, i now have the wavenumber and i need to transform it to the frequencyspectrum.

Thanks a lot for your help so far Im making progress.
Falu is offline   Reply With Quote

Old   January 9, 2021, 09:25
Default
  #6
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,252
Rep Power: 67
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
You don't have to padd the velocity vector with zeros but you have to use the correct dimension. The nyquist frequency (dimensional) is pi/dt and corresponds to the max wavenumber (a number) of your Fourier representation, don't make confusion. In your case you can write pi/dt=n_max*2*pi/T -> n_max=0.5*T/dt is the max wavenumber.



Then, T is not necessarily the total time of your simulation but that characteristic physical period of your flow. If you use the total time of the computation, you are implicitly assuming that the flow repeat yourself every T time. This way you need of N*T times samples to perform the statistical averaging in time.
FMDenaro is offline   Reply With Quote

Old   January 10, 2021, 07:40
Default
  #7
New Member
 
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4
Falu is on a distinguished road
Okay, I'm really confused right now I dont get why i would want to know T or the wavenumbers. Im interested in the frequencies. Im doing an fft in the time-domain, so I would say im getting my result in the frequencydomain. I dont need any transformationen to the wavenumberdomain, do I?

And Im still confused about the nyquist frequency defintion. Is my nyquist definition wrong with 5000Hz? [Delta_t=0.0001 -> f_sample=1/Delta_t=10.000Hz -> f_nyquist=f_sample/2).
Because if I use your definiton, im getting a f_nyquist_2=pi/Delta_t=31416Hz.
Falu is offline   Reply With Quote

Old   January 10, 2021, 08:15
Default
  #8
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,252
Rep Power: 67
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by Falu View Post
Okay, I'm really confused right now I dont get why i would want to know T or the wavenumbers. Im interested in the frequencies. Im doing an fft in the time-domain, so I would say im getting my result in the frequencydomain. I dont need any transformationen to the wavenumberdomain, do I?

And Im still confused about the nyquist frequency defintion. Is my nyquist definition wrong with 5000Hz? [Delta_t=0.0001 -> f_sample=1/Delta_t=10.000Hz -> f_nyquist=f_sample/2).
Because if I use your definiton, im getting a f_nyquist_2=pi/Delta_t=31416Hz.

The Nyquist wavenumber can be defined by the ratio between the max wavelenght T and the smallest one representable on a finite grid of size dt. This is

n_max=T/2dt = N*dt/2*dt=N/2.

The Nyquist frequency is

f_max=n_max*2*pi/T= pi/dt
Falu likes this.
FMDenaro is offline   Reply With Quote

Old   January 10, 2021, 11:00
Default
  #9
New Member
 
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4
Falu is on a distinguished road
Thank you very much FMDenaro, your derivation makes perfectly sense. I'm still wondering, why my way of defining the nyqust frequency is wrong? I've tried out both definitions and both lead to the same graph but with a different scaling (see attachments). This is due to my definition of my function, which has an dependency on f_sample:

P1Y =Y_f(1:N/2+1);
PY=(1/(f_sample*N)) * abs(P1Y).^2;
PY(2:end-1)=2*PY(2:end-1);

Does anyone know why my first way of calculating f_nyquist is wrong? I really want to know my error.
Attached Images
File Type: jpg EnergySpectrumVelU6_Denaro.jpg (20.6 KB, 14 views)
File Type: jpg EnergySpectrumVelU6.jpg (22.7 KB, 9 views)
File Type: jpg VelocityU6_Frequenzbereich_Hann-Window.jpg (21.2 KB, 9 views)
File Type: jpg VelocityU6_Frequenzbereich_Hann-Window_Denaro.jpg (20.9 KB, 7 views)
Falu is offline   Reply With Quote

Old   January 10, 2021, 11:04
Default
  #10
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,252
Rep Power: 67
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by Falu View Post
Thank you very much FMDenaro, your derivation makes perfectly sense. I'm still wondering, why my way of defining the nyqust frequency is wrong? I've tried out both definitions and both lead to the same graph but with a different scaling (see attachments). This is due to my definition of my function, which has an dependency on f_sample:

P1Y =Y_f(1:N/2+1);
PY=(1/(f_sample*N)) * abs(P1Y).^2;
PY(2:end-1)=2*PY(2:end-1);

Does anyone know why my first way of calculating f_nyquist is wrong? I really want to know my error.



As you see, your f_nyquist is actually n_max, that is a number, not a frequency. From your f_nyquist you have to evaluate the cut-off frequency.
FMDenaro is offline   Reply With Quote

Old   January 10, 2021, 11:45
Default
  #11
New Member
 
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4
Falu is on a distinguished road
Aaah yes now i get it. Sorry, thank you very much

Another topic, you said earlier i dont need to add zeros to my velocity vector. I had done this, because i want to get 8192 entries (2^13), so that matlab uses the FFT algorithm and not the DFT one. Im wondering, if i corrupted my data with this approach and now the results are wrong and I maybe shouldn't do the addition off the zeros and instead use the DFT algrotihm?

Im havent seen so many spectrum diagramms until now, so I'm not entirely sure my diagramms are totally correct. Especially the higher frequencies in my diagramms are irritating me, but i guess this is due numerical problems with the simulation.
Falu is offline   Reply With Quote

Old   January 10, 2021, 14:20
Default
  #12
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,252
Rep Power: 67
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by Falu View Post
Aaah yes now i get it. Sorry, thank you very much

Another topic, you said earlier i dont need to add zeros to my velocity vector. I had done this, because i want to get 8192 entries (2^13), so that matlab uses the FFT algorithm and not the DFT one. Im wondering, if i corrupted my data with this approach and now the results are wrong and I maybe shouldn't do the addition off the zeros and instead use the DFT algrotihm?

Im havent seen so many spectrum diagramms until now, so I'm not entirely sure my diagramms are totally correct. Especially the higher frequencies in my diagramms are irritating me, but i guess this is due numerical problems with the simulation.



As reported in mathworks:


"You can interpolate the DFT by zero padding. Zero padding enables you to obtain more accurate amplitude estimates of resolvable signal components. On the other hand, zero padding does not improve the spectral (frequency) resolution of the DFT. The resolution is determined by the number of samples and the sample rate."


https://it.mathworks.com/help/signal...o-padding.html
FMDenaro is offline   Reply With Quote

Old   January 11, 2021, 07:17
Default
  #13
New Member
 
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4
Falu is on a distinguished road
Again, thank you very much FMDenaro

Finally I would like to come back to the last open question in this thread regarding the spatial fft analysis. In my understanding the situation is like this:

I have 24 spatial monitoring Points, so my N_sp=24. So n_max_sp=N_sp/2=12. The corresponding nyquist frequency is f_ny_sp=pi/dx =31.4150Hz because the distance between the points is 0.1m. So to improve my resolution in the space domain i would need to install more monitor points. If im doing everything right, i should be able to see the -5/3 slope in the result diagramm(in the wavenumber domain). Finally, when i want to know my total kinetic energy for every time step E is evaluated like this: E_sp=sqrt(ku^2+ky^2+kz^2). ku,ky and kz are hereby my fft transformed velocity signals. With this i would get for all my 8192 time steps a spatial energy plot, which i can finally sum to the total energy: E_sp_total=Sum_i(E_sp_i*dt)/t_total with i={1,...,8192} which results to the total energy density spectrum, avareged over time. Here, you should be able to see -5/3 slope, too (again, wavenumber domain).

Am i getting something wrong here?
Falu is offline   Reply With Quote

Old   January 11, 2021, 09:53
Default
  #14
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,252
Rep Power: 67
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
You need to have at least one direction of homogeneity, for example a spanwise periodic direction. Then You collect all the grid points of tour solution for the 1d FFT
FMDenaro is offline   Reply With Quote

Old   January 11, 2021, 10:38
Default
  #15
New Member
 
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4
Falu is on a distinguished road
Okay, assuming my x direction is homogenous. I still can do a spatial fft with the monitoring points, or am Im wrong? Then my example would be correct if i edit it like this:

I have 24 spatial monitoring Points, so my N_sp=24. So n_max_sp=N_sp/2=12. The corresponding nyquist frequency is f_ny_sp=pi/dx =31.4150Hz because the distance between the points is 0.1m. So to improve my resolution in the space domain i would need to install more monitor points. If im doing everything right, i should be able to see the -5/3 slope in the result diagramm(in the wavenumber domain). Finally, when i want to know my total kinetic energy for every time step E is evaluated like this: E_sp=ku^2. ku is hereby my fft transformed velocity signal. With this i would get for all my 8192 time steps a spatial energy plot, which i can finally sum up to the total energy: E_sp_total=Sum_i(E_sp_i*dt)/t_total with i={1,...,8192} which results to the total energy density spectrum, avareged over time. Here, you should be able to see -5/3 slope, too (again, wavenumber domain).
Falu is offline   Reply With Quote

Old   January 11, 2021, 12:15
Default
  #16
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,252
Rep Power: 67
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
It makes no sense using the probe points, you have to use the whole nodes of the numerical solution.
FMDenaro is offline   Reply With Quote

Old   January 11, 2021, 12:19
Default
  #17
New Member
 
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4
Falu is on a distinguished road
The problem is that my grid resolution changes, the further the flow is from the opening the coarser the grid becomes.
Falu is offline   Reply With Quote

Old   January 11, 2021, 12:25
Default
  #18
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,252
Rep Power: 67
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by Falu View Post
The problem is that my grid resolution changes, the further the flow is from the opening the coarser the grid becomes.

If your flow problem has inflow/outflow, you have walls, you have no homogeneity and it makes no sense a spectral analysis. I think you can only do the temporal analysis.
FMDenaro is offline   Reply With Quote

Old   January 12, 2021, 06:44
Default
  #19
New Member
 
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4
Falu is on a distinguished road
Thank you very much for all your help FMDenaro

Perhaps I should have mentioned that I am simulating a free jet into a large volume. But Im gonna stay now with the temporal analysis.
Falu is offline   Reply With Quote

Old   January 12, 2021, 06:46
Default
  #20
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,252
Rep Power: 67
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by Falu View Post
Thank you very much for all your help FMDenaro

Perhaps I should have mentioned that I am simulating a free jet into a large volume. But Im gonna stay now with the temporal analysis.
Ok, but did you set periodic conditions in the other two directions? If yes you can do the 1d spectra
FMDenaro is offline   Reply With Quote

Reply

Tags
energy, matlab, spectrum

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
[Other] Contribution a new utility: refine wall layer mesh based on yPlus field lakeat OpenFOAM Community Contributions 58 December 23, 2021 03:36
pressure in incompressible solvers e.g. simpleFoam chrizzl OpenFOAM Running, Solving & CFD 13 March 28, 2017 06:49
Stuck in a Rut- interDyMFoam! xoitx OpenFOAM Running, Solving & CFD 14 March 25, 2016 08:09
Moving mesh Niklas Wikstrom (Wikstrom) OpenFOAM Running, Solving & CFD 122 June 15, 2014 07:20
Upgraded from Karmic Koala 9.10 to Lucid Lynx10.04.3 bookie56 OpenFOAM Installation 8 August 13, 2011 05:03


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