
[Sponsors] 
January 31, 2007, 05:25 
Hi, Rasmus
Thank you for your

#41 
New Member
Guanghao Wu
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 15
Rep Power: 8 
Hi, Rasmus
Thank you for your quick reply. In my case, time step=1e4, cell size=1e3m order. As you said, I should take shorter time step when I use the particleparticle model. In alphaEqn.H, // Correct the particleparticle force magnitude ppMagf = g0*rUaAf*min(exp(preAlphaExp*(alphaf  alphaMax)), expMax); So, I think g0, expMax, etc. are to calculate ppMagf which play a diffusion effect like viscosity. Right? As to g0, is there a recommendable value of g0? or g0 is case by case? Best regards, Guanghao 

January 31, 2007, 05:42 
That is correct, the particle

#42 
Senior Member
Rasmus Hemph
Join Date: Mar 2009
Location: Sweden
Posts: 108
Rep Power: 8 
That is correct, the particle collision pressure in twoPhaseEuelerFOAM is implemented, in the classical model, as a diffusive term in the alpha equation. If you have access to, for instance, Enwald et al, Int J. Multiphase Flow, vol22, 1996, there is a review over three of these models, with experimentally found values. They all have a g0 of 1 (with a unit of Pa).
Take a look at thread http://www.cfdonline.com/OpenFOAM_D...tml?1169210617 for more discussions regarding this term. Cheers, Rasmus 

January 31, 2007, 06:26 
Thank you, Rasmus.
Your comm

#43 
New Member
Guanghao Wu
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 15
Rep Power: 8 
Thank you, Rasmus.
Your comment and recommended materials are very helpful for me. Thank you again. 

January 31, 2007, 17:56 
Hello Guanghao,
the problems

#44 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Hello Guanghao,
the problems you noticed are due to the normal stress modulus (ppMagf) or, if you use the kinetic theory, to the granular pressure gradient in the momentum equation of the particle phase. When you're close to the packing limit, these quantities becomes suddenly big and this causes divergence. Of course the problem can be limited, but this requires quite a deep change to the solver. For example FLUENT treats the particle phase as incompressible when the packing limit is reached, and solves a pressure correction equation for it in the packed cells. MFIX uses a volume fraction correction equation instead of directly solving for the volume fraction of the dispersed phase, and then underrelax the volume fraction where particles are close to the packing limit. You can see the details in the MFIX manual. The implementation of one of these solutions (based on the SIMPLE algorithm) would help a lot. Btw, sooner or later MFIX should be rewritten using OpenFOAM, according to what reported by Hrvoje. But on the official site of MFIX there's no mention to that. Has anyone some information on this? Regards, Alberto 

February 6, 2007, 06:46 
Thank you Rasmus and Alberto.

#45 
New Member
Guanghao Wu
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 15
Rep Power: 8 
Thank you Rasmus and Alberto.
After many test calculations, I found the diameter of particle is a sensitive factor for the void fraction. Whene the diameter of particle is too big, the packing limit is reached very quickly. Now, I set the diameter a small value, then somehow, the calculation runs smoothly without divergence. And there is another question. In UEqns.H + (fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca) I think it should be as follows. + (fvc::grad(alpha)/( alpha + scalar(0.001)) & Rca) Why take the average of alpha? Just to keep the conservation? Btw, MFIX documents is very helpful to me, thank you Alberto. Best Regards, Guanghao 

February 6, 2007, 10:21 
I've an additional question to

#46 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
I've an additional question too about the alpha equation. If the particleparticle force model is enabled (g0 != 0), an additional term is added to the phase continuity equation to implicitly make it sensitive to the particle phase normal stress modulus.
If the kinetic theory is enabled, nothing equivalent is done. Why? The derivative of the granular pressure against the volume fraction is the normal stress modulus of the phase, so a similar implementation should be possible. Am I right? Regards, Alberto 

March 1, 2007, 09:50 
Hi Alberto. I looked (briefly)

#47 
Senior Member
Rasmus Hemph
Join Date: Mar 2009
Location: Sweden
Posts: 108
Rep Power: 8 
Hi Alberto. I looked (briefly) at the implementation of the Pfterm kinetic theory model. I believe the approach you are suggesting is possible. A seemingly simple way is to just comment out the pfterm in the kineticTheory model and instead set g0 to 1 and use ppMagf for frictional pressure. This should make it more stable since the ppMagf term is used semiimplicitly in the alphaequation. Observe that you in that case would need to take special care to the frictional viscosity muf, since it normally gets it pressure from pf!
I have a slightly different question. If I increase the time step above 1e4, I get problems, not with particle pressure, but with checkerboarding of the alphafield. Are any of you seeing this behavior? It dissapears if the time step is lowered. This could indicate a problem with discretization or something similar. I think it might be relate to the use of slipboundary condition for alpha, but I am not certain. (the checkerboard field would appear when slip is used). 

