Courant Number become bigger and bigger!!!
 

I just begin to compute my case with icoFoam.

But the Courant Number becomes bigger and bigger even I use very small time step such as 1.0e-7. And I also assigned the maxCo with only 0.1. But it doesn't work. I really don't know the reason. If anyone has some experience, please share me with it. Or if you have some good advice, please inform me.

Thanks!!

Hi,

in order to use adpative timestep with icoFoam you need to follow this thread:
http://www.cfd-online.com/Forums/ope...co-number.html

and concerning instability of your simulation, from my experience this might be some mesh issue or messed boundaries.
you should plot here your boundaries settings.
and do checkMesh.

Basicaly there's no general answer for your question, you must clear it out a bit.

Regards,
Chris

i had same situation with rhoSonicFoam .. I did follow the thread chris posted .... i included the "setDelta.H " "readTimeControls.H" files. This hasnt solved the problem entirely , even though the solution runs for quite a long time.. Maybe Initial timestep is a key .. As im new to both CFD & Openfoam [taking a class in CFD this sem] , thats what i guess
@ flying .... keep posting about your progress .im sure we can share info

In fact, I also try it to fixed it as Problems with adjustable timestep control and maxCo Number, but it doesn't work, and I really don't know how to do it. Then I change my dt is very small about 1e -6. Then it works. In fact, maxCo is only 0.03. So it is still not well. Do you have some other way to solve this problem?

Thanks and I would like to share information with you!

Hi,

as chris_sev already mentioned: the problems of your small timestep maybe related to either mesh issues or messed up boundary conditions (or both). Run a checkMesh and post the output. Some information on your BC wouldn't be so bad either ;)

so long

as i posted in my thread:

1. check your boundary conditions. fixedValue for U and p at the inlet lets the firtst cells "explode" (U=4e25 after some iterations)

2. try to use bounded discretisations. rhoSonicFoam needs bounded schemes for div, like upwind or limitedLinear.

/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.4.1 |
| \\ / A nd | Web: http://www.openfoam.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/

Exec : checkMesh . .
Date : Apr 22 2009
Time : 09:04:16
Host : bonito
PID : 24378
Root : /home/msamel/OpenFOAM/OpenFOAM-1.4.1/tutorials/rhoSonicFoam/more_fine
Case : .
Nprocs : 1
Create time

Create polyMesh for time = constant

Time = constant

Mesh stats
points: 15328
edges: 37947
faces: 30140
internal faces: 14813
cells: 7520
boundary patches: 7
point zones: 0
face zones: 0
cell zones: 0

Number of cells of each type:
hexahedra: 7353
prisms: 167
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 0
polyhedra: 0

Checking topology...
Boundary definition OK.
Point usage OK.
Upper triangular ordering OK.
Topological cell zip-up check OK.
Face vertices OK.
Face-face connectivity OK.
Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
Patch Faces Points Surface
inlet 10 21 ok (not multiply connected)
wall 217 436 ok (not multiply connected)
outlet2 60 121 ok (not multiply connected)
axis 0 0 ok (empty)
frontAndBackPlanes 0 0 ok (empty)
front 7520 7748 ok (not multiply connected)
back 7520 7748 ok (not multiply connected)

Checking geometry...
Domain bounding box: (0 0 -0.763521) (150 35 0.763521)
Boundary openness (-5.74591e-19 6.68058e-16 8.53495e-16) OK.
Max cell openness = 1.8787e-16 OK.
Max aspect ratio = 82.5128 OK.
Minumum face area = 0.0054537. Maximum face area = 1.37128. Face area magni tudes OK.
Min volume = 0.00489197. Max volume = 0.815718. Total volume = 2830.48. Ce ll volumes OK.
Mesh non-orthogonality Max: 8.53774e-07 average: 0
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.332699 OK.
Min/max edge length = 0.0218148 1.52704 OK.
All angles in faces OK.
Face flatness (1 = flat, 0 = butterfly) : average = 1 min = 1
All face flatness OK.

