[guanghaowu]

Hi, Rasmus
Thank you for your quick reply.
In my case, time step=1e-4, cell size=1e-3m order. As you said, I should take shorter time step when I use the particle-particle model.

In alphaEqn.H,
// Correct the particle-particle 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

[hemph]

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.cfd-online.com/OpenFOAM_D...tml?1169210617
for more discussions regarding this term.

Cheers,
Rasmus

[guanghaowu]

Thank you, Rasmus.
Your comment and recommended materials are very helpful for me. Thank you again.

[alberto]

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 under-relax 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

[guanghaowu]

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

[alberto]

I've an additional question too about the alpha equation. If the particle-particle 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

[hemph]

Hi Alberto. I looked (briefly) at the implementation of the Pf-term kinetic theory model. I believe the approach you are suggesting is possible. A seemingly simple way is to just comment out the pf-term 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 semi-implicitly in the alpha-equation. 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 1e-4, I get problems, not with particle pressure, but with checkerboarding of the alpha-field. 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 slip-boundary condition for alpha, but I am not certain. (the checkerboard field would appear when slip is used).

[alberto]

Hello Rasmus,
the solution you proposed definetly works. I did the simulations I talked about in the bug-report 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.0e-5s). I'll do some experiment on this too.

Btw, do you solve for alpha and beta, or for alpha only and then find beta = 1-alpha? I found that the second approach seems to be more stable.

Regards,
Alberto

[hemph]

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 bug-report. It would be useful to clean up (remove) the Pf and muf-part of the kinetic theory and instead use it in the way ppMagf is used right now. I have made the particle-particle 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 ppMag-term improved, it should be possible to increase the time step. Currently I need a Courant number of 0.001 (corresponding to a deltaT of 2e-5s) 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 7e-4 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 slip-boundary condition for Ua, not alpha, sorry about the unclearness!
Thanks for the tip about soving the beta-equation. I will try it.

Cheers,
Rasmus

[alberto]

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.0e-5s) while on FLUENT a time step of 1.0e-3 / 5.0e-4s 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 PC-SIMPLE instead of PISO, with under-relaxation.

I thought to the drag term management myself as a source of problems because I had some non-physical 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

[juho]

Hi,(again)

What did you end up doing with the pf and muf? Did you implement it somewhere else?

Regards,
Juho

[danielle]

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 k-epsilon 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!

[alberto]

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

[danielle]

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 !

[alberto]

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 p-equation
- fvm::Sp(beta/rhoa*K, Ua)
//+ beta/rhoa*K*Ub // Explicit drag transfered to p-equation
- 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

[danielle]

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

[niklas]

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.

[danielle]

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?

[alberto]

Dear Danielle,

what kind of system are you trying to simulate?

A.

[su_junwei]

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