pressure/velocity stability problem

Submitted by evansmuts on Mon, 08/15/2011 - 16:59

Hi

I am having trouble with stability in my model and I was wondering if anyone can give me some idea as to what is causing the problem.

I have simple laminar flow between two concentric cylinders, initially starting from rest. A moving wall boundary on the inner wall starts moving the fluid, and eventually the whole domain gets "up to speed". The particles in the flow are initially evenly distributed, but as the flow speeds up, the particles start moving towards the outer wall.

After a reasonably short amount of time (in terms of the whole simulation time) the Courant number, which is usually quite low (around 0.04), suddenly explodes within 3-10 timesteps and the CFD solution crashes.

I suspect it is due to the particles bunching at the outer wall as the particle concentration goes up. However, the lowest void fraction is around 0.9994, which means that the particle concentration is still very low in these regions. Reducing the number of particles allows the simulation to run for longer before crashing, but I am still far away from the number of particles I want to model.

I have also noticed that during the PISO loops, the pressure loops are taking MUCH longer to converge than the velocity ones. I have used the same settings in the fvSchemes and fvSolution files as the tutorial examples that came with the CFDEM coupling. Occasionally I have changed the number of nCorrectors and nOrthogonalCorrectors and the relTol of the pressure scheme to see their effect, but nothing significantly affected stability.

The error message when OpenFOAM crashes indicates that the crash happens in either the PBiCG solver or the DILU preconditioner (depending on the run). Both these setting are for the U field and not the pressure. Or does this just mean that the velocity field is crashing due to a bad pressure field?

I should mention that I am still using version 1 of the coupling software, version 1.2.7 of LIGGGHTS and OpenFOAM version 1.6. I have developed some code that will require some effort to transfer to the new versions (OpenFOAM and LIGGGHTS), so I want to see if I can get my model to work with these older versions first.

Any comments would be appreciated.

Thanks

Evan

Forums:

cgoniva's picture

cgoniva | Mon, 08/15/2011 - 18:36

Hi Evan,

It's nice to hear from you in this forum!

you could try a multi-stage approach to tackle the problem:
- does the simulation converge if you run a pure CFD simulation. (otherwise BCs might be wrong)
- does the simulation converge if you neglect the volume of the particles (e.g. set min voidfraction to 1, or change the eqns. to be solved - otherwise momentum exchange might cause the trouble)
- what happens if you limit the momentum exchange terms?
- is the coupling interfal fine enough so that the exchange fields do not vary too much between two coupling steps?
- how big are the particles compared to the computational cells (you checked that allready)
- does it help to use smaller couplinga and fluid time steps?
- are the DEM time steps small enough? particle velocity after wall contact might be un-physically high?

Hope that these hints can help you a little further.
btw it might be good to set up a routine which allows you to "merge" you code with the current liggghts and cfdem versions (not necessarily a "GIT merge") to profit from the development and bug-fixes.

Kind regards,
Chris

evansmuts | Tue, 08/16/2011 - 16:45

Hi Christoph

Yes, it is great to finally be working on a coupled model. After thinking about what you said, I have some more questions/comments.

- the pure CFD model works perfectly and converges easily. The DEM code also works well (I tested both the DEM and CFD codes independently).

- what do you mean by "set min voidfraction to 1"? Where can I do this?
- what do you mean by "change the eqns. to be solved "? Would I have to edit the equations in the source code to not include momentum transfer?
- "neglect the volume of the particles" I think granular particles in LIGGGHTS need to have volume. Or am I mis-understanding you?
- "limit momentum exchange terms" Do you mean limit the information (i.e. number of particles) that is transfered each timestep?

- I usually run with a CFD timestep of 0.001s and DEM step of 0.00001s. I tried running the simulation with a smaller CFD timestep (0.0001), but that did not seem to make a difference. I deleted the results by mistake, so I will re-run and check it again. I will also try with a smaller DEM timestep.

- I have now also run the coupled simulation with no particles (actually 1, which got deleted after the first timestep) to see if the coupled simulation can re-create my pure CFD results. The simulation ran as expected for a while (third of the way), but then simulation started giving strange behaviour in the pressure field. Pockets of high/low pressure appeared and stayed until the end of the simulation. The pockets appeared to be periodic versions of each other (i.e. they looked identical to one another). I will investigate this further as the result should be the same as the pure CFD case, or is it "legal" to run the coupled simulation without any particles (i.e. zero momentum transfer)?

Thanks
Evan

cgoniva's picture

cgoniva | Tue, 08/16/2011 - 19:06

Hi Evan

- what do you mean by "set min voidfraction to 1"? Where can I do this?

