[mAlletto]

Originally Posted by HappyS5

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

[HappyS5]

Originally Posted by mAlletto

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.

[mAlletto]

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]

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.

#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.

[HappyS5]

Originally Posted by dlahaye

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.

Originally Posted by dlahaye

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.

[HappyS5]

Originally Posted by HappyS5

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

[dlahaye]

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]

Originally Posted by dlahaye

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.

[HappyS5]

Originally Posted by dlahaye

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]

Wonderful! Keep us updated!

Domenico.

[HappyS5]

Originally Posted by mAlletto

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.

[B_R_Khan]

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!