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

Bubble Foam momentum equation

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

Reply
 
LinkBack Thread Tools Display Modes
Old   September 17, 2010, 00:56
Default Bubble Foam momentum equation
  #1
Senior Member
 
Balkrishna Patankar
Join Date: Mar 2009
Location: Pune
Posts: 119
Rep Power: 7
balkrishna is on a distinguished road
I want to use the standard momentum equation in bubble Foam rather than the phase intensive one as suggested by Henrik Rusche in his PhD thesis... I got the following problems :

Code:
UEqns.H

volVectorField Uamod = alpha*Ua ;
volVectorField Ubmod = beta*Ub;
fvVectorMatrix UaEqn(Ua, Ua.dimensions()*dimVol/dimTime);
fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);

{
    volTensorField Rca = -nuEffa*(fvc::grad(Ua)().T());
    Rca = Rca + (2.0/3.0)*sqr(Ct)*I*k - (2.0/3.0)*I*tr(Rca);

    Rca = alpha*Rca ;
    surfaceScalarField phiRa =
      - fvc::interpolate(nuEffa)
       *mesh.magSf()*fvc::snGrad(alpha);


    UaEqn =
    (
        (scalar(1) + Cvm*rhob*beta/rhoa)*
        (
            fvm::ddt(Uamod)
            + fvm::div(phia,Uamod,"div(phia,Ua,alpha)")
            - fvm::Sp(fvc::div(phia,alpha), Ua)
        )

        - fvm::laplacian(nuEffa,Uamod)
        + fvc::div(Rca)

      + fvm::div(phiRa,Ua,"div(phia,Ua)")
      - fvm::Sp(fvc::div(phiRa), Ua)
      + (fvc::grad(alpha)&Rca)
     ==
    //  g                        // Buoyancy term transfered to p-equation
      - fvm::Sp(alpha*beta/rhoa*dragCoef, Ua)
    //    + beta/rhoa*dragCoef*Ub    // Explicit drag transfered to p-equation
      - alpha*beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb)
    );

    UaEqn.relax();
This compiles fine but I am not able to run it . The fvSchemes file is modified as follows :
Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.7.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phia,Uamod)    Gauss limitedLinearV 1;
    div(phib,Ubmod)    Gauss limitedLinearV 1;
    div(phia,Ua)    Gauss limitedLinearV 1;
    div(phib,Ub)    Gauss limitedLinearV 1;
    div(phib,k)     Gauss limitedLinear 1;
    div(phib,epsilon) Gauss limitedLinear 1;
    div(phi,alpha)  Gauss limitedLinear01 1;
    div((-nuEffa*grad(Ua).T())) Gauss linear;
    div((-nuEffb*grad(Ub).T())) Gauss linear;
}

laplacianSchemes
{
    default         none;
    laplacian(nuEffa,(alpha*Ua)) Gauss linear corrected;
    laplacian(nuEffb,(beta*Ub)) Gauss linear corrected;
    laplacian(nuEffa,Ua) Gauss linear corrected;
    laplacian(nuEffb,Ub) Gauss linear corrected;
    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
    laplacian(DkEff,k) Gauss linear corrected;
    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
    laplacian(DREff,R) Gauss linear corrected;
    laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p               ;
}
This gives the following error .... what can be done to correct this ???
Code:
    gradientInternalCoeffs cannot be called for a calculatedFvPatchField
    on patch inlet of field (alpha*Ua) in file "/home/ifmg/MySIM/CHECK/balBubbleFoam/bubbleColumn/0/(alpha*Ua)"
    You are probably trying to solve for a field with a default boundary condition.

    From function calculatedFvPatchField<Type>::gradientInternalCoeffs() const
    in file fields/fvPatchFields/basic/calculated/calculatedFvPatchField.C at line 186.
Any suggestion to overcoming this problem ...
balkrishna is offline   Reply With Quote

Old   September 17, 2010, 04:29
Default
  #2
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 367
Rep Power: 13
l_r_mcglashan will become famous soon enough
Are you sure about doing this? There are some very good reasons for dividing the momentum equation by the phase fraction! Alberto Passalacqua was working on this, have a read through this thread:

Conservative Two Phase Euler

