[vainilreb]

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

[chegdan]

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]

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?

[chegdan]

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

http://www.cfd-online.com/Forums/ope...tml#post193470


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

[vainilreb]

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]

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]

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]

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]

What version of OpenFOAM are you using?

[vainilreb]

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]

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.

Originally Posted by vainilreb

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]

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]

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.

[vainilreb]

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]

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]

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

[s.m]

Originally Posted by yhaomin2007

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]

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]

Originally Posted by yhaomin2007

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?

[yhaomin2007]

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.