CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

How to calculate predominant shedding frequency in Strouhal Number

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

Like Tree25Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   December 8, 2020, 02:22
Default
  #21
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
Quote:
Originally Posted by HappyS5 View Post
Michael,

Thanks for the excellent contribution. It has helped me progress.

Where did you get the 2500 for N? What is N? What is Nev?



Can you explain a bit better. I do not understand what you mean
mAlletto is offline   Reply With Quote

Old   December 8, 2020, 11:28
Default
  #22
Member
 
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9
HappyS5 is on a distinguished road
Quote:
Originally Posted by mAlletto View Post
Can you explain a bit better. I do not understand what you mean

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
HappyS5 is offline   Reply With Quote

Old   December 9, 2020, 03:21
Default
  #23
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
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
HappyS5 likes this.
mAlletto is offline   Reply With Quote

Old   December 9, 2020, 11:45
Default Working Dominate Frequency Extraction Snippet in Python
  #24
Member
 
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9
HappyS5 is on a distinguished road
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)
I know this code works because it matched the Strouhal numbers I mentioned previously. As usual, Python estimates better than I, but while the simple sinusoidal method of getting the period between two Coefficient of Lift maximum's for the predominate shedding frequency exists, a senior member mentioned the versatile benefit of FFT. I can verify the usefulness of Python.
raj kumar saini likes this.

Last edited by HappyS5; December 11, 2020 at 11:28. Reason: Added Python version
HappyS5 is offline   Reply With Quote

Old   December 9, 2020, 17:10
Default
  #25
Member
 
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9
HappyS5 is on a distinguished road
Quote:
Originally Posted by dlahaye View Post
Very happy to read about the progress you are making.

I am slightly confused about the following two issues.

Issue 1: my (always limited) understanding is that the oscillation of lift vs. time (and thus the Strouhal number) should be *independent* of the initial guess for the transient solver that you impose. That is, the initial guess only affects the value for the lift in the initial stages of the solver output. After some time, the curve of output (lift) vs. time is independent of the initial guess that you impose. It is true that initial transients disappear slower or faster depending on the initial guess imposed. After sufficiently long time however, the initial guess is completely removed from the simulation results.
From my recent experience, I got different Coefficient of Lift values, and nearly the same Strouhal number, for transient solver (pimpleFoam) seeded with steady state(simpleFoam) solver latest time folder and transient(icofoam) latest time folder as a seed(used as 0 folder). With that said, you are much more knowledgeable than I. I have not done an analysis of the effect of "seed" and final results, but I was looking for a difference.

Quote:
Originally Posted by dlahaye View Post
Issue 2: I am unfamiliar with the differences between icoFoam and pimpleFoam. I am thus wondering whether the differences you observe are due to solver settings (mesh, outer correction, relative residuals, etc)

I am keen to see your FFT results. Cheers.
Although IcoFOAM gave me the best results, I did not originally choose IcoFoam. It was suggested by a tutorial. It is a laminar only solver and pimpleFoam uses PISO and is a laminar or turbulence solver.

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.
HappyS5 is offline   Reply With Quote

Old   December 18, 2020, 17:17
Default
  #26
Member
 
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9
HappyS5 is on a distinguished road
Quote:
Originally Posted by HappyS5 View Post
I have completed all the above steps, and will now concentrate on the FFT data and dominant frequency extraction from the FFT data. I am interested to see the results.


I have new and better preliminary results, from steady-state to transient suggestion by dlahaye, that more closely match literature that I have seen published. My Coefficient of Lift = 1.27 which approximates the trend for data that covers Re=50 to Re=200. That data had a slight downward trend and the Coefficient of Lift for Re = 200 was 1.30-1.34. My estimated Strouhal number (St) = 0.156 and about 0.02 off from calculated value. Literature also shows that the St number should be near 0.2 and greater than 0.195. I will have a better estimate after I perform FFT. The article linked in the above post suggests that Strouhal Number should be near ST = 0.2 which I get from IcoFoam. Note, the calculated Strouhal number is 0.18 at Re=250.

IcoFoam gave me results of 1.356 Coefficient of Drag, and a Strouhal number = St = 0.2.
Correction! I forgot the Coefficient of Drag (Cd) oscillates too so the Cd has to be an average.

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
HappyS5 is offline   Reply With Quote

Old   December 20, 2020, 14:45
Default
  #27
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 722
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
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.
HappyS5 likes this.
dlahaye is offline   Reply With Quote

Old   December 20, 2020, 19:12
Default
  #28
Member
 
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9
HappyS5 is on a distinguished road
Quote:
Originally Posted by dlahaye View Post
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.

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.
dlahaye likes this.

Last edited by HappyS5; December 22, 2020 at 22:22.
HappyS5 is offline   Reply With Quote

Old   January 18, 2021, 17:35
Default
  #29
Member
 
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9
HappyS5 is on a distinguished road
Quote:
Originally Posted by dlahaye View Post
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.
I am working on this project. Sadly, my Linux system crashed, because the SanDisk hard drive failed, and I can't even log on. I recovered my project data, and working on deleting unneeded data.




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.
dlahaye likes this.

Last edited by HappyS5; January 25, 2021 at 15:28. Reason: Hard drive failure
HappyS5 is offline   Reply With Quote

Old   January 19, 2021, 02:51
Default
  #30
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 722
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
Wonderful! Keep us updated!

Domenico.
HappyS5 likes this.
dlahaye is offline   Reply With Quote

Old   February 26, 2021, 13:34
Default
  #31
Member
 
Chris Harding
Join Date: Dec 2016
Posts: 76
Rep Power: 9
HappyS5 is on a distinguished road
Quote:
Originally Posted by mAlletto View Post
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

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.
HappyS5 is offline   Reply With Quote

Old   April 4, 2022, 04:33
Default
  #32
Member
 
Bushra Rasheed
Join Date: Dec 2020
Posts: 97
Rep Power: 5
B_R_Khan is on a distinguished road
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!
Attached Images
File Type: png fft.png (19.1 KB, 12 views)
B_R_Khan is offline   Reply With Quote

Reply

Tags
frequency, strouhal number

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
decomposePar no field transfert Jeanp OpenFOAM Pre-Processing 3 June 18, 2022 12:01
how to calculate vortex shedding frequency mohammad FLUENT 1 December 1, 2020 18:05
decomposePar -allRegions stru OpenFOAM Pre-Processing 2 August 25, 2015 03:58
Cluster ID's not contiguous in compute-nodes domain. ??? Shogan FLUENT 1 May 28, 2014 15:03
AMI interDyMFoam for mixer danny123 OpenFOAM Running, Solving & CFD 4 June 19, 2013 04:49


All times are GMT -4. The time now is 17:19.