[xpqiu]

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"

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"

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"

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:

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:

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 :

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:

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.

       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

[LassiL]

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.

[sharonyue]

Originally Posted by xpqiu

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
http://www.cfd-online.com/Forums/ope...ty-result.html

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

Best,

[sharonyue]

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.

[xpqiu]

Hi Dongyue，

Originally Posted by sharonyue

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]

Hi, Dongyue

Originally Posted by sharonyue

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

[LassiL]

Dear experts,

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

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?

"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,

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.

[sharonyue]

Originally Posted by xpqiu

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,

[xpqiu]

Dear LassiL,

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

But in

"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,

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.