The error is probably due to the boundary condition not being explicitly stated. You can either setup the Uamod field in the 0/folder, use the == operator instead of = (not sure if that works... try it) or

Code:
       volVectorField Uamod
        (
            IOobject
            (
                "Uamod",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            1.0,
            Ua.boundaryField().types()
        );
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   September 17, 2010, 04:33
Default
  #3
Senior Member
 
Balkrishna Patankar
Join Date: Mar 2009
Location: Pune
Posts: 119
Rep Power: 7
balkrishna is on a distinguished road
thanks for the link ....
balkrishna is offline   Reply With Quote

Old   September 17, 2010, 04:40
Default
  #4
Senior Member
 
Balkrishna Patankar
Join Date: Mar 2009
Location: Pune
Posts: 119
Rep Power: 7
balkrishna is on a distinguished road
I read the same arguments in Dr. Rusche's thesis .... The doubt I had was what is exactly meant by equations becoming singular ... I didnt quite understand the effect it would have on the flow ... For the same purpose , I decided to implement momentum eqn in the conservative form and check out how it affected the solution .....

Another doubt I face is :
How is the non conservative form of the momentum equation derived ? i.e. how is it possible to divide by alpha when alpha is a variable in itself ....
balkrishna is offline   Reply With Quote

Old   September 17, 2010, 04:42
Default
  #5
Senior Member
 
Balkrishna Patankar
Join Date: Mar 2009
Location: Pune
Posts: 119
Rep Power: 7
balkrishna is on a distinguished road
Presently I am getting the following error on running the program :

Code:
 FOAM FATAL ERROR: 
incompatible fields for operation 
    [Ua] - [Uamod]

    From function checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)
    in file /home/ifmg/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude/fvMatrix.C at line 1181.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) in "/home/ifmg/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libOpenFOAM.so"
#1  Foam::error::abort() in "/home/ifmg/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libOpenFOAM.so"
#2  

#3  

#4  

#5  __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#6  

Aborted
balkrishna is offline   Reply With Quote

Old   September 17, 2010, 04:48
Default
  #6
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 367
Rep Power: 13
l_r_mcglashan will become famous soon enough
Another link in case you haven't spotted it is the bubbleFoam wiki article. Again, thanks to Alberto.

I haven't got first hand experience of the difficulties you'll encounter, but effectively in areas where the phase fraction is zero, your momentum becomes zero and the velocity becomes undefined. The matrix probably becomes ill conditioned.

I have a full derivation. If you want to see how to get it for yourself, just expand the derivatives in the momentum equation and then divide the whole equation by its phase fraction.
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   September 17, 2010, 04:50
Default
  #7
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 367
Rep Power: 13
l_r_mcglashan will become famous soon enough
I just spotted that your UaEqn/UbEqn should have Uamod and Ubmod as the variable to solve for, not Ua and Ub. There may be other mistakes.
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   September 17, 2010, 04:51
Default
  #8
Senior Member
 
Balkrishna Patankar
Join Date: Mar 2009
Location: Pune
Posts: 119
Rep Power: 7
balkrishna is on a distinguished road
can u send me a link to the derivation .... thanks ...
balkrishna is offline   Reply With Quote

Old   September 17, 2010, 05:24
Default
  #9
Senior Member
 
Balkrishna Patankar
Join Date: Mar 2009
Location: Pune
Posts: 119
Rep Power: 7
balkrishna is on a distinguished road
got it!!! derived it myself ... but are there cases when using a non conservative form for multiphase flow may result in spurious results ...
balkrishna is offline   Reply With Quote

Old   September 18, 2010, 05:03
Default
  #10
Senior Member
 
Balkrishna Patankar
Join Date: Mar 2009
Location: Pune
Posts: 119
Rep Power: 7
balkrishna is on a distinguished road
Have you shared (made it open source ) the PEA source code ???
balkrishna is offline   Reply With Quote

Old   September 20, 2010, 11:36
Default
  #11
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 367
Rep Power: 13
l_r_mcglashan will become famous soon enough
Not yet, maybe soon!
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   September 20, 2010, 12:04
Default
  #12
Senior Member
 
Balkrishna Patankar
Join Date: Mar 2009
Location: Pune
Posts: 119
Rep Power: 7
balkrishna is on a distinguished road
I have started working on the same ..... right now working on two things ...
A spalding approach to calculate the volume fractions and the second one being the PEA ... Have completed the first one ... and will start the second part soon ....

