33 #include "surfaceTensionModel.H" 39 #include "phaseCompressibleTurbulenceModel.H" 43 void Foam::diameterModels::populationBalanceModel::registerVelocityGroups()
47 if (isA<velocityGroup>(fluid_.
phases()[
phasei].dPtr()()))
50 refCast<const velocityGroup>(fluid_.
phases()[
phasei].dPtr()());
52 if (velGroup.popBalName() == this->
name())
54 velocityGroups_.resize(velocityGroups_.size() + 1);
58 velocityGroups_.size() - 1,
62 forAll(velGroup.sizeGroups(), i)
64 this->registerSizeGroups
66 const_cast<sizeGroup&>(velGroup.sizeGroups()[i])
75 void Foam::diameterModels::populationBalanceModel::registerSizeGroups
82 sizeGroups_.size() != 0
84 group.x().value() <= sizeGroups_.last().x().value()
88 <<
"Size groups must be entered according to their representative" 93 sizeGroups_.resize(sizeGroups_.size() + 1);
94 sizeGroups_.set(sizeGroups_.size() - 1, &
group);
97 if (sizeGroups_.size() == 1)
106 sizeGroups_.last().x()
117 sizeGroups_.last().x()
128 sizeGroups_[sizeGroups_.size()-2].x()
129 + sizeGroups_.last().x()
139 sizeGroups_.last().x()
144 delta_.append(
new PtrList<dimensionedScalar>());
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::correct()
222 forAll(velocityGroups_, v)
224 velocityGroups_[
v].preSolve();
227 forAll(coalescence_, model)
229 coalescence_[model].correct();
234 breakup_[model].correct();
236 breakup_[model].dsdPtr()().correct();
239 forAll(binaryBreakup_, model)
241 binaryBreakup_[model].correct();
246 drift_[model].correct();
249 forAll(nucleation_, model)
251 nucleation_[model].correct();
256 void Foam::diameterModels::populationBalanceModel::
269 for (
label i = j; i < sizeGroups_.size(); i++)
274 if (Gamma.value() == 0)
continue;
282 0.5*fi.x()*coalescenceRate_()*Gamma
283 *fj*fj.phase()/fj.x()
284 *fk*fk.phase()/fk.x();
289 fi.x()*coalescenceRate_()*Gamma
290 *fj*fj.phase()/fj.x()
291 *fk*fk.phase()/fk.x();
296 const phasePairKey pairij
302 if (pDmdt_.
found(pairij))
304 const scalar dmdtSign
309 pDmdt_[pairij]->ref() += dmdtSign*fj.x()/v*Sui_*fi.phase().rho();
312 const phasePairKey pairik
318 if (pDmdt_.
found(pairik))
320 const scalar dmdtSign
325 pDmdt_[pairik]->ref() += dmdtSign*fk.x()/v*Sui_*fi.phase().rho();
331 void Foam::diameterModels::populationBalanceModel::
341 SuSp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x();
345 SuSp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x();
350 void Foam::diameterModels::populationBalanceModel::
359 for (
label i = 0; i <=
k; i++)
364 fi.x()*breakupRate_()*breakup_[model].dsdPtr()().nik(i, k)
365 *fk*fk.phase()/fk.x();
369 const phasePairKey pair
375 if (pDmdt_.
found(pair))
377 const scalar dmdtSign
382 pDmdt_[pair]->ref() += dmdtSign*Sui_*fi.phase().rho();
388 void Foam::diameterModels::populationBalanceModel::deathByBreakup(
const label i)
392 SuSp_[i] += breakupRate_()*fi.phase();
396 void Foam::diameterModels::populationBalanceModel::calcDeltas()
400 if (delta_[i].empty())
402 for (
label j = 0; j <= sizeGroups_.size() - 1; j++)
410 v_[i+1].value() - v_[i].value()
416 v_[i].value() < 0.5*sizeGroups_[j].
x().value()
418 0.5*sizeGroups_[j].
x().value() < v_[i+1].value()
421 delta_[i][j] =
mag(0.5*sizeGroups_[j].
x() - v_[i]);
423 else if (0.5*sizeGroups_[j].
x().value() <= v_[i].value())
425 delta_[i][j].value() = 0;
433 void Foam::diameterModels::populationBalanceModel::
443 Sui_ = fi.x()*binaryBreakupRate_()*delta_[i][j]*fj*fj.phase()/fj.x();
447 const phasePairKey pairij
453 if (pDmdt_.
found(pairij))
455 const scalar dmdtSign
460 pDmdt_[pairij]->ref() += dmdtSign*Sui_*fi.phase().rho();
466 for (
label k = 0; k <= j; k++)
471 if (Gamma.value() == 0)
continue;
478 sizeGroups_[
k].x()*binaryBreakupRate_()*delta_[i][j]*Gamma
479 *fj*fj.phase()/fj.x();
483 const phasePairKey pairkj
489 if (pDmdt_.
found(pairkj))
491 const scalar dmdtSign
500 pDmdt_[pairkj]->ref() += dmdtSign*Suk*fi.phase().rho();
506 void Foam::diameterModels::populationBalanceModel::
515 SuSp_[i] += alphai*binaryBreakupRate_()*delta_[j][i];
519 void Foam::diameterModels::populationBalanceModel::drift(
const label i)
525 rx_() =
pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x();
527 else if (i == sizeGroups_.size() - 1)
529 rx_() =
neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
534 pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x()
535 +
neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
539 (
neg(1 - rx_()) +
neg(rx_() - rx_()/(1 - rx_())))*driftRate_()
540 *fp.phase()/((rx_() - 1)*fp.x());
545 if (i == sizeGroups_.size() - 2)
547 rx_() =
pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x();
551 *(sizeGroups_[i+1].x() - sizeGroups_[i].x())
552 /(sizeGroups_[i].
x() - sizeGroups_[i-1].x());
554 else if (i < sizeGroups_.size() - 2)
556 rx_() =
pos(driftRate_())*sizeGroups_[i+2].x()/sizeGroups_[i+1].x();
560 *(sizeGroups_[i+2].x() - sizeGroups_[i+1].x())
561 /(sizeGroups_[i+1].
x() - sizeGroups_[i].x());
566 rx_() +=
neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
570 *(sizeGroups_[i].x() - sizeGroups_[i-1].x())
571 /(sizeGroups_[i+1].
x() - sizeGroups_[i].x());
575 rx_() +=
neg(driftRate_())*sizeGroups_[i-2].x()/sizeGroups_[i-1].x();
579 *(sizeGroups_[i-1].x() - sizeGroups_[i-2].x())
580 /(sizeGroups_[i].
x() - sizeGroups_[i-1].x());
583 if (i != sizeGroups_.size() - 1)
589 pos(driftRate_())*driftRate_()*rdx_()
590 *fp*fp.phase()/fp.x()
595 const phasePairKey pairij
601 if (pDmdt_.
found(pairij))
603 const scalar dmdtSign
608 pDmdt_[pairij]->ref() -= dmdtSign*Sue*fp.phase().rho();
618 neg(driftRate_())*driftRate_()*rdx_()
619 *fp*fp.phase()/fp.x()
624 const phasePairKey pairih
630 if (pDmdt_.
found(pairih))
632 const scalar dmdtSign
637 pDmdt_[pairih]->ref() -= dmdtSign*Suw*fp.phase().rho();
643 void Foam::diameterModels::populationBalanceModel::nucleation(
const label i)
645 Su_[i] += sizeGroups_[i].x()*nucleationRate_();
649 void Foam::diameterModels::populationBalanceModel::sources()
664 pDmdt_(phasePairIter())->ref() =
Zero;
674 if (coalescence_.size() != 0)
676 for (
label j = 0; j <= i; j++)
680 if (fi.x() + fj.x() > sizeGroups_.last().x())
break;
682 coalescenceRate_() =
Zero;
684 forAll(coalescence_, model)
686 coalescence_[model].addToCoalescenceRate
694 birthByCoalescence(i, j);
696 deathByCoalescence(i, j);
700 if (breakup_.size() != 0)
704 breakup_[model].setBreakupRate(breakupRate_(), i);
706 birthByBreakup(i, model);
712 if (binaryBreakup_.size() != 0)
716 while (delta_[j][i].value() != 0)
718 binaryBreakupRate_() =
Zero;
720 forAll(binaryBreakup_, model)
722 binaryBreakup_[model].addToBinaryBreakupRate
724 binaryBreakupRate_(),
730 birthByBinaryBreakup(j, i);
732 deathByBinaryBreakup(j, i);
738 if (drift_.size() != 0)
744 drift_[model].addToDriftRate(driftRate_(), i);
750 if (nucleation_.size() != 0)
752 nucleationRate_() =
Zero;
754 forAll(nucleation_, model)
756 nucleation_[model].addToNucleationRate(nucleationRate_(), i);
765 void Foam::diameterModels::populationBalanceModel::dmdt()
767 forAll(velocityGroups_, v)
771 velGroup.dmdtRef() =
Zero;
775 if (&sizeGroups_[i].phase() == &velGroup.phase())
779 velGroup.dmdtRef() += fi.phase().rho()*(Su_[i] - SuSp_[i]*fi);
786 void Foam::diameterModels::populationBalanceModel::calcAlphas()
790 forAll(velocityGroups_, v)
792 const phaseModel& phase = velocityGroups_[
v].phase();
794 alphas_() +=
max(phase, phase.residualAlpha());
800 Foam::diameterModels::populationBalanceModel::calcDsm()
802 tmp<volScalarField> tInvDsm
814 forAll(velocityGroups_, v)
816 const phaseModel& phase = velocityGroups_[
v].phase();
818 invDsm +=
max(phase, phase.residualAlpha())/(phase.d()*alphas_());
825 void Foam::diameterModels::populationBalanceModel::calcVelocity()
829 forAll(velocityGroups_, v)
831 const phaseModel& phase = velocityGroups_[
v].phase();
833 U_() +=
max(phase, phase.residualAlpha())*phase.U()/alphas_();
837 bool Foam::diameterModels::populationBalanceModel::updateSources()
839 const bool result = sourceUpdateCounter_ % sourceUpdateInterval() == 0;
841 ++ sourceUpdateCounter_;
851 const phaseSystem& fluid,
853 HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>& pDmdt
861 fluid.
time().constant(),
867 mesh_(fluid_.
mesh()),
871 fluid_.subDict(
"populationBalanceCoeffs").subDict(name_)
873 pimple_(mesh_.lookupObject<pimpleControl>(
"solutionControl")),
876 mesh_.lookupObject<phaseModel>
900 dict_.
lookup(
"coalescenceModels"),
901 coalescenceModel::iNew(*this)
906 dict_.
lookup(
"breakupModels"),
907 breakupModel::iNew(*this)
912 dict_.
lookup(
"binaryBreakupModels"),
913 binaryBreakupModel::iNew(*this)
915 binaryBreakupRate_(),
918 dict_.
lookup(
"driftModels"),
919 driftModel::iNew(*this)
926 dict_.
lookup(
"nucleationModels"),
927 nucleationModel::iNew(*this)
938 this->registerVelocityGroups();
940 this->createPhasePairs();
942 if (sizeGroups_.size() < 3)
945 <<
"The populationBalance " << name_
946 <<
" requires a minimum number of three sizeGroups to be specified." 951 if (coalescence_.size() != 0)
960 mesh_.time().timeName(),
969 if (breakup_.size() != 0)
978 fluid_.time().timeName(),
987 if (binaryBreakup_.size() != 0)
989 binaryBreakupRate_.set
996 fluid_.time().timeName(),
1002 "binaryBreakupRate",
1010 if (drift_.size() != 0)
1019 fluid_.time().timeName(),
1034 fluid_.time().timeName(),
1049 fluid_.time().timeName(),
1058 if (nucleation_.size() != 0)
1067 fluid_.time().timeName(),
1081 if (velocityGroups_.size() > 1)
1090 fluid_.time().timeName(),
1107 fluid_.time().timeName(),
1124 fluid_.time().timeName(),
1147 return autoPtr<populationBalanceModel>(
nullptr);
1174 lowerBoundary = sizeGroups_[i-1].x();
1177 if (i == sizeGroups_.size() - 1)
1183 upperBoundary = sizeGroups_[i+1].x();
1186 if (v < lowerBoundary || v > upperBoundary)
1190 else if (v.value() <= xi.value())
1192 return (v - lowerBoundary)/(xi - lowerBoundary);
1196 return (upperBoundary - v)/(upperBoundary - xi);
1204 const phaseModel& dispersedPhase
1208 fluid_.lookupSubModel<surfaceTensionModel>
1210 phasePair(dispersedPhase, continuousPhase_)
1224 continuousPhase_.name()
1232 const dictionary& solutionControls = mesh_.solverDict(name_);
1233 bool solveOnFinalIterOnly =
1234 solutionControls.lookupOrDefault<
bool>(
"solveOnFinalIterOnly",
false);
1236 if (!solveOnFinalIterOnly || pimple_.finalPimpleIter())
1239 const scalar tolerance =
1240 readScalar(solutionControls.lookup(
"tolerance"));
1248 scalar maxInitialResidual = 1;
1250 while (++iCorr <= nCorr && maxInitialResidual > tolerance)
1252 Info<<
"populationBalance " 1258 if (updateSources())
1265 maxInitialResidual = 0;
1269 sizeGroup& fi = sizeGroups_[i];
1270 const phaseModel& phase = fi.phase();
1278 + fi.VelocityGroup().mvConvection()->fvmDiv
1280 phase.alphaRhoPhi(),
1286 - fi.VelocityGroup().dmdt(),
1296 sizeGroupEqn.relax();
1298 maxInitialResidual =
max 1300 sizeGroupEqn.solve().initialResidual(),
1308 forAll(velocityGroups_, i)
1310 velocityGroups_[i].postSolve();
1314 if (velocityGroups_.size() > 1)
1323 sizeGroups_.first()*sizeGroups_.first().phase()
1328 sizeGroups_.last()*sizeGroups_.last().phase()
1331 Info<< this->
name() <<
" sizeGroup phase fraction first, last = " 1332 << fAlpha0.weightedAverage(this->
mesh().V()).value()
1333 <<
' ' << 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/m^2/K^4].
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.
word group() const
Return group (extension part of name)
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)
regIOobject(const IOobject &, const bool isTime=false)
Construct from IOobject. Optional flag for if IOobject is the.
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Mesh &, const dimensionSet &, const word &patchFieldType=fvPatchField< scalar >::calculatedType())
Return a temporary field constructed from name, mesh, dimensionSet.
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.
Calculate the field for explicit evaluation of implicit and explicit sources.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
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(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-7/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
Internal & ref()
Return a reference to the dimensioned internal field.
#define notImplemented(functionName)
Issue a FatalErrorIn for a function not currently implemented.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
const dimensionSet dimLength(0, 1, 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.
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
Calculate the matrix for implicit and explicit sources.
const dimensionSet dimVelocity