What is the best "way" to get vortex shedding ?
 

Hi,

I am super stuck with my problem for a while now and I am running out of ideas.

Here is what I am doing:

I have 3D cases: sphere (D=1m and 2m) in a rectangular domain.

I am able to reproduce the drag coefficients for Re 50, 100, 150 and 200 using icoFoam.

So the tricky part is now, that for Re greater than 250 the flowfield should become asymmetric and for Re greater 300 there should be vortex shedding.

Yeah - you guessed it - that's not happening. The flow stays symmetric and no vortex shedding occurs (ever :( ).

I tried asymmetric initialization of the domain different Re. (up to 500) nothing happens.

I am open for any Ideas. One of my meshes is a snappyHexMesh - so if someone wants to play with the case I could easily zip my case and send it to you.

On the other hand if someone has a working sphere case with vortex shedding I would love to take a look at it maybe I can borrow some ideas.

If you want more Information on the case - I would be happy to post everything you need.

Thank you

- harly

I have a question about the use of terminology:

Doesn't the onset of asymmetry in the velocity field downstream of the sphere beyond a critical Re imply that the flow field is unsteady and that vortex shedding has begun? Or is there a range of Reynolds numbers after the critical Re for which there is no vortex shedding as harly puts it.

my reference
 

Hi,

so from what I am understanding - it seems there are in fact (at least) 3 different "areas".

A symmetric in the region up to Re 200

An asymmetric in the region 200 < Re < 250 in which there is supposed to be a asymmetric field distribution but no shedding - which results in constant Drag and lift coefficients (lift not zero).

A shedding domain at Re > 300

So my reference on this is T.A. Johnson and V. C. Patel

"Flow past a sphere up to a Reynolds number of 300"

yours,
Daniel

P.S. If you want I can send you my reference paper.

A year or so ago I did some vortex shedding simulations with a flat plate.

My experience was that with a poor tetrahedral mesh the vortex shedding started on its own. With a perfectly symmetrical hexahedral mesh this didn't happen.

I kickstarted the shedding by giving a small wall normal velocity component on the top boundary of the channel for a couple of timesteps. ie turning the top boundary temporarely into an inlet. After this kickstart the shedding continued indefinately.

Hey,

thanx for the input - maybe I will give your method a shot. Right now I try to kickstart it with initializing the Field lets say with 0.305 for +y and 0.295 for -y while having an inlet speed of 0.3 (which by the way works like a charm in 2D but somehow it is not working for my 3D{I have to add, that the 2D case is a Re 20000 of a square cylinder computed with turbFoam - which means constant timesteps}).

For your Idea with the mesh - I had the exact same idea and therefore I've created my snappyHexMesh - which is far away form being symmetric in any direction :) - my hexahedrical mesh was created in Gambit.

Both meshs refuse to work, though.

-harly

One question, was the plate 3D or 2D - if it was 3D I would love to see your fvScheme File as well as your init files for U,p,... - of course only if you have the files by your hand if not don't worry - somehow I can't imagine, that there is a mistake in my files that would prevent vortex shedding ...

First a 2D case and then a 3D with oodles. This was a closed channel with the plate not quite touching the sides of the channel. So really a 3D geometry.The 2D definately needed the kickstart, but I can't remember what was the case with the 3D.

Unfortunately I no longer have the case files but I don't think there was anything remarkable in them. I believe you should get vortices with the default settings in icoFoam.

Hey,

Yeah I would have thought so too, but I just don't get it to work - if someone has the time and the motivation to give it a try I would be more than happy to send you my case. Meanwhile I have the following two ideas running:

Calculated a Sphere with Re = 300 with icoFoam until it is 3/4 to convergence than I switched to turbFoam without the turbulencemodels(only for the timestep).

Additional to that I am running a version with a really "big" kickstart and will see how that goes.

I am really desperate right now - I am working for over two months on this problem now ....

-harly

Hi Daniel,

How do your fvSchemes and fvSolution files look like?

If you want, you can post your case here in the forum, and I may give it a try.

Regards,
Jose Santos

link to my case
 

Hi,

here is the link to my case with the variables set to Re = 300 (but feel free to doublecheck):

http://harly.de.vu/www/spherecase.zip

-harly

Used fvSchemes and fvSolution files
 

Hi - here are the config files I am using:

FoamFile
{
class dictionary;
object fvSchemes;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
default CrankNicholson 0.5;
}

