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

Calculate center position and velocity of two bubbles

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

Like Tree2Likes
  • 1 Post By piu58
  • 1 Post By Taataa

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 25, 2016, 10:32
Default Calculate center position and velocity of two bubbles
  #1
Member
 
Arsalan
Join Date: Jul 2014
Posts: 74
Rep Power: 10
arsalan.dryi is on a distinguished road
Hi Foamers,

I'm doing a 3D simulation of two and three bubble rising using a modified interFoam solver and I need to center position, velocity and surface area of the bubbles.

For a single bubble rising I used swak4Foam expressions for example for center position in Y as follows :

Code:
    bubbleCentreY
    {
        type swakExpression;
        valueType internalField;
        verbose true;
    variables (
    "Vol= sum (alpha1 < 0.5 ? vol() : 0);"
    "VolY= sum (alpha1 < 0.5 ? pos().y*vol() : 0);"    
    );        
    expression "VolY/Vol";
        accumulations (
        min        
        );

     
    }
Is there a way to compute these for two or three bubbles in this manner?

Thanks in advance,
Best Regards,
Arsalan.
arsalan.dryi is offline   Reply With Quote

Old   May 26, 2016, 12:01
Default
  #2
Member
 
Arsalan
Join Date: Jul 2014
Posts: 74
Rep Power: 10
arsalan.dryi is on a distinguished road
Does anybody have an idea?

Any help and comment will be greatly appreciated.

Regards,
Arsalan.
arsalan.dryi is offline   Reply With Quote

Old   June 8, 2016, 12:06
Default
  #3
Member
 
Arsalan
Join Date: Jul 2014
Posts: 74
Rep Power: 10
arsalan.dryi is on a distinguished road
I still didn't solve this problem, does anyone have an idea?!

Thanks in advance,
Regards.
arsalan.dryi is offline   Reply With Quote

Old   March 29, 2018, 17:32
Default
  #4
New Member
 
Lionel GAMET
Join Date: Nov 2013
Location: Lyon
Posts: 17
Rep Power: 11
x86_64leon is on a distinguished road
Hi Arsalan

Did you solve this problem of tracking more than one bubble and computing the center of mass, velocity, surface, volume ... etc of each of the bubbles ?

Regards
x86_64leon is offline   Reply With Quote

Old   March 30, 2018, 00:02
Default
  #5
Senior Member
 
Taher Chegini
Join Date: Nov 2014
Location: Houston, Texas
Posts: 125
Rep Power: 11
Taataa is on a distinguished road
One way is to use a multiphase flow solver like multiphaseInterFoam and for each bubble use a different name but same physical parameters. This way you can track each phase using the method OP mentioned or use this function that I have written for a wedge mesh which works in parallel as well:

