33 #include "surfaceTensionModel.H" 38 #include "phaseCompressibleTurbulenceModel.H" 43 Foam::diameterModels::populationBalanceModel::registerVelocityAndSizeGroups()
47 if (isA<velocityGroup>(fluid_.
phases()[
phasei].dPtr()()))
50 refCast<const velocityGroup>(fluid_.
phases()[
phasei].dPtr()());
52 if (velGroup.popBalName() == this->
name())
54 velocityGroups_.append(&const_cast<velocityGroup&>(velGroup));
56 forAll(velGroup.sizeGroups(), i)
60 &const_cast<sizeGroup&>(velGroup.sizeGroups()[i])
69 void Foam::diameterModels::populationBalanceModel::add(
sizeGroup* group)
73 sizeGroups_.size() != 0
75 group->x().value() <= sizeGroups_.last()->x().value()
80 "populationBalance::add" 82 ) <<
"Size groups must be entered according to their representative" 88 sizeGroups_.append(group);
91 if (sizeGroups_.size() == 1)
100 sizeGroups_.last()->x()
111 sizeGroups_.last()->x()
122 sizeGroups_[sizeGroups_.size()-2]->x()
123 + sizeGroups_.last()->x()
133 sizeGroups_.last()->x()
138 delta_.append(
new PtrList<dimensionedScalar>());
182 void Foam::diameterModels::populationBalanceModel::createPhasePairs()
184 forAll(velocityGroups_, i)
186 const phaseModel& phasei = velocityGroups_[i]->phase();
188 forAll(velocityGroups_, j)
190 const phaseModel& phasej = velocityGroups_[j]->phase();
192 if (&phasei != &phasej)
194 const phasePairKey key
201 if (!phasePairs_.
found(key))
222 void Foam::diameterModels::populationBalanceModel::correct()
226 forAll(velocityGroups_, v)
228 velocityGroups_[
v]->preSolve();
231 forAll(coalescence_, model)
233 coalescence_[model].correct();
238 breakup_[model].correct();
240 breakup_[model].dsdPtr()().correct();
243 forAll(binaryBreakup_, model)
245 binaryBreakup_[model].correct();
250 drift_[model].correct();
253 forAll(nucleation_, model)
255 nucleation_[model].correct();
260 void Foam::diameterModels::populationBalanceModel::
273 for (
label i = j; i < sizeGroups_.size(); i++)
278 if (Gamma.value() == 0)
continue;
286 0.5*fi.x()*coalescenceRate_()*Gamma
287 *fj*fj.phase()/fj.x()
288 *fk*fk.phase()/fk.x();
293 fi.x()*coalescenceRate_()*Gamma
294 *fj*fj.phase()/fj.x()
295 *fk*fk.phase()/fk.x();
302 const phasePairKey pairij
308 if (pDmdt_.
found(pairij))
310 const scalar dmdtSign
315 pDmdt_[pairij]->ref() += dmdtSign*ratio*Sui_*fi.phase().rho();
318 const phasePairKey pairik
324 if (pDmdt_.
found(pairik))
326 const scalar dmdtSign
331 pDmdt_[pairik]->ref() += dmdtSign*(1 - ratio)*Sui_*fi.phase().rho();
337 void Foam::diameterModels::populationBalanceModel::
347 SuSp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x();
351 SuSp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x();
356 void Foam::diameterModels::populationBalanceModel::
365 for (
label i = 0; i <=
k; i++)
370 fi.x()*breakupRate_()*breakup_[model].dsdPtr()().nik(i, k)
371 *fk*fk.phase()/fk.x();
375 const phasePairKey pair
381 if (pDmdt_.
found(pair))
383 const scalar dmdtSign
388 pDmdt_[pair]->ref() += dmdtSign*Sui_*fi.phase().rho();
394 void Foam::diameterModels::populationBalanceModel::deathByBreakup(
const label i)
398 SuSp_[i] += breakupRate_()*fi.phase();
402 void Foam::diameterModels::populationBalanceModel::calcDeltas()
406 if (delta_[i].empty())
408 for (
label j = 0; j <= sizeGroups_.size() - 1; j++)
416 v_[i+1].value() - v_[i].value()
422 v_[i].value() < 0.5*sizeGroups_[j]->
x().value()
424 0.5*sizeGroups_[j]->
x().value() < v_[i+1].value()
427 delta_[i][j] =
mag(0.5*sizeGroups_[j]->
x() - v_[i]);
429 else if (0.5*sizeGroups_[j]->
x().value() <= v_[i].value())
439 void Foam::diameterModels::populationBalanceModel::
449 Sui_ = fi.x()*binaryBreakupRate_()*delta_[i][j]*fj*fj.phase()/fj.x();
453 const phasePairKey pairij
459 if (pDmdt_.
found(pairij))
461 const scalar dmdtSign
466 pDmdt_[pairij]->ref() += dmdtSign*Sui_*fi.phase().rho();
472 for (
label k = 0; k <= j; k++)
477 if (Gamma.value() == 0)
continue;
484 sizeGroups_[
k]->x()*binaryBreakupRate_()*delta_[i][j]*Gamma
485 *fj*fj.phase()/fj.x();
489 const phasePairKey pairkj
495 if (pDmdt_.
found(pairkj))
497 const scalar dmdtSign
506 pDmdt_[pairkj]->ref() += dmdtSign*Suk*fi.phase().rho();
512 void Foam::diameterModels::populationBalanceModel::
521 SuSp_[i] += alphai*binaryBreakupRate_()*delta_[j][i];
525 void Foam::diameterModels::populationBalanceModel::drift(
const label i)
531 rx_() =
pos(driftRate_())*sizeGroups_[i+1]->x()/sizeGroups_[i]->x();
533 else if (i == sizeGroups_.size() - 1)
535 rx_() =
neg(driftRate_())*sizeGroups_[i-1]->x()/sizeGroups_[i]->x();
540 pos(driftRate_())*sizeGroups_[i+1]->x()/sizeGroups_[i]->x()
541 +
neg(driftRate_())*sizeGroups_[i-1]->x()/sizeGroups_[i]->x();
545 (
neg(1 - rx_()) +
neg(1 - rx_()/(1 - rx_())))*driftRate_()
546 *fi.phase()/((rx_() - 1)*sizeGroups_[i]->
x());
551 if (i < sizeGroups_.size() - 2)
553 rx_() +=
pos(driftRate_())*sizeGroups_[i+2]->x()/sizeGroups_[i+1]->x();
557 *(sizeGroups_[i+2]->x() - sizeGroups_[i+1]->x())
558 /(sizeGroups_[i+1]->
x() - sizeGroups_[i]->x());
560 else if (i == sizeGroups_.size() - 2)
562 rx_() +=
pos(driftRate_())*sizeGroups_[i+1]->x()/sizeGroups_[i]->x();
566 *(sizeGroups_[i+1]->x() - sizeGroups_[i]->x())
567 /(sizeGroups_[i]->
x() - sizeGroups_[i-1]->x());
573 neg(driftRate_())*sizeGroups_[i-1]->x()
574 /sizeGroups_[i]->x();
578 *(sizeGroups_[i]->x() - sizeGroups_[i-1]->x())
579 /(sizeGroups_[i+1]->
x() - sizeGroups_[i]->x());
583 rx_() +=
neg(driftRate_())*sizeGroups_[i-2]->x()/sizeGroups_[i-1]->x();
587 *(sizeGroups_[i-1]->x() - sizeGroups_[i-2]->x())
588 /(sizeGroups_[i]->
x() - sizeGroups_[i-1]->x());
591 if (i != sizeGroups_.size() - 1)
597 pos(driftRate_())*driftRate_()*rdx_()
598 *fi*fi.phase()/fi.x()
603 const phasePairKey pairij
609 if (pDmdt_.
found(pairij))
611 const scalar dmdtSign
616 pDmdt_[pairij]->ref() -= dmdtSign*Suj*fi.phase().rho();
626 neg(driftRate_())*driftRate_()*rdx_()
627 *fi*fi.phase()/fi.x()
632 const phasePairKey pairih
638 if (pDmdt_.
found(pairih))
640 const scalar dmdtSign
645 pDmdt_[pairih]->ref() -= dmdtSign*Suh*fi.phase().rho();
651 void Foam::diameterModels::populationBalanceModel::nucleation(
const label i)
655 Su_[i] += sizeGroups_[i]->x()*nucleationRate_();
659 void Foam::diameterModels::populationBalanceModel::sources()
674 pDmdt_(phasePairIter())->ref() *= 0.0;
684 if (coalescence_.size() != 0)
686 for (
label j = 0; j <= i; j++)
690 if (fi.x() + fj.x() > sizeGroups_.last()->x())
break;
692 coalescenceRate_() *= 0.0;
694 forAll(coalescence_, model)
696 coalescence_[model].addToCoalescenceRate
704 birthByCoalescence(i, j);
706 deathByCoalescence(i, j);
710 if (breakup_.size() != 0)
714 breakup_[model].setBreakupRate(breakupRate_(), i);
716 birthByBreakup(i, model);
722 if (binaryBreakup_.size() != 0)
726 while (delta_[j][i].value() != 0.0)
728 binaryBreakupRate_() *= 0.0;
730 forAll(binaryBreakup_, model)
732 binaryBreakup_[model].addToBinaryBreakupRate
734 binaryBreakupRate_(),
740 birthByBinaryBreakup(j, i);
742 deathByBinaryBreakup(j, i);
748 if (drift_.size() != 0)
754 drift_[model].addToDriftRate(driftRate_(), i);
760 if (nucleation_.size() != 0)
762 nucleationRate_() *= 0.0;
764 forAll(nucleation_, model)
766 nucleation_[model].addToNucleationRate(nucleationRate_(), i);
775 void Foam::diameterModels::populationBalanceModel::dmdt()
777 forAll(velocityGroups_, v)
781 velocityGroups_[
v]->dmdt() *= 0.0;
785 if (&sizeGroups_[i]->phase() == &VelocityGroup.phase())
789 VelocityGroup.dmdt() += fi.phase().rho()*(Su_[i] - SuSp_[i]*fi);
796 void Foam::diameterModels::populationBalanceModel::calcAlphas()
800 forAll(velocityGroups_, v)
802 alphas_ += velocityGroups_[
v]->phase();
807 void Foam::diameterModels::populationBalanceModel::calcVelocity()
811 forAll(velocityGroups_, v)
813 const phaseModel& phase = velocityGroups_[
v]->phase();
815 U_ += phase*phase.U()/
max(alphas_, phase.residualAlpha());
824 const phaseSystem& fluid,
826 HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>& pDmdt
844 fluid.subDict(
"populationBalanceCoeffs").subDict(name_)
846 pimple_(mesh_.lookupObject<pimpleControl>(
"solutionControl")),
849 mesh_.lookupObject<phaseModel>
873 dict_.
lookup(
"coalescenceModels"),
874 coalescenceModel::iNew(*this)
879 dict_.
lookup(
"breakupModels"),
880 breakupModel::iNew(*this)
885 dict_.
lookup(
"binaryBreakupModels"),
886 binaryBreakupModel::iNew(*this)
888 binaryBreakupRate_(),
891 dict_.
lookup(
"driftModels"),
892 driftModel::iNew(*this)
899 dict_.
lookup(
"nucleationModels"),
900 nucleationModel::iNew(*this)
940 this->registerVelocityAndSizeGroups();
942 this->createPhasePairs();
944 if (sizeGroups_.size() < 3)
947 <<
"The populationBalance " << name_
948 <<
" requires a minimum number of three sizeGroups to be" 954 if (coalescence_.size() != 0)
956 coalescenceRate_.reset
963 fluid.time().timeName(),
972 if (breakup_.size() != 0)
981 fluid.time().timeName(),
990 if (binaryBreakup_.size() != 0)
992 binaryBreakupRate_.reset
999 fluid.time().timeName(),
1005 "binaryBreakupRate",
1013 if (drift_.size() != 0)
1022 fluid.time().timeName(),
1037 fluid.time().timeName(),
1052 fluid.time().timeName(),
1061 if (nucleation_.size() != 0)
1063 nucleationRate_.reset
1070 fluid.time().timeName(),
1096 return autoPtr<populationBalanceModel>(
nullptr);
1123 lowerBoundary = sizeGroups_[i-1]->x();
1126 if (i == sizeGroups_.size() - 1)
1132 upperBoundary = sizeGroups_[i+1]->x();
1135 if (v < lowerBoundary || v > upperBoundary)
1139 else if (v.value() <= xi.value())
1141 return (v - lowerBoundary)/(xi - lowerBoundary);
1145 return (upperBoundary - v)/(upperBoundary - xi);
1153 const phaseModel& dispersedPhase
1157 fluid_.lookupSubModel<surfaceTensionModel>
1159 phasePair(dispersedPhase, continuousPhase_)
1173 continuousPhase_.name()
1181 const dictionary& solutionControls = mesh_.solverDict(name_);
1182 bool solveOnFinalIterOnly
1184 solutionControls.lookupOrDefault<
bool>
1186 "solveOnFinalIterOnly",
1194 if (!solveOnFinalIterOnly || pimple_.finalIter())
1199 readScalar(solutionControls.lookup(
"tolerance"))
1208 scalar maxInitialResidual = 1;
1210 while (++iCorr <= nCorr && maxInitialResidual > tolerance)
1212 Info<<
"populationBalance " 1222 maxInitialResidual = 0;
1226 sizeGroup& fi = *sizeGroups_[i];
1227 const phaseModel& phase = fi.phase();
1235 + fi.VelocityGroup().mvConvection()->fvmDiv
1237 phase.alphaRhoPhi(),
1243 - fi.VelocityGroup().dmdt(),
1253 sizeGroupEqn.relax();
1255 maxInitialResidual =
max 1257 sizeGroupEqn.solve().initialResidual(),
1265 forAll(velocityGroups_, i)
1267 velocityGroups_[i]->postSolve();
1273 *sizeGroups_.first()*sizeGroups_.first()->phase()
1278 *sizeGroups_.last()*sizeGroups_.last()->phase()
1281 Info<< this->
name() <<
" sizeGroup phase fraction first, last = " 1282 << fAlpha0.weightedAverage(this->
mesh().V()).value()
1283 <<
' ' << fAlphaN.weightedAverage(this->
mesh().V()).value()
fvMatrix< scalar > fvScalarMatrix
#define forAll(list, i)
Loop across all elements in list.
ThermalDiffusivity< PhaseCompressibleTurbulenceModel< phaseModel > > phaseCompressibleTurbulenceModel
Typedef for phaseCompressibleTurbulenceModel.
tmp< fvMatrix< Type > > SuSp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
const word & name() const
Return name.
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)
const PtrList< dimensionedScalar > & v() const
Return the sizeGroup boundaries.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
static int compare(const Pair< word > &a, const Pair< word > &b)
Compare Pairs.
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.
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
const phasePairTable & phasePairs() const
Return list of unordered phasePairs in this populationBalance.
const dimensionedScalar gamma(const label i, const dimensionedScalar &v) const
Return allocation coefficient.
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
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 > &)
dimensionedScalar neg(const dimensionedScalar &ds)
regIOobject(const IOobject &, const bool isTime=false)
Construct from IOobject. Optional flag for if IOobject is the.
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
autoPtr< populationBalanceModel > clone() const
Return clone.
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Calculate the first temporal derivative.
const phaseModelList & phases() const
Return the phase models.
dimensionedScalar pos(const dimensionedScalar &ds)
const dimensionSet dimVolume(pow3(dimLength))
const Key & key() const
Return the Key corresponding to the iterator.
stressControl lookup("compactNormalStress") >> compactNormalStress
bool found(const Key &) const
Return true if hashedEntry is found in table.
static const word propertiesName
Default name of the turbulence properties dictionary.
static word groupName(Name name, const word &group)
const phaseModel & phase() const
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Calculate the matrix for the first temporal derivative.
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
label readLabel(Istream &is)
Calculate the divergence of the given field.
const phaseCompressibleTurbulenceModel & continuousTurbulence() const
Return reference to turbulence model of the continuous phase.
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/ubuntu/OpenFOAM-6/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
word name(const complex &)
Return a string representation of a complex.
friend class velocityGroup
#define notImplemented(functionName)
Issue a FatalErrorIn for a function not currently implemented.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const Time & time() const
Return time.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
dimensioned< scalar > mag(const dimensioned< Type > &)
const fvMesh & mesh() const
Return reference to the mesh.
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.
A class for managing temporary objects.
void solve()
Solve the population balance equation.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Calculate the matrix for implicit and explicit sources.
const dimensionSet dimVelocity