# Calculate center position and velocity of two bubbles

 Register Blogs Members List Search Today's Posts Mark Forums Read May 25, 2016, 10:32 Calculate center position and velocity of two bubbles #1 Member   Arsalan Join Date: Jul 2014 Posts: 74 Rep Power: 10 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.   May 26, 2016, 12:01 #2 Member   Arsalan Join Date: Jul 2014 Posts: 74 Rep Power: 10 Does anybody have an idea? Any help and comment will be greatly appreciated. Regards, Arsalan.   June 8, 2016, 12:06 #3 Member   Arsalan Join Date: Jul 2014 Posts: 74 Rep Power: 10 I still didn't solve this problem, does anyone have an idea?! Thanks in advance, Regards.   March 29, 2018, 17:32 #4 New Member   Lionel GAMET Join Date: Nov 2013 Location: Lyon Posts: 17 Rep Power: 11 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   March 30, 2018, 00:02 #5 Senior Member   Taher Chegini Join Date: Nov 2014 Location: Houston, Texas Posts: 125 Rep Power: 11 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("alpha.water"); const volVectorField& U = mesh().lookupObject("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 ( "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()); reduce(sumX, sumOp()); reduce(sumMass, sumOp()); reduce(sumUY, sumOp()); reduce(sumVol, sumOp()); reduce(vl, sumOp()); // 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("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()); 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..   March 30, 2018, 02:41 #6 Senior Member   Uwe Pilz Join Date: Feb 2017 Location: Leipzig, Germany Posts: 727 Rep Power: 14 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)   April 4, 2018, 05:10 #7
Member

Arsalan
Join Date: Jul 2014
Posts: 74
Rep Power: 10 Quote:
 Originally Posted by Taataa 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("alpha.water"); const volVectorField& U = mesh().lookupObject("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 ( "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()); reduce(sumX, sumOp()); reduce(sumMass, sumOp()); reduce(sumUY, sumOp()); reduce(sumVol, sumOp()); reduce(vl, sumOp()); // 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("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()); 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.   April 4, 2018, 05:13 #8
Member

Arsalan
Join Date: Jul 2014
Posts: 74
Rep Power: 10 Quote:
 Originally Posted by piu58 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.
Thanks.   April 4, 2018, 07:32 #9
Senior Member

Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 727
Rep Power: 14 Quote:
 Originally Posted by arsalan.dryi 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.
__________________
Uwe Pilz
--
Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950)   April 5, 2018, 00:05 #10
Senior Member

Taher Chegini
Join Date: Nov 2014
Location: Houston, Texas
Posts: 125
Rep Power: 11 Quote:
 Originally Posted by arsalan.dryi 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!   March 28, 2022, 12:52 question
#11
Senior Member

Join Date: Nov 2019
Posts: 157
Rep Power: 5 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,

Quote:
 Originally Posted by arsalan.dryi 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.  Tags bubble centre, bubble velocity, swak4foam, vof Thread Tools Search this Thread Show Printable Version Email this Page Search this Thread: Advanced Search Display Modes Linear Mode Switch to Hybrid Mode Switch to Threaded Mode 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 OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules Similar Threads Thread Thread Starter Forum Replies Last Post arsalan.dryi OpenFOAM Post-Processing 2 May 21, 2016 06:28 souria Fluent Multiphase 1 January 30, 2015 06:58 dm2747 FLUENT 0 April 17, 2009 02:29 Abbas FLUENT 1 May 26, 2003 08:23 Art Stretton Phoenics 5 April 2, 2002 06:59

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