With the help of Ram, a student in the group where I work, we started validating the code with basic test cases, using a commercial code and MFIX as reference.

Here there is just one picture of a bubbling bed with a central jet, and the typical first bubble.



As you can notice, the interface does not show wrinkles during the expansion, even in this case where the flow rate is pretty high (the bubbles degenerate in a slug after a while), which leads to a high particle pressure right below the interface, as you can see from the higher particle concentration there.

Additionally profiles across the interface are smooth (trust me... I don't know how to attach more than one picture!), the solution is stable until the end of the simulation (2s).

These results have been obtained using the kinetic theory model with algebraic equation for the granular temperature with Gidaspow's closures and drag, Schaeffer's frictional model.
The PEA algorithm is used for the drag term. Convective schemes are all linearUpwind, while linear has been used elsewhere.

More to come...

P.S. Before everyone asks. The code is *not* the one in my public git repository. It will be there soon, after tests are passed, with the test-cases as tutorials.

Other picture from one run of a uniform bubbling bed:

Hi Alberto,

These are interesting results. Is the maximum alpha reached the same one set in the dict files? I'm mostly interested how this code will behave in packing situations such as a pressure stagnation point for instance.

Ziad

Ziad Boutanios wrote:
Hi Alberto,

These are interesting results. Is the maximum alpha reached the same one set in the dict files? I'm mostly interested how this code will behave in packing situations such as a pressure stagnation point for instance.

Ziad

Hi Ziad,

thanks for your question. The volume fraction superior limit is enforced by the frictional model in the current implementation of the code.

In other words, there are two regions of volume fraction:

1) alpha < alphaMinFriction

2) alphaMinFriction < alpha < alphaMax

The particle pressure is computed from the kinetic theory closures in region 1, and it depends on the radial distribution function g0.

In region 2, the particle pressure given by the kinetic theory is used in the equation for the granular temperature, and in the determination of the kinetic stresses.
However the normal stress term that appears in the modified equation for alpha (ppMagf) is computed using the expression for frictional pressure given by Schaeffer's model. The reason of this choice is both numerical and physical:

- Numerically, the expression of the pressure given by this model is continuous, and, even if it grows quickly, it does not diverge.

- Physically, the kinetic theory closures are not really valid in the packed regions, since they assume binary collisions. As a consequence, it makes sense to use the frictional model in that range of volume fractions.

Notice that, in this way, you obtain a bounded value of the phase fraction, which can be used to safely compute the radial distribution g0 where you need it.

Additionally, if you do not want to switch to a frictional model, you can use exactly the same approach, by dealing with the discontinuity of g0 in some way (an expansion).

Now, back to your question. In the case set-up, you specify both

alphaMinFriction , set to 0.61 in Schaeffer's model
alphaMax, set to 0.63 in the results shown above

The value of alphaMax is used to compute the radial distribution function g0. However the actual maximum value reached by the dispersed phase fraction depends only on the equation of state that gives you the frictional pressure. What happens in the calculations is that you have alphaMax between 0.61 and 0.63.

In the case of a completely quiescent system, like a settled fluidised bed, the algorithm is able to enforce the packing limit at approximately the maximum value (it won't be exactly 0.63, but something around it, depending on the equation of state), as described above. In this case the force balance is quite delicate, and you need to use a relatively small time step to have a converged solution. We run some cases to test this behaviour. I will dig them out and show you what happens.

On this point, I have to say that in these years I tried a large variety of algorithms suggested in the literature as "the" solution to this problem. In particular, I spent a good amount of time working on the implementation of a pressure-based solver, where the particle pressure is treated as unknown, and the dispersed-phase continuity equation is transformed into an equation for the particle pressure.
These algorithms look very good on paper, but their efficiency and stability did not seem great to me. In particular, it is true that the particle pressure equation degenerates into an incompressibility condition, independent from the equation of state of the particle pressure, when alpha tends to alphaMax, but it is also true that the equation is undefined when alpha -> 0, and it is troublesome when alpha is small. Additionally, you have to deal with the non-linearity of the equation of state, and its inversion, which might introduce conservation errors.
So, in the end, after many attempts, I switched back to an equation for the volume fraction, sensitized to the particle pressure. This does not have limitations for alpha -> 0, but requires some care when alpha -> alphaMax (read: a continuous function for the particle pressure), which I solved following what is done in MFIX, as explained above.

If you have suggestions on this topic, I am very open to them

Well thank you for a nicely detailed reply. I guess you answered my question and any other ones I might have

So far I've been using conventional 2-fluid formulations for bubble and droplet systems. Now that I work on flows with ice particles I thought it could prove more appropriate to move into granular dynamics. Still fairly new to the subject though so I don't really have any suggestions so far.

Playing with twoPhaseEulerFoam I noticed some "issues", such as very small negative alpha for instance (replace with log(alpha) in the solver to enforce positivity?), and unrealistic alpha values in packed situations.

I lately got MFIX and it's nice to see all the documentation that comes with it. That can only help!

Looking forward to try your solver...

Ziad Boutanios wrote:

Playing with twoPhaseEulerFoam I noticed some "issues", such as very small negative alpha for instance (replace with log(alpha) in the solver to enforce positivity?), and unrealistic alpha values in packed situations.


Yes... "some" . These are two separate issues that appear in multi-fluid models.

The small negative values of alpha are often coming out from linear solvers (evidence is the fact that even with the upwind scheme, you still have them). If you dig into the MFIX code, you notice that there are are good set of "ADJUST" functions, whose role is exactly to take care of issues like these for the phase fraction, the granular temperature and so on. As long as those values are small (about machine precision), they should not give problems.

The high values of phase fraction in the packed regions have been a major pain caused by twoPhaseEulerFoam. The cause is mixed: on one hand the PISO algorithm is not suitable to deal with the strong non-linearity of the particle pressure term (you need under-relaxation), on the other hand the particle pressure term is treated fully explicitly if the kinetic theory is used. This is well known to give problems (MFIX does the same for n > 2 phases, but they have an aggressive adaption of the time step and a solution-recovery mechanism that keep things going).

Best,
Alberto