gradSchemes
{
default fourth;
grad(p) fourth;

//snGradCorr(cellMotionU) Gauss linear;
}

divSchemes
{
default none;
div(phi,U) Gauss limitedLinear 1.0 phi;
}

laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian((1|A(U)),p) Gauss linear corrected;

laplacian(diffusivity,cellMotionU) Gauss linear corrected;
}

interpolationSchemes
{
default linear;
interpolate(HbyA) linear corrected;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p;
}


// ************************************************** ******* //


FoamFile
{
class dictionary;
object fvSolution;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p PCG
{
preconditioner DIC;
tolerance 1e-07;
relTol 0;
};

U PBiCG
{
preconditioner DILU;
tolerance 1e-06;
relTol 0;
};
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************** *********//

-harly

Hi,

did you try with the linear scheme for the divergence term and backward as time scheme? This would ensure second order accuracy.

Regards,

Hi,

so far I only tried the linear with the CrankNicholson(not backward) - but - isn't that second order too ?

Do you think that would help with the shedding ?

I can give it a try as soon as one of my clients frees up.

-harly

Hi Harly,
I asked because your fvScheme contains CrankNicholson 0.5, while the pure CrankNicholson scheme should be specified with "1", to give second order accuracy (see user's manual, p.U-112).
In my experience backward gives better results in LES for examples, that's why I suggested it. However you should see the shedding also with CranckNicholson.

Best,

Are you getting a very elongated attached wake with time?

Using a second order scheme to trigger shedding does not seem to be the most direct and may be effective way to obtain results. What triggers shedding is the asymmetry due to perturbations from the outside environment or imperfect geometry in real life and numerical inaccuracies in CFD. However, for coarse grids one may be dissipating the small numerical errors/asymmetries, so refining the mesh might help in these cases. If you are going this track, also rotate the mesh geometry by a very small angle to make sure it is no longer symmetric about center plane lined up with the free stream. I think the easiest way to start is to perturb the instability aggressively by introducing an error in a variable like u or p and keep patient. I had a similar situation with my PhD code for flow over cylinder and this approach worked. I think someone above suggested this and I stand behind it. You say you switched to a turbulent solver but with turbulence switched off? What use is this, I could not sense?

Cem

Hi,

I am not able to download your file.

Regards,
Jose Santos

That's exactly the point of using more accurate schemes too: reducing numerical dissipation.

Do you mean introducing "non-orthogonality" in the direction of the flow? This should not be necessary to have the formation of the structures. Anyway, with a sphere inside the domain, you do not have a perfectly symmetric and orthogonal mesh anyway around the obstacle.

Btw, did you see Frank Bos work on the flow around a cylinder? It is discussed here

http://www.cfd-online.com/Forums/ope...-re-3d150.html

and his case files can be downloaded from his website here

http://www.aero.lr.tudelft.nl/~frank/index.php?id=downloads/files

Regards,

@santos

I am sorry about the link, but it works for me if you could send me an email to openfoam.messageboard at gmail.com - I will reply to you and attach my case.


@all

Thanks for the suggestions - in fact I saw the work Frank Bos did with cylinders and I will have a closer look at it now, because I started a 2D case of a circular cylinder and again was not able to obtain vortex shedding. I will report back as soon one of my cases gives me any usable result. But for now the 2D square cylinder is the only one which gives me vortex shedding.

Oh and I checked the snappyHexMesh version of my mesh is not even symmetric to the symmetry planes - so it should be the perfect mesh ...

-harly

Hi Daniel,

I had a look at your case, and I think your mesh is not fine enough in the region around the sphere. The critical region may be the wake downstream the sphere.

Regards,
Jose Santos

Thx for checking my case, so you think the problem is in the mesh? - I will create a finer mesh and will report back before I let it run. Do you have any suggestions regarding the schemes and solver I use ? Once the Calculation is running it will probably take 5 days depending on how much bigger the mesh gets.

-harly

@albcem:

I am very sorry somehow I overread your first post. Can you explain to me how you would preturbate the flow ? And to your question of the turbulence model switched off - generated a question below.

@all

So from what I gathered:

- a switch to complete 2nd order aka backward in time and Gauss linear in space

- a refinement of the mesh

- preturbation of flow

could do the trick, right ?

Which solver would you use ? - so far I was using icoFoam, because I was under the impression I don't need a turbulence model like for the 2D case Frank Bos presents. That explains also, why I used turbulence model switched off in turbFoam. - But maybe that could be the mistake?
-harly

IIRC at least in some version of OpenFOAM it was:

Turbulence off:
Turbulence is not solved but the turbulent viscosity is added to the effective viscosity. Meaning that if there are non-zero k and epsilon fields present the effective viscosity is higher than the laminar viscosity.

"Laminar" turbulence model:
Returns a zero turbulent viscosity -> correct laminar viscosity.

Just a note, as Hrv suggested, to have second order on all meshes, use leastSquares for gradient evaluation. However the use of Gauss should not prevent the vortex shedding.

Imho this is not necessary. As you can see in the case of the cylinder done by Frank Bos, he doesn't perturbate the flow, and the structures start to form immediately.

Yes, icoFoam should do the job if your Re is similar to those considered by Frank Bos. Btw, what is your Re, defined as U*d/nu, where d is the cylinder/sphere diameter?

Regards,

So,

3D simulation is on its way - once it is finished I will post pictures (pressure, velocity plots) and then we can discuss the matter further.

I use the "standard" formulation of the Reynoldsnumber: Re = U*d/nu and I am using air properties for T=20(where the heck is the degree symbol on an American keyboard?) Celcius.

In the meantime I will set up a 2D case for a circular and a square cylinder in low Re(around 100 - that vortex shedding occurs) and will then try to verify the results with a reference. I will probably be able to post the first results in a couple of hours plus pictures of the meshs I used.

-harly

... couple of hours later ...
 

Hi everyone,

so here are some interesting results I'd like to discuss with you:

I took the case Frank Bos offers on his homepage and ran the case with a force configuration in the controldict[1] as presented in the Thread Forces in OF-1.5 and a different viscosity (0.01 instead of 0.0066667) to realize Re = 100.

You can see the C_L / C_D - time - plot for the original(Re=150) here.
http://www.harly.de.vu/www/cylinder_Re150.jpg
You can see the C_L / C_D - time - plot for the Re=100 here.
http://www.harly.de.vu/www/circular_...nk_Bos_50k.jpg
Both were created on the 50k mesh

I took this as my reference for how my cases should look like. The first thing I wanted to change was the parameters Frank used (he sets U to 1 and adapts nu to get the Re) and replace them with "my" standard values for air:

density: 1.204
dynamic viscosity(eta): 1.83692474747684e-5
kinematic viscosity(nu): 1.525684971e-5

with d=1 that would result in U=0.001525....

From here I wanted to reproduce a result I have in a reference[2] for a circular cylinder in 2D at Re = 100.
The results should look as following:
http://www.harly.de.vu/www/reference_Re100.jpg
(I will try to scan the picture, unfortunately I only have a hard copy of my reference - so I took a picture with my camera)

Here is a (bad) picture of the original mesh the results were obtained with:
http://www.harly.de.vu/www/reference_mesh.jpg

Unfortunately I couldn't get any useful results with the "real" air configuration and my question would be WHY? - I figured they would just take longer to get some results but even running them overnight didn't bring useful results, !I could not obtain shedding! - this is very important can someone please comment if she/he had the same experience? Maybe that can solve the problems I am having with my sphere? - What do you think?

I took the cases I had and switched to U=1 and calculated the viscosity with Re = U*d/nu (with U and d always being 1 that means nu = 1/Re).


and here are the competing meshs:

Frank Bos(50k):
http://www.harly.de.vu/www/50k.png
Frank Bos(100k):
http://www.harly.de.vu/www/100k.png
A mesh from me, created in Gambit:
http://www.harly.de.vu/www/mymesh.png
and here are the results:

50k :
http://www.harly.de.vu/www/circular_...nk_Bos_50k.jpg
100k :
http://www.harly.de.vu/www/circular_...k_Bos_100k.jpg
my mesh:
http://www.harly.de.vu/www/circular_Re100_own_mesh.jpg
I packed my case for you to have a look at:

here

I didn't want to repack Frank's stuff so please visit his website and download his testcases:
here

As for the conigurations which were used in all three cases I summed up the most important bits below:

-----------------------------------------------------------------------------------------------------
ddtSchemes
{
default CrankNicholson 0.5;
}

gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}

divSchemes
{
default none;
div(phi,U) Gauss linear;
}

laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian((1|A(U)),p) Gauss linear corrected;
}

interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p;
}
--------------------------------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------------------------------
solvers
{
p ICCG 1e-06 0;
U BICCG 1e-05 0;
}

PISO
{
momentumPredictor yes;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
---------------------------------------------------------------------------------------------------------------

Any suggestions in how I would get my results closer to what it is supposed to be? Please check my controlDict if I am making a stupid mistake while calculating the forces. And Frank if you read this maybe you could tell us if your results were of the same magnitude. Looking forward to your comments.

-harly

[1]: link to the used controlDict
[2]: Computational Modelling of vortex shedding flows by Vlado Przulj

Hi Daniel,

I encounter the same problems as you with respect to the absence of any vortex shedding behind the sphere. Did you already solve the problem with the hints of your last message?

At the moment, I am running a case with a perturbated inital velocity field, but I am pessimistic about the result...

Many thanks.

Regards,
Bob

I had similar problems when running a LES with wall-functions.
Using a wall-resolved grid helped to get rid of this problem.

Markus

I am in the same point as you with a cube centered in the domain by using RANS models...
did you find the solution???

Antonio

Hi Harly and All,

Thanks for the posting, it's a good discussion. I'm also doing about the same simulation, flow over a square, 2d at Re= 22k. I have three grid resolution, course:30cells x 4 edges, medium:50 cells x 4 edges and fine: 80 cells x 4 edges. The course and medium have wall function and the fine without wall function.

The results for coarse and medium are about the same for experiments, but for fine grid is not.

Any comment why the fine grid gives the 'wrong' results? The y+ for fine grid is less than 3, so I don't use wall function.

Thanks,
Mali

Dear harly,
Did you figured out why you sometimes cannot obtain shedding?

I'm doing simulation around rectangular prism with Re 1e5 and I can't get oscillations of Cd and Cl and vortex shedding using RANS.
The results of LES calculations are ok.

Regards,
Paulina

Paulina, can you post some of your log output?

I'm a little bit out of topic, my calculations were made in Fluent.

I've tried to recreate simulation from an article (http://www.sciencedirect.com/science...43974X06001593) were they used same grid for LES and RANS, and enhanced wall funcion for Kato Launder k-epsilon.
Using above mentioned configuration my Cd and Cl were completely flat (as if it was laminar flow).

Yesterday I run simulation for standard wall function for a grid with higher y+. I've got oscillations of Cl, but after 2s it went flat.

I still don't know what is missing in the calculations (grid? set up?).

How long did you let it running after the decay of the initial oscillations? It sometimes took a really long time in Fluent until the actual von Karman street started.

The longest time was about 4 seconds. I assumed that nothing new will happen later. The period of one oscillation (from the article and my LES simulation) is about 0,13s.

While using same grid as for LES initial oscillations were quite small and ended even before 1s.

Do you think that 10s would be enough for vortex shedding to occur?

I'll run simulations few more hours and let you know about the results.

How long is your time step and did you switch all the numerical settings to 2nd order?

Yes, all numerical settings are in 2nd order (And Fluent 15 run in double precision). BC are velocity inlet and outflow.

The time step is 0.002s or 0.001s. I used even 0.0005s but haven't seen difference in the results.
As a convergence criteria I set up epslion equal to 1e-7, so now it is 5 iterations per 1 time step.

You can see settings and how Cd and Cl look after few more seconds of calculations.

I can not guarantee that, but I once had this behavior and after a long time of waiting the vortex street started. If you have time, give it a try ;)

10 seconds passed and nothing happened. The flow field looks like after steady flow calculations.

I think this is a matter of bad lack, or rather bad mesh.

I'm out of ideas what to change to get right results.

I've tried few more things but the situation is still the same.
Initial oscillations are gone after few seconds and nothing else happens.

If you would like to try and check here you can find my mesh and inlet profiles:
https://www.dropbox.com/sh/5hui9t2vk...toMth-SJa?dl=0

Hi all,
Did you solve your problems? I was actually solving flow past acylinder with D/H=0.3 and Re=100. I am getting correct strouhal number but getting very low drag co-efficient. Can you tell me why? I am really stuck at this for relly long time. I am using openFoam 3.0.1 and ubuntu 16.04. i have attached my case.

Thanks in advance...

regards
Subhankar

Subhankar,

Have u try higher order scheme?