CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

mass conservation of solid phase violated when using groovyBC with twoPhaseEulerFoam

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

Like Tree1Likes
  • 1 Post By xpqiu

Reply
 
LinkBack Thread Tools Display Modes
Old   June 1, 2015, 23:24
Default mass conservation of solid phase violated when using groovyBC with twoPhaseEulerFoam
  #1
New Member
 
Qiu Xiaoping
Join Date: Apr 2013
Location: IPE CAS China
Posts: 8
Rep Power: 5
xpqiu is on a distinguished road
Hi, every one

I am working on OpenFOAM-2.1.1, I need to simulate a gas-solid two phase flow with twoPhaseEulerFoam, the geometry of my simulation is a 2d-rectangular, bottom is inlet, top is outlet, and side patches are walls.
gas phase velocity BC is set as follows:

"Ub"
Code:
boundaryField
{
    walls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    outlet
    {
       type             zeroGradient;
    }
    inlet
    {
        type               fixedValue;
        value              uniform (0 0.20 0);
    }
    frontAndBackPlanes
    {
        type               empty;
    }
}
For solid phase, I need the solid exitting the top outlet recirculate back to the bottom inlet, so I used groovyBC for inlet of solid phase velocity, as follows:

"Ua"
Code:
boundaryField
{
    walls
    {
        type               fixedValue;
        value              uniform (0 0 0);
    }
    outlet
    {
        type               zeroGradient;
    //type                 inletOutlet;
    //inletValue           uniform (0 0 0 );
    //value                uniform (0 0 0 );
    }
    inlet
    {
        //type               fixedValue;
        //value              uniform (0 0 0);
    type                 groovyBC;
    valueExpression      "-inVel*normal()"
    value                uniform ( 0 0 0 );
    variables (
        "A=sum(area());"
        "outFlow{outlet}=sum(Ua&normal()*area()*alpha);"
        "myFlow=outFlow/alpha;"
        "inVel=myFlow/A;"
    );

    }
    frontAndBackPlanes
    {
        type               empty;
    }
}
where "alpha" is the volume fraction of solid phase.

Boundary condition for alpha are set as:


"alpha"
Code:
boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 0.2;
    }
    outlet
    {
        type              zeroGradient;
    }
    walls
    {
        type            zeroGradient;
    }
    frontAndBackPlanes
    {
        type            empty;
    }
Note that I used a fixed inlet solid phase volume fraction.


The simulation worked fine at the beginning, but as the solid phase approached the top outlet, and started to recirculated back to bottom inlet, things got strange. It seems that the mass conservation of solid phase is violated. Here are some of my logs:

At the beginning , total solid phase volume fraction is:
Code:
Dispersed phase volume fraction = 0.2925  Min(alpha) = 0  Max(alpha) = 0.6
But at Time = 8.65s, solid phase approached at the top outlet, and total solid phase volume fraction started to increase:

Code:
Time = 8.65

DILUPBiCG:  Solving for alpha, Initial residual = 0.00712221, Final residual = 6.46064e-11, No Iterations 3
Dispersed phase volume fraction = 0.292501  Min(alpha) = 0.0141432  Max(alpha) = 0.565962
from that on, total solid phase volume fraction increased monotonically, and reached 0.293174 at time 20.4s :

Code:
Time = 20.4

DILUPBiCG:  Solving for alpha, Initial residual = 0.00496807, Final residual = 1.18516e-11, No Iterations 3
Dispersed phase volume fraction = 0.293174  Min(alpha) = 0.0697732  Max(alpha) = 0.513855


I monitored the flux of solid phase at outlet and inlet with "simpleSwakFunctionObjects" , here is the code in the controlDict file:

Code:
massConservationTest
    {
        type patchExpression;
        outputControlMode   outputTime;
    patches (
        inlet
        outlet
        );
    verbose true;
    expression "(Ua & normal())*alpha*area()";
    accumulations (
        sum
    );
    }


