CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   Sudden crash caused by k-epsilon (

vainilreb November 9, 2012 05:59

Sudden crash caused by k-epsilon

I'm running a calculation of ammonia flowing into a room.
The BC for this is inletFlowRateVelocity with a massflow of 4e-04 for the first 0.1 s, linearly decreasing to 0 during the next 59.9 s .

My problem is the sudden explosion of k-epsilon calculation after about 1.09 seconds. Everything works fine and all of a sudden "Bounding epsilon" occurs. A few steps later "Bounding k" follows and soon the simulation crashes due to floating point exception. During these few timesteps Epsilon and k reach maximum values of x*e20 and higher, starting from physically reasonable values.

I hope anyone can give me a hint if it's a mesh problem or more related to BC or perhaps solver settings?

Thanks a lot.


chegdan November 9, 2012 06:27

This could be a whole list of things, but I will start with a list of things that I would check:

0. what solver are you using, what are your settings, etc.? Others may be able to help if you provide this info.

1. run checkMesh to see if you have some bad cells in there. Search the forum if you receive multiple failed mesh checks
2. What are your discretization schemes? you might be using one that is unbounded. try using limitedLinear or try linearUpwind with a cellLimited leastSquares gradient.
3. Boundary conditions. Inlet, exit, walls, etc. e.g. If there is rotation at the exit or a zero defined at the wall, this may be the cause.
4. Is your initial condition nice? Look in the user's guide for how to initialize the k and epsilon fields.

hope this helps and good luck.

vainilreb November 9, 2012 06:45

2 Attachment(s)

thanks for the quick response.

I've attached my /0 and /system files

I used reactingFoam with combustion off and chemistry off. CheckMesh says my mesh is OK. Max Skewness is around 1, Av. non-orthogonality is 19, max. 65.

What puzzles me is that the calculation seems to run rockstable and then explodes in only a few time steps. I just tried to run it with upwind schemes for the div() terms and up to now it seems to work. The magical 1.09 seconds have already been passed.

But I've got some worries about numerical diffusion, which might strongly influence the solution, because there is no air flow besides convective flow. Is an upwind scheme as bad for my case as i imagine?

chegdan November 9, 2012 07:06

First of all, where is your outlet? This is probably why everything is blowing up.

For everything else:

For your inlet k, i would try turbulentIntensityKineticEnergyInlet BC
for epsilon, try turbulentMixingLengthDissipationRateInlet. Example located at

For discretization schemes. Linear upwind is as you noted, diffuse. linearUpwind is unbounded if you use an unbounded gradient scheme. You could use something like

div(phi,k)    Gauss linearUpwind grad(k);
and for gradient


grad(k) cellMDLimited leastSquares 1;
this will ensure a second order scheme. you can change the limiter from cellMDLimited -> cellLimited -> faceMDlimited->faceLimited from least diffuse to most diffuse limiters. However, you may see stagnation in your convergence using limiters (noted somewhere in Fluent documentation) be aware. You could also try limitedLinear for your divergence as some have noted that it is less diffuse than linearUpwind in certain cases.

vainilreb November 9, 2012 08:12


I really hope that the missing outlet is not the problem. I simulate this case to validate the simulation accuracy by comparing my CFD results to a similar experiment. It's not that much mass being added to the system so I thought a compressible solver should be able to handle the pressure increase?

I'll try to change the BC for k and epsilon, but I don't see why it should be the problem as my calculation won't blow up before 1.09 s? My time step is around 6e-05 at Co<0.3 so there's a lot of time steps being calculated without any problems before the crash?

EDIT: My OpenFOAM version doesn't accept the linearUpwind entry. Is "Gauss linearUpwind grad(k)" right?

Thanks again for your help!

vainilreb November 12, 2012 03:17

Hey again,

I started the simulation using "Gauss upwind grad(X); grad(X) cellMDLimited leastSquares 1;" - with X here being placeholder for U, Yi_h, h, K, p, epsilon and k.
So far it is running rockstable, but only 43 secs were calculated since friday. I wonder if the accuracy will come up to my expectations...

I'll keep you up to date.

chegdan November 12, 2012 06:04

Just note that since upwind does not use a gradient term in its discretization, adding grad(X) after the upwind statement does not do anything.

vainilreb November 12, 2012 06:33

Okay, I thought simply "upwind" was the same as "linearUpwind", since my OpenFOAM version strictly refuses to accept "linearUpwind".

So I'm afraid my simulation results won't be satisfying... :(

Which is the right FVSchemes entry for what I really wanted to implement?

EDIT: Okay, I was able to fix this myself. ;) Now all that's left, is the question whether using linearUpwind for all div() discretizations is really a good decision...?

chegdan November 12, 2012 07:50

What version of OpenFOAM are you using?

vainilreb November 12, 2012 07:55

It's okay, I already fixed it. Is it a good decision to use linearUpwind for all div() schemes?

Thanks a lot for your help, Chegdan.

chegdan November 12, 2012 08:05

Well, for practical could use upwind for fields like k and epsilon to see some better behavior (i.e. boundedness). I use linearUpwind or limitedLinear for most all fields. Its a matter of what works and what is stable. I tend to run simulations in stages (like many others), changing from lower order schemes to higher order to produce more realistic and iterative "renditions" of the solution...mapping the previous solution if its steady-state.

No problem, glad its working out.


Originally Posted by vainilreb (Post 391629)
It's okay, I already fixed it. Is it a good decision to use linearUpwind for all div() schemes?

Thanks a lot for your help, Chegdan.

vainilreb November 15, 2012 03:47

I varied my mesh settings a bit and last night my workstation managed to compute 60 secs without floating point exception or any other crash scenario. :)