Mesh OK.

End
^^^^ my checkMech ... btw this is an axisymetric mesh ... prepared using fluentmeshToFoam & makeAxialmesh

p
internalField uniform 50000;

boundaryField
{
inlet
{
type zeroGradient;

}

outlet2
{
type fixedValue;
value uniform 1e5;
}

front
{
type wedge;
}


back
{
type wedge;
}

wall
{
type zeroGradient;
}

defaultFaces
{
type empty;
}
}

U
internalField uniform (280 0 0);

boundaryField
{
inlet
{
type fixedValue;
value uniform (360 0 0);
}

outlet2
{
type zeroGradient;
}

front
{
type wedge;
}

back
{
type wedge;
}

wall
{
type fixedValue;
value uniform (0 0 0);
}

please let me know if you neeed any more info ..

thanx a lot

Create time

Create polyMesh for time = constant

Time = constant

Mesh stats
points: 56392
faces: 163161
internal faces: 157599
cells: 53460
boundary patches: 6
point zones: 0
face zones: 1
cell zones: 1

Number of cells of each type:
hexahedra: 53460
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 0
polyhedra: 0

Checking topology...
Boundary definition OK.
Point usage OK.
Upper triangular ordering OK.
Topological cell zip-up check OK.
Face vertices OK.
Face-face connectivity OK.
Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
Patch Faces Points Surface
WALL1 864 900 ok (not multiply connected)
WALL2 1188 1224 ok (not multiply connected)
WALL3 1548 1584 ok (not multiply connected)
WALL4 1152 1188 ok (not multiply connected)
INLET 405 424 ok (not multiply connected)
OUTLET 405 424 ok (not multiply connected)

Checking geometry...
Domain bounding box: (-0.088514 -0.164031 -0.0121122) (0 0.0210195 0.0121122)
Boundary openness (-5.56024e-16 -6.2903e-17 -1.37923e-16) OK.
Max cell openness = 2.94802e-15 OK.
Max aspect ratio = 79.7319 OK.
Minumum face area = 1.15377e-08. Maximum face area = 9.02275e-06. Face area magnitudes OK.
Min volume = 2.20732e-11. Max volume = 1.19129e-08. Total volume = 4.45843e-05. Cell volumes OK.
Mesh non-orthogonality Max: 56.8885 average: 11.598
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 1.04553 OK.
All angles in faces OK.
Face flatness (1 = flat, 0 = butterfly) : average = 0.999978 min = 0.991932
All face flatness OK.

Mesh OK.

End

did you check your boundary types ???/ that is also an issue

In the forwardstep tutorial , the top wall of the domain is designated as symmetryPlane , & not as a patch . I had been designating the upper wall of my axisymmetric nozzle as 'patch' . i changed it to symmetryPlane & now its working fine ....

habe to see the final solution though

I really don't think it will be work:

The problem is that we can't set the Maxco. And the function for setting doesn't work well.

Although the Courant number has relation with the boundary condition and mesh quantity, it should be set the maximum value.

It is very pity for me to know few how it implement this function in the openfoam.

implemeting maxCo is a easy thing ... you just have to add a few lines to your code & compile it ...
just like how its explained in the thread ..

I think there is no problem about my boundary condition and with small time step, I have got some result about my case, and it seems reasonable. If there is some problem about the boundary, it should not be reasonable

If I didn't misremember it, there are lines that introduced in the thread. But it doesn't work, maybe I need to try again and compile it again. I will try again as the thread. Thanks for your reply.

check your PM !!!

I am sorry I don't know the meaning of "PM". Thanks!

check you personal inbox of this forum .... ive given my gmail id.. i can tel ll what you have to do for variable timestep for rhoSonic.. i guess i know why its not working for you

if the coNum gets bigger and bigger it has nothing to do with variable timestepping.

if the solver isnt stable with very low deltaT then it wouldn't be stable for variable timestepping, either. CoNum is just a nessesary condition for stability.
the icoFoam solver is just a simple solver for incompressible N-S-Equations with a PISO correction.
For most FV solvers it is enough if you use a CoNum of <0.5 . The CoNum is U*deltaT/deltaX. I dont know the maximum of U in the N-S Equations but for the Euler Equations (rhoSonic) this is U+c, cause this is an eigenvalue of this system and so its a speed of error waves for this equations.
if you know the maximum speed of error waves, you can set your CoNum small enough and you have the nessesary condition for stability.

if the solver is still not stable, then its not a problem of the CoNum. then its a problem of your BC or of your schemes.

For example: if you use a linear scheme for the div term in rhoSonic then your solver will diverge no matter what you do. Even for deltaT=1e-100 it will crash cause it is an unbounded scheme what isnt stable for these kind of equations.
so try some more stable schemes like upwind or minmod.

hey thats true!! sorry for not mentioning that ,!!!! thanx

Hi,

I used sonicFoam for solving the flow through a nozzle. I encountered the same problem of diverging Courant number after a few iterations.

Changing the numerical scheme for the divergence terms did solve the problem. After some time, I found that my specific heat capacity (c_v) was specified incorrectly. When I changed c_v from 1.78571 to 717.51 and started the iteration process, the Courant number did also diverge, however it is stable over many time step before it suddenly diverges.

Can anyone explain this to me?

did you just change Cv and nothing else?

if you just set up Cv and don't change R and U, then you get a completely different case.

with Cv the speed of sound is defined as:
c=sqr(R*p/rho*Cv + p/rho)

if you just change Cv, you change c.

check all other variables, possibly there is a problem.

I changed the Cv value, only. Now, I'm sure that I defined R and Cv correctly.

After some investigation of the flow, I encountered another possible source for the diverging Courant-number. The following happens:

I'm simulating a supersonic convergent-divergent nozzle. There is a shock in the divergent part of the nozzle causing the boundary layer to separate, which creates a vortice. This vortice is transported towards the outlet. When it hit the outlet, there is a sudden and dramatic increase in velocity in this region, which is pointing into the nozzle. Thus the Courant-number increases.
I'm quiet sure that this increase in velocity is artificial. I suppose, it occurs because of reversed flow due to the vortice. The BC at the outlet obivously does not let the vortice pass. I already used the 'totalPressure' BC and a 'fixedValue' BC for the outlet, however, the effect occurs for both types.

Can anyone suggest a BC which is more appropriate?

ok, i never used 'totalpressure', but 'fixedValue' is wrong.

for mu=0 the sonicFoam solver solves the Euler-Equations.
These Equations have the eigenvalues v-c, v and v+c.
For a supersonic flow all these eigenvalues are >0 so error shocks run all in your flow-direction.

So have to set 3 fixedValues at inlet and no value at the outlet.
try 'zeroGradient' at the outlet.

ok here is my update which i also posted in another thread :

Hey ive been following this discussion

I am trying to model barrell shocks in an axisymmetric model. Inlet air is M=1. The exit is a subsonic outlet.

I am using rhoSonicFoam [modified to read p , rho, T, U fields , so that the solver can accept derived BCs ] , Im using non-reflective BCs at the exit for rho, since for subsonic outlet , the eigenvalue correspding to rho is -ve. However i am not getting the inlet BCs correct .


I ll summarize the Bc i have tried :--

p
Inlet : totalPressure outlet : fixedValue [ static]

rho
inlet: fixedValue outlet : nonReflective.

U
inlet : 350 outlet : zeroGradient

T
Inlet : 250 outlet : zerogradient

U & T are fine..
kindly suggest me a better combination ....

I used the BC 'waveTransmissive' in order to have a non-reflecting BC. Find more infos here: http://www.openfoamwiki.net/index.ph...dary_condition

and here:

http://www.cfd-online.com/Forums/ope...-velocity.html

hey thank you ..

ok , they have said that the "waveTransmissive" BC is more general , but just for confirmation . can it be used for U & p [subsonic outlet ] both ???

You can use is for U, as well as for p. You have to set the 'field' variable in the definition of the outlet, accordingly:

I don't know how adequate this BC is for subsonic outlets. All I know is that it simulates a farfield. This way, waves should be transmitted and not reflected.

hey thanx..
t
inlet fixed
outlet zeroGradient

these r my BC .are they correct ? . ill proceed with my new simulation & see where it leads me

Well, I think your BCs are correctly defined. Let me know if your results are okay.

CFL criterion
 

I am currently simulating vortex breakdown in a conical diffuser and am doing unsteady simulations using transientSimpleFoam and LAunderGibsonRSTM model.

I have a doubt. The generally suggested criteria that CFL<1, is it the Courant Number mean value or the max courant number value???

I have a mean Courant Number of close to 0.13 and Maximum Courant Number of close to 7.

(Courant Number mean: 0.130184 max: 7.10512 velocity magnitude: 3.59437)

And I am able to capture the vortex. (URANS using conventional k-e and k-w models dont work!!!!)

Is it the correct solution that I have got or is the solution unphysical as the CFL condition states??

So you know im trying to simulate barrel shock & mach disc in a free underexpanded jet . These BCs are a much better improvement over the previous results . However , the reflections arent completely eliminated , just damped. So for initial runtimes the results are very satisfactory , with very sharp shockfronts, & satisfactory comparison with literature data. However, as time progresses , the reflections still prevail. If you want i can post some results, [just tell me how to ] . FOr now im pasting the fluxes across the each boundary field which may give a clue as to whats going on :--
P.S. I ve been playing around with the lInf value ... varying it by an order of a decimal

Dear Statesman,

for my simulations (critical CD-nozzle) I'm using the 'rhoCentralFoam' solver. I read, that it should be more precise than the other incompressible solvers http://openfoamwiki.net/index.php/TestLucaG

Maybe, this will help.

I'm currently using the rhoSonicFoam solver (OpenFoam 1.6) and my Courant number is exploding as well, although I'm sure my bc are correct, since my case is very similar to the forewardStep tutorial.
Therefore I went back to the forewardStep tutorial and changed the mesh. I just moved the vertice of the obstacle corner a little right (from (0.6 0.2 z) to (1 0.2 z)). But even with this small change the courant number explodes after a while. Since the changing of the mesh causes some cells to be deformed and some getting a bit smaller I decreased the time step to 0.001 (from 0.002). But it didn't help.

For me rhoSonicFoam only worked for the turorial cases. What can I do?

sandrak, plese give some more info.

what are your BCs?
what is your U?
how is your start CoNum?
what are your initial values?
what schemes do you use?

1. rhoSonic is just stable for Ma>2. If you use a lower Ma, you need a far smaller timestep.

2. you have to use bounded schemes for div. rhoSonic solves the Euler Equations, they are convection dominated and if the convection is not bounded, the solver is unstable.

Hi Sandrak,

maybe you should try to use 'rhoCentralFoam'. According to the literature I gave in my last post, this solvers is more accurate than 'rhoSonicFoam'. Furthermore, from my point of view it seems to be even more stable, as well, at least in my case.

Generally, my problem is that pressure waves are not transmitted across the output. They are reflected, even tough, I use the 'waveTransmissive' BC for p at the outlet (I am currently testing different values for lInf, so maybe I will get better results in some time). With 'rhoSonicFoam' my simulations always crashed after the pressure wave has been reflected. With 'rhoCentralFoam' the pressure wave is reflected as well, however, the simulation does not crash. It's just that the pressure wave will travel towards the inlet of my nozzle. Thus, I suppose that 'rhoCentralFoam' is more stable. Furthermore, the shock is resolved more distinct with 'rhoCentralFoam' than with 'rhoSonicFoam'.

Sandrack

So you are where each of us [ me , Julian , Joern] was a few months ago..

Please post your BCs that you are applying so that we can know better what you are actually doing.

• So yes the problem you are facing is more than just of timesteps. As Julian said , study the Eigenstructure of the Euler equations for inviscid flow & you shall understand what he means by ""pressure waves are not transmitted across the output""about & how that affects the BCs . You may need nonReflective boundary conditions & rhoSonic isnt the best code for application of these. Try sonicFoam which applies waveTransmissive [which is a non-reflective BC] .
• & Yes you need bounded schemes , for divSchemes. try the TVD family of schemes to know which one is best suited for your case . You may wanna read how TVD schemes behave w.r.t. shock-resolution.

Just read CFD texts to understand the above mentioned items


Julian
I ve been using sonicFoam with Minmod schemes for divSchemes. & backward for ddtSchemes. SO far I'm very happy in terms of both stability & results.
However would you mind sharing with me the boundary conditions of rhoCentralFoam ? I used the same as i posted earlier , but my code wouldn't run.

Thanks for the last post Mihir, you have stated the problem quite well.
Here are my boundary conditions I use for my simulation and a brief explanation of what I am doing:

My domain is 2D, axi-Symmetric and consists of a CD noozle,only. That means, I do not have a far field before or after the nozzle. That's probably why I get problems with the outlet BC, because, if I had a far field, the pressure waves could leave the outlet and due to increasing cell size, the waves would be damped, so that at the farfield outlet, there would be a weak BC interaction, only. Anyway, I induce the flow with a pressure difference at dp=300mbar. At the inlet I have atmospheric pressure and thus at the outlet the low pressure.
I am using rhoCentralFoam.
I have uploaded a little video of the pressure contours:http://www.youtube.com/watch?v=-hhriqus2-8

The BCs are:

p
U
T
fvSchemes
fvSolution
Do you have any suggestions? Especially, for the settings in the fvSchemes and fvSolutions.
If you need more information, let me know.

Mihir, could you also post your fvSchemes? I'd be very interested. Maybe also you rBC setup for P U and T. Thanks.

Thanks for your replies, they helped. As I said, my case is very similar to the forewardStep case, just a slightly different mesh.
U: inlet fixedValue: (2 0 0), outlet: zeroGradient, wall/obstacle: fixedValue (0 0 0), top: symmetryPlane
p: inlet fixedValue: 1, outlet: zeroGradient, wall/obstacle: zeroGradient, top: symmetryPlane

The shock waves were a really good hint. In my case a pressure wave was reflected at the obstacle and travelled back to my inlet, where I had stated a constant pressure of 1. It happend that a region of very high pressure and very low pressure were at a small distance and thus I got very high velocity behind my inlet and that was, what caused the crash.
After I changed the p bc at the inlet to zeroGradient as well, the program runs stable, even when I change the velocity at the inlet to Mach smaller than 2.

So thanks. I understand the problem much more now.

sandrak,

what you discribe is a problem of rhosonicfoam. the solver is a "direct" solver for the decoupled euler-equations. if you use Ma 2 as speed the shockwaves all point into your domain, thats ok.
but the speed is not fast enough to transport reflected shockwaves. if you use Ma 3 it should work.

the real problem is, that the decoupled equations produce an error in U. If the speed is fast enough the solve of the mass equation (for rho) corrects this error.
if the speed is slower this error is not corrected. An solution for that gives the sonicFoam solver. This solver does an pressure-correction (PISO) to correct U and p. there you should be able to use all Ma speeds correct with a possible CoNum of ~1. (The BCs have to be set right)
For sonicFoam make shure that you do outer an inner corrections (fvSolution-PISOControls). If one of them is 0 you never solve the pressure equation and do no correction.