March 4, 2007, 11:18 
Hello Rasmus,
the solution yo

#48 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Hello Rasmus,
the solution you proposed definetly works. I did the simulations I talked about in the bugreport and they're OK. The solution is more stable and I have no divergence problem. Moreover, as expected, the change to the algorithm doesn't seem to influence the results. I'm working right now on extending this to the kinetic theory. I didn't notice the checkerboarding of alpha using zeroGradient conditions (equivalent to slip, being alpha a scalar), but I'm using small time steps to get convergence (~ 1.0e5s). I'll do some experiment on this too. Btw, do you solve for alpha and beta, or for alpha only and then find beta = 1alpha? I found that the second approach seems to be more stable. Regards, Alberto 

March 5, 2007, 05:57 
Hi,
Glad to hear that the ch

#49 
Senior Member
Rasmus Hemph
Join Date: Mar 2009
Location: Sweden
Posts: 108
Rep Power: 8 
Hi,
Glad to hear that the changed interpolation improved stability in your case as well. I wrote a little illustration to the problem on the bugreport. It would be useful to clean up (remove) the Pf and mufpart of the kinetic theory and instead use it in the way ppMagf is used right now. I have made the particleparticle forces runTime selectable and removed muf out of kinetic, I could send it to you too look at if you wish. To the next point. I believe that the issue of checkerboarding is an important one. With the stability of the ppMagterm improved, it should be possible to increase the time step. Currently I need a Courant number of 0.001 (corresponding to a deltaT of 2e5s) to get a stable solution. This seems very low. I recently saw a Fluent simulation (with kinetic theory) and high volume fraction particles running stable at 7e4 s. To investigate the issue, I tried setting drag forces to zero, only accelerating the particles by gravity. This seemed to remove the checkerboard behavior. Perhaps is it related to the semi explicit treatment of drag, as discussed in Henrik Rusche's thesis? Regarding the boundary condition, I meant setting a slipboundary condition for Ua, not alpha, sorry about the unclearness! Thanks for the tip about soving the betaequation. I will try it. Cheers, Rasmus 

March 5, 2007, 13:57 
Hello,
I read your comment on

#50 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Hello,
I read your comment on the bug report. The same observation is perfectly suitable for bubbling beds with a gas superficial velocity close to the minimum fluidization one. About the frictional model, I agree. The frictional model is not necessary in the kineticTheoryModel class (the equation for Theta has to be solved with pa(Theta), without Pf). Moreover Pf and muf can be added easily somewhere else, making the code more versatile. Are you using Pf instead of ppMagf? It would be interesting to see how you implemented it. About the time step. I have to use the same time step size for diluite flows (dt = 1.0e5s) while on FLUENT a time step of 1.0e3 / 5.0e4s is possible. I think the obtimisation obtained by FLUENT is due to a combination of factors. They manage the particle pressure fully implicitly close to the packing limit (no details available, the only reference I found is the FLUENT 4 manual), and they manage the drag term through the partial elimination algorithm. Plus they use PCSIMPLE instead of PISO, with underrelaxation. I thought to the drag term management myself as a source of problems because I had some nonphysical results in diluite flows simulation in a vertical pipe (i.e. particles accumulating at the bottom, with too high volume fractions). Reducing the time step, the problem disappeared. Actually, when the drag is treated as a source term, the time step should be small enough to grant a small K in the K*Ur product. This might explain the issues we meet. See you, Alberto 

June 5, 2008, 09:34 
Hi,(again)
What did you end

#51 
Member
Juho Peltola
Join Date: Mar 2009
Location: Finland
Posts: 80
Rep Power: 8 
Hi,(again)
What did you end up doing with the pf and muf? Did you implement it somewhere else? Regards, Juho 

October 13, 2008, 00:35 
Hello every body,
I need to

#52 
Member
Danielle PRL
Join Date: Mar 2009
Posts: 42
Rep Power: 8 
Hello every body,
I need to understand twoPhaseEulerFoam solver, for this I count on your help !! 1. In the momentum Eq I don't found grad(p) term (neither in UaEqn nor in UbEqn) ? 2. In the momentum Eq I found an additional term (fvm::Sp(fvc::div(phiRb), Ub) in UbEqn). Why you added it ?? If I added it, my case diverge!! 3.In kepsilon file why you write G like: G= 2*nutb*(tgradUb() && dev(symm(tgradUb()))); I think that the correct form is: G = nutb*( (symm(tgradU())skew(tgradU())) && tgradU() );  Thanks for your help! 