in CFD/constant/couplingProperties you have the sub-dictionary for the voidfraction model, which wants "alphaMin" as an input parameter - set it to 1 and voidfraction should be ignored.

- what do you mean by "change the eqns. to be solved "? Would I have to edit the equations in the source code to not include momentum transfer?

yes

- "neglect the volume of the particles" I think granular particles in LIGGGHTS need to have volume. Or am I mis-understanding you?

above measures will ignore the particles volume in the CFD mass balance

- "limit momentum exchange terms" Do you mean limit the information (i.e. number of particles) that is transfered each timestep?

no, limit the value of Ksl (momentum exchange field) respectively to a maximum value. you can do so by adding KslLimit 10; (or other value) to the momentum coupling subdict in CFD/constant/couplinProperties (see constructor of implicitCouple.C) This will ignore the momentum (due to drag etc.) transferred from particles to fluid, but not the other way round.

- I usually run with a CFD timestep of 0.001s and DEM step of 0.00001s. I tried running the simulation with a smaller CFD timestep (0.0001), but that did not seem to make a difference. I deleted the results by mistake, so I will re-run and check it again. I will also try with a smaller DEM timestep.

The smaller the TS and the closer the coupling the more stable is the simulation (in general).

- I have now also run the coupled simulation with no particles (actually 1, which got deleted after the first timestep) to see if the coupled simulation can re-create my pure CFD results. The simulation ran as expected for a while (third of the way), but then simulation started giving strange behaviour in the pressure field. Pockets of high/low pressure appeared and stayed until the end of the simulation. The pockets appeared to be periodic versions of each other (i.e. they looked identical to one another). I will investigate this further as the result should be the same as the pure CFD case, or is it "legal" to run the coupled simulation without any particles (i.e. zero momentum transfer)?

when the exchange fields voidfraction and Ksl approach 1 and 0 the results should be the same as for pure CFD as the eqns should reduce to transient incompressible NSE. To be honest I did not investigate that in detail. Please keep me posted on that topic.

Cheers,
Chris

evansmuts | Fri, 08/19/2011 - 17:14

Hi Chris

Thanks for that detailed feedback. I am still trying a few things out, but here are some observations in the mean time:

- I changed alphaMin to 1, but I did not get an improvement in stability. In fact it crashed earlier than before. The pressure values were still acting as if there were particles in the flow, yet the voidFraction stayed at a value of 1. So I think I will have to edit out the momentum exchange terms to model the effect of removing the particles.

- I don't think I can change Ksl as version 1 of the coupling still uses explicit coupling, not implicit like version 2. The explicit coupling doesn't seem to take Ksl into account. Changing the KslLimit value to 1 did nothing, which confirms this.

- A smaller DEM timestep does improve stability, but not by much.

- I have run a sim with the small DEM timstep (smaller by factor of 10), and a smaller CFD timestep (smaller by factor of 2). It was actually more unstable. This could be due to the increase in DEM timsteps between updates. It kept giving me an error when I reduced my CFD timestep to smaller than 0.0005s. I will look into this.

- I have noticed that I lose a significant amount of particles (from 14000 to 100) during the timesteps that the Courant number is "blowing up". I think this is more a result of the high Courant number, and not a cause of the problem.

Cheers,
Evan

cgoniva's picture

cgoniva | Sun, 08/21/2011 - 22:26

Hi Evan,

you can of course also limit f (explicit momentum exchange) but there is no predefined handle for that - just go into the code of expMomSource.C and insert a limit.

you might lose particles as the flow solution diverges and then they are simply blown away...

some problems need really small timesteps...

good luck,
chris

evansmuts | Fri, 09/16/2011 - 10:50

Hi Christoph

I wanted to give some feedback on my stability problem.

After trying a number of things, I found that just reducing the timestep for both DEM and CFD did not help with stability. A combination of the smaller timesteps and a slightly smaller domain (narrower gap between concentric cylinders) helped with stability. While the bigger domain was fine, in hindsight, the smaller domain gave "better" CFD results, so I should have stuck with that geometry anyway. The resulting higher particle loading with the smaller domain has not badly affected the stability.

I still need to do some testing with the latest version of the coupling software, which should hopefully improve stability (implicit vs. explicit schemes). I will let you know how that goes.

Cheers,
Evan

vinaym | Mon, 03/19/2018 - 10:14

Dear Chris,

I am facing stability issues similar to Evan. I have fluidized bed setup for which a few simulations are unstable. while for most inlet flow rates they run fine.
How do i make a good estimate of KslLimit? Like the maximum value?

Thanks and kind regards,

Vinay