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

TwoPhaseEulerFoam convergence problems

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

Reply
 
LinkBack Thread Tools Display Modes
Old   August 25, 2005, 17:13
Default I'm trying to simulate a bubbl
  #1
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
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 Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   August 26, 2005, 09:21
Default I tried with another test case
  #2
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
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
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   August 26, 2005, 11:01
Default Here's what I've learned...
  #3
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21
niklas will become famous soon enoughniklas will become famous soon enough
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
niklas is offline   Reply With Quote

Old   August 26, 2005, 12:14
Default Thanks Niklas. Actually my in
  #4
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
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
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   August 26, 2005, 12:33
Default Hi, Starting with a void fr
  #5
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21
niklas will become famous soon enoughniklas will become famous soon enough
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 is offline   Reply With Quote

Old   August 26, 2005, 13:24
Default May I just add that this code
  #6
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21
niklas will become famous soon enoughniklas will become famous soon enough
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
niklas is offline   Reply With Quote

Old   August 29, 2005, 02:51
Default I checked the kinetic theory a
  #7
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
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 Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 3, 2005, 03:23
Default Hello, why in kineticTheoryMo
  #8
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
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
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 4, 2005, 03:29
Default Hi, mistake, should be pa
  #9
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21
niklas will become famous soon enoughniklas will become famous soon enough
Hi,

mistake, should be
pa_ = PsCoeff*Theta;

N
niklas is offline   Reply With Quote

Old   September 18, 2005, 10:19
Default Hello, I completed my tests o
  #10
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
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
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 19, 2005, 03:01
Default Hi, The radial distribution
  #11
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21
niklas will become famous soon enoughniklas will become famous soon enough
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
niklas is offline   Reply With Quote

Old   September 19, 2005, 17:09
Default Sure, I'll try. However th
  #12
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
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
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 20, 2005, 02:26
Default I'm guessing the problem is th
  #13
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21
niklas will become famous soon enoughniklas will become famous soon enough
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
niklas is offline   Reply With Quote

Old   September 23, 2005, 11:45
Default Hello, I tried to disable pf
  #14
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
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
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 23, 2005, 12:13
Default OK, so there's a division by z
  #15
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21
niklas will become famous soon enoughniklas will become famous soon enough
OK, so there's a division by zero somewhere (probably).
What does the line above the Foam::divide() say?
niklas is offline   Reply With Quote

Old   September 23, 2005, 12:59
Default If you mean in the code, I don
  #16
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
If you mean in the code, I don't know how to see it.

Alberto
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 23, 2005, 13:02
Default Type 'where' in gdb. Have a se
  #17
Super Moderator
 
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,416
Rep Power: 17
mattijs is on a distinguished road
Type 'where' in gdb. Have a search on this board for e.g. FOAM_DEBUG.
mattijs is offline   Reply With Quote

Old   September 23, 2005, 13:20
Default Thanks Mattijs. #0 0x405eb
  #18
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
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 Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 23, 2005, 18:17
Default I replaced the Gidspow formula
  #19
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
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
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 25, 2005, 02:31
Default You say that you dont have a s
  #20
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21
niklas will become famous soon enoughniklas will become famous soon enough
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
niklas 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
Some problems with twoPhaseEulerFoam su_junwei OpenFOAM Running, Solving & CFD 2 November 2, 2012 02:12
Convergence problems!! Elleana FLUENT 9 June 10, 2008 04:39
Convergence problems - please help M Liddell FLUENT 3 February 8, 2005 20:06
convergence problems jeremy FLUENT 7 May 30, 2002 06:41
Convergence Problems Prateep Chatterjee FLUENT 7 October 9, 2001 09:29


All times are GMT -4. The time now is 06:37.