34 const scalar Foam::multiphaseMixture::convertToRad =
40 void Foam::multiphaseMixture::calcAlphas()
47 alphas_ += level*iter();
64 phases_(lookup(
"phases"),
phase::iNew(U, phi)),
65 refPhase_(*phases_.lookup(
word(lookup(
"refPhase")))),
97 zeroGradientFvPatchScalarField::typeName
100 sigmas_(lookup(
"sigmas")),
101 dimSigma_(1, 0, -2, 0, 0),
113 alphaTable_.add(iter());
126 for(++iter; iter != phases_.
end(); ++iter)
128 trho() += iter()*iter().rho();
141 for(++iter; iter != phases_.
end(); ++iter)
143 tmu() += iter()*iter().rho()*iter().nu();
157 for(++iter; iter != phases_.
end(); ++iter)
188 "surfaceTensionForce",
189 mesh_.time().timeName(),
195 "surfaceTensionForce",
206 const phase& alpha1 = iter1();
211 for(; iter2 != phases_.
end(); ++iter2)
213 const phase& alpha2 = iter2();
215 sigmaTable::const_iterator sigma =
216 sigmas_.find(interfacePair(alpha1, alpha2));
218 if (sigma == sigmas_.end())
220 FatalErrorIn(
"multiphaseMixture::surfaceTensionForce() const")
221 <<
"Cannot find interface " << interfacePair(alpha1, alpha2)
222 <<
" in list of sigma values"
252 mesh_.solutionDict().subDict(
"PISO").lookup(
"nAlphaSubCycles")
258 readLabel(mesh_.solutionDict().subDict(
"PISO").lookup(
"nAlphaCorr"))
263 Switch(mesh_.solutionDict().subDict(
"PISO").lookup(
"cycleAlpha"))
268 readScalar(mesh_.solutionDict().subDict(
"PISO").lookup(
"cAlpha"))
282 !(++alphaSubCycle).end();
286 rhoPhiSum += (runTime.
deltaT()/totalDeltaT)*rhoPhi_;
322 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
333 return nHatfv(alpha1, alpha2) & mesh_.Sf();
343 void Foam::multiphaseMixture::correctContactAngle
350 const volScalarField::GeometricBoundaryField& gbf
351 = refPhase_.boundaryField();
357 if (isA<alphaContactAngleFvPatchScalarField>(gbf[
patchi]))
360 refCast<const alphaContactAngleFvPatchScalarField>(gbf[
patchi]);
365 mesh_.Sf().boundaryField()[
patchi]
366 /mesh_.magSf().boundaryField()[
patchi];
376 "multiphaseMixture::correctContactAngle"
377 "(const phase& alpha1, const phase& alpha2, "
378 "fvPatchVectorFieldField& nHatb) const"
379 ) <<
"Cannot find interface " << interfacePair(alpha1, alpha2)
380 <<
"\n in table of theta properties for patch "
381 << acap.patch().name()
385 bool matched = (tp.key().first() == alpha1.
name());
387 scalar theta0 = convertToRad*tp().theta0(matched);
388 scalarField theta(boundary[patchi].size(), theta0);
390 scalar uTheta = tp().uTheta();
395 scalar thetaA = convertToRad*tp().thetaA(matched);
396 scalar thetaR = convertToRad*tp().thetaR(matched);
400 U_.boundaryField()[
patchi].patchInternalField()
401 - U_.boundaryField()[
patchi];
402 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
406 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch;
409 nWall /= (
mag(nWall) + SMALL);
415 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
429 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
437 nHatPatch = a*AfHatPatch + b*nHatPatch;
439 nHatPatch /= (
mag(nHatPatch) + deltaN_.
value());
453 correctContactAngle(alpha1, alpha2, tnHatfv().
boundaryField());
456 return -
fvc::div(tnHatfv & mesh_.Sf());
470 mesh_.time().timeName(),
481 tnearInt() =
max(tnearInt(),
pos(alphaf - 0.01)*
pos(0.99 - alphaf));
488 void Foam::multiphaseMixture::solveAlphas
491 const bool cycleAlpha,
495 static label nSolves=-1;
498 word alphaScheme(
"div(phi,alpha)");
499 word alphacScheme(
"div(phirb,alpha)");
508 mesh_.divScheme(alphaScheme)
513 phic =
min(cAlpha*phic,
max(phic));
517 phase* refPhasePtr = &refPhase_;
522 for(label i=0; i<nSolves%phases_.size(); i++)
526 refPhasePtr = &refPhaseIter();
529 phase& refPhase = *refPhasePtr;
534 rhoPhi_ = phi_*refPhase.
rho();
538 phase& alpha = iter();
540 if (&alpha == &refPhase)
continue;
550 phase& alpha2 = iter2();
552 if (&alpha2 == &alpha)
continue;
558 alphaEqn +=
fvm::div(phirb12, alpha, alphacScheme);
561 alphaEqn.solve(mesh_.solver(
"alpha"));
563 rhoPhi_ += alphaEqn.flux()*(alpha.
rho() - refPhase.
rho());
565 Info<< alpha.
name() <<
" volume fraction, min, max = "
567 <<
' ' <<
min(alpha).value()
568 <<
' ' <<
max(alpha).value()
571 refPhaseNew == refPhaseNew - alpha;
574 refPhase == refPhaseNew;
576 Info<< refPhase.
name() <<
" volume fraction, min, max = "
578 <<
' ' <<
min(refPhase).value()
579 <<
' ' <<
max(refPhase).value()
598 readOK &= iter().read(phaseData[phasei++].dict());
601 lookup(
"sigmas") >> sigmas_;