CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   TwoPhaseEulerFoam convergence problems (

alberto August 25, 2005 17:13

I'm trying to simulate a bubbl
I'm trying to simulate a bubbling bed with twoPhaseEulerFoam, to which I added the drag model of Syamlal, Rogers and OBrien.

The bed (from Boeder's paper) is a 0,57x1m square box filled up to 0.5m with solid particles at minimum fluidization. The gas phase is uniformally fed to the bottom (0,676 m/s). Both the maximum packing limit and the particle particle restitution coefficient are set to 0.6. Gidaspow's models are used for viscosity and granular conductivity. The granular temperature is calculated with the algebraic equation (equilibrium = on),
The simulation has to be performed for 2 seconds.

The problem is the following:

- Using limitedLinear 1 scheme for convection, I get a solution close to that given by FLUENT with second order upwind until 0.08s. After this time the solution becomes unstable.

- If I use the upwind scheme I obtain the same behaviour just with a worse form of the initial bubble.

I also noticed that the solid volume fraction can become greater than the maximum packing limit. This should be impossible.

What do you suggest?

Thanks in advance,

alberto August 26, 2005 09:21

I tried with another test case
I tried with another test case (the one used by Guenther and Syamlal to test numerical schemes): a bed with a central jet.
The problem is almost the same, but appears almost immediatly, just a few iterations from the beginning.
Underrelaxation just delays the problem. :-(


niklas August 26, 2005 11:01

Here's what I've learned...
Here's what I've learned...

You have to break down the 'symmetry' very fast to avoid over-packing, either by a small asymmetry in inlet velocity or by small differences in initial void-fraction.
anything else than upwind is asking for trouble.
As the flow develops you can switch to something else, but if you start with very symmetric and smooth fields, use upwind.

using a 'large' time step initially works better than reducing it.

play around with the number of PISO and alpha correction loops, it really makes a difference

...and last, use the particle particle interaction force to prevent too high packing values.
The particle-particle interaction force treats the packing implictly, while the kinetic theory treats is explicit, and is much more stable.


alberto August 26, 2005 12:14

Thanks Niklas. Actually my in
Thanks Niklas.
Actually my initial condition sets the solid volume fraction at the maximum packing limit, which may be the main cause of the problem. I'll try to change this a bit to see what happens.

I'm a bit concerned about the first order upwind scheme because it has been proven to give bubbles with a wrong shape, but maybe it's not that critical if used at the beginning of the calculation.

The particle-partilce interaction force is switched on if g0 is greater than zero in the ppProperties dictionary. Is this right?

Best regards,

niklas August 26, 2005 12:33

Hi, Starting with a void fr

Starting with a void fraction at maximum packing is not good. This is most likely your problem.
Try to modify the code that sets the initial conditions so that
alpha is something like
2.0*alpha0*(0.5 - epsilon*rnd(0 to 1))

there's no 'switch' for the particle-particle force
so yes, set g0 to a value greater than zero to use it and disable it by setting g0 to zero.


niklas August 26, 2005 13:24

May I just add that this code
May I just add that this code was released in an experimental state hoping that the users will improve it, the kinetic theory stuff is relatively untested and may contain bugs, (Henry found some fortunately and did a great job cleaning it out), but I thought it would be good to release it anyway so we dont have to invent the wheel twice.


alberto August 29, 2005 02:51

I checked the kinetic theory a
I checked the kinetic theory and drag models and I didn't find errors :-)

The only doubt is about the formulation of the frictional viscosity. If in kineticTheoryModel.C, line 228:

volScalarField muf = 0.5*pf*sin(phi_)*time1;

there should be the Schaeffer correlation for the frictional viscosity, then it should be divided by the square root of the second invariant of the solid strain rate tensor.

I added the Syamlal and the Gidaspow (Ergun + Wen&Yu) drag models. The code is on the wiki, under code snippets. Tell me if you need the files.
In my implementation of Gidaspow's drag probably there's a better workaround to return a tmp volScalarField from a volScalarField than multipling it by 1.0 on the return line :-)


alberto September 3, 2005 03:23

Hello, why in kineticTheoryMo
why in kineticTheoryModel.C, lines 217 - 224 pf is summed to pa_ twice, first by adding it to the PsCoeff and then in the pa_ calculation?

// add frictional stress for alpha > alphaMin
PsCoeff += pf/(Theta_+Tsmall);


// update particle pressure
pa_ = PsCoeff*Theta_ + pf;

Also, has someone simulated a bubble formation/eruption in a bubbling bed with twoPhaseEulerFoam?

Best regards,

niklas September 4, 2005 03:29

Hi, mistake, should be pa

mistake, should be
pa_ = PsCoeff*Theta;


alberto September 18, 2005 10:19

Hello, I completed my tests o
I completed my tests of the code on the bubbling bed. I mainly considered bubble formations and eruption starting from a packed or almost packed bed.
Unfortunately results are far from beeing good and far from what it's possible to obtain with FLUENT.

In my opinion the problems depend on the management of the packing limit.

If I use the particle particle interaction force by setting g0 != 0 in the ppProperties dictionary, I get convergence, but the solid volume fraction is wrong. The bubble has a wrong form and in some case it just breaks instead of closing on itself and going up.

If I set g0 = 0 in order to manage the packing limit with the kinetic theory approach, the result seems to be better, but it's impossible to complete a calculation due to convergence problems. After a few iteration I get a "nan" in the continuity equation residual. This should be related to the radial distribution function g0 which becomes infinite if the packing limit is reached.

I examined the implementation of the kinetic theory in MFIX. It adopts a different approach where the granular pressure is under-relaxed in the packed zones, but I have no idea if this is possible in OpenFoam.

In the kineticTheory.C file there's a comment about a linear extrapolation for g0. What does it mean?

Best regards,

niklas September 19, 2005 03:01

Hi, The radial distribution

The radial distribution is one of the things that
'has room for improvement'.
As you probably have seen, g0 is clipped when
alpha goes above maximum packing, so it shouldnt become infinite, it also behaves very bad for alpha > alphaMax. Instead of clipping it you can have it increase linearly with alpha when alpha goes above a certain value (close to alphaMax).
That is what I mean with linear extrapolation.

When you get NaN's could you please check where it fails by setting

setenv FOAM_SIGFPE 1
setenv FOAM_ABORT 1
and recompiling the code and libraries with these flags added to the Make/options file


and then run it through gdb.


alberto September 19, 2005 17:09

Sure, I'll try. However th
Sure, I'll try.

However the problem happens systematically when alpha becomes greater than alpha_max.
The solid volume fraction reaches alpha_max and at the following iteration nua is inf. The nan follows immediatly.


niklas September 20, 2005 02:26

I'm guessing the problem is th
I'm guessing the problem is the fricitional normal stress term pf.
You could try and set Fr to 0.0 in kineticTheoryProperties to disable it and/or

have a look in kineticTheoryModel.C

Instead of setting PsCoeff.max and min, above PsCoeff += ... set pf.max and min instead.


alberto September 23, 2005 11:45

Hello, I tried to disable pf
I tried to disable pf both setting Fr to zero and commenting it in the kineticTheoryModel.C source file, but I still get the nan problem. Also the limitation seems to have no effect.

I replaced the pf formulation with the one proposed by Syamlal (pf = A(alpha - alphaMin)^n, wih A = 10^25 and n = 10), but nothing changed.

gdb gives:

Program received signal SIGFPE, Arithmetic exception.
[Switching to Thread 1088750368 (LWP 13242)]
0x405eb246 in Foam::divide () from /home/alberto/OpenFOAM/OpenFOAM-1.2/lib/linuxGcc4Opt/

but I'm not familiar at all with this tool, so I don't know how to locate the corresponding code line.


niklas September 23, 2005 12:13

OK, so there's a division by z
OK, so there's a division by zero somewhere (probably).
What does the line above the Foam::divide() say?

alberto September 23, 2005 12:59

If you mean in the code, I don
If you mean in the code, I don't know how to see it.


mattijs September 23, 2005 13:02

Type 'where' in gdb. Have a se
Type 'where' in gdb. Have a search on this board for e.g. FOAM_DEBUG.

alberto September 23, 2005 13:20

Thanks Mattijs. #0 0x405eb
Thanks Mattijs.

#0 0x405eb246 in Foam::divide () from /home/alberto/OpenFOAM/OpenFOAM-1.2/lib/linuxGcc4Opt/
#1 0x0809e169 in Foam::divide<foam::fvpatchfield,> (result=@0x8315d90, ds=@0xbfff4e2c, gsf=@0x8315d90) at GeometricScalarField.C:123
#2 0x080ab9dc in Foam::operator/<foam::fvpatchfield,> (ds=@0xbfff4e2c, tgsf=@0xbfff4e70) at GeometricScalarField.C:183
#3 0x40285110 in Foam::kineticTheoryModel::g0 () from /home/alberto/OpenFOAM/alberto-1.2/lib/linuxGcc4Opt/
#4 0x40285b62 in Foam::kineticTheoryModel::solve () from /home/alberto/OpenFOAM/alberto-1.2/lib/linuxGcc4Opt/
#5 0x08083b5e in main (argc=3, argv=0xbfffde34) at twoPhaseEulerFoam.C:95

The radial distribution function g0 seems to be the guilty.


alberto September 23, 2005 18:17

I replaced the Gidspow formula
I replaced the Gidspow formulation of g0 with the Carnahan and Starling (from Van Wachem phd thesis) one:

g0 = 1.0/(1.0 - alpha) + 3.0*alpha/(2.0*pow(1.0 - alpha, 2)) + pow(alpha, 2)/(2.0*pow(1.0 - alpha,3));

but at alpha slightly greater than 0.65 (my alphaMax) I still have the problem (gdb keeps telling me it's in g0), even if the function has no singularity there.

I tried to make the solver write also the granular pressure and the viscosity by changing the NO_WRITE flag to AUTO_WRITE in kineticTheoryModel.C, but it doesn't work. Why?

Also, if I use the

Info << "Text"

in the code I get no output text.
I always do an ./Allwclean and then an ./Allwmake to recompile the code. Is it OK?


niklas September 25, 2005 02:31

You say that you dont have a s
You say that you dont have a singularity,
but what happens if you use max(1.0-alpa, SMALL) instead...

AUTO_WRITE should work and also Info...
Are you running the 'correct' code.
Are you on a NFS system, i.e you're compiling on
one computer and running it on another one.
Then you might have a timing problem and it can take a while before the other compuer sees the update.

Yes it's ok, but Allwclean shouldnt be necessary.


All times are GMT -4. The time now is 05:02.