October 13, 2008, 14:45 
Hi Danielle,
you can find t

#53 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Hi Danielle,
you can find the details in H. Rusche thesis, which can be downloaded from http://powerlab.fsb.hr/ped/kturbo/Op...chePhD2002.pdf Only a comment about your first question. What you read in the code is not the full momentum equation, but the momentum predictor. The effect of the pressure gradient is accounted for at a later stage. Regards, Alberto 

October 13, 2008, 20:46 
Thanks Alberto,
I don't found

#54 
Member
Danielle PRL
Join Date: Mar 2009
Posts: 42
Rep Power: 8 
Thanks Alberto,
I don't found any think in this thesis about the additional term (fvm::Sp(fvc::div(phiRb), Ub) ! For the first question, I will reread this thesis. I will added the MRF in this solver if you have an Idea? I added the Coriolis Force in the both equations UaEqn and UbEqn, but my problem is the correction flux ! 

October 15, 2008, 15:49 
Hi,
when the UaEqn is defin

#55 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Hi,
when the UaEqn is defined you have: UaEqn = ( (scalar(1) + Cvm*rhob*beta/rhoa)* ( fvm::ddt(Ua) + fvm::div(phia, Ua, "div(phia,Ua)")  fvm::Sp(fvc::div(phia), Ua) )  fvm::laplacian(nuEffa, Ua) + fvc::div(Rca) + fvm::div(phiRa, Ua, "div(phia,Ua)")  fvm::Sp(fvc::div(phiRa), Ua) + (fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca) == // g // Buoyancy term transfered to pequation  fvm::Sp(beta/rhoa*K, Ua) //+ beta/rhoa*K*Ub // Explicit drag transfered to pequation  beta/rhoa*(liftCoeff  Cvm*rhob*DDtUb) ); Pay attention that you need to calculate the _total_convective_flux_: phi_a^# = phi_a  nu_eff (S_f sn_grad(alpha)_f)/(alpha_f + small_number) with the second term of the sum and the gradient contained in it calculated on faces. If you look at the first part of the equation, you have: fvm::div(phia, Ua, "div(phia,Ua)")  fvm::Sp(fvc::div(phia), Ua) To have the complete equation, you still need: + fvm::div(phiRa, Ua, "div(phia,Ua)")  fvm::Sp(fvc::div(phiRa), Ua) which are considered separately because they are not affected by the virtual mass coefficient Cvm. Did this help? Let me know if you need further details. It is a bit tricky to write equations here. :) I have never coded MRF procedures. Regards, Alberto 

October 23, 2008, 20:48 
Hi,
I'm back for asking help

#56 
Member
Danielle PRL
Join Date: Mar 2009
Posts: 42
Rep Power: 8 
Hi,
I'm back for asking help. In order to calculate the drag coefficient we need the diameter db for the phase b. how can I give a diameter for a continuous phase? Thanks 

October 24, 2008, 02:34 
A phase is defined by its dens

#57 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 20 
A phase is defined by its density, viscosity and characteristic size.
dragPhase is used to define from which phase the size information is to be taken, ie. which one is to be used to calculate the drag. So if you have particles in air and phase a is the particles, set dragPhase to a and it will not use the size information from phase b. Hence you can set it to 0 or 1, it doesnt matter. only use dragPhase blended if both phases have discrete properties and the size in phase b is sensible. note that the dragModels assumes phase a to be the discrete one. 

October 24, 2008, 09:00 
Thanks Niklas,
In my case I h

#58 
Member
Danielle PRL
Join Date: Mar 2009
Posts: 42
Rep Power: 8 
Thanks Niklas,
In my case I have water like a continuous phase, the air is the dispersed one. In this case we can also use the bubbleFoam Solver. maybe we should gives a diametre for this phase (water) which strongly infuence the drag phase (dragPhase should be set to blended ). how can I give db? 

October 24, 2008, 16:56 
Dear Danielle,
what kind of

#59 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Dear Danielle,
what kind of system are you trying to simulate? A. 

October 25, 2008, 06:11 
Hi Danielle
You can set

#60 
Senior Member

Hi Danielle
You can set the diameter for both phase though the dict of transportProperties in the dir of constant The d in the subdict of phase(phasea or phaseb) is the diameter of that phase. Junwei 

Thread Tools  
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 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 