arsalan.dryi |
April 4, 2018 05:10 |
Quote:
Originally Posted by Taataa
(Post 687089)
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.
|