Code:
functions
{
  riseVelocity
  {
    libs            ("libutilityFunctionObjects.so");
    
    type coded;
    
    name riseVelocity;
    
    codeExecute
    #{
        const volScalarField& alpha =
            mesh().lookupObject<volScalarField>("alpha.water");
        const volVectorField& U =
            mesh().lookupObject<volVectorField>("U");
        const volVectorField& centers = mesh().C();
        const scalarField& vols = mesh().V();
        
        const objectRegistry& db = mesh().thisDb();
        scalar ts = db.time().value();
        const dictionary& transportProperties =
            db.lookupObject<IOdictionary>
            (
                "transportProperties"
            );

        dictionary alpha2
        (
            transportProperties.subDict("air")
        );
        
        dimensionedScalar rhoBubble
        (
            "rhoBubble",
            dimDensity,
            alpha2.lookup("rho")
        );
        
        // threshold for interface
        scalar th = 0.5;
        // water level in the well
        scalar ht = 0.3;
        // sum of mass
        scalar sumMass = 0;
        // sum of volume
        scalar sumVol = 0;
        // sum of y by mass
        scalar sumX = 0;
        scalar sumY = 0;
        // sum of Uy by volume
        scalar sumUY = 0;
        //volume of liquid in well
        scalar vl = 0;
        
        forAll(centers, I)
        {
            if ( alpha[I] < th && centers[I].y() < ht )
            {
                //Mass center
                sumY    += (1 - alpha[I])*rhoBubble.value()*vols[I]*centers[I].y();
                sumX    += (1 - alpha[I])*rhoBubble.value()*vols[I]*centers[I].x();
                sumMass += (1 - alpha[I])*rhoBubble.value()*vols[I];
	              
	        //Rising Velocity
                sumUY   += (1 - alpha[I])*vols[I]*U[I].y();
                sumVol  += (1 - alpha[I])*vols[I];
            } else if ( centers[I].y() < ht ) {
                vl += alpha[I]*vols[I];
            }
        }
        
        reduce(sumY, sumOp<scalar>());
        reduce(sumX, sumOp<scalar>());
        reduce(sumMass, sumOp<scalar>());
        reduce(sumUY, sumOp<scalar>());
        reduce(sumVol, sumOp<scalar>());
        reduce(vl, sumOp<scalar>());
        
        // center of gravity
        scalar cgX = sumX/max(sumMass, SMALL);
        scalar cgY = sumY/max(sumMass, SMALL);
        
        // rising velocity
        scalar vr = sumUY/max(sumVol, SMALL);
        
        const volScalarField& p =
            mesh().lookupObject<volScalarField>("p");
        const faceList& ff = mesh().faces();
        const pointField& pp = mesh().points();
        // Pressure inside the bubble
        scalar pBubble = 0;
        
        forAll(centers, I)
        {
            const cell& cc = mesh().cells()[I];
            labelList pLabels(cc.labels(ff));
            pointField pLocal(pLabels.size(), vector::zero);

            forAll (pLabels, J)
            {
                pLocal[J] = pp[pLabels[J]];
            }
            
            if
            (
              sign(Foam::max(pLocal& vector(1,0,0)) - cgX) == 1 &&
              sign(cgX - Foam::min(pLocal& vector(1,0,0))) == 1 &&
              sign(Foam::max(pLocal& vector(0,1,0)) - cgY) == 1 &&
              sign(cgY - Foam::min(pLocal& vector(0,1,0))) == 1
            )
            {
                pBubble = mag(p[I]);
                break;
            }
        }
        reduce(pBubble, sumOp<scalar>());
        
        if (Pstream::master())
        {
            const fileName& outputFile = "plots/riseUCG";
            OFstream os
            (
                outputFile,
                IOstream::ASCII,
                IOstream::currentVersion,
                IOstream::UNCOMPRESSED,
                true
            );
            
            scalar tot = 360.0/5.0;
            Info<< "Writing rising velocity, CG, "
                << "volume of bubble, volume of liquid "
                << "and pressure inside the bubble for current time" << endl;
            os  << ts << "  " << vr << "  " << cgX  << "  " << cgY 
                << "  " << sumVol*tot << "  " << vl*tot << "  " << pBubble << endl;
        }
    #};
  }
}
It can easily be extended to multiphase e.g. for each phase add a new lookup, say alpha1 for alpha.water, alpha2 for alpha.air1 etc..
Taataa is offline   Reply With Quote

Old   March 30, 2018, 02:41
Default
  #6
Senior Member
 
piu58's Avatar
 
Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 727
Rep Power: 14
piu58 is on a distinguished road
I solved this problem by a small external program which analyses the alpha value for each cell. It is very easy to calculate the values you need form that.
__________________
Uwe Pilz
--
Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950)
piu58 is offline   Reply With Quote

Old   April 4, 2018, 05:10
Default
  #7
Member
 
Arsalan
Join Date: Jul 2014
Posts: 74
Rep Power: 10
arsalan.dryi is on a distinguished road
Quote:
Originally Posted by Taataa View Post
One way is to use a multiphase flow solver like multiphaseInterFoam and for each bubble use a different name but same physical parameters. This way you can track each phase using the method OP mentioned or use this function that I have written for a wedge mesh which works in parallel as well:

Code:
functions
{
  riseVelocity
  {
    libs            ("libutilityFunctionObjects.so");
    
    type coded;
    
    name riseVelocity;
    
    codeExecute
    #{
        const volScalarField& alpha =
            mesh().lookupObject<volScalarField>("alpha.water");
        const volVectorField& U =
            mesh().lookupObject<volVectorField>("U");
        const volVectorField& centers = mesh().C();
        const scalarField& vols = mesh().V();
        
        const objectRegistry& db = mesh().thisDb();
        scalar ts = db.time().value();
        const dictionary& transportProperties =
            db.lookupObject<IOdictionary>
            (
                "transportProperties"
            );

        dictionary alpha2
        (
            transportProperties.subDict("air")
        );
        
        dimensionedScalar rhoBubble
        (
            "rhoBubble",
            dimDensity,
            alpha2.lookup("rho")
        );
        
        // threshold for interface
        scalar th = 0.5;
        // water level in the well
        scalar ht = 0.3;
        // sum of mass
        scalar sumMass = 0;
        // sum of volume
        scalar sumVol = 0;
        // sum of y by mass
        scalar sumX = 0;
        scalar sumY = 0;
        // sum of Uy by volume
        scalar sumUY = 0;
        //volume of liquid in well
        scalar vl = 0;
        
        forAll(centers, I)
        {
            if ( alpha[I] < th && centers[I].y() < ht )
            {
                //Mass center
                sumY    += (1 - alpha[I])*rhoBubble.value()*vols[I]*centers[I].y();
                sumX    += (1 - alpha[I])*rhoBubble.value()*vols[I]*centers[I].x();
                sumMass += (1 - alpha[I])*rhoBubble.value()*vols[I];
	              
	        //Rising Velocity
                sumUY   += (1 - alpha[I])*vols[I]*U[I].y();
                sumVol  += (1 - alpha[I])*vols[I];
            } else if ( centers[I].y() < ht ) {
                vl += alpha[I]*vols[I];
            }
        }
        
        reduce(sumY, sumOp<scalar>());
        reduce(sumX, sumOp<scalar>());
        reduce(sumMass, sumOp<scalar>());
        reduce(sumUY, sumOp<scalar>());
        reduce(sumVol, sumOp<scalar>());
        reduce(vl, sumOp<scalar>());
        
        // center of gravity
        scalar cgX = sumX/max(sumMass, SMALL);
        scalar cgY = sumY/max(sumMass, SMALL);
        
        // rising velocity
        scalar vr = sumUY/max(sumVol, SMALL);
        
        const volScalarField& p =
            mesh().lookupObject<volScalarField>("p");
        const faceList& ff = mesh().faces();
        const pointField& pp = mesh().points();
        // Pressure inside the bubble
        scalar pBubble = 0;
        
        forAll(centers, I)
        {
            const cell& cc = mesh().cells()[I];
            labelList pLabels(cc.labels(ff));
            pointField pLocal(pLabels.size(), vector::zero);

            forAll (pLabels, J)
            {
                pLocal[J] = pp[pLabels[J]];
            }
            
            if
            (
              sign(Foam::max(pLocal& vector(1,0,0)) - cgX) == 1 &&
              sign(cgX - Foam::min(pLocal& vector(1,0,0))) == 1 &&
              sign(Foam::max(pLocal& vector(0,1,0)) - cgY) == 1 &&
              sign(cgY - Foam::min(pLocal& vector(0,1,0))) == 1
            )
            {
                pBubble = mag(p[I]);
                break;
            }
        }
        reduce(pBubble, sumOp<scalar>());
        
        if (Pstream::master())
        {
            const fileName& outputFile = "plots/riseUCG";
            OFstream os
            (
                outputFile,
                IOstream::ASCII,
                IOstream::currentVersion,
                IOstream::UNCOMPRESSED,
                true
            );
            
            scalar tot = 360.0/5.0;
            Info<< "Writing rising velocity, CG, "
                << "volume of bubble, volume of liquid "
                << "and pressure inside the bubble for current time" << endl;
            os  << ts << "  " << vr << "  " << cgX  << "  " << cgY 
                << "  " << sumVol*tot << "  " << vl*tot << "  " << pBubble << endl;
        }
    #};
  }
}
It can easily be extended to multiphase e.g. for each phase add a new lookup, say alpha1 for alpha.water, alpha2 for alpha.air1 etc..
Multiphase solvers need a surface tension between each two different phases and the results can be completely wrong especially when bubbles get close to each other, for example in bubbles coalescence.
arsalan.dryi is offline   Reply With Quote

