Alberto Passalacqua wrote:

A first version of twoPhaseEulerPimpleFoam is available here: git://github.com/AlbertoPa/twoPhaseEulerPimpleFoam.git


Moving the twoPhaseEulerPimpleFoam discussion from the "wall" to this discussion forum:

Alberto Passalacqua wrote:

1) I store alpha at the beginning of the outer corrector loop, not at the beginning of the time step.


What I meant was that why do you store the alpha.storePrevIter() at all, when the alpha is relaxed with alphaEqn.relax() instead of alpha.relax()? I've thought the storePrevIter() is only needed when the field (not the matrix) is relaxed.


Alberto Passalacqua wrote:

2) It seems to work fine in many cases, as long as your time step is not extremely large. What has a very important role in stabilizing the computation is the alpha-U-p coupling: if you ensure simultaneous convergence of these equations, you are in a good position to have a robust algorithm. That's why I feel PIMPLE (or unsteady SIMPLE) are a better alternative to PISO when the particle pressure term is introduced, since it really changes too quickly to have an acceptable time-step and, at the same time, a strict enforcement of the maximum packing limit.

Notice that in the released code nothing has been changed in the kinetic theory / frictional stress models. To further stabilize the solution, I did something similar to what done in MFIX (and if I remember correctly, in your thesis too): it is possible to either write the radial distribution function g0 as a power series around the packing limit, or replacing it with a continuous expression for the particle pressure taken from the frictional models.


I did test different continuous formulations for the frictional pressure and g0, but they didn't solve the packing limit problems with the PISO + underrelaxed alphaEqn formulation I used. In the final version with the frictional pressure flux limiter, the standard Johnson & Jackson formulation was used with less strict limiters.

The standard limiters in the twoPhaseEulerFoam are very strict: alphaMax-0.01 for the g0 and alphaMax-0.05 for the denominator of the Johnson & Jackson frictional pressure model. With commonly used values for alphaMinFriction (~0.6) and alphaMax (~0.63) the whole frictional stress model is drastically affected by the limiter.

JohnsonJacksonFrictionalStress::
frictionalPressure

Alberto,

does this new code suffer from the same parallel error/instability as twoPhaseEulerFoam ?

Juho Peltola wrote:

What I meant was that why do you store the alpha.storePrevIter() at all, when the alpha is relaxed with alphaEqn.relax() instead of alpha.relax()? I've thought the storePrevIter() is only needed when the field (not the matrix) is relaxed.

I was experimenting with the local relaxation of alpha, as in MFIX. You can drop alpha storage: it is not used anywhere. The local relaxation in packed regions was not satisfactory (conservation error became definitely too large, without improvements on the stability).


I did test different continuous formulations for the frictional pressure and g0, but they didn't solve the packing limit problems with the PISO + underrelaxed alphaEqn formulation I used. In the final version with the frictional pressure flux limiter, the standard Johnson & Jackson formulation was used with less strict limiters.


Yes, I tried the same with PISO, and I did not stabilize the solution. Sooner or later the equation for alpha does not converge exactly, independently from the number of correctors performed.


The standard limiters in the twoPhaseEulerFoam are very strict: alphaMax-0.01 for the g0 and alphaMax-0.05 for the denominator of the Johnson & Jackson frictional pressure model. With commonly used values for alphaMinFriction (~0.6) and alphaMax (~0.63) the whole frictional stress model is drastically affected by the limiter.


Yes, I agree. All the tests I did were done with Schaeffer's model. For Johnson and Jackson, you can use a power series expansion to remove the discontinuity (recently implemented in MFIX too), so you do not actually have to enforce any strict limit. Similarly for g0.

Short addition. In the code I released, one change I made is missing since I made it after switching to the conservative form of the momentum equation.
A similar treatment to the one for ppMagf is used for the particle pressure, by writing

grad(pa) = (dpa/dalpha)*grad(alpha)

and including its effect in the alphaEqn, UEqn and pEqn, as done currently for ppMagf. This explains why the released code is still a bit troublesome when the kinetic theory is used (even if stabler than twoPhaseEulerFoam).

Best,
A.

David Huckaby wrote:
Alberto,

does this new code suffer from the same parallel error/instability as twoPhaseEulerFoam ?

Hi David,

it seems to me the code runs fine in parallel. I added a test case on 3 CPU's (same as jetBed but without kinetic theory).

If you run with the kinetic theory on, it is likely you have more difficulties, since that part of the code is not final yet. Also, if you run the tutorial, you will see the interfaces after a while breaks down, which demonstrate what I wrote before.

I should clarify that this code is really a starting point for further improvements, and *not* a final version of a solver (so take it "as is", which does not mean "stable, accurate and bugless", at least for now).

Best,
Alberto

I updated code and tutorials as follows:

- The particle pressure term is not managed "implicitly" in alphaEqn also for the kinetic theory. This requires the evaluation of the first derivative of g0.
I added a function to kineticTheoryModel so that ppMagf is easier to update, and it is computed using pa from KTGF if alpha < alphaMinFriction, and using pf if alpha >= alphaMinFriction (same as in MFIX).

- To avoid problems with discontinuities, I switched the jetBed tutorials to the CarnahanStarling distribution. In this way the frictional model becomes responsible of enforcing the packing limit. This means that you must set alphaMinFriction < alphaMax.
If you use the Schaeffer model (recommended, since the Johnson & Jackson model still needs some work to remove the discontinuity), alphaMinFriction = alphaMax - 0.01 or alphaMinFriction = alphaMax - 0.02 do the job.

Notice that, in this way, whatever model is used for ppMagf, the momentum predictor and the flux in pEqn do not need any modification.

I also removed the storage of the RHS in momentum predictors (not necessary at this point).

Interfaces are still terribly wavy, but that will come in another episode

Alberto Passalacqua wrote:
I updated code and tutorials as follows:

- The particle pressure term is not managed "implicitly" in alphaEqn also for the kinetic theory. This requires the evaluation of the first derivative of g0.
I added a function to kineticTheoryModel so that ppMagf is easier to update, and it is computed using pa from KTGF if alpha < alphaMinFriction, and using pf if alpha >= alphaMinFriction (same as in MFIX).

- To avoid problems with discontinuities, I switched the jetBed tutorials to the CarnahanStarling distribution. In this way the frictional model becomes responsible of enforcing the packing limit. This means that you must set alphaMinFriction < alphaMax.
If you use the Schaeffer model (recommended, since the Johnson & Jackson model still needs some work to remove the discontinuity), alphaMinFriction = alphaMax - 0.01 or alphaMinFriction = alphaMax - 0.02 do the job.

Notice that, in this way, whatever model is used for ppMagf, the momentum predictor and the flux in pEqn do not need any modification.

I also removed the storage of the RHS in momentum predictors (not necessary at this point).

Interfaces are still terribly wavy, but that will come in another episode

In the current implementation, the granular and frictional viscosities are summed up and the alpha used to calculate the gs0 is limited with

min(max(alpha_, 1e-6), alphaMax_ - 0.01)

On the other hand, the solids pressure is switched at the alphaMinFriction between the granular and frictional solids pressure models.

I believe the treatment of the viscosities and pressures should be consistent. Wouldn't it be better to just use the viscosity approach also for the pressure? That is, just sum up the granular and frictional ppMagf contributions and limit the alpha values in the gs0 calculation to avoid singularities at the packing limit? This would avoid the discontinuity at alphaMinFriction caused by the sharp switch between the models. The frictional pressure should still be responsible for enforcing the packing limit, because the granular temperature - and therefore the granular pressure - should approach zero at the packing limit (in theory).