
[Sponsors] 
How to calculate predominant shedding frequency in Strouhal Number 

LinkBack  Thread Tools  Search this Thread  Display Modes 
December 8, 2020, 02:22 

#21 
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 611
Rep Power: 14 

December 8, 2020, 11:28 

#22  
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 8 
Quote:
N, Nev are declared in fftGridConvergence.py: https://gitlab.com/mAlletto/openfoam...7573b3f94fe187 How did you decide on these values? I am trying to calculate the sample spacing or time step. I have followed your code and extracted my time values: <times = data[:,0]>. I am now trying to create the linspace code. You use N in that code<t = np.linspace(times[len(times)/2],times[1],N)>, and I am wondering if N is number of samples. What are N and Nev at the beginning of the above link. Also, I get an error when I try your code for <t = np.linspace(times[len(times)/2],times[1],N)>. It tells me that "IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis(`None`) and integer or boolean arrays are valid indices" After debugging a little, I think the problem is that <len(times)/2> produces a double. When I entered times[int(len(times)/2)], the code ran. Finally, why did you divide <len(times)/2> by the two? I ask because I think linspace is my problem. I tried using my times[1]  times[0] because I thought the time between sample was the timestep. I am getting small values in my output so I think my problem is with specifying linspace for sample spacing. Can you explain how you came up with your linspace criteria? As an example, why did you start in the middle of your data, go to end, and decide that it should be divided into 2500 samples. Last edited by HappyS5; December 8, 2020 at 23:28. Reason: clarify 

December 9, 2020, 03:21 

#23 
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 611
Rep Power: 14 
N is the number of interpolation points. Since the output of the force signal is not at equidistant time steps, we need to interpolate the signal to equidistant time steps.
The number N should not be too small to have high computational costs to do the fft later and not too large to resolve the signal. N = 2500 was a reasonable value in my mind. You should maybe resolve on period with 10 points. Nev is when to start to the fft. You want to cut off the initial transition period where the flow transits from the initial state to some periodically oscillating state. So for Nev = 1000 i cut off a bit less than the first part of the signal. I obtained this be visually inspecting the signal. If you want some more theoretical approach how to detect initial transient stats you can read this: https://www.researchgate.net/publica...LENT_FLOW_DATA 

December 9, 2020, 11:45 
Working Dominate Frequency Extraction Snippet in Python

#24 
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 8 
Hello,
After the help from the senior members, I have developed the following working code to do a Fast Fourier Transform and then to extract the dominate frequency. I used Python3.8. Please use Python manual to look up the use of each. Code:
#1/usr/bin/env python3.8 import sys import numpy as np import scipy.fftpack as fftpack import matplotlib.pyplot as plt print("Hello World!") N = 2500 Nev = 1000 data = np.loadtxt('CoefficientLiftData.dat', usecols= (0,3)) times = data[:,0] length = int(len(times)/2) forcez= data[:,1] t = np.linspace(times[length], times[1], 2500) forcezint = np.interp(t, times, forcez) fourier = fftpack.fft(forcezint[Nev1:N1]) frequencies = fftpack.fftfreq(forcezint[Nev1:N1].size, d=t[1]t[0]) #print(frequencies) freq = frequencies[np.argmax(np.abs(fourier))] print(freq) Last edited by HappyS5; December 11, 2020 at 11:28. Reason: Added Python version 

December 9, 2020, 17:10 

#25  
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 8 
Quote:
Quote:
Based off my training, I set up fvSchemes and fvSolutions. I also created the geometry and mesh. The Mesh could be better because it was my first attempt. I think I should have made my area around the cylinder a little larger to have consistent Coefficient of Lift calculations. As for the fvSchemes, it is second order accuracy. I used 2 nOuterCorrectors. Also, I could build a better mesh with less aspect ratio. One more refinement would probably be better. See past post for results. I am running a pimpleFoam only run next. Last edited by HappyS5; December 13, 2020 at 17:48. 

December 18, 2020, 17:17 

#26  
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 8 
Quote:
RE = 250 IcoFoam: Cd = 1.3912; Strouhal number (St) = 0.20618. PimpleFoam w/ SimpleFoam seed: Cd = 1.2723; St = 0.15528. PimpleFoam w/ IcoFoam seed: Cd = 1.2723; St = 0.15320. PimpleFoam: Cd = 1.2738; St = 0.15566 

December 20, 2020, 14:45 

#27 
Senior Member

Seems like wonderful progress!
Refining the mesh seems like a very good idea. In case of value to you (and in that case only), I suggest you place your results in a report. I will be happy to give your report a look. 

December 20, 2020, 19:12 

#28  
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 8 
Quote:
Thanks for the offer! I am retired (BS Chemical Engineering from Oregon State University) so I have no professional value. With that said, I could work on a report. Your support would be evidence that I am learning OpenFOAM. Your support will make me a better owner and administrator at Beginning OpenFOAM on Facebook. Last edited by HappyS5; December 22, 2020 at 22:22. 

January 18, 2021, 17:35 

#29  
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 8 
Quote:
I started getting better results with pimpleFoam after I configured it properly. The Courant number seems to play a major factor in pimpleFoam flow around a cylinder accuracy. Last edited by HappyS5; January 25, 2021 at 15:28. Reason: Hard drive failure 

February 26, 2021, 13:34 

#31  
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 8 
Quote:
Dr. Alletto, I am using the python code that I posted on this page. I have 50,000 datapoints. When I use N=2500 and Nev = 1000, I get 0.0 as a dominate frequency. I changed N to 2500, 7500, 10,000, 12,500, 15,000, 20,000 and still got 0.0. Can you advise me on adjusting Nev and N? I am lost. Last edited by HappyS5; March 5, 2021 at 15:10. 

April 4, 2022, 04:33 

#32 
Member
Bushra Rasheed
Join Date: Dec 2020
Posts: 97
Rep Power: 4 
Hi everyone!
I found this discussion very informative and thought maybe someone can provide me some insights. The latest paraview version now allows to take FFT of a table ie if you calculate lift and drag coefficients through function objects in openfoam and convert .dat file of force coefficients into csv, you can open it in paraview and apply table FFT over it to get shedding frequency. However, I don't know which one is the shedding frequency from the FFT graph. I am attaching the graph here The graph is of FFT of coefficient of lift magnitude. Paraview automatically sets the sampling frequency to 10 Hz (I don't know why it can't detect my time column). Any lead would be appreciated. Thanks! 

Tags 
frequency, strouhal number 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
decomposePar no field transfert  Jeanp  OpenFOAM PreProcessing  3  June 18, 2022 12:01 
how to calculate vortex shedding frequency  mohammad  FLUENT  1  December 1, 2020 18:05 
decomposePar allRegions  stru  OpenFOAM PreProcessing  2  August 25, 2015 03:58 
Cluster ID's not contiguous in computenodes domain. ???  Shogan  FLUENT  1  May 28, 2014 15:03 
AMI interDyMFoam for mixer  danny123  OpenFOAM Running, Solving & CFD  4  June 19, 2013 04:49 