Old   April 4, 2018, 05:13
Default
  #8
Member
 
Arsalan
Join Date: Jul 2014
Posts: 74
Rep Power: 10
arsalan.dryi is on a distinguished road
Quote:
Originally Posted by piu58 View Post
I solved this problem by a small external program which analyses the alpha value for each cell. It is very easy to calculate the values you need form that.
Could you explain more about that, please? And also post your written code here if possible.
Thanks.
arsalan.dryi is offline   Reply With Quote

Old   April 4, 2018, 07:32
Default
  #9
Senior Member
 
piu58's Avatar
 
Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 727
Rep Power: 14
piu58 is on a distinguished road
Quote:
Originally Posted by arsalan.dryi View Post
Could you explain more about that, please? And also post your written code here if possible.
Thanks.
That is quite easy. You go row wise through the alpha field, calculate the x and y coordinates according the setting of your region. If the alpha value is below a limiting value ( I used 0.2 for discriminating air / water), you add the x and y coordinates to get finally the average of them. This average should be close to the center of mass.

I don't believe that my C program is of any use for you. It is tailored to my setting and comments are in German.
arsalan.dryi likes this.
__________________
Uwe Pilz
--
Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950)
piu58 is offline   Reply With Quote

Old   April 5, 2018, 00:05
Default
  #10
Senior Member
 
Taher Chegini
Join Date: Nov 2014
Location: Houston, Texas
Posts: 125
Rep Power: 11
Taataa is on a distinguished road
Quote:
Originally Posted by arsalan.dryi View Post
Multiphase solvers need a surface tension between each two different phases and the results can be completely wrong especially when bubbles get close to each other, for example in bubbles coalescence.
Oh, I didn't know interFoam is not a multiphase solver!
cdritselis78 likes this.
Taataa is offline   Reply With Quote

Old   March 28, 2022, 12:52
Default question
  #11
Senior Member
 
Farzad Faraji
Join Date: Nov 2019
Posts: 157
Rep Power: 5
farzadmech is on a distinguished road
Dear Arslan
Have you added these lines to controlDict? what other lines have you added? which library have you added in the libs section?
for exmample;
Code:
libs (
    "libgroovyBC.so"
) ;
Thanks,
Farzad

Quote:
Originally Posted by arsalan.dryi View Post
Hi Foamers,

I'm doing a 3D simulation of two and three bubble rising using a modified interFoam solver and I need to center position, velocity and surface area of the bubbles.

For a single bubble rising I used swak4Foam expressions for example for center position in Y as follows :

Code:
    bubbleCentreY
    {
        type swakExpression;
        valueType internalField;
        verbose true;
    variables (
    "Vol= sum (alpha1 < 0.5 ? vol() : 0);"
    "VolY= sum (alpha1 < 0.5 ? pos().y*vol() : 0);"    
    );        
    expression "VolY/Vol";
        accumulations (
        min        
        );

     
    }
Is there a way to compute these for two or three bubbles in this manner?

Thanks in advance,
Best Regards,
Arsalan.
farzadmech is offline   Reply With Quote

Reply

Tags
bubble centre, bubble velocity, swak4foam, vof

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Calculate Bubbles Centre position and Velocity arsalan.dryi OpenFOAM Post-Processing 2 May 21, 2016 06:28
DPM : Fluid velocity at particle position ? souria Fluent Multiphase 1 January 30, 2015 06:58
DPM UDF particle position using the macro P_POS(p)[i] dm2747 FLUENT 0 April 17, 2009 02:29
Wall velocity changing with position Abbas FLUENT 1 May 26, 2003 08:23
Combustion Convergence problems Art Stretton Phoenics 5 April 2, 2002 06:59


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