[alberto]

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

[alberto]

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. :-(

Regards,
Alberto

[niklas]

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.

N

[alberto]

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,
Alberto

[niklas]

Hi,

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.

N

[niklas]

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.

N

[alberto]

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

Regards,
Alberto

[alberto]

Hello,
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);

PsCoeff.min(1.0e+10);
PsCoeff.max(-1.0e+10);

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

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

Best regards,
Alberto

[niklas]

Hi,

mistake, should be
pa_ = PsCoeff*Theta;

N

[alberto]

Hello,
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,
Alberto

[niklas]

Hi,

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

-DFULLDEBUG -g -O0

and then run it through gdb.

N

[alberto]

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.

AP

[niklas]

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.

N

[alberto]

Hello,
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/libOpenFOAM.so

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

Regards,
Alberto

[niklas]

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

[alberto]

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

Alberto

[mattijs]

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

[alberto]

Thanks Mattijs.

#0 0x405eb246 in Foam::divide () from /home/alberto/OpenFOAM/OpenFOAM-1.2/lib/linuxGcc4Opt/libOpenFOAM.so
#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/libkineticTheoryModel.so
#4 0x40285b62 in Foam::kineticTheoryModel::solve () from /home/alberto/OpenFOAM/alberto-1.2/lib/linuxGcc4Opt/libkineticTheoryModel.so
#5 0x08083b5e in main (argc=3, argv=0xbfffde34) at twoPhaseEulerFoam.C:95

The radial distribution function g0 seems to be the guilty.

Alberto

[alberto]

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?

Regards,
Alberto

[niklas]

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.

N