How are you programming the PEA ? As a general elimination algorithm or hard coding it ...

Last edited by balkrishna; September 20, 2010 at 12:04. Reason: appropriate question
balkrishna is offline   Reply With Quote

Old   September 20, 2010, 12:09
Default
  #13
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 367
Rep Power: 13
l_r_mcglashan will become famous soon enough
The PEA is Spalding's approach, no? That's the method where you substitute one momentum equation into the other and rearrange.

I created a multiphase PISO abstract class and have semiImplicit/partialElimination/multiphasePartialElimination derived classes.

I don't think you need the PEA for bubbles. It's only necessary when the drag force is large, which is true for small particles in fluidised beds.
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   September 20, 2010, 12:21
Default
  #14
Senior Member
 
Balkrishna Patankar
Join Date: Mar 2009
Location: Pune
Posts: 119
Rep Power: 7
balkrishna is on a distinguished road
I meant the spalding approach to solving the volume fraction ... discretise the volume fraction eqns and then calculate alpha = alphadis/(alphadis + betadis) where alphadis and betadis are the discretized volume fraction eqns ...
balkrishna is offline   Reply With Quote

Old   September 23, 2010, 02:17
Default
  #15
Senior Member
 
Balkrishna Patankar
Join Date: Mar 2009
Location: Pune
Posts: 119
Rep Power: 7
balkrishna is on a distinguished road
I have come up with an implementation of the PEA .. .....
The equations i have used are on this page ... https://sites.google.com/site/balkrishnanitt/eq

Can anyone tell me whether the implementation is correct ....
Code:
fvVectorMatrix UaEqn(Ua, Ua.dimensions()*dimVol/dimTime);
fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);

{
    volTensorField Rca = -nuEffa*(fvc::grad(Ua)().T());
    Rca = Rca + (2.0/3.0)*sqr(Ct)*I*k - (2.0/3.0)*I*tr(Rca);

    surfaceScalarField phiRa =
      - fvc::interpolate(nuEffa)
       *mesh.magSf()*fvc::snGrad(alpha)/fvc::interpolate(alpha + scalar(0.001));

    UaEqn =
    (

        (scalar(1) + Cvm*rhob*beta/rhoa)*
        (
            fvm::ddt(Ua)
          + fvm::div(phia, Ua, "div(phia,Ua)")
          - fvm::Sp(fvc::div(phia), Ua)
        )

      - fvm::laplacian(nuEffa, Ua)
      
      + fvm::div(phiRa, Ua, "div(phia,Ua)")
      - fvm::Sp(fvc::div(phiRa), Ua)
      ==
      - beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb)


    );

    UaEqn.relax();


    volTensorField Rcb = -nuEffb*fvc::grad(Ub)().T();
    Rcb = Rcb + (2.0/3.0)*I*k - (2.0/3.0)*I*tr(Rcb);

    surfaceScalarField phiRb =
       - fvc::interpolate(nuEffb)
        *mesh.magSf()*fvc::snGrad(beta)/fvc::interpolate(beta + scalar(0.001));

    UbEqn =
    (
        (scalar(1) + Cvm*rhob*alpha/rhob)*
        (
            fvm::ddt(Ub)
          + fvm::div(phib, Ub, "div(phib,Ub)")
          - fvm::Sp(fvc::div(phib), Ub)
        )

      - fvm::laplacian(nuEffb, Ub)
      
      + fvm::div(phiRb, Ub, "div(phib,Ub)")
      - fvm::Sp(fvc::div(phiRb), Ub)

      ==
      
      
      + alpha/rhob*(liftCoeff + Cvm*rhob*DDtUa)
    );

    UbEqn.relax();

    volScalarField beta12 = alpha/rhoa*dragCoef ;
    volScalarField beta21 = alpha/rhob*dragCoef ;
    volScalarField denomA = UaEqn.A()*UbEqn.A() + UbEqn.A()*beta12 +UaEqn.A()*beta21 ;
    volScalarField denomB = UbEqn.A()*UaEqn.A() + UbEqn.A()*beta12 +UaEqn.A()*beta21 ;

    Ua = ((UbEqn.A()+beta21)*(UaEqn.H()- fvc::div(Rca)-(fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca) + g) 
          + beta12 * (UbEqn.H() -  fvc::div(Rcb)-(fvc::grad(beta)/(fvc::average(beta) + scalar(0.001)) & Rcb) + g) ) / denomA ;

    Ub = ((beta21)*(UaEqn.H()- fvc::div(Rca)-(fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca)+g) 
          + (UaEqn.A() + beta12) * (UbEqn.H() -  fvc::div(Rcb)-(fvc::grad(beta)/(fvc::average(beta) + scalar(0.001)) & Rcb)+g) ) / denomB ;


}
Its still premature ... hence a bit messy ... sry for that ...

