
[Sponsors] 
August 25, 2005, 17:13 
I'm trying to simulate a bubbl

#1 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36 
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 (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

August 26, 2005, 09:21 
I tried with another test case

#2 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36 
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 (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

August 26, 2005, 11:01 
Here's what I've learned...

#3 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 29 
Here's what I've learned...
You have to break down the 'symmetry' very fast to avoid overpacking, either by a small asymmetry in inlet velocity or by small differences in initial voidfraction. 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 particleparticle interaction force treats the packing implictly, while the kinetic theory treats is explicit, and is much more stable. N 

August 26, 2005, 12:14 
Thanks Niklas.
Actually my in

#4 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36 
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 particlepartilce 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 (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

August 26, 2005, 12:33 
Hi,
Starting with a void fr

#5 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 29 
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 particleparticle force so yes, set g0 to a value greater than zero to use it and disable it by setting g0 to zero. N 

August 26, 2005, 13:24 
May I just add that this code

#6 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 29 
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 

August 29, 2005, 02:51 
I checked the kinetic theory a

#7 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36 
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 (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 3, 2005, 03:23 
Hello,
why in kineticTheoryMo

#8 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36 
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 (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 4, 2005, 03:29 
Hi,
mistake, should be
pa

#9 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 29 
Hi,
mistake, should be pa_ = PsCoeff*Theta; N 

September 18, 2005, 10:19 
Hello,
I completed my tests o

#10 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36 
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 underrelaxed 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 (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 19, 2005, 03:01 
Hi,
The radial distribution

#11 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 29 
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 

September 19, 2005, 17:09 
Sure, I'll try.
However th

#12 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36 
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 (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 20, 2005, 02:26 
I'm guessing the problem is th

#13 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 29 
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 

September 23, 2005, 11:45 
Hello,
I tried to disable pf

#14 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36 
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/OpenFOAM1.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 (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 23, 2005, 12:13 
OK, so there's a division by z

#15 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 29 
OK, so there's a division by zero somewhere (probably).
What does the line above the Foam::divide() say? 

September 23, 2005, 12:59 
If you mean in the code, I don

#16 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36 
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 (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 23, 2005, 13:02 
Type 'where' in gdb. Have a se

#17 
Senior Member
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,419
Rep Power: 26 
Type 'where' in gdb. Have a search on this board for e.g. FOAM_DEBUG.


September 23, 2005, 13:20 
Thanks Mattijs.
#0 0x405eb

#18 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36 
Thanks Mattijs.
#0 0x405eb246 in Foam::divide () from /home/alberto/OpenFOAM/OpenFOAM1.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/alberto1.2/lib/linuxGcc4Opt/libkineticTheoryModel.so #4 0x40285b62 in Foam::kineticTheoryModel::solve () from /home/alberto/OpenFOAM/alberto1.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 (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 23, 2005, 18:17 
I replaced the Gidspow formula

#19 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36 
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 (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 25, 2005, 02:31 
You say that you dont have a s

#20 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 29 
You say that you dont have a singularity,
but what happens if you use max(1.0alpa, 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 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Some problems with twoPhaseEulerFoam  su_junwei  OpenFOAM Running, Solving & CFD  2  November 2, 2012 01:12 
Convergence problems!!  Elleana  FLUENT  9  June 10, 2008 04:39 
Convergence problems  please help  M Liddell  FLUENT  3  February 8, 2005 19:06 
convergence problems  jeremy  FLUENT  7  May 30, 2002 06:41 
Convergence Problems  Prateep Chatterjee  FLUENT  7  October 9, 2001 09:29 