36#include "surfaceTensionModel.H"
42#include "phaseCompressibleTurbulenceModel.H"
46void 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])
78void 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(),
181void 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))
221void 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();
259void 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
309 Pair<word>::compare(pDmdt_.find(pairij).key(), pairij)
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
325 Pair<word>::compare(pDmdt_.find(pairik).key(), pairik)
328 pDmdt_[pairik]->ref() += dmdtSign*fk.x()/v*Sui_*fi.phase().rho();
334void 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();
353void 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
382 Pair<word>::compare(pDmdt_.find(pair).key(), pair)
385 pDmdt_[pair]->ref() += dmdtSign*Sui_*fi.phase().rho();
391void Foam::diameterModels::populationBalanceModel::deathByBreakup(
const label i)
393 const sizeGroup& fi = sizeGroups_[i];
395 SuSp_[i] += breakupRate_()*fi.phase();
399void 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;
436void 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
460 Pair<word>::compare(pDmdt_.find(pairij).key(), pairij)
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();
509void Foam::diameterModels::populationBalanceModel::
518 SuSp_[i] += alphai*binaryBreakupRate_()*delta_[j][i];
522void 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
608 Pair<word>::compare(pDmdt_.find(pairij).key(), pairij)
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
637 Pair<word>::compare(pDmdt_.find(pairih).key(), pairih)
640 pDmdt_[pairih]->ref() -= dmdtSign*Suw*fp.phase().rho();
646void Foam::diameterModels::populationBalanceModel::nucleation(
const label i)
648 Su_[i] += sizeGroups_[i].x()*nucleationRate_();
652void 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);
768void 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);
789void Foam::diameterModels::populationBalanceModel::calcAlphas()
793 forAll(velocityGroups_, v)
795 const phaseModel& phase = velocityGroups_[v].phase();
797 alphas_() +=
max(phase, phase.residualAlpha());
803Foam::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_());
828void Foam::diameterModels::populationBalanceModel::calcVelocity()
832 forAll(velocityGroups_, v)
834 const phaseModel& phase = velocityGroups_[v].phase();
836 U_() +=
max(phase, phase.residualAlpha())*phase.U()/alphas_();
840bool 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")),
881 IOobject::groupName(
"alpha", dict_.get<
word>(
"continuousPhase"))
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;
1283 phase.alphaRhoPhi(),
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 = "
dimensioned< Type > weightedAverage(const DimensionedField< scalar, GeoMesh > &weightField) const
Calculate and return weighted average.
Internal & ref(const bool updateAccessTime=true)
Return a reference to the dimensioned internal field.
A HashTable of pointers to objects of type <T>, with deallocation management of the pointers.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
const word & name() const noexcept
Return the object name.
const Time & time() const
Return Time associated with the objectRegistry.
static word groupName(StringType base, const word &group)
Create dot-delimited name.group string.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
static word timeName(const scalar t, const int precision=precision_)
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Base class for binary breakup models which give the breakup rate between a sizeGroup pair directly,...
Base class for breakup models which give a total breakup rate and a separate daughter size distributi...
Base class for coalescence models.
Constant dispersed-phase particle diameter model.
Base class for drift models.
Base class for nucleation models.
Return a pointer to a new populationBalanceModel object created on.
Class that solves the univariate population balance equation by means of a class method (also called ...
bool writeData(Ostream &) const
Dummy write for regIOobject.
const dimensionedScalar gamma(const label i, const dimensionedScalar &v) const
Return allocation coefficient.
autoPtr< populationBalanceModel > clone() const
Return clone.
friend class velocityGroup
const phaseCompressibleTurbulenceModel & continuousTurbulence() const
Return reference to turbulence model of the continuous phase.
const tmp< volScalarField > sigmaWithContinuousPhase(const phaseModel &dispersedPhase) const
Return the surface tension coefficient between a given dispersed.
virtual ~populationBalanceModel()
Destructor.
void solve()
Solve the population balance equation.
This class represents a single sizeGroup belonging to a velocityGroup. The main property of a sizeGro...
const velocityGroup & VelocityGroup() const
Return const-reference to the velocityGroup.
const phaseModel & phase() const
Return const-reference to the phase.
const volScalarField & dmdt() const
Return const-reference to the mass transfer rate.
const tmp< fv::convectionScheme< scalar > > & mvConvection() const
Return const-reference to multivariate convectionScheme.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
const Type & value() const
Return const reference to value.
A special matrix type and solver, designed for finite volume solutions of scalar equations....
void relax(const scalar alpha)
Relax matrix (for steady-state solution).
SolverPerformance< Type > solve(const dictionary &)
Solve returning the solution statistics.
const Time & time() const
Return the top-level database.
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...
Class to represent a system of phases and model interfacial transfers between them.
const phaseModelList & phases() const
Return the phase models.
Single incompressible phase derived from the phase-fraction. Used as part of the multiPhaseMixture fo...
const dimensionedScalar & rho() const
Return const-access to phase1 density.
PIMPLE control class to supply convergence information/checks for the PIMPLE loop.
Lookup type of boundary radiation properties.
regIOobject is an abstract class derived from IOobject to handle automatic object registration with t...
A class for managing temporary objects.
static const word propertiesName
Default name of the turbulence properties dictionary.
A class for handling words, derived from Foam::string.
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
Calculate the first temporal derivative.
Calculate the divergence of the given field.
Calculate the field for explicit evaluation of implicit and explicit sources.
Calculate the matrix for the first temporal derivative.
Calculate the finiteVolume matrix for implicit and explicit sources.
constexpr const char *const group
Group name for atomic constants.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
zeroField SuSp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
zeroField Sp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
dimensionedScalar pos(const dimensionedScalar &ds)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
const dimensionSet dimless
Dimensionless.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
GeometricField< scalar, fvPatchField, volMesh > volScalarField
const dimensionSet dimVelocity
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
dimensionedScalar neg(const dimensionedScalar &ds)
static constexpr const zero Zero
Global zero (0)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
const dimensionSet dimVolume(pow3(dimLength))
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
errorManipArg< error, int > exit(error &err, const int errNo=1)
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object.