33 #include "surfaceTensionModel.H" 43 namespace diameterModels
52 void Foam::diameterModels::populationBalanceModel::registerVelocityGroups()
56 if (isA<velocityGroup>(fluid_.phases()[
phasei].dPtr()()))
58 const velocityGroup& velGroup =
59 refCast<const velocityGroup>(fluid_.phases()[
phasei].dPtr()());
61 if (velGroup.popBalName() == this->
name())
63 velocityGroups_.resize(velocityGroups_.size() + 1);
67 velocityGroups_.size() - 1,
68 &
const_cast<velocityGroup&
>(velGroup)
71 forAll(velGroup.sizeGroups(), i)
73 this->registerSizeGroups
75 const_cast<sizeGroup&>(velGroup.sizeGroups()[i])
84 void Foam::diameterModels::populationBalanceModel::registerSizeGroups
91 sizeGroups().size() != 0
93 group.x().value() <= sizeGroups().last().x().value()
97 <<
"Size groups must be entered according to their representative" 102 sizeGroups_.resize(sizeGroups().size() + 1);
103 sizeGroups_.set(sizeGroups().size() - 1, &group);
105 if (sizeGroups().size() == 1)
112 sizeGroups().last().
x()
121 sizeGroups().last().
x()
130 sizeGroups()[sizeGroups().size()-2].x()
131 + sizeGroups().last().x()
139 sizeGroups().last().
x()
144 delta_.append(
new PtrList<dimensionedScalar>());
153 fluid_.time().timeName(),
168 fluid_.time().timeName(),
178 void Foam::diameterModels::populationBalanceModel::createPhasePairs()
180 forAll(velocityGroups_, i)
182 const phaseModel& phasei = velocityGroups_[i].phase();
184 forAll(velocityGroups_, j)
186 const phaseModel& phasej = velocityGroups_[j].phase();
188 if (&phasei != &phasej)
190 const phasePairKey key
197 if (!phasePairs_.found(key))
218 void Foam::diameterModels::populationBalanceModel::precompute()
220 forAll(coalescenceModels_, model)
222 coalescenceModels_[model].precompute();
225 forAll(breakupModels_, model)
227 breakupModels_[model].precompute();
229 breakupModels_[model].dsdPtr()->precompute();
232 forAll(binaryBreakupModels_, model)
234 binaryBreakupModels_[model].precompute();
239 drift_[model].precompute();
242 forAll(nucleation_, model)
244 nucleation_[model].precompute();
249 void Foam::diameterModels::populationBalanceModel::birthByCoalescence
255 const sizeGroup& fj = sizeGroups()[j];
256 const sizeGroup& fk = sizeGroups()[
k];
261 for (
label i = j; i < sizeGroups().size(); i++)
265 if (Eta.value() == 0)
continue;
267 const sizeGroup& fi = sizeGroups()[i];
272 0.5*fi.x()/(fj.x()*fk.x())*Eta
273 *coalescenceRate_()*fj*fj.phase()*fk*fk.phase();
278 fi.x()/(fj.x()*fk.x())*Eta
279 *coalescenceRate_()*fj*fj.phase()*fk*fk.phase();
284 const phasePairKey pairij
290 if (pDmdt_.found(pairij))
292 const scalar dmdtSign
297 *pDmdt_[pairij] += dmdtSign*fj.x()/v*Sui_*fj.phase().rho();
300 const phasePairKey pairik
306 if (pDmdt_.found(pairik))
308 const scalar dmdtSign
313 *pDmdt_[pairik] += dmdtSign*fk.x()/v*Sui_*fk.phase().rho();
316 sizeGroups_[i].shapeModelPtr()->addCoalescence(Sui_, fj, fk);
321 void Foam::diameterModels::populationBalanceModel::deathByCoalescence
327 const sizeGroup& fi = sizeGroups()[i];
328 const sizeGroup& fj = sizeGroups()[j];
330 Sp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x();
334 Sp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x();
339 void Foam::diameterModels::populationBalanceModel::birthByBreakup
345 const sizeGroup& fk = sizeGroups()[
k];
347 for (
label i = 0; i <=
k; i++)
349 const sizeGroup& fi = sizeGroups()[i];
352 fi.x()*breakupModels_[model].dsdPtr()().nik(i, k)/fk.x()
353 *breakupRate_()*fk*fk.phase();
357 const phasePairKey pair
363 if (pDmdt_.found(pair))
365 const scalar dmdtSign
370 *pDmdt_[pair] += dmdtSign*Sui_*fk.phase().rho();
373 sizeGroups_[i].shapeModelPtr()->addBreakup(Sui_, fk);
378 void Foam::diameterModels::populationBalanceModel::deathByBreakup(
const label i)
380 Sp_[i] += breakupRate_()*sizeGroups()[i].phase();
384 void Foam::diameterModels::populationBalanceModel::calcDeltas()
388 if (delta_[i].empty())
390 for (
label j = 0; j <= sizeGroups().size() - 1; j++)
392 const sizeGroup& fj = sizeGroups()[j];
400 v_[i+1].value() - v_[i].value()
406 v_[i].value() < 0.5*fj.x().value()
408 0.5*fj.x().value() < v_[i+1].value()
411 delta_[i][j] =
mag(0.5*fj.x() - v_[i]);
413 else if (0.5*fj.x().value() <= v_[i].value())
415 delta_[i][j].
value() = 0;
423 void Foam::diameterModels::populationBalanceModel::birthByBinaryBreakup
429 const sizeGroup& fi = sizeGroups()[i];
430 const sizeGroup& fj = sizeGroups()[j];
434 Sui_ = fi.x()*delta_[i][j]/fj.x()*
Su;
438 sizeGroups_[i].shapeModelPtr()->addBreakup(Sui_, fj);
440 const phasePairKey pairij
446 if (pDmdt_.found(pairij))
448 const scalar dmdtSign
453 *pDmdt_[pairij] += dmdtSign*Sui_*fj.phase().rho();
459 for (
label k = 0; k <= j; k++)
463 if (Eta.value() == 0)
continue;
465 const sizeGroup& fk = sizeGroups()[
k];
469 Suk = fk.x()*delta_[i][j]*Eta/fj.x()*
Su;
473 const phasePairKey pairkj
479 if (pDmdt_.found(pairkj))
481 const scalar dmdtSign
485 pDmdt_.find(pairkj).key(),
490 *pDmdt_[pairkj] += 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_()*fp.phase();
526 fe.x()/(fp.x()*(fe.x() - fp.x()))
527 *
pos(driftRate_())*driftRate_()*fp*fp.phase();
531 const phasePairKey pairij
537 if (pDmdt_.found(pairij))
539 const scalar dmdtSign
544 *pDmdt_[pairij] -= dmdtSign*Sue*fp.phase().rho();
547 sizeGroups_[i+1].shapeModelPtr()->addDrift(Sue, fp, model);
550 if (i == sizeGroups().size() - 1)
552 Sp_[i] -=
pos(driftRate_())*driftRate_()*fp.phase()/fp.x();
557 const sizeGroup& fw = sizeGroups()[i-1];
560 Sp_[i] += 1/(fw.x() - fp.x())*
neg(driftRate_())*driftRate_()*fp.phase();
563 fw.x()/(fp.x()*(fw.x() - fp.x()))
564 *
neg(driftRate_())*driftRate_()*fp*fp.phase();
568 const phasePairKey pairih
574 if (pDmdt_.found(pairih))
576 const scalar dmdtSign
581 *pDmdt_[pairih] -= dmdtSign*Suw*fp.phase().rho();
584 sizeGroups_[i-1].shapeModelPtr()->addDrift(Suw, fp, model);
589 Sp_[i] -=
neg(driftRate_())*driftRate_()*fp.phase()/fp.x();
594 void Foam::diameterModels::populationBalanceModel::nucleation
597 nucleationModel& model
600 const sizeGroup& fi = sizeGroups()[i];
602 model.addToNucleationRate(nucleationRate_(), i);
604 Sui_ = fi.x()*nucleationRate_();
608 sizeGroups_[i].shapeModelPtr()->addNucleation(Sui_, fi, model);
612 void Foam::diameterModels::populationBalanceModel::sources()
616 sizeGroups_[i].shapeModelPtr()->reset();
628 *pDmdt_(phasePairIter()) =
Zero;
631 forAll(coalescencePairs_, coalescencePairi)
633 label i = coalescencePairs_[coalescencePairi].first();
634 label j = coalescencePairs_[coalescencePairi].second();
636 coalescenceRate_() =
Zero;
638 forAll(coalescenceModels_, model)
640 coalescenceModels_[model].addToCoalescenceRate
648 birthByCoalescence(i, j);
650 deathByCoalescence(i, j);
653 forAll(binaryBreakupPairs_, binaryBreakupPairi)
655 label i = binaryBreakupPairs_[binaryBreakupPairi].first();
656 label j = binaryBreakupPairs_[binaryBreakupPairi].second();
658 binaryBreakupRate_() =
Zero;
660 forAll(binaryBreakupModels_, model)
662 binaryBreakupModels_[model].addToBinaryBreakupRate
664 binaryBreakupRate_(),
670 birthByBinaryBreakup(j, i);
672 deathByBinaryBreakup(j, i);
677 forAll(breakupModels_, model)
679 breakupModels_[model].setBreakupRate(breakupRate_(), i);
681 birthByBreakup(i, model);
690 drift(i, drift_[model]);
693 forAll(nucleation_, model)
695 nucleationRate_() =
Zero;
697 nucleation(i, nucleation_[model]);
703 void Foam::diameterModels::populationBalanceModel::calcAlphas()
707 forAll(velocityGroups_, v)
709 const phaseModel& phase = velocityGroups_[v].phase();
711 alphas_() +=
max(phase, phase.residualAlpha());
717 Foam::diameterModels::populationBalanceModel::calcDsm()
719 tmp<volScalarField> tInvDsm
731 forAll(velocityGroups_, v)
733 const phaseModel& phase = velocityGroups_[v].phase();
735 invDsm +=
max(phase, phase.residualAlpha())/(phase.d()*alphas_());
742 void Foam::diameterModels::populationBalanceModel::calcVelocity()
746 forAll(velocityGroups_, v)
748 const phaseModel& phase = velocityGroups_[v].phase();
750 U_() +=
max(phase, phase.residualAlpha())*phase.U()/alphas_();
754 bool Foam::diameterModels::populationBalanceModel::updateSources()
756 const bool result = sourceUpdateCounter_ % sourceUpdateInterval() == 0;
758 ++ sourceUpdateCounter_;
768 const phaseSystem& fluid,
770 HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>& pDmdt
778 fluid.time().constant(),
784 mesh_(fluid_.
mesh()),
788 fluid_.subDict(
"populationBalanceCoeffs").subDict(name_)
790 pimple_(mesh_.lookupObject<pimpleControl>(
"solutionControl")),
793 mesh_.lookupObject<phaseModel>
795 IOobject::groupName(
"alpha", dict_.
lookup(
"continuousPhase"))
817 dict_.
lookup(
"coalescenceModels"),
818 coalescenceModel::iNew(*this)
824 dict_.
lookup(
"breakupModels"),
825 breakupModel::iNew(*this)
830 dict_.
lookup(
"binaryBreakupModels"),
831 binaryBreakupModel::iNew(*this)
833 binaryBreakupRate_(),
834 binaryBreakupPairs_(),
837 dict_.
lookup(
"driftModels"),
838 driftModel::iNew(*this)
843 dict_.
lookup(
"nucleationModels"),
844 nucleationModel::iNew(*this)
850 sourceUpdateCounter_(0)
852 this->registerVelocityGroups();
854 this->createPhasePairs();
856 if (sizeGroups().size() < 3)
859 <<
"The populationBalance " << name_
860 <<
" requires a minimum number of three sizeGroups to be specified." 865 if (coalescenceModels_.size() != 0)
874 mesh_.time().timeName(),
884 for (
label j = 0; j <= i; j++)
886 coalescencePairs_.append(
labelPair(i, j));
891 if (breakupModels_.size() != 0)
900 fluid_.time().timeName(),
909 if (binaryBreakupModels_.size() != 0)
911 binaryBreakupRate_.set
918 fluid_.time().timeName(),
937 while (delta_[j][i].value() != 0)
939 binaryBreakupPairs_.append(
labelPair(i, j));
945 if (drift_.size() != 0)
954 fluid_.time().timeName(),
963 if (nucleation_.size() != 0)
972 fluid_.time().timeName(),
986 if (velocityGroups_.size() > 1)
995 fluid_.time().timeName(),
1012 fluid_.time().timeName(),
1029 fluid_.time().timeName(),
1052 return autoPtr<populationBalanceModel>(
nullptr);
1075 if (i != 0) lowerBoundary = sizeGroups()[i-1].x();
1077 if (i != sizeGroups().size() - 1) upperBoundary = sizeGroups()[i+1].x();
1079 if ((i == 0 && v < x0) || (i == sizeGroups().size() - 1 && v > xm))
1083 else if (v < lowerBoundary || v > upperBoundary)
1087 else if (v.value() == xi.value())
1093 return (upperBoundary - v)/(upperBoundary - xi);
1097 return (v - lowerBoundary)/(xi - lowerBoundary);
1105 const phaseModel& dispersedPhase
1109 fluid_.lookupSubModel<surfaceTensionModel>
1111 phasePair(dispersedPhase, continuousPhase_)
1124 momentumTransportModel::typeName,
1125 continuousPhase_.name()
1133 const dictionary& solutionControls = mesh_.solverDict(name_);
1134 bool solveOnFinalIterOnly =
1135 solutionControls.lookupOrDefault<
bool>(
"solveOnFinalIterOnly",
false);
1137 if (!solveOnFinalIterOnly || pimple_.finalPimpleIter())
1140 const scalar tolerance =
1141 solutionControls.lookup<scalar>(
"tolerance");
1149 scalar maxInitialResidual = 1;
1151 while (++iCorr <= nCorr && maxInitialResidual > tolerance)
1153 Info<<
"populationBalance " 1159 if (updateSources())
1164 maxInitialResidual = 0;
1168 sizeGroup& fi = sizeGroups_[i];
1169 const phaseModel& phase = fi.phase();
1181 + fluid_.fvModels().source(alpha, rho, fi)/rho
1186 sizeGroupEqn.relax();
1187 fluid_.fvConstraints().constrain(sizeGroupEqn);
1189 maxInitialResidual =
max 1191 sizeGroupEqn.solve().initialResidual(),
1195 fluid_.fvConstraints().constrain(fi);
1201 sizeGroups().first()*sizeGroups().first().phase()
1206 sizeGroups().last()*sizeGroups().last().phase()
1209 Info<< this->
name() <<
" sizeGroup phase fraction first, last = " 1210 << fAlpha0.weightedAverage(this->
mesh().V()).value()
1211 <<
' ' << fAlphaN.weightedAverage(this->
mesh().V()).value()
1219 if (velocityGroups_.size() > 1)
Templated abstract base class for multiphase compressible turbulence models.
fvMatrix< scalar > fvScalarMatrix
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
virtual ~populationBalanceModel()
Destructor.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
populationBalanceModel(const phaseSystem &fluid, const word &name, HashPtrTable< volScalarField, phasePairKey, phasePairKey::hash > &pDmdt)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
static int compare(const Pair< word > &a, const Pair< word > &b)
Compare Pairs.
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< scalar >> &)
Return a temporary field constructed from name,.
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
const phaseCompressible::momentumTransportModel & continuousTurbulence() const
Return reference to turbulence model of the continuous phase.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const dimensionSet dimless
label k
Boltzmann constant.
Generic dimensioned Type class.
bool writeData(Ostream &) const
Dummy write for regIOobject.
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
GeometricField< vector, fvPatchField, volMesh > volVectorField
dimensionedScalar neg(const dimensionedScalar &ds)
const dimensionedScalar eta(const label i, const dimensionedScalar &v) const
Return allocation coefficient.
const dimensionSet dimLength
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
autoPtr< populationBalanceModel > clone() const
Return clone.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
Calculate the first temporal derivative.
dimensionedScalar pos(const dimensionedScalar &ds)
stressControl lookup("compactNormalStress") >> compactNormalStress
static word groupName(Name name, const word &group)
virtual void precompute()
Precompute diameter independent expressions.
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
const Type & value() const
Return const reference to value.
void correct()
Correct derived quantities.
Pair< label > labelPair
Label pair.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
const dimensionSet dimVelocity
defineTypeNameAndDebug(combustionModel, 0)
word name(const complex &)
Return a string representation of a complex.
Internal & ref()
Return a reference to the dimensioned internal field.
#define notImplemented(functionName)
Issue a FatalErrorIn for a function not currently implemented.
PhaseCompressibleMomentumTransportModel< dynamicTransportModel > momentumTransportModel
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const dimensionSet dimVolume
dimensioned< scalar > mag(const dimensioned< Type > &)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
const tmp< volScalarField > sigmaWithContinuousPhase(const phaseModel &dispersedPhase) const
Return the surface tension coefficient between a given dispersed.
const tmp< volScalarField::Internal > & Su
A class for managing temporary objects.
void solve()
Solve the population balance equation.