[hemph]

Description:
The application gives erroneous solution for the Ua-field when the particle-particle force is included (g0 > 0). If settling particles are simulated, the velocity of the particles does not approach zero as they reach the bottom of the domain.

Solution. The ppMagf-term in alphaEqn.H: should be divided by alphaf on row 26.

Old:
ppMagf = g0*rUaAf*min(exp(preAlphaExp*(alphaf - alphaMax)), expMax);

New:
ppMagf = 1.0/(rhoa*max(alphaf,SMALL))*g0*rUaAf*min(exp(preAlpha Exp*(alphaf - alphaMax)), expMax);

Also row 32 in alphaEqn.H needs to be changed.

Old:
alphaEqn -= fvm::laplacian(ppMagf, alpha);

New:
alphaEqn -= fvm::laplacian(alphaf*ppMagf, alpha);

Solver/Application:
twoPhaseEulerFoam

Source file:
alphaEqn.H

Testcase:


Platform:
All

Version:
All

Notes:
The dimension of g0 is in Pascal in most references. In twoPhaseEulerFoam the dimensions of g0 is Pa/(kg/m3). To aid comparison between models, the ppMagf term should be divided by rho, as
ppMagf = 1.0/(rhoa*max(alphaf,SMALL))*g0*rUaAf*min(exp(preAlpha Exp*(alphaf - alphaMax)), expMax);

with a corresonding change to the dimensions of g0 in constant/ppProperties

[niklas]

Not a bug, they were left out deliberately.

ppMagf is only important around and below alphaMax and the g0 constant can be modified to account for the division/multiplication of alpha.

Niklas

[hemph]

Sorry, but why is this term left out? According to the equations they should be there. With the proposed fix included I get expected results. Without it, physics are wrong and results suspicious. Try the supplied test case for an example. I also disagree about ppMagf. It is certainly important over a wide range of volume fractions in packing of for instance a soft material.

//Rasmus

[hemph]

To expand the discussion. Although you could exclude the division of alpha in the ppMagf-calculation (for hard materials), the "bug", I believe, is in the laplacian-term. There ppMagf needs to be multiplied by alpha to give correct results.

[alberto]

Does the term in pEqn.H stay unchanged?

Regards,
Alberto

[hemph]

Sorry, missed your message Alberto. I believe that pEqn.H remains unchanged. The change in the code applies to the ppMagf-term, which is calculated in alphaEqn.H.

[alberto]

OK. Thanks for your confirmation. :-)

[juho]

It seems to me the equations that Rasmus reports as erroneous are still used in the twoPhaseEulerFoam in the 1.4.1 release.

Has there been chenges to the solver after the bug report? Is the bug report still valid?

[henry]

Are we all agreed that the proposed change is correct? If so I will put it in 1.5. Do we have a reference or derivation which covers this?

Henry

[alberto]

Expressions for the solids normal stress modulus are empirical. You can find a summary and a list of references here

http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0100-73862001000200006&l ng=&nrm=iso&tlng=#tab01

together with a short explanation of how they're used inside the particle phase momentum equation.

I have almost all the original references, if necessary.

Alberto

[henry]

From these do you agree with Rasmus' change?

Henry

[hemph]

Which change is under discussion here?

I have reported five different issues. I do not have the 1.41 source heres, so I can not tell which ones have been corrected.

First, there is a missing alpha in the laplacian-expression of the alphaEqn. The inclusion of alpha can be shown from the equations and was found to be important to the physics.

Second, ppMagf should be divided by alpha*rho. Niklas left this out and left it to the user to modify the expression for G. However, for soft materials, there is a particle-particle force over a wide range of volume fractions, which falsifies the assumption of alpha at packing being a constant. (Also, as alpha has no dimensions, the user would need to read through the source code to understand that she needed to multiply her constant by alpha..)

Thirdly, there is (was?) a grave bug in the interpolation of particle-particle forces to the faces. It has been reported here http://www.cfd-online.com/OpenFOAM_D...tml?1173086743

Fourthly, there is a problem with certain boundary conditions for particles. It was reported here http://www.cfd-online.com/OpenFOAM_D...tml?1170065229

Rasmus

[henry]

Hi Rasmus,

I assumed that it was the change at the top of this list being discussed. I am happy to include any changes you think appropriate but as I am not currently developing or running this solver I can't test or verify it. If you could send me the latest code I would be happy to merge it into the 1.5 release if you would like me to.

Thanks

Henry

[juho]

As far as I can see at least issues 1-3 exist in the 1.4.1.

[hemph]

Hi Henry and Juho, I think that it would be good if the changes came into 1.5. They have been found to increase stability quite a lot. The code I have been running is largely modified for a possible article and is hard to submit for inclusion as is.

However, here is a suggestion of a fix for the bug in interpolation (third bug in my post) and particle forces calculation (first and second). The changes needs to go into alphaEqn.H. Warning, I do not have access to a linux computer at the moment, so there is possibly some syntatic errors. Note that the dimensions of g0 need to change in the tutorial test case as well. I would be very happy if someone was willing to test it on the tutorials.

Old code:
if (g0.value() > 0.0)
{
surfaceScalarField alphaf = fvc::interpolate(alpha);
// Correct the particle-particle force magnitude
ppMagf = g0*rUaAf*min(exp(preAlphaExp*(alphaf - alphaMax)), expMax);
alphaEqn -= fvm::laplacian(ppMagf, alpha);
}

New code:
if (g0.value() > 0.0)
{
volScalarField G = (1.0/(rhoa*(alpha+scalar(0.0001))))*g0*min(exp(preAlpha Exp*(alpha - alphaMax)), expMax));

ppMagf = rUaAf*fvc::interpolate(G);
surfaceScalarField alphaf = fvc::interpolate(alpha);

alphaEqn -= fvm::laplacian(((alphaf+scalar(0.0001))*ppMagf), alpha);
}

Rasmus

[juho]

Thank you. It compiled fine. Requires an ugly entry to fvScemes.

I can run tests on the tutorials. At least it seems to run fine on another test case I have.

[hemph]

Good. Thanks for trying it out. To get rid of the ugly entry, I believe you can give the term a name, which is then called in fvSchemes. In Ueqn.H for example, it is fvm::div(phia, Ua, "div(phia,Ua)").

[juho]

Is there a reason not to use

max(alphaf,SMALL)
max(alpha,SMALL)

instead of

(alphaf+scalar(0.0001))
(alpha+scalar(0.0001))

?

The latter seems a lot less stable.

[juho]

False alarm. The instability wasn't related to the max(alphaf,SMALL) vs. (alphaf+scalar(0.0001)).

But is there a reason to choose on over the other?


Does max() compare and replace values of each cell in of any field?

A difference between the approaches is of course that the max() filters out even negative values of magnitude larger than 0.0001 and for debugging reasons it migth be nice that such values cause an error.

[hemph]

The reason for going with 0.0001 here was an attempt to make the division a little better numerically, as SMALL is 1e-15. The same term is multiplied in the laplacian equation later. I do not know if it made any practical difference, but I agree that is a bit ugly with a free scalar like that.

If I recall correctly, you can not do max(alphaf, SMALL) unless both terms are of the same dimension (because of the replacement), so you need to create a new dimensioned scalar, but principally it should be ok.