|
[Sponsors] |
How to calculate predominant shedding frequency in Strouhal Number |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
December 8, 2020, 03:22 |
|
#21 |
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 616
Rep Power: 16 |
||
December 8, 2020, 12:28 |
|
#22 | |
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9 |
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 9, 2020 at 00:28. Reason: clarify |
||
December 9, 2020, 04:21 |
|
#23 |
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 616
Rep Power: 16 |
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, 12:45 |
Working Dominate Frequency Extraction Snippet in Python
|
#24 |
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9 |
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[Nev-1:N-1]) frequencies = fftpack.fftfreq(forcezint[Nev-1:N-1].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 12:28. Reason: Added Python version |
|
December 9, 2020, 18:10 |
|
#25 | ||
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9 |
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 18:48. |
|||
December 18, 2020, 18:17 |
|
#26 | |
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9 |
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, 15: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, 20:12 |
|
#28 | |
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9 |
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 23:22. |
||
January 18, 2021, 18:35 |
|
#29 | |
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9 |
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 16:28. Reason: Hard drive failure |
||
February 26, 2021, 14:34 |
|
#31 | |
Member
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9 |
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 16:10. |
||
April 4, 2022, 05:33 |
|
#32 |
Member
Bushra Rasheed
Join Date: Dec 2020
Posts: 97
Rep Power: 5 |
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 i-e 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 Pre-Processing | 3 | June 18, 2022 13:01 |
how to calculate vortex shedding frequency | mohammad | FLUENT | 1 | December 1, 2020 19:05 |
decomposePar -allRegions | stru | OpenFOAM Pre-Processing | 2 | August 25, 2015 04:58 |
Cluster ID's not contiguous in compute-nodes domain. ??? | Shogan | FLUENT | 1 | May 28, 2014 16:03 |
AMI interDyMFoam for mixer | danny123 | OpenFOAM Running, Solving & CFD | 4 | June 19, 2013 05:49 |