46 namespace diameterModels
55 void Foam::diameterModels::populationBalanceModel::registerVelocityGroups()
59 if (isA<velocityGroup>(fluid_.
phases()[phasei].dPtr()()))
61 const velocityGroup& velGroup =
62 refCast<const velocityGroup>(fluid_.
phases()[phasei].dPtr()());
64 if (velGroup.popBalName() == this->name())
66 velocityGroupPtrs_.insert(velGroup.phase().name(), &velGroup);
68 dilatationErrors_.insert
70 velGroup.phase().name(),
78 velGroup.phase().name()
88 forAll(velGroup.sizeGroups(), i)
90 this->registerSizeGroups
92 const_cast<sizeGroup&
>(velGroup.sizeGroups()[i])
101 void Foam::diameterModels::populationBalanceModel::registerSizeGroups
108 sizeGroups().size() != 0
110 group.x().value() <= sizeGroups().last().
x().value()
114 <<
"Size groups must be entered according to their representative"
119 sizeGroups_.resize(sizeGroups().size() + 1);
120 sizeGroups_.set(sizeGroups().size() - 1, &
group);
122 if (sizeGroups().size() == 1)
129 sizeGroups().last().
x()
138 sizeGroups().last().
x()
147 sizeGroups()[sizeGroups().size()-2].x()
148 + sizeGroups().last().x()
156 sizeGroups().last().
x()
161 delta_.append(
new PtrList<dimensionedScalar>());
170 fluid_.time().name(),
185 fluid_.time().name(),
195 void Foam::diameterModels::populationBalanceModel::initialiseDmdtfs()
199 HashTable<const diameterModels::velocityGroup*>,
204 const diameterModels::velocityGroup& velGrp1 = *iter1();
208 HashTable<const diameterModels::velocityGroup*>,
213 const diameterModels::velocityGroup& velGrp2 = *iter2();
215 const phaseInterface interface(velGrp1.phase(), velGrp2.phase());
221 !dmdtfs_.found(interface)
224 fluid_.validateMassTransfer
225 <diameterModels::populationBalanceModel>(interface);
239 mesh().time().
name(),
252 void Foam::diameterModels::populationBalanceModel::precompute()
254 forAll(coalescenceModels_, model)
256 coalescenceModels_[model].precompute();
259 forAll(breakupModels_, model)
261 breakupModels_[model].precompute();
263 breakupModels_[model].dsdPtr()->precompute();
266 forAll(binaryBreakupModels_, model)
268 binaryBreakupModels_[model].precompute();
273 drift_[model].precompute();
276 forAll(nucleation_, model)
278 nucleation_[model].precompute();
283 void Foam::diameterModels::populationBalanceModel::birthByCoalescence
289 const sizeGroup& fj = sizeGroups()[j];
290 const sizeGroup& fk = sizeGroups()[
k];
295 for (
label i = j; i < sizeGroups().size(); i++)
299 if (Eta.value() == 0)
continue;
301 const sizeGroup& fi = sizeGroups()[i];
306 0.5*fi.x()/(fj.x()*fk.x())*Eta
307 *coalescenceRate_()*fj*fj.phase()*fk*fk.phase();
312 fi.x()/(fj.x()*fk.x())*Eta
313 *coalescenceRate_()*fj*fj.phase()*fk*fk.phase();
318 const phaseInterface interfaceij(fi.phase(), fj.phase());
320 if (dmdtfs_.found(interfaceij))
322 const scalar dmdtSign =
323 interfaceij.index(fi.phase()) == 0 ? +1 : -1;
325 *dmdtfs_[interfaceij] += dmdtSign*fj.x()/v*Sui_*fj.phase().rho();
328 const phaseInterface interfaceik(fi.phase(), fk.phase());
330 if (dmdtfs_.found(interfaceik))
332 const scalar dmdtSign =
333 interfaceik.index(fi.phase()) == 0 ? +1 : -1;
335 *dmdtfs_[interfaceik] += dmdtSign*fk.x()/v*Sui_*fk.phase().rho();
338 sizeGroups_[i].shapeModelPtr()->addCoalescence(Sui_, fj, fk);
343 void Foam::diameterModels::populationBalanceModel::deathByCoalescence
349 const sizeGroup& fi = sizeGroups()[i];
350 const sizeGroup& fj = sizeGroups()[j];
352 Sp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x();
356 Sp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x();
361 void Foam::diameterModels::populationBalanceModel::birthByBreakup
367 const sizeGroup& fk = sizeGroups()[
k];
369 for (
label i = 0; i <=
k; i++)
371 const sizeGroup& fi = sizeGroups()[i];
374 fi.x()*breakupModels_[model].dsdPtr()().nik(i,
k)/fk.x()
375 *breakupRate_()*fk*fk.phase();
379 const phaseInterface interface(fi.phase(), fk.phase());
381 if (dmdtfs_.found(interface))
383 const scalar dmdtSign =
384 interface.index(fi.phase()) == 0 ? +1 : -1;
386 *dmdtfs_[interface] += dmdtSign*Sui_*fk.phase().rho();
389 sizeGroups_[i].shapeModelPtr()->addBreakup(Sui_, fk);
394 void Foam::diameterModels::populationBalanceModel::deathByBreakup(
const label i)
396 Sp_[i] += breakupRate_()*sizeGroups()[i].phase();
400 void Foam::diameterModels::populationBalanceModel::calcDeltas()
404 if (delta_[i].empty())
406 for (
label j = 0; j <= sizeGroups().size() - 1; j++)
408 const sizeGroup& fj = sizeGroups()[j];
416 v_[i+1].value() - v_[i].value()
422 v_[i].value() < 0.5*fj.x().value()
424 0.5*fj.x().value() < v_[i+1].value()
427 delta_[i][j] =
mag(0.5*fj.x() - v_[i]);
429 else if (0.5*fj.x().value() <= v_[i].value())
431 delta_[i][j].value() = 0;
439 void Foam::diameterModels::populationBalanceModel::birthByBinaryBreakup
445 const sizeGroup& fi = sizeGroups()[i];
446 const sizeGroup& fj = sizeGroups()[j];
450 Sui_ = fi.x()*delta_[i][j]/fj.x()*
Su;
454 sizeGroups_[i].shapeModelPtr()->addBreakup(Sui_, fj);
456 const phaseInterface interfaceij(fi.phase(), fj.phase());
458 if (dmdtfs_.found(interfaceij))
460 const scalar dmdtSign =
461 interfaceij.index(fi.phase()) == 0 ? +1 : -1;
463 *dmdtfs_[interfaceij] += dmdtSign*Sui_*fj.phase().rho();
473 if (Eta.value() == 0)
continue;
475 const sizeGroup& fk = sizeGroups()[
k];
479 Suk = fk.x()*delta_[i][j]*Eta/fj.x()*
Su;
483 const phaseInterface interfacekj(fk.phase(), fj.phase());
485 if (dmdtfs_.found(interfacekj))
487 const scalar dmdtSign =
488 interfacekj.index(fk.phase()) == 0 ? +1 : -1;
490 *dmdtfs_[interfacekj] += dmdtSign*Suk*fj.phase().rho();
493 sizeGroups_[
k].shapeModelPtr()->addBreakup(Suk, fj);
498 void Foam::diameterModels::populationBalanceModel::deathByBinaryBreakup
504 Sp_[i] += sizeGroups()[i].phase()*binaryBreakupRate_()*delta_[j][i];
508 void Foam::diameterModels::populationBalanceModel::drift
514 model.addToDriftRate(driftRate_(), i);
516 const sizeGroup& fp = sizeGroups()[i];
518 if (i < sizeGroups().size() - 1)
520 const sizeGroup& fe = sizeGroups()[i+1];
523 Sp_[i] += 1/(fe.x() - fp.x())*
pos(driftRate_())*driftRate_();
526 fe.x()/(fp.x()*(fe.x() - fp.x()))*
pos(driftRate_())*driftRate_()*fp;
530 const phaseInterface interfacepe(fp.phase(), fe.phase());
532 if (dmdtfs_.found(interfacepe))
534 const scalar dmdtSign =
535 interfacepe.index(fp.phase()) == 0 ? +1 : -1;
537 *dmdtfs_[interfacepe] -= dmdtSign*Sue*fp.phase().rho();
540 sizeGroups_[i+1].shapeModelPtr()->addDrift(Sue, fp, model);
543 if (i == sizeGroups().size() - 1)
545 Sp_[i] -=
pos(driftRate_())*driftRate_()/fp.x();
550 const sizeGroup& fw = sizeGroups()[i-1];
553 Sp_[i] += 1/(fw.x() - fp.x())*
neg(driftRate_())*driftRate_();
556 fw.x()/(fp.x()*(fw.x() - fp.x()))*
neg(driftRate_())*driftRate_()*fp;
560 const phaseInterface interfacepw(fp.phase(), fw.phase());
562 if (dmdtfs_.found(interfacepw))
564 const scalar dmdtSign =
565 interfacepw.index(fp.phase()) == 0 ? +1 : -1;
567 *dmdtfs_[interfacepw] -= dmdtSign*Suw*fp.phase().rho();
570 sizeGroups_[i-1].shapeModelPtr()->addDrift(Suw, fp, model);
575 Sp_[i] -=
neg(driftRate_())*driftRate_()/fp.x();
580 void Foam::diameterModels::populationBalanceModel::nucleation
583 nucleationModel& model
586 const sizeGroup& fi = sizeGroups()[i];
588 model.addToNucleationRate(nucleationRate_(), i);
590 Sui_ = fi.x()*nucleationRate_();
594 sizeGroups_[i].shapeModelPtr()->addNucleation(Sui_, fi, model);
598 void Foam::diameterModels::populationBalanceModel::sources()
602 sizeGroups_[i].shapeModelPtr()->reset();
612 forAll(coalescencePairs_, coalescencePairi)
614 label i = coalescencePairs_[coalescencePairi].first();
615 label j = coalescencePairs_[coalescencePairi].second();
617 coalescenceRate_() =
Zero;
619 forAll(coalescenceModels_, model)
621 coalescenceModels_[model].addToCoalescenceRate
629 birthByCoalescence(i, j);
631 deathByCoalescence(i, j);
634 forAll(binaryBreakupPairs_, binaryBreakupPairi)
636 label i = binaryBreakupPairs_[binaryBreakupPairi].first();
637 label j = binaryBreakupPairs_[binaryBreakupPairi].second();
639 binaryBreakupRate_() =
Zero;
641 forAll(binaryBreakupModels_, model)
643 binaryBreakupModels_[model].addToBinaryBreakupRate
645 binaryBreakupRate_(),
651 birthByBinaryBreakup(j, i);
653 deathByBinaryBreakup(j, i);
658 forAll(breakupModels_, model)
660 breakupModels_[model].setBreakupRate(breakupRate_(), i);
662 birthByBreakup(i, model);
671 drift(i, drift_[model]);
674 forAll(nucleation_, model)
676 nucleationRate_() =
Zero;
678 nucleation(i, nucleation_[model]);
684 void Foam::diameterModels::populationBalanceModel::correctDilatationError()
688 HashTable<volScalarField>,
694 const word& phaseName = iter.key();
695 const phaseModel& phase = fluid_.phases()[phaseName];
696 const velocityGroup& velGroup = *velocityGroupPtrs_[phaseName];
704 forAll(velGroup.sizeGroups(), i)
706 const sizeGroup& fi = velGroup.sizeGroups()[i];
708 dilatationError -= Su_[fi.i()] -
fvc::Sp(Sp_[fi.i()], fi);
714 void Foam::diameterModels::populationBalanceModel::calcAlphas()
718 forAllConstIter(HashTable<const velocityGroup*>, velocityGroupPtrs_, iter)
720 const phaseModel& phase = iter()->phase();
722 alphas_() +=
max(phase, phase.residualAlpha());
728 Foam::diameterModels::populationBalanceModel::calcDsm()
730 tmp<volScalarField> tInvDsm
742 forAllConstIter(HashTable<const velocityGroup*>, velocityGroupPtrs_, iter)
744 const phaseModel& phase = iter()->phase();
746 invDsm +=
max(phase, phase.residualAlpha())/(phase.d()*alphas_());
753 void Foam::diameterModels::populationBalanceModel::calcVelocity()
757 forAllConstIter(HashTable<const velocityGroup*>, velocityGroupPtrs_, iter)
759 const phaseModel& phase = iter()->phase();
761 U_() +=
max(phase, phase.residualAlpha())*phase.U()/alphas_();
765 bool Foam::diameterModels::populationBalanceModel::updateSources()
767 const bool result = sourceUpdateCounter_ % sourceUpdateInterval() == 0;
769 ++ sourceUpdateCounter_;
794 mesh_(fluid_.mesh()),
798 fluid_.subDict(
"populationBalanceCoeffs").subDict(name_)
804 IOobject::groupName(
"alpha", dict_.lookup(
"continuousPhase"))
825 dict_.lookup(
"coalescenceModels"),
832 dict_.lookup(
"breakupModels"),
838 dict_.lookup(
"binaryBreakupModels"),
841 binaryBreakupRate_(),
842 binaryBreakupPairs_(),
845 dict_.lookup(
"driftModels"),
851 dict_.lookup(
"nucleationModels"),
858 sourceUpdateCounter_(0)
860 this->registerVelocityGroups();
865 <<
"The populationBalance " << name_
866 <<
" requires a minimum number of three sizeGroups to be specified."
870 this->initialiseDmdtfs();
872 if (coalescenceModels_.size() != 0)
891 for (
label j = 0; j <= i; j++)
893 coalescencePairs_.append(
labelPair(i, j));
898 if (breakupModels_.size() != 0)
916 if (binaryBreakupModels_.size() != 0)
918 binaryBreakupRate_.set
944 while (delta_[j][i].value() != 0)
946 binaryBreakupPairs_.append(
labelPair(i, j));
952 if (drift_.size() != 0)
970 if (nucleation_.size() != 0)
993 if (velocityGroupPtrs_.size() > 1)
1082 if (i != 0) lowerBoundary = sizeGroups()[i-1].x();
1084 if (i != sizeGroups().size() - 1) upperBoundary = sizeGroups()[i+1].x();
1086 if ((i == 0 && v < x0) || (i == sizeGroups().size() - 1 && v > xm))
1090 else if (v < lowerBoundary || v > upperBoundary)
1100 return (upperBoundary - v)/(upperBoundary - xi);
1104 return (v - lowerBoundary)/(xi - lowerBoundary);
1125 continuousPhase_.name()
1132 if (!solveOnFinalIterOnly() || fluid_.pimple().finalIter())
1134 const label nCorr = this->nCorr();
1135 const scalar tolerance =
1136 mesh_.solution().solverDict(name_).lookup<scalar>(
"tolerance");
1138 const bool updateSrc = updateSources();
1140 if (nCorr > 0 && updateSrc)
1146 scalar maxInitialResidual = 1;
1148 while (++iCorr <= nCorr && maxInitialResidual > tolerance)
1150 Info<<
"populationBalance "
1161 correctDilatationError();
1163 maxInitialResidual = 0;
1172 dilatationErrors_[phase.
name()];
1178 +
fvm::Sp(-(1 - small)*dilatationError, fi)
1189 /this->mesh().time().deltaT(),
1195 sizeGroupEqn.
relax();
1197 fluid_.fvConstraints().constrain(sizeGroupEqn);
1199 maxInitialResidual =
max
1201 sizeGroupEqn.
solve().initialResidual(),
1205 fluid_.fvConstraints().constrain(fi);
1211 sizeGroups().
first()*sizeGroups().
first().phase()
1216 sizeGroups().last()*sizeGroups().last().phase()
1219 Info<< this->
name() <<
" sizeGroup phase fraction first, last = "
1229 if (velocityGroupPtrs_.size() > 1)
#define forAll(list, i)
Loop across all elements in list.
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
dimensioned< Type > weightedAverage(const DimensionedField< scalar, GeoMesh > &) const
Calculate and return weighted average.
Generic GeometricField class.
static tmp< GeometricField< Type, PatchField, GeoMesh > > New(const word &name, const Internal &, const PtrList< PatchField< Type >> &)
Return a temporary field constructed from name,.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
const Time & time() const
Return time.
IOobject(const word &name, const fileName &instance, const objectRegistry ®istry, readOption r=NO_READ, writeOption w=NO_WRITE, bool registerObject=true)
Construct from name, instance, registry, io options.
const word & name() const
Return name.
static word groupName(Name name, const word &group)
bool good() const
Return true if next operation might succeed.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Base class for binary breakup models that provide a breakup rate between a size class pair directly,...
Base class for breakup models which provide a total breakup rate and a separate daughter size distrib...
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.
Model for tracking the evolution of a dispersed phase size distribution due to coalescence (synonymou...
const dimensionedScalar eta(const label i, const dimensionedScalar &v) const
Return allocation coefficient.
bool writeData(Ostream &) const
Dummy write for regIOobject.
void correct()
Correct derived quantities.
autoPtr< populationBalanceModel > clone() const
Return clone.
populationBalanceModel(const phaseSystem &fluid, const word &name)
Construct for a fluid.
const UPtrList< sizeGroup > & sizeGroups() const
Return the size groups belonging to this populationBalance.
const tmp< volScalarField > sigmaWithContinuousPhase(const phaseModel &dispersedPhase) const
Return the surface tension coefficient between a given dispersed.
const phaseCompressible::momentumTransportModel & continuousTurbulence() const
Return reference to momentumTransport model of the continuous phase.
virtual ~populationBalanceModel()
Destructor.
void solve()
Solve the population balance equation.
Single size class fraction field representing a fixed particle volume as defined by the user through ...
const phaseModel & phase() const
Return const-reference to the phase.
const Type & value() const
Return const reference to value.
const word & name() const
Return const reference to name.
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 segregated or coupled returning the solution statistics.
const Time & time() const
Return the top-level database.
Templated abstract base class for multiphase compressible turbulence models.
Class to represent an interface between phases. Derivations can further specify the configuration of ...
tmp< volScalarField > sigma() const
Surface tension coefficient.
const dimensionedScalar & residualAlpha() const
Return the residual phase-fraction for given phase.
const word & name() const
Return the name of this phase.
virtual const volScalarField & rho() const =0
Return the density field.
virtual tmp< surfaceScalarField > alphaPhi() const =0
Return the volumetric flux of the phase.
Class to represent a system of phases and model interfacial transfers between them.
const phaseModelList & phases() const
Return the phase models.
HashPtrTable< volScalarField, phaseInterfaceKey, phaseInterfaceKey::hash > dmdtfTable
regIOobject is an abstract class derived from IOobject to handle automatic object registration with t...
A class for managing temporary objects.
A class for handling words, derived from string.
#define notImplemented(functionName)
Issue a FatalErrorIn for a function not currently implemented.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
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 matrix for the divergence of the given field and flux.
Calculate the matrix for implicit and explicit sources.
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
const char *const group
Group name for atomic constants.
defineTypeNameAndDebug(constant, 0)
tmp< VolField< Type > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
tmp< VolField< Type > > Su(const VolField< Type > &su, const VolField< Type > &vf)
tmp< VolField< Type > > Sp(const volScalarField &sp, const VolField< Type > &vf)
tmp< VolField< Type > > div(const SurfaceField< Type > &ssf)
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const VolField< Type > &)
tmp< fvMatrix< Type > > SuSp(const volScalarField::Internal &, const VolField< Type > &)
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const VolField< Type > &vf, const word &name)
tmp< fvMatrix< Type > > ddt(const VolField< Type > &vf)
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
dimensionedScalar pos(const dimensionedScalar &ds)
Pair< label > labelPair
Label pair.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const dimensionSet dimless
const dimensionSet dimLength
labelList first(const UList< labelPair > &p)
const dimensionSet dimTime
const dimensionSet dimDensity
const dimensionSet dimVolume
VolField< scalar > volScalarField
dimensioned< scalar > mag(const dimensioned< Type > &)
word typedName(Name name)
Return the name of the object within the given type.
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix.
dimensionedScalar neg(const dimensionedScalar &ds)
const dimensionSet dimVelocity
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
word name(const complex &)
Return a string representation of a complex.