And the flux at inlet and outlet are always equal in value and opposite in direction.
Code:
       Time       flux of inlet       flux of outlet
          8.65     -0.000376597          0.000376597
          8.7     -0.000382609          0.000382609
         8.75     -0.000377901          0.000377901
          8.8     -0.000374998          0.000374998
         8.85     -0.000373855          0.000373855
          8.9     -0.000377711          0.000377711
         8.95     -0.000387632          0.000387632
            9     -0.000393048          0.000393048
         9.05      -0.00039903           0.00039903
          9.1     -0.000403242          0.000403242
         9.15     -0.000376861          0.000376861
          9.2     -0.000323375          0.000323375
         9.25     -0.000147949          0.000147949
          9.3      -0.00010581           0.00010581
         9.35     -0.000138232          0.000138232
          9.4     -0.000138196          0.000138196
         9.45     -4.15751e-05          4.15751e-05
          9.5     -2.36378e-05          2.36378e-05
         9.55      -9.4021e-05           9.4021e-05
          9.6     -8.63874e-05          8.63874e-05
         9.65     -8.27537e-05          8.27537e-05
          9.7     -9.60432e-05          9.60432e-05
         9.75     -0.000152863          0.000152863
          9.8     -0.000159683          0.000159683
         9.85     -0.000191661          0.000191661
          9.9      -0.00022118           0.00022118
         9.95     -0.000152557          0.000152557
           10     -0.000189653          0.000189653
        10.05     -0.000233514          0.000233514
         10.1     -0.000197505          0.000197505
        10.15     -0.000214339          0.000214339
         10.2     -0.000235228          0.000235228
        10.25     -0.000221006          0.000221006
         10.3     -0.000226268          0.000226268
        10.35      -0.00025188           0.00025188
         10.4     -0.000264869          0.000264869
        10.45     -0.000206712          0.000206712
         10.5     -0.000236892          0.000236892
        10.55     -0.000204415          0.000204415
         10.6     -6.89294e-05          6.89294e-05
        10.65     -7.52851e-05          7.52851e-05
         10.7     -9.99905e-05          9.99905e-05
        10.75     -6.23968e-05          6.23968e-05
         10.8      8.60961e-05         -8.60961e-05
        10.85     -4.74735e-05          4.74735e-05
         10.9     -7.44399e-05          7.44399e-05
        10.95     -3.77192e-05          3.77192e-05
           11     -1.52147e-05          1.52147e-05


It seems that the flux on both inlet and outlet are reasonable, and flux on
inlet is exactly conserved with flux on outlet, but why the total volume fraction of solid phase is not conserved ?

Thanks a lot for any hints

Last edited by xpqiu; June 4, 2015 at 21:15.
xpqiu is offline   Reply With Quote

Old   June 3, 2015, 03:01
Default
  #2
New Member
 
Join Date: Apr 2015
Posts: 4
Rep Power: 3
LassiL is on a distinguished road
Dear experts,

I am struggling with the same kind of issue and thus I want to bring out that your help will be appreciated by me, too.

Regards,
Lassi L.
LassiL is offline   Reply With Quote

Old   June 4, 2015, 10:52
Default
  #3
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Torino, Italy
Posts: 743
Rep Power: 9
sharonyue is on a distinguished road
Quote:
Originally Posted by xpqiu View Post
Hi, every one


It seems that the flux on both inlet and outlet are reasonable, and flux on [/COLOR][/FONT] inlet is exactly conserved with flux on outlet, but why the total volume fraction of solid phase in not conserved ?

Thanks a lot for any hints
Hi,

this happens on tons of cases, here are some threads maybe u are interested into it. now its not solved out. Anyway u can force it to be constant.

BTW, this is not a serious problem. But yeah, it weird. This happens in 22x,23x,211(by your test).

http://www.openfoam.org/mantisbt/view.php?id=1700
http://www.openfoam.org/mantisbt/view.php?id=1237#c3117
conservative form twoPhaseEulerFoam23x gives an unsatisfactory velocity result.

I expect in 23x this should be solved. but its not. even with an conservative equation implemented.