But I need some further advice: My time step is 1e-5 s at Co<=0.3, because the velocity at the gas outlet is 5 m/s and the outlet diameter is only 6 mm. Thus the cell size at the outlet is very small. How can I get rid of this problem? I thought about making the mesh more coarse in direction of the gas flow and as fine as needed in the perpendicular directions, but I can't do it with gmsh... :/

chegdan November 15, 2012 12:14

If you can provide a STL file of your surfaces, why not give snappyHexMesh a try? I personally have never used i can't help you there.

vainilreb November 19, 2012 03:05

I've already tried to save my models as STL. But they are built with gmsh and it doesn't really work with this tool.

vainilreb November 20, 2012 12:23

If anyone is interested in how I finally solved my stability problem: I simply use kOmegaSST instead of kEpsilon. No more crashes since I changed the turbulence model. :)

yhaomin2007 December 25, 2012 23:35

I think you can change the relaxation factor for k and epsilon equation

s.m August 20, 2013 07:56


Originally Posted by yhaomin2007 (Post 399009)
I think you can change the relaxation factor for k and epsilon equation

Hi Haomin
how should we change the under-relaxation-factors to have a more stable and also accurate result?

e.g p=0.3, and the other ones(U, nuTilda,omega,k) 0.7, now, we should decrease their value of URF or increase them?
Thank you very much:)

yhaomin2007 August 20, 2013 10:45

I assume you are doing steady state problem. This URF should decrease to insure stability. But lower URF will require more steps to get to your final answer.

s.m August 20, 2013 11:06


Originally Posted by yhaomin2007 (Post 446912)
I assume you are doing steady state problem. This URF should decrease to insure stability. But lower URF will require more steps to get to your final answer.

Thank you. :) yes i am doing the steay state problem,
this question may very stupid to you,:o but what is stability?, what is up when a simulation don't have stability and when has it?

yhaomin2007 August 20, 2013 11:10

the simulation will blows up, produce some extreme large number(NAN), if the simulation is not stable.~~
You can just try some lower URF to see it the simulation do not stop in the middle.

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