CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Bugs (
-   -   TwoPhaseEulerFOAM application (

hemph November 16, 2006 08:25

Description: The application
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.

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

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.

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

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


Source file:




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 November 16, 2006 08:33

Not a bug, they were left out
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.


hemph November 16, 2006 13:09

Sorry, but why is this term le
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.


hemph November 16, 2006 14:18

To expand the discussion. Alth
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 November 22, 2006 16:45

Does the term in pEqn.H stay u
Does the term in pEqn.H stay unchanged?


hemph January 19, 2007 08:01

Sorry, missed your message Alb
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 January 19, 2007 08:43

OK. Thanks for your confirmati
OK. Thanks for your confirmation. :-)

juho May 26, 2008 03:12

It seems to me the equations t
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 May 26, 2008 05:59

Are we all agreed that the pro
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?


alberto May 26, 2008 23:09

Expressions for the solids nor
Expressions for the solids normal stress modulus are empirical. You can find a summary and a list of references here 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.


henry May 27, 2008 03:04

From these do you agree with R
From these do you agree with Rasmus' change?


hemph May 27, 2008 07:26

Which change is under discussi
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

Fourthly, there is a problem with certain boundary conditions for particles. It was reported here


henry May 27, 2008 07:50

Hi Rasmus, I assumed that i
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.



juho May 27, 2008 07:51

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

hemph May 27, 2008 09:00

Hi Henry and Juho, I think tha
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);


juho May 27, 2008 09:17

Thank you. It compiled fine. R
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 May 27, 2008 09:55

Good. Thanks for trying it out
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 May 27, 2008 10:26

Is there a reason not to use
Is there a reason not to use


instead of



The latter seems a lot less stable.

juho May 27, 2008 11:09

False alarm. The instability w
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 May 27, 2008 14:32

The reason for going with 0.00
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.

All times are GMT -4. The time now is 21:56.