Best,
__________________
Im the translator of OpenFOAM User Guide Chinese Edition! But always newbie on CFD.

http://dyfluid.com/en.html
sharonyue is offline   Reply With Quote

Old   June 4, 2015, 10:56
Default
  #4
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Torino, Italy
Posts: 743
Rep Power: 9
sharonyue is on a distinguished road
Have u ever tried a case without swakFoam? I think this is not related with swak. This dispersed phase volume not constant problem also happens with the tutorial cases.
__________________
Im the translator of OpenFOAM User Guide Chinese Edition! But always newbie on CFD.

http://dyfluid.com/en.html
sharonyue is offline   Reply With Quote

Old   June 4, 2015, 21:23
Default
  #5
New Member
 
Qiu Xiaoping
Join Date: Apr 2013
Location: IPE CAS China
Posts: 8
Rep Power: 5
xpqiu is on a distinguished road
Hi Dongyue,
Quote:
Originally Posted by sharonyue View Post
Have u ever tried a case without swakFoam? I think this is not related with swak.
Thanks for the reply,I have tried bubbling fluidized bed without swak4Foam, in that case, Ua set as "uniform (0 0 0)" on both inlet and outlet, and alpha is set as "uniform 0" on both inlet and outlet. In my bubbling fluidized bed case, mass of solid phase is exactly conserved.
__________________
xpqiu
q.giskard@gmail.com
xpqiu is offline   Reply With Quote

Old   June 4, 2015, 21:46
Default
  #6
New Member
 
Qiu Xiaoping
Join Date: Apr 2013
Location: IPE CAS China
Posts: 8
Rep Power: 5
xpqiu is on a distinguished road
Hi, Dongyue
Quote:
Originally Posted by sharonyue View Post
this happens on tons of cases, here are some threads maybe u are interested into it.
Thanks for your reply. Yesterday, I searched this forum for "twoPhaseEulerFoam" related topic, and I have read those threads posted by you, thank you for the information.

I have tried some bubbling bed cases and circulating fluidized bed(CFB) cases, and violation of mass only happens in my CFB cases.

As I know, in OpenFOAM-2.1.1, alphaEqn is solved with iteration method, but since OpenFOAM-2.2.x, alphaEqn is solved with MULES. I didn't try to run my case on OpenFOAM-2.2.x and 2.3.x.

I happened to known that MULES is an effective way to limit the maximum value of alpha, to avoid an overpacking of solid phase, so I tried to modify the "twoPhaseEulerFoam" in OpenFOAM-2.1.1, using MULES to solve the alphaEqn. I found that for the same circulating fluidized bed case, violation of mass conversation is more serious in my modified solver, so I don't known if to some extent MULES is responsible for the mass loss(addition). By the way, do you have some reference about MULES?

Best
__________________
xpqiu
q.giskard@gmail.com

Last edited by xpqiu; June 17, 2015 at 01:55.
xpqiu is offline   Reply With Quote

Old   June 15, 2015, 04:04
Default
  #7
New Member
 
Join Date: Apr 2015
Posts: 4
Rep Power: 3
LassiL is on a distinguished road
Dear experts,

I have a question related to the code presented in this thread.
Ua:

Code:
inlet     {         
//type               fixedValue;         
//value              uniform (0 0 0);     
type                 groovyBC;     
valueExpression      "-inVel*normal()"     
value                uniform ( 0 0 0 );     
variables (         
          "A=sum(area());"         
          "outFlow{outlet}=sum(Ua&normal()*area()*alpha);"
          "myFlow=outFlow/alpha;"         
          "inVel=myFlow/A;"     
           );
}
In the phrase "A=sum(area());" the word area() obviously stands for the area of inlet. But does the world area() mean the area of inlet or the area of outlet on the following row?
Code:
"outFlow{outlet}=sum(Ua&normal()*area()*alpha);"
I ask because in my case the areas of inlet and outlet are not equal and thus it makes a difference. Also, in the code of the controlDict file,
Code:
massConservationTest
    {
        type patchExpression;
        outputControlMode   outputTime;
    patches (
        inlet
        outlet
        );
    verbose true;
    expression "(Ua & normal())*alpha*area()";
    accumulations (
        sum
    );
    }
