36 #include "surfaceTensionModel.H" 42 #include "phaseCompressibleTurbulenceModel.H" 46 void Foam::diameterModels::populationBalanceModel::registerVelocityGroups()
50 if (isA<velocityGroup>(fluid_.
phases()[
phasei].dPtr()()))
53 refCast<const velocityGroup>(fluid_.
phases()[
phasei].dPtr()());
55 if (velGroup.popBalName() == this->
name())
57 velocityGroups_.resize(velocityGroups_.size() + 1);
61 velocityGroups_.size() - 1,
65 forAll(velGroup.sizeGroups(), i)
67 this->registerSizeGroups
69 const_cast<sizeGroup&>(velGroup.sizeGroups()[i])
78 void Foam::diameterModels::populationBalanceModel::registerSizeGroups
85 sizeGroups_.size() != 0
87 group.x().value() <= sizeGroups_.last().x().value()
91 <<
"Size groups must be entered according to their representative" 96 sizeGroups_.resize(sizeGroups_.size() + 1);
97 sizeGroups_.set(sizeGroups_.size() - 1, &
group);
100 if (sizeGroups_.size() == 1)
109 sizeGroups_.last().x()
120 sizeGroups_.last().x()
131 sizeGroups_[sizeGroups_.size()-2].x()
132 + sizeGroups_.last().x()
142 sizeGroups_.last().x()
147 delta_.append(
new PtrList<dimensionedScalar>());
156 fluid_.time().timeName(),
171 fluid_.time().timeName(),
181 void Foam::diameterModels::populationBalanceModel::createPhasePairs()
183 forAll(velocityGroups_, i)
185 const phaseModel&
phasei = velocityGroups_[i].phase();
187 forAll(velocityGroups_, j)
189 const phaseModel& phasej = velocityGroups_[j].phase();
193 const phasePairKey
key 200 if (!phasePairs_.found(
key))
221 void Foam::diameterModels::populationBalanceModel::correct()
225 forAll(velocityGroups_, v)
227 velocityGroups_[v].preSolve();
230 forAll(coalescence_, model)
232 coalescence_[model].correct();
237 breakup_[model].correct();
239 breakup_[model].dsdPtr()().
correct();
242 forAll(binaryBreakup_, model)
244 binaryBreakup_[model].correct();
249 drift_[model].correct();
252 forAll(nucleation_, model)
254 nucleation_[model].correct();
259 void Foam::diameterModels::populationBalanceModel::
266 const sizeGroup& fj = sizeGroups_[j];
267 const sizeGroup& fk = sizeGroups_[
k];
272 for (label i = j; i < sizeGroups_.size(); i++)
277 if (Gamma.value() == 0)
continue;
279 const sizeGroup& fi = sizeGroups_[i];
285 0.5*fi.x()*coalescenceRate_()*Gamma
286 *fj*fj.phase()/fj.x()
287 *fk*fk.phase()/fk.x();
292 fi.x()*coalescenceRate_()*Gamma
293 *fj*fj.phase()/fj.x()
294 *fk*fk.phase()/fk.x();
299 const phasePairKey pairij
305 if (pDmdt_.found(pairij))
307 const scalar dmdtSign
312 pDmdt_[pairij]->ref() += dmdtSign*fj.x()/v*Sui_*fi.phase().rho();
315 const phasePairKey pairik
321 if (pDmdt_.found(pairik))
323 const scalar dmdtSign
328 pDmdt_[pairik]->ref() += dmdtSign*fk.x()/v*Sui_*fi.phase().rho();
334 void Foam::diameterModels::populationBalanceModel::
341 const sizeGroup& fi = sizeGroups_[i];
342 const sizeGroup& fj = sizeGroups_[j];
344 SuSp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x();
348 SuSp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x();
353 void Foam::diameterModels::populationBalanceModel::
360 const sizeGroup& fk = sizeGroups_[
k];
362 for (label i = 0; i <=
k; i++)
364 const sizeGroup& fi = sizeGroups_[i];
367 fi.x()*breakupRate_()*breakup_[model].dsdPtr()().nik(i,
k)
368 *fk*fk.phase()/fk.x();
372 const phasePairKey pair
378 if (pDmdt_.found(pair))
380 const scalar dmdtSign
385 pDmdt_[pair]->ref() += dmdtSign*Sui_*fi.phase().rho();
391 void Foam::diameterModels::populationBalanceModel::deathByBreakup(
const label i)
393 const sizeGroup& fi = sizeGroups_[i];
395 SuSp_[i] += breakupRate_()*fi.phase();
399 void Foam::diameterModels::populationBalanceModel::calcDeltas()
403 if (delta_[i].empty())
405 for (label j = 0; j <= sizeGroups_.size() - 1; j++)
413 v_[i+1].value() - v_[i].value()
419 v_[i].value() < 0.5*sizeGroups_[j].
x().value()
421 0.5*sizeGroups_[j].
x().value() < v_[i+1].value()
424 delta_[i][j] =
mag(0.5*sizeGroups_[j].
x() - v_[i]);
426 else if (0.5*sizeGroups_[j].
x().value() <= v_[i].value())
428 delta_[i][j].value() = 0;
436 void Foam::diameterModels::populationBalanceModel::
443 const sizeGroup& fj = sizeGroups_[j];
444 const sizeGroup& fi = sizeGroups_[i];
446 Sui_ = fi.x()*binaryBreakupRate_()*delta_[i][j]*fj*fj.phase()/fj.x();
450 const phasePairKey pairij
456 if (pDmdt_.found(pairij))
458 const scalar dmdtSign
463 pDmdt_[pairij]->ref() += dmdtSign*Sui_*fi.phase().rho();
469 for (label
k = 0;
k <= j;
k++)
474 if (Gamma.value() == 0)
continue;
476 const sizeGroup& fk = sizeGroups_[
k];
481 sizeGroups_[
k].x()*binaryBreakupRate_()*delta_[i][j]*Gamma
482 *fj*fj.phase()/fj.x();
486 const phasePairKey pairkj
492 if (pDmdt_.found(pairkj))
494 const scalar dmdtSign
498 pDmdt_.find(pairkj).key(),
503 pDmdt_[pairkj]->ref() += dmdtSign*Suk*fi.phase().rho();
509 void Foam::diameterModels::populationBalanceModel::
518 SuSp_[i] += alphai*binaryBreakupRate_()*delta_[j][i];
522 void Foam::diameterModels::populationBalanceModel::drift(
const label i)
524 const sizeGroup& fp = sizeGroups_[i];
528 rx_() =
pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x();
530 else if (i == sizeGroups_.size() - 1)
532 rx_() =
neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
537 pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x()
538 +
neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
542 (
neg(1 - rx_()) +
neg(rx_() - rx_()/(1 - rx_())))*driftRate_()
543 *fp.phase()/((rx_() - 1)*fp.x());
548 if (i == sizeGroups_.size() - 2)
550 rx_() =
pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x();
554 *(sizeGroups_[i+1].x() - sizeGroups_[i].x())
555 /(sizeGroups_[i].
x() - sizeGroups_[i-1].x());
557 else if (i < sizeGroups_.size() - 2)
559 rx_() =
pos(driftRate_())*sizeGroups_[i+2].x()/sizeGroups_[i+1].x();
563 *(sizeGroups_[i+2].x() - sizeGroups_[i+1].x())
564 /(sizeGroups_[i+1].
x() - sizeGroups_[i].x());
569 rx_() +=
neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
573 *(sizeGroups_[i].x() - sizeGroups_[i-1].x())
574 /(sizeGroups_[i+1].
x() - sizeGroups_[i].x());
578 rx_() +=
neg(driftRate_())*sizeGroups_[i-2].x()/sizeGroups_[i-1].x();
582 *(sizeGroups_[i-1].x() - sizeGroups_[i-2].x())
583 /(sizeGroups_[i].
x() - sizeGroups_[i-1].x());
586 if (i != sizeGroups_.size() - 1)
588 const sizeGroup& fe = sizeGroups_[i+1];
592 pos(driftRate_())*driftRate_()*rdx_()
593 *fp*fp.phase()/fp.x()
598 const phasePairKey pairij
604 if (pDmdt_.found(pairij))
606 const scalar dmdtSign
611 pDmdt_[pairij]->ref() -= dmdtSign*Sue*fp.phase().rho();
617 const sizeGroup& fw = sizeGroups_[i-1];
621 neg(driftRate_())*driftRate_()*rdx_()
622 *fp*fp.phase()/fp.x()
627 const phasePairKey pairih
633 if (pDmdt_.found(pairih))
635 const scalar dmdtSign
640 pDmdt_[pairih]->ref() -= dmdtSign*Suw*fp.phase().rho();
646 void Foam::diameterModels::populationBalanceModel::nucleation(
const label i)
648 Su_[i] += sizeGroups_[i].x()*nucleationRate_();
652 void Foam::diameterModels::populationBalanceModel::sources()
667 pDmdt_(phasePairIter())->ref() *= 0.0;
675 const sizeGroup& fi = sizeGroups_[i];
677 if (coalescence_.size() != 0)
679 for (label j = 0; j <= i; j++)
681 const sizeGroup& fj = sizeGroups_[j];
683 if (fi.x() + fj.x() > sizeGroups_.last().x())
break;
685 coalescenceRate_() *= 0.0;
687 forAll(coalescence_, model)
689 coalescence_[model].addToCoalescenceRate
697 birthByCoalescence(i, j);
699 deathByCoalescence(i, j);
703 if (breakup_.size() != 0)
707 breakup_[model].setBreakupRate(breakupRate_(), i);
709 birthByBreakup(i, model);
715 if (binaryBreakup_.size() != 0)
719 while (delta_[j][i].value() != 0)
721 binaryBreakupRate_() *= 0.0;
723 forAll(binaryBreakup_, model)
725 binaryBreakup_[model].addToBinaryBreakupRate
727 binaryBreakupRate_(),
733 birthByBinaryBreakup(j, i);
735 deathByBinaryBreakup(j, i);
741 if (drift_.size() != 0)
747 drift_[model].addToDriftRate(driftRate_(), i);
753 if (nucleation_.size() != 0)
755 nucleationRate_() *= 0.0;
757 forAll(nucleation_, model)
759 nucleation_[model].addToNucleationRate(nucleationRate_(), i);
768 void Foam::diameterModels::populationBalanceModel::dmdt()
770 forAll(velocityGroups_, v)
772 velocityGroup& velGroup = velocityGroups_[v];
774 velGroup.dmdtRef() *= 0.0;
778 if (&sizeGroups_[i].phase() == &velGroup.phase())
780 sizeGroup& fi = sizeGroups_[i];
782 velGroup.dmdtRef() += fi.phase().rho()*(Su_[i] - SuSp_[i]*fi);
789 void Foam::diameterModels::populationBalanceModel::calcAlphas()
793 forAll(velocityGroups_, v)
795 const phaseModel& phase = velocityGroups_[v].phase();
797 alphas_() +=
max(phase, phase.residualAlpha());
803 Foam::diameterModels::populationBalanceModel::calcDsm()
805 tmp<volScalarField> tInvDsm
817 forAll(velocityGroups_, v)
819 const phaseModel& phase = velocityGroups_[v].phase();
821 invDsm +=
max(phase, phase.residualAlpha())/(phase.d()*alphas_());
828 void Foam::diameterModels::populationBalanceModel::calcVelocity()
832 forAll(velocityGroups_, v)
834 const phaseModel& phase = velocityGroups_[v].phase();
836 U_() +=
max(phase, phase.residualAlpha())*phase.U()/alphas_();
840 bool Foam::diameterModels::populationBalanceModel::updateSources()
842 const bool result = sourceUpdateCounter_ % sourceUpdateInterval() == 0;
844 ++ sourceUpdateCounter_;
870 mesh_(fluid_.
mesh()),
874 fluid_.subDict(
"populationBalanceCoeffs").subDict(name_)
876 pimple_(mesh_.lookupObject<
pimpleControl>(
"solutionControl")),
903 dict_.
lookup(
"coalescenceModels"),
909 dict_.
lookup(
"breakupModels"),
915 dict_.
lookup(
"binaryBreakupModels"),
918 binaryBreakupRate_(),
921 dict_.
lookup(
"driftModels"),
929 dict_.
lookup(
"nucleationModels"),
941 this->registerVelocityGroups();
943 this->createPhasePairs();
945 if (sizeGroups_.size() < 3)
948 <<
"The populationBalance " << name_
949 <<
" requires a minimum number of three sizeGroups to be specified." 954 if (coalescence_.size() != 0)
956 coalescenceRate_.reset
972 if (breakup_.size() != 0)
990 if (binaryBreakup_.size() != 0)
992 binaryBreakupRate_.reset
1005 "binaryBreakupRate",
1013 if (drift_.size() != 0)
1061 if (nucleation_.size() != 0)
1063 nucleationRate_.reset
1084 if (velocityGroups_.size() > 1)
1178 lowerBoundary = sizeGroups_[i-1].x();
1181 if (i == sizeGroups_.size() - 1)
1187 upperBoundary = sizeGroups_[i+1].x();
1190 if (v < lowerBoundary || v > upperBoundary)
1196 return (v - lowerBoundary)/(xi - lowerBoundary);
1200 return (upperBoundary - v)/(upperBoundary - xi);
1214 phasePair(dispersedPhase, continuousPhase_)
1228 continuousPhase_.name()
1236 const dictionary& solutionControls = mesh_.solverDict(name_);
1237 const bool solveOnFinalIterOnly =
1238 solutionControls.
getOrDefault(
"solveOnFinalIterOnly",
false);
1240 if (!solveOnFinalIterOnly || pimple_.finalIter())
1243 const scalar tolerance = solutionControls.
get<scalar>(
"tolerance");
1251 scalar maxInitialResidual = 1;
1253 while (++iCorr <= nCorr && maxInitialResidual > tolerance)
1255 Info<<
"populationBalance " 1261 if (updateSources())
1268 maxInitialResidual = 0;
1272 sizeGroup& fi = sizeGroups_[i];
1273 const phaseModel& phase = fi.phase();
1281 + fi.VelocityGroup().mvConvection()->fvmDiv
1283 phase.alphaRhoPhi(),
1289 - fi.VelocityGroup().dmdt(),
1299 sizeGroupEqn.relax();
1301 maxInitialResidual =
max 1303 sizeGroupEqn.solve().initialResidual(),
1311 forAll(velocityGroups_, i)
1313 velocityGroups_[i].postSolve();
1317 if (velocityGroups_.size() > 1)
1326 sizeGroups_.first()*sizeGroups_.first().phase()
1331 sizeGroups_.last()*sizeGroups_.last().phase()
1334 Info<< this->
name() <<
" sizeGroup phase fraction first, last = " 1335 << fAlpha0.weightedAverage(this->
mesh().V()).value()
1336 <<
' ' << fAlphaN.weightedAverage(this->
mesh().V()).value()
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
const Type & value() const noexcept
Return const reference to value.
fvMatrix< scalar > fvScalarMatrix
Base class for drift models.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
populationBalanceModel(const phaseSystem &fluid, const word &name, HashPtrTable< volScalarField, phasePairKey, phasePairKey::hash > &pDmdt)
const word & name() const noexcept
Return the object name.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const tmp< volScalarField > sigmaWithContinuousPhase(const phaseModel &dispersedPhase) const
Return the surface tension coefficient between a given dispersed.
label k
Boltzmann constant.
const dimensionSet dimless
Dimensionless.
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T. FatalIOError if not found, or if the number of tokens is incorrect.
bool writeData(Ostream &) const
Dummy write for regIOobject.
const Time & time() const
Return the top-level database.
Lookup type of boundary radiation properties.
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object.
GeometricField< vector, fvPatchField, volMesh > volVectorField
dimensionedScalar sigma("sigma", dimMass/sqr(dimTime), transportProperties)
dimensionedScalar neg(const dimensionedScalar &ds)
autoPtr< populationBalanceModel > clone() const
Return clone.
Base class for coalescence models.
#define forAll(list, i)
Loop across all elements in list.
constexpr const char *const group
Group name for atomic constants.
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Base class for breakup models which give a total breakup rate and a separate daughter size distributi...
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Calculate the first temporal derivative.
Base class for binary breakup models which give the breakup rate between a sizeGroup pair directly...
const phaseModelList & phases() const
Return the phase models.
dimensionedScalar pos(const dimensionedScalar &ds)
const dimensionSet dimVolume(pow3(dimLength))
static word groupName(StringType base, const word &group)
Create dot-delimited name.group string.
Base class for nucleation models.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
static int compare(const Pair< word > &a, const Pair< word > &b)
Compare Pairs.
Class to represent a system of phases and model interfacial transfers between them.
static const word propertiesName
Default name of the turbulence properties dictionary.
A class for handling words, derived from Foam::string.
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Calculate the matrix for the first temporal derivative.
zeroField SuSp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
Calculate the field for explicit evaluation of implicit and explicit sources.
Constant dispersed-phase particle diameter model.
zeroField Sp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
const Time & time() const noexcept
Return Time associated with the objectRegistry.
Calculate the divergence of the given field.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
static word timeName(const scalar t, const int precision=precision_)
Return time name of given scalar time formatted with the given precision.
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(psi *p+alphal *rhol0+((alphav *psiv+alphal *psil) - psi) *pSat, rhoMin);# 1 "/home/chef2/andy/OpenFOAM/release/v2306/OpenFOAM-v2306/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=clamp((rho - rholSat)/(rhovSat - rholSat), zero_one{});alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
OBJstream os(runTime.globalPath()/outputName)
friend class velocityGroup
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
PIMPLE control class to supply convergence information/checks for the PIMPLE loop.
Automatically write from objectRegistry::writeObject()
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
regIOobject is an abstract class derived from IOobject to handle automatic object registration with t...
messageStream Info
Information stream (stdout output on master, null elsewhere)
const phaseCompressibleTurbulenceModel & continuousTurbulence() const
Return reference to turbulence model of the continuous phase.
Internal & ref(const bool updateAccessTime=true)
Same as internalFieldRef()
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T, or return the given default value. FatalIOError if it is found and the number of...
A class for managing temporary objects.
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Mesh &mesh, const dimensionSet &dims, const word &patchFieldType=fvPatchField< scalar >::calculatedType())
Return tmp field from name, mesh, dimensions and patch type.
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
void solve()
Solve the population balance equation.
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Defines the attributes of an object for which implicit objectRegistry management is supported...
const dimensionedScalar gamma(const label i, const dimensionedScalar &v) const
Return allocation coefficient.
Single incompressible phase derived from the phase-fraction. Used as part of the multiPhaseMixture fo...
Description for mass transfer between a pair of phases. The direction of the mass transfer is from th...
Calculate the finiteVolume matrix for implicit and explicit sources.
virtual ~populationBalanceModel()
Destructor.
static constexpr const zero Zero
Global zero (0)
const dimensionSet dimVelocity