34 #include "phaseSystem.H"
35 #include "surfaceTensionModel.H"
41 #include "phaseCompressibleTurbulenceModel.H"
45 void Foam::diameterModels::populationBalanceModel::registerVelocityGroups()
49 if (isA<velocityGroup>(fluid_.
phases()[
phasei].dPtr()()))
52 refCast<const velocityGroup>(fluid_.
phases()[
phasei].dPtr()());
54 if (velGroup.popBalName() == this->
name())
56 velocityGroups_.resize(velocityGroups_.size() + 1);
60 velocityGroups_.size() - 1,
64 forAll(velGroup.sizeGroups(), i)
66 this->registerSizeGroups
68 const_cast<sizeGroup&
>(velGroup.sizeGroups()[i])
77 void Foam::diameterModels::populationBalanceModel::registerSizeGroups
84 sizeGroups_.size() != 0
86 group.x().value() <= sizeGroups_.last().x().value()
90 <<
"Size groups must be entered according to their representative"
95 sizeGroups_.resize(sizeGroups_.size() + 1);
96 sizeGroups_.set(sizeGroups_.size() - 1, &
group);
99 if (sizeGroups_.size() == 1)
108 sizeGroups_.last().x()
119 sizeGroups_.last().x()
130 sizeGroups_[sizeGroups_.size()-2].x()
131 + sizeGroups_.last().x()
141 sizeGroups_.last().x()
146 delta_.append(
new PtrList<dimensionedScalar>());
155 fluid_.time().timeName(),
170 fluid_.time().timeName(),
180 void Foam::diameterModels::populationBalanceModel::createPhasePairs()
182 forAll(velocityGroups_, i)
184 const phaseModel&
phasei = velocityGroups_[i].phase();
186 forAll(velocityGroups_, j)
188 const phaseModel& phasej = velocityGroups_[j].phase();
192 const phasePairKey key
199 if (!phasePairs_.found(key))
220 void Foam::diameterModels::populationBalanceModel::correct()
224 forAll(velocityGroups_, v)
226 velocityGroups_[v].preSolve();
229 forAll(coalescence_, model)
231 coalescence_[model].correct();
236 breakup_[model].correct();
238 breakup_[model].dsdPtr()().
correct();
241 forAll(binaryBreakup_, model)
243 binaryBreakup_[model].correct();
248 drift_[model].correct();
251 forAll(nucleation_, model)
253 nucleation_[model].correct();
258 void Foam::diameterModels::populationBalanceModel::
265 const sizeGroup& fj = sizeGroups_[j];
266 const sizeGroup& fk = sizeGroups_[
k];
271 for (
label i = j; i < sizeGroups_.size(); i++)
276 if (Gamma.value() == 0)
continue;
278 const sizeGroup& fi = sizeGroups_[i];
284 0.5*fi.x()*coalescenceRate_()*Gamma
285 *fj*fj.phase()/fj.x()
286 *fk*fk.phase()/fk.x();
291 fi.x()*coalescenceRate_()*Gamma
292 *fj*fj.phase()/fj.x()
293 *fk*fk.phase()/fk.x();
298 const phasePairKey pairij
304 if (pDmdt_.found(pairij))
306 const scalar dmdtSign
311 pDmdt_[pairij]->ref() += dmdtSign*fj.x()/v*Sui_*fi.phase().rho();
314 const phasePairKey pairik
320 if (pDmdt_.found(pairik))
322 const scalar dmdtSign
327 pDmdt_[pairik]->ref() += dmdtSign*fk.x()/v*Sui_*fi.phase().rho();
333 void Foam::diameterModels::populationBalanceModel::
340 const sizeGroup& fi = sizeGroups_[i];
341 const sizeGroup& fj = sizeGroups_[j];
343 SuSp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x();
347 SuSp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x();
352 void Foam::diameterModels::populationBalanceModel::
359 const sizeGroup& fk = sizeGroups_[
k];
361 for (
label i = 0; i <=
k; i++)
363 const sizeGroup& fi = sizeGroups_[i];
366 fi.x()*breakupRate_()*breakup_[model].dsdPtr()().nik(i,
k)
367 *fk*fk.phase()/fk.x();
371 const phasePairKey pair
377 if (pDmdt_.found(pair))
379 const scalar dmdtSign
384 pDmdt_[pair]->ref() += dmdtSign*Sui_*fi.phase().rho();
390 void Foam::diameterModels::populationBalanceModel::deathByBreakup(
const label i)
392 const sizeGroup& fi = sizeGroups_[i];
394 SuSp_[i] += breakupRate_()*fi.phase();
398 void Foam::diameterModels::populationBalanceModel::calcDeltas()
402 if (delta_[i].empty())
404 for (
label j = 0; j <= sizeGroups_.size() - 1; j++)
412 v_[i+1].value() - v_[i].value()
418 v_[i].value() < 0.5*sizeGroups_[j].
x().value()
420 0.5*sizeGroups_[j].
x().value() < v_[i+1].value()
423 delta_[i][j] =
mag(0.5*sizeGroups_[j].
x() - v_[i]);
425 else if (0.5*sizeGroups_[j].
x().value() <= v_[i].value())
427 delta_[i][j].value() = 0;
435 void Foam::diameterModels::populationBalanceModel::
442 const sizeGroup& fj = sizeGroups_[j];
443 const sizeGroup& fi = sizeGroups_[i];
445 Sui_ = fi.x()*binaryBreakupRate_()*delta_[i][j]*fj*fj.phase()/fj.x();
449 const phasePairKey pairij
455 if (pDmdt_.found(pairij))
457 const scalar dmdtSign
462 pDmdt_[pairij]->ref() += dmdtSign*Sui_*fi.phase().rho();
473 if (Gamma.value() == 0)
continue;
475 const sizeGroup& fk = sizeGroups_[
k];
480 sizeGroups_[
k].x()*binaryBreakupRate_()*delta_[i][j]*Gamma
481 *fj*fj.phase()/fj.x();
485 const phasePairKey pairkj
491 if (pDmdt_.found(pairkj))
493 const scalar dmdtSign
497 pDmdt_.find(pairkj).key(),
502 pDmdt_[pairkj]->ref() += dmdtSign*Suk*fi.phase().rho();
508 void Foam::diameterModels::populationBalanceModel::
517 SuSp_[i] += alphai*binaryBreakupRate_()*delta_[j][i];
521 void Foam::diameterModels::populationBalanceModel::drift(
const label i)
523 const sizeGroup& fp = sizeGroups_[i];
527 rx_() =
pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x();
529 else if (i == sizeGroups_.size() - 1)
531 rx_() =
neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
536 pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x()
537 +
neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
541 (
neg(1 - rx_()) +
neg(rx_() - rx_()/(1 - rx_())))*driftRate_()
542 *fp.phase()/((rx_() - 1)*fp.x());
547 if (i == sizeGroups_.size() - 2)
549 rx_() =
pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x();
553 *(sizeGroups_[i+1].x() - sizeGroups_[i].x())
554 /(sizeGroups_[i].
x() - sizeGroups_[i-1].x());
556 else if (i < sizeGroups_.size() - 2)
558 rx_() =
pos(driftRate_())*sizeGroups_[i+2].x()/sizeGroups_[i+1].x();
562 *(sizeGroups_[i+2].x() - sizeGroups_[i+1].x())
563 /(sizeGroups_[i+1].
x() - sizeGroups_[i].x());
568 rx_() +=
neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
572 *(sizeGroups_[i].x() - sizeGroups_[i-1].x())
573 /(sizeGroups_[i+1].
x() - sizeGroups_[i].x());
577 rx_() +=
neg(driftRate_())*sizeGroups_[i-2].x()/sizeGroups_[i-1].x();
581 *(sizeGroups_[i-1].x() - sizeGroups_[i-2].x())
582 /(sizeGroups_[i].
x() - sizeGroups_[i-1].x());
585 if (i != sizeGroups_.size() - 1)
587 const sizeGroup& fe = sizeGroups_[i+1];
591 pos(driftRate_())*driftRate_()*rdx_()
592 *fp*fp.phase()/fp.x()
597 const phasePairKey pairij
603 if (pDmdt_.found(pairij))
605 const scalar dmdtSign
610 pDmdt_[pairij]->ref() -= dmdtSign*Sue*fp.phase().rho();
616 const sizeGroup& fw = sizeGroups_[i-1];
620 neg(driftRate_())*driftRate_()*rdx_()
621 *fp*fp.phase()/fp.x()
626 const phasePairKey pairih
632 if (pDmdt_.found(pairih))
634 const scalar dmdtSign
639 pDmdt_[pairih]->ref() -= dmdtSign*Suw*fp.phase().rho();
645 void Foam::diameterModels::populationBalanceModel::nucleation(
const label i)
647 Su_[i] += sizeGroups_[i].x()*nucleationRate_();
651 void Foam::diameterModels::populationBalanceModel::sources()
666 pDmdt_(phasePairIter())->ref() *= 0.0;
674 const sizeGroup& fi = sizeGroups_[i];
676 if (coalescence_.size() != 0)
678 for (
label j = 0; j <= i; j++)
680 const sizeGroup& fj = sizeGroups_[j];
682 if (fi.x() + fj.x() > sizeGroups_.last().x())
break;
684 coalescenceRate_() *= 0.0;
686 forAll(coalescence_, model)
688 coalescence_[model].addToCoalescenceRate
696 birthByCoalescence(i, j);
698 deathByCoalescence(i, j);
702 if (breakup_.size() != 0)
706 breakup_[model].setBreakupRate(breakupRate_(), i);
708 birthByBreakup(i, model);
714 if (binaryBreakup_.size() != 0)
718 while (delta_[j][i].value() != 0)
720 binaryBreakupRate_() *= 0.0;
722 forAll(binaryBreakup_, model)
724 binaryBreakup_[model].addToBinaryBreakupRate
726 binaryBreakupRate_(),
732 birthByBinaryBreakup(j, i);
734 deathByBinaryBreakup(j, i);
740 if (drift_.size() != 0)
746 drift_[model].addToDriftRate(driftRate_(), i);
752 if (nucleation_.size() != 0)
754 nucleationRate_() *= 0.0;
756 forAll(nucleation_, model)
758 nucleation_[model].addToNucleationRate(nucleationRate_(), i);
767 void Foam::diameterModels::populationBalanceModel::dmdt()
769 forAll(velocityGroups_, v)
771 velocityGroup& velGroup = velocityGroups_[v];
773 velGroup.dmdtRef() *= 0.0;
777 if (&sizeGroups_[i].phase() == &velGroup.phase())
779 sizeGroup& fi = sizeGroups_[i];
781 velGroup.dmdtRef() += fi.phase().rho()*(Su_[i] - SuSp_[i]*fi);
788 void Foam::diameterModels::populationBalanceModel::calcAlphas()
792 forAll(velocityGroups_, v)
794 const phaseModel& phase = velocityGroups_[v].phase();
796 alphas_() +=
max(phase, phase.residualAlpha());
802 Foam::diameterModels::populationBalanceModel::calcDsm()
804 tmp<volScalarField> tInvDsm
816 forAll(velocityGroups_, v)
818 const phaseModel& phase = velocityGroups_[v].phase();
820 invDsm +=
max(phase, phase.residualAlpha())/(phase.d()*alphas_());
827 void Foam::diameterModels::populationBalanceModel::calcVelocity()
831 forAll(velocityGroups_, v)
833 const phaseModel& phase = velocityGroups_[v].phase();
835 U_() +=
max(phase, phase.residualAlpha())*phase.U()/alphas_();
839 bool Foam::diameterModels::populationBalanceModel::updateSources()
841 const bool result = sourceUpdateCounter_ % sourceUpdateInterval() == 0;
843 ++ sourceUpdateCounter_;
863 fluid.time().constant(),
869 mesh_(fluid_.mesh()),
873 fluid_.subDict(
"populationBalanceCoeffs").subDict(name_)
875 pimple_(mesh_.lookupObject<
pimpleControl>(
"solutionControl")),
894 mesh_.time().timeName(),
902 dict_.lookup(
"coalescenceModels"),
908 dict_.lookup(
"breakupModels"),
914 dict_.lookup(
"binaryBreakupModels"),
917 binaryBreakupRate_(),
920 dict_.lookup(
"driftModels"),
928 dict_.lookup(
"nucleationModels"),
937 (mesh_.time().timeIndex()*
nCorr())%sourceUpdateInterval()
940 this->registerVelocityGroups();
942 this->createPhasePairs();
944 if (sizeGroups_.size() < 3)
947 <<
"The populationBalance " << name_
948 <<
" requires a minimum number of three sizeGroups to be specified."
953 if (coalescence_.size() != 0)
962 mesh_.time().timeName(),
971 if (breakup_.size() != 0)
980 fluid_.time().timeName(),
989 if (binaryBreakup_.size() != 0)
991 binaryBreakupRate_.set
998 fluid_.time().timeName(),
1004 "binaryBreakupRate",
1012 if (drift_.size() != 0)
1021 fluid_.time().timeName(),
1036 fluid_.time().timeName(),
1051 fluid_.time().timeName(),
1060 if (nucleation_.size() != 0)
1069 fluid_.time().timeName(),
1083 if (velocityGroups_.size() > 1)
1092 fluid_.time().timeName(),
1109 fluid_.time().timeName(),
1126 fluid_.time().timeName(),
1176 lowerBoundary = sizeGroups_[i-1].x();
1179 if (i == sizeGroups_.size() - 1)
1185 upperBoundary = sizeGroups_[i+1].x();
1188 if (v < lowerBoundary || v > upperBoundary)
1194 return (v - lowerBoundary)/(xi - lowerBoundary);
1198 return (upperBoundary - v)/(upperBoundary - xi);
1212 phasePair(dispersedPhase, continuousPhase_)
1226 continuousPhase_.name()
1234 const dictionary& solutionControls = mesh_.solverDict(name_);
1235 bool solveOnFinalIterOnly =
1236 solutionControls.
lookupOrDefault<
bool>(
"solveOnFinalIterOnly",
false);
1238 if (!solveOnFinalIterOnly || pimple_.finalIter())
1241 const scalar tolerance =
1242 readScalar(solutionControls.
lookup(
"tolerance"));
1250 scalar maxInitialResidual = 1;
1252 while (++iCorr <= nCorr && maxInitialResidual > tolerance)
1254 Info<<
"populationBalance "
1260 if (updateSources())
1267 maxInitialResidual = 0;
1282 phase.alphaRhoPhi(),
1298 sizeGroupEqn.
relax();
1300 maxInitialResidual =
max
1302 sizeGroupEqn.solve().initialResidual(),
1310 forAll(velocityGroups_, i)
1312 velocityGroups_[i].postSolve();
1316 if (velocityGroups_.size() > 1)
1325 sizeGroups_.first()*sizeGroups_.first().phase()
1330 sizeGroups_.last()*sizeGroups_.last().phase()
1333 Info<< this->
name() <<
" sizeGroup phase fraction first, last = "
1334 << fAlpha0.weightedAverage(this->
mesh().V()).value()
1335 <<
' ' << fAlphaN.weightedAverage(this->
mesh().V()).value()