Does the word area() stand for area in general (different area in different patch), area of inlet or area of outlet?

Thank you.
LassiL is offline   Reply With Quote

Old   June 15, 2015, 04:16
Default
  #8
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Torino, Italy
Posts: 743
Rep Power: 9
sharonyue is on a distinguished road
Quote:
Originally Posted by xpqiu View Post
Hi, Dongyue


Thanks for your reply. Yesterday, I searched this forum for "twoPhaseEulerFoam" related topic, and I have read those threads posted by you, thank you for the information.

I have tried some bubbling bed cases and circulating fluidized bed(CFB) cases, and violation of mass only happens in my CFB cases.

As I know, in OpenFOAM-2.1.1, alphaEqn is solved with iteration method, but since OpenFOAM-2.2.x, alphaEqn is solved with MULES. I didn't try to run my case on OpenFOAM-2.2.x and 2.3.x.

I happened to known that MULES is an effect way to limit the maximum value of alpha, to avoid an overpacking of solid phase, so I tried to modify the "twoPhaseEulerFoam" in OpenFOAM-2.1.1, using MULES to solve the alphaEqn. I found that for the same circulating fluidized bed case, violation of mass conversation is more serious in my modified solver, so I don't known if to some extent MULES is responsible for the mass loss(addition). By the way, do you have some reference about MULES?

Best
Hello, xiaoping,

That has been a really long time since I touched the OpenFOAM-2.1.1 last time. From your depiction, If MULES violate the mass conservation, thats interesting. I never tried with solve alphaEqn by conventional method. But maybe solving alphaEqn(without MULES) in OpenFOAM 2.2.x or 2.3.x can provides a mass conservation result?

Maybe someday I can make a test, for now I just force it to be mass conservation.

Also maybe u can report a bug~

Best,
__________________
Im the translator of OpenFOAM User Guide Chinese Edition! But always newbie on CFD.

http://dyfluid.com/en.html
sharonyue is offline   Reply With Quote

Old   June 17, 2015, 02:08
Default
  #9
New Member
 
Qiu Xiaoping
Join Date: Apr 2013
Location: IPE CAS China
Posts: 8
Rep Power: 5
xpqiu is on a distinguished road
Dear LassiL,

In the phrase "A=sum(area())" , "area()" means area of current patch.

But in
Quote:
"outFlow{outlet}=sum(Ua&normal()*area()*alpha) ;"
"area()" means area of outlet.

The trick is that when you use "{patch_name}" after a variable, then all the terms on the right hand side are properties of the specified patch "patch_name".

And in the following code,
Quote:
Code:
massConservationTest    
 {         
type patchExpression;        
 outputControlMode   outputTime;   
  patches (  
       inlet    
     outlet        
 );     
verbose true;     
expression "(Ua & normal())*alpha*area()";     
accumulations (  
       sum  
   );     }
The code above calculate the flux on both inlet and outlet, and when flux on patch "inlet" is calculated, "area()" means area of inlet, when flux on patch "outlet" is calculated, "area()" means area of outlet.
LassiL likes this.
__________________
xpqiu
q.giskard@gmail.com
xpqiu is offline   Reply With Quote

Reply

Tags
groovybc, twophaseeulerfoam

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
Mass Transfer without Phase Change in InterPhaseChangeFoam Parisa_Khiabani OpenFOAM Running, Solving & CFD 5 May 26, 2015 09:40
twoPhaseEulerFoam, mass loss and velocity profile problems mwaqas OpenFOAM Running, Solving & CFD 0 November 14, 2014 18:44
How to calculate mass flow of the discrete phase in Euler-lagrange model juliom CFX 0 November 8, 2011 20:50
Mass conservation Diego Nogueira Main CFD Forum 1 July 30, 2004 15:50
How to calculate density of solid phase zhou Main CFD Forum 0 December 17, 1999 20:06


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