CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Bugs

TwoPhaseEulerFOAM application

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   November 16, 2006, 08:25
Default Description: The application
  #1
Senior Member
 
Rasmus Hemph
Join Date: Mar 2009
Location: Sweden
Posts: 108
Rep Power: 8
hemph is on a distinguished road
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
hemph is offline   Reply With Quote

Old   November 16, 2006, 08:33
Default Not a bug, they were left out
  #2
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 19
niklas will become famous soon enough
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
niklas is offline   Reply With Quote

Old   November 16, 2006, 13:09
Default Sorry, but why is this term le
  #3
Senior Member
 
Rasmus Hemph
Join Date: Mar 2009
Location: Sweden
Posts: 108
Rep Power: 8
hemph is on a distinguished road
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 is offline   Reply With Quote

Old   November 16, 2006, 14:18
Default To expand the discussion. Alth
  #4
Senior Member
 
Rasmus Hemph
Join Date: Mar 2009
Location: Sweden
Posts: 108
Rep Power: 8
hemph is on a distinguished road
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.
hemph is offline   Reply With Quote

Old   November 22, 2006, 16:45
Default Does the term in pEqn.H stay u
  #5
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Does the term in pEqn.H stay unchanged?

Regards,
Alberto
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   January 19, 2007, 08:01
Default Sorry, missed your message Alb
  #6
Senior Member
 
Rasmus Hemph
Join Date: Mar 2009
Location: Sweden
Posts: 108
Rep Power: 8
hemph is on a distinguished road
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.
hemph is offline   Reply With Quote

Old   January 19, 2007, 08:43
Default OK. Thanks for your confirmati
  #7
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
OK. Thanks for your confirmation. :-)
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   May 26, 2008, 03:12
Default It seems to me the equations t
  #8
Member
 
Juho Peltola
Join Date: Mar 2009
Location: Finland
Posts: 80
Rep Power: 8
juho is on a distinguished road
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?
juho is offline   Reply With Quote

Old   May 26, 2008, 05:59
Default Are we all agreed that the pro
  #9
Senior Member
 
Join Date: Mar 2009
Posts: 854
Rep Power: 13
henry is on a distinguished road
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
henry is offline   Reply With Quote

Old   May 26, 2008, 23:09
Default Expressions for the solids nor
  #10
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
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
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   May 27, 2008, 03:04
Default From these do you agree with R
  #11
Senior Member
 
Join Date: Mar 2009
Posts: 854
Rep Power: 13
henry is on a distinguished road
From these do you agree with Rasmus' change?

Henry
henry is offline   Reply With Quote

Old   May 27, 2008, 07:26
Default Which change is under discussi
  #12
Senior Member
 
Rasmus Hemph
Join Date: Mar 2009
Location: Sweden
Posts: 108
Rep Power: 8
hemph is on a distinguished road
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
hemph is offline   Reply With Quote

Old   May 27, 2008, 07:50
Default Hi Rasmus, I assumed that i
  #13
Senior Member
 
Join Date: Mar 2009
Posts: 854
Rep Power: 13
henry is on a distinguished road
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
henry is offline   Reply With Quote

Old   May 27, 2008, 07:51
Default As far as I can see at least i
  #14
Member
 
Juho Peltola
Join Date: Mar 2009
Location: Finland
Posts: 80
Rep Power: 8
juho is on a distinguished road
As far as I can see at least issues 1-3 exist in the 1.4.1.
juho is offline   Reply With Quote

Old   May 27, 2008, 09:00
Default Hi Henry and Juho, I think tha
  #15
Senior Member
 
Rasmus Hemph
Join Date: Mar 2009
Location: Sweden
Posts: 108
Rep Power: 8
hemph is on a distinguished road
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
hemph is offline   Reply With Quote

Old   May 27, 2008, 09:17
Default Thank you. It compiled fine. R
  #16
Member
 
Juho Peltola
Join Date: Mar 2009
Location: Finland
Posts: 80
Rep Power: 8
juho is on a distinguished road
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.
juho is offline   Reply With Quote

Old   May 27, 2008, 09:55
Default Good. Thanks for trying it out
  #17
Senior Member
 
Rasmus Hemph
Join Date: Mar 2009
Location: Sweden
Posts: 108
Rep Power: 8
hemph is on a distinguished road
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)").
hemph is offline   Reply With Quote

Old   May 27, 2008, 10:26
Default Is there a reason not to use
  #18
Member
 
Juho Peltola
Join Date: Mar 2009
Location: Finland
Posts: 80
Rep Power: 8
juho is on a distinguished road
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 is offline   Reply With Quote

Old   May 27, 2008, 11:09
Default False alarm. The instability w
  #19
Member
 
Juho Peltola
Join Date: Mar 2009
Location: Finland
Posts: 80
Rep Power: 8
juho is on a distinguished road
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.
juho is offline   Reply With Quote

Old   May 27, 2008, 14:32
Default The reason for going with 0.00
  #20
Senior Member
 
Rasmus Hemph
Join Date: Mar 2009
Location: Sweden
Posts: 108
Rep Power: 8
hemph is on a distinguished road
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.
hemph is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
TwoPhaseEulerFoam sara OpenFOAM Running, Solving & CFD 2 November 6, 2008 20:26
Bug in twoPhaseEulerFoam alberto OpenFOAM Bugs 2 May 20, 2008 21:25
TwoPhaseEulerFoam Bug alondono OpenFOAM Bugs 1 February 19, 2008 21:01
TwoPhaseEulerFOAM application hemph OpenFOAM Bugs 0 November 16, 2006 08:27
TwoPhaseEulerFoam newbee OpenFOAM 0 March 27, 2006 08:41


All times are GMT -4. The time now is 12:45.