
[Sponsors] 
January 8, 2021, 13:16 
How to plot time dependent energy spectrum

#1 
New Member
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4 
Hey CFDCommunity,
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. 

January 8, 2021, 14:26 

#2 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,599
Rep Power: 70 
First of all, use the log scale but it is also mandatory that you plot the range of frequencies only up to Nyquist cutoff (pi/dt). In a timefrequency analsysis you should adopt a timewindowing technique, why don't you start first by doing a spatial analysis?


January 9, 2021, 05:01 

#3 
New Member
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4 
Thanks FMDenaro for your reply.
I still have a few questions. I have implemented the HanWindow 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 timeenergy 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) 

January 9, 2021, 05:20 

#4 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,599
Rep Power: 70 
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 timesamples. The relation between frequency (f) and wavenumber (n) is f=n*2*pi/T. Therefore you need to evaluate the period T. 

January 9, 2021, 08:03 

#5  
New Member
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4 
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:
Quote:
Quote:
Thanks a lot for your help so far Im making progress. 

January 9, 2021, 08:25 

#6 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,599
Rep Power: 70 
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. 

January 10, 2021, 06:40 

#7 
New Member
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4 
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 timedomain, 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. 

January 10, 2021, 07:15 

#8  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,599
Rep Power: 70 
Quote:
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 

January 10, 2021, 10:00 

#9 
New Member
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4 
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:end1)=2*PY(2:end1); Does anyone know why my first way of calculating f_nyquist is wrong? I really want to know my error. 

January 10, 2021, 10:04 

#10  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,599
Rep Power: 70 
Quote:
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 cutoff frequency. 

January 10, 2021, 10:45 

#11 
New Member
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4 
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. 

January 10, 2021, 13:20 

#12  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,599
Rep Power: 70 
Quote:
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...opadding.html 

January 11, 2021, 06:17 

#13 
New Member
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4 
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? 

January 11, 2021, 08:53 

#14 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,599
Rep Power: 70 
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


January 11, 2021, 09:38 

#15 
New Member
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4 
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). 

January 11, 2021, 11:15 

#16 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,599
Rep Power: 70 
It makes no sense using the probe points, you have to use the whole nodes of the numerical solution.


January 11, 2021, 11:19 

#17 
New Member
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4 
The problem is that my grid resolution changes, the further the flow is from the opening the coarser the grid becomes.


January 11, 2021, 11:25 

#18  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,599
Rep Power: 70 
Quote:
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. 

January 12, 2021, 05:44 

#19 
New Member
Fab
Join Date: Nov 2020
Posts: 25
Rep Power: 4 
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. 

January 12, 2021, 05:46 

#20 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,599
Rep Power: 70 
Ok, but did you set periodic conditions in the other two directions? If yes you can do the 1d spectra


Tags 
energy, matlab, spectrum 
Thread Tools  Search this Thread 
Display Modes  


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 02:36 
pressure in incompressible solvers e.g. simpleFoam  chrizzl  OpenFOAM Running, Solving & CFD  13  March 28, 2017 05:49 
Stuck in a Rut interDyMFoam!  xoitx  OpenFOAM Running, Solving & CFD  14  March 25, 2016 07:09 
Moving mesh  Niklas Wikstrom (Wikstrom)  OpenFOAM Running, Solving & CFD  122  June 15, 2014 06:20 
Upgraded from Karmic Koala 9.10 to Lucid Lynx10.04.3  bookie56  OpenFOAM Installation  8  August 13, 2011 04:03 