Last edited by balkrishna; September 23, 2010 at 04:08. Reason: Improper description of solution method ...
balkrishna is offline   Reply With Quote

Old   September 23, 2010, 04:01
Default
  #16
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 367
Rep Power: 13
l_r_mcglashan will become famous soon enough
The equations are not in that link. Best to check those before code.

You're not solving the equations simultaneously, you're eliminating the other phase's velocity from each velocity equation.
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   September 23, 2010, 04:06
Default
  #17
Senior Member
 
Balkrishna Patankar
Join Date: Mar 2009
Location: Pune
Posts: 119
Rep Power: 7
balkrishna is on a distinguished road
I have updated the link ....
balkrishna is offline   Reply With Quote

Old   September 27, 2010, 04:11
Default
  #18
Senior Member
 
Balkrishna Patankar
Join Date: Mar 2009
Location: Pune
Posts: 119
Rep Power: 7
balkrishna is on a distinguished road
I am getting the following error ... How can it be resolved ??
<code>
--> FOAM FATAL ERROR:
incompatible fields for operation
[Ubmod] == [Ub]

From function checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)
in file /home/ifmg/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude/fvMatrix.C at line 1181.

FOAM aborting

#0 Foam::error:rintStack(Foam::Ostream&) in "/home/ifmg/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/home/ifmg/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libOpenFOAM.so"
#2

#3

#4 __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#5

Aborted
</code>
When is the function checkMethod executed ??
balkrishna is offline   Reply With Quote

Old   October 19, 2012, 04:08
Default
  #19
Member
 
Join Date: Jun 2011
Posts: 35
Rep Power: 5
mikeP is on a distinguished road
Not very relevant to the original poster's problem, but I'd like to have make sure about something in bubblefoam momentum equation.

Code:
 surfaceScalarField phiRa =       - fvc::interpolate(nuEffa)        *mesh.magSf()*fvc::snGrad(alpha)/fvc::interpolate(alpha +         scalar(0.001));
As I understand, scalar(0.001) is there to avoid division by zero. Is that correct?

If yes, then why not VSMALL but a random scalar number?
mikeP is offline   Reply With Quote

Old   October 19, 2012, 04:08
Default
  #20
Member
 
Join Date: Jun 2011
Posts: 35
Rep Power: 5
mikeP is on a distinguished road
Not very relevant to the original poster's problem, but I'd like to have make sure about something in bubblefoam momentum equation.

Code:
 surfaceScalarField phiRa =
       - fvc::interpolate(nuEffa)
       *mesh.magSf()*fvc::snGrad(alpha)/fvc::interpolate(alpha + scalar(0.001));
As I understand, scalar(0.001) is there to avoid division by zero. Is that correct?

If yes, then why not VSMALL but a random scalar number?
mikeP is offline   Reply With Quote

Reply

Tags
bubblefoam, equations, momentum

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
BlockMesh FOAM warning gaottino OpenFOAM Native Meshers: blockMesh 7 July 19, 2010 14:11
Import problem ARC Open Source Meshers: Gmsh, Netgen, CGNS, ... 0 February 27, 2010 10:56
Derivation of Momentum Equation in Integral Form Demonwolf Main CFD Forum 2 October 29, 2009 20:53
Axisymmetrical mesh Rasmus Gjesing (Gjesing) OpenFOAM Native Meshers: blockMesh 10 April 2, 2007 14:00
Import gmsh msh to Foam adorean Open Source Meshers: Gmsh, Netgen, CGNS, ... 24 April 27, 2005 08:19


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