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

Sudden crash caused by k-epsilon

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

Like Tree3Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   November 9, 2012, 05:59
Default Sudden crash caused by k-epsilon
  #1
Member
 
Robert
Join Date: Aug 2012
Location: Berlin
Posts: 74
Rep Power: 4
vainilreb is on a distinguished road
Hi,

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.

Robert
vainilreb is offline   Reply With Quote

Old   November 9, 2012, 06:27
Default
  #2
Senior Member
 
chegdan's Avatar
 
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 543
Rep Power: 18
chegdan will become famous soon enough
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.
__________________
Dan

Find me on twitter @dancombest and LinkedIn
chegdan is offline   Reply With Quote

Old   November 9, 2012, 06:45
Default
  #3
Member
 
Robert
Join Date: Aug 2012
Location: Berlin
Posts: 74
Rep Power: 4
vainilreb is on a distinguished road
Hey,

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?
Attached Files
File Type: zip system.zip (4.4 KB, 25 views)
File Type: zip 0.zip (11.3 KB, 9 views)
vainilreb is offline   Reply With Quote

Old   November 9, 2012, 07:06
Default
  #4
Senior Member
 
chegdan's Avatar
 
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 543
Rep Power: 18
chegdan will become famous soon enough
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

About kEpsilon turbulence model


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
Code:
div(phi,k)    Gauss linearUpwind grad(k);
and for gradient

Code:
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)....so be aware. You could also try limitedLinear for your divergence as some have noted that it is less diffuse than linearUpwind in certain cases.
s.m likes this.
__________________
Dan

Find me on twitter @dancombest and LinkedIn

Last edited by chegdan; November 9, 2012 at 09:40. Reason: Gauss
chegdan is offline   Reply With Quote

Old   November 9, 2012, 08:12
Default
  #5
Member
 
Robert
Join Date: Aug 2012
Location: Berlin
Posts: 74
Rep Power: 4
vainilreb is on a distinguished road
Hi,

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

Old   November 12, 2012, 03:17
Default
  #6
Member
 
Robert
Join Date: Aug 2012
Location: Berlin
Posts: 74
Rep Power: 4
vainilreb is on a distinguished road
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.
vainilreb is offline   Reply With Quote

Old   November 12, 2012, 06:04
Default
  #7
Senior Member
 
chegdan's Avatar
 
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 543
Rep Power: 18
chegdan will become famous soon enough
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.
s.m likes this.
__________________
Dan

Find me on twitter @dancombest and LinkedIn
chegdan is offline   Reply With Quote

Old   November 12, 2012, 06:33
Default
  #8
Member
 
Robert
Join Date: Aug 2012
Location: Berlin
Posts: 74
Rep Power: 4
vainilreb is on a distinguished road
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...?

Last edited by vainilreb; November 12, 2012 at 07:52.
vainilreb is offline   Reply With Quote

Old   November 12, 2012, 07:50
Default
  #9
Senior Member
 
chegdan's Avatar
 
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 543
Rep Power: 18
chegdan will become famous soon enough
What version of OpenFOAM are you using?
__________________
Dan

Find me on twitter @dancombest and LinkedIn
chegdan is offline   Reply With Quote

Old   November 12, 2012, 07:55
Default
  #10
Member
 
Robert
Join Date: Aug 2012
Location: Berlin
Posts: 74
Rep Power: 4
vainilreb is on a distinguished road
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 is offline   Reply With Quote

Old   November 12, 2012, 08:05
Default
  #11
Senior Member
 
chegdan's Avatar
 
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 543
Rep Power: 18
chegdan will become famous soon enough
Well, for practical cases...you 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.

Quote:
Originally Posted by vainilreb View Post
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.
__________________
Dan

Find me on twitter @dancombest and LinkedIn
chegdan is offline   Reply With Quote

Old   November 15, 2012, 03:47
Default
  #12
Member
 
Robert
Join Date: Aug 2012
Location: Berlin
Posts: 74
Rep Power: 4
vainilreb is on a distinguished road
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... :/
vainilreb is offline   Reply With Quote

Old   November 15, 2012, 12:14
Default
  #13
Senior Member
 
chegdan's Avatar
 
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 543
Rep Power: 18
chegdan will become famous soon enough
If you can provide a STL file of your surfaces, why not give snappyHexMesh a try? I personally have never used gmsh...so i can't help you there.
__________________
Dan

Find me on twitter @dancombest and LinkedIn
chegdan is offline   Reply With Quote

Old   November 19, 2012, 03:05
Default
  #14
Member
 
Robert
Join Date: Aug 2012
Location: Berlin
Posts: 74
Rep Power: 4
vainilreb is on a distinguished road
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 is offline   Reply With Quote

Old   November 20, 2012, 12:23
Default
  #15
Member
 
Robert
Join Date: Aug 2012
Location: Berlin
Posts: 74
Rep Power: 4
vainilreb is on a distinguished road
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.
vainilreb is offline   Reply With Quote

Old   December 25, 2012, 23:35
Default
  #16
Member
 
Haomin Yuan
Join Date: Jan 2012
Location: Madison, Wisconsin, USA
Posts: 54
Rep Power: 5
yhaomin2007 is on a distinguished road
I think you can change the relaxation factor for k and epsilon equation
yhaomin2007 is offline   Reply With Quote

Old   August 20, 2013, 07:56
Default
  #17
s.m
Senior Member
 
saeideh mohamadi
Join Date: Aug 2012
Posts: 229
Rep Power: 5
s.m is on a distinguished road
Quote:
Originally Posted by yhaomin2007 View Post
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
s.m is offline   Reply With Quote

Old   August 20, 2013, 10:45
Default
  #18
Member
 
Haomin Yuan
Join Date: Jan 2012
Location: Madison, Wisconsin, USA
Posts: 54
Rep Power: 5
yhaomin2007 is on a distinguished road
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.
yhaomin2007 is offline   Reply With Quote

Old   August 20, 2013, 11:06
Default
  #19
s.m
Senior Member
 
saeideh mohamadi
Join Date: Aug 2012
Posts: 229
Rep Power: 5
s.m is on a distinguished road
Quote:
Originally Posted by yhaomin2007 View Post
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, but what is stability?, what is up when a simulation don't have stability and when has it?
s.m is offline   Reply With Quote

Old   August 20, 2013, 11:10
Default
  #20
Member
 
Haomin Yuan
Join Date: Jan 2012
Location: Madison, Wisconsin, USA
Posts: 54
Rep Power: 5
yhaomin2007 is on a distinguished road
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.
yhaomin2007 is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
SimpleFoam k and epsilon bounded nedved OpenFOAM Running, Solving & CFD 13 November 4, 2013 15:13
epsilon and K blowing up. sivakumar OpenFOAM Running, Solving & CFD 1 October 25, 2012 04:50
Orifice Plate with a fully developed flow - Problems with convergence jonmec OpenFOAM Running, Solving & CFD 3 July 28, 2011 05:24
SimpleFoam k and epsilon bounded nedved OpenFOAM Running, Solving & CFD 1 November 25, 2008 21:21
CFX Solver : Sudden crash Hervé CFX 2 June 16, 2008 06:40


All times are GMT -4. The time now is 21:13.