36 #include "surfaceTensionModel.H" 42 #include "phaseCompressibleTurbulenceModel.H" 46 void Foam::diameterModels::populationBalanceModel::registerVelocityGroups()
50 if (isA<velocityGroup>(fluid_.
phases()[
phasei].dPtr()()))
53 refCast<const velocityGroup>(fluid_.
phases()[
phasei].dPtr()());
55 if (velGroup.popBalName() == this->
name())
57 velocityGroups_.resize(velocityGroups_.size() + 1);
61 velocityGroups_.size() - 1,
65 forAll(velGroup.sizeGroups(), i)
67 this->registerSizeGroups
69 const_cast<sizeGroup&>(velGroup.sizeGroups()[i])
78 void Foam::diameterModels::populationBalanceModel::registerSizeGroups
85 sizeGroups_.size() != 0
87 group.x().value() <= sizeGroups_.last().x().value()
91 <<
"Size groups must be entered according to their representative" 96 sizeGroups_.resize(sizeGroups_.size() + 1);
97 sizeGroups_.set(sizeGroups_.size() - 1, &
group);
100 if (sizeGroups_.size() == 1)
109 sizeGroups_.last().x()
120 sizeGroups_.last().x()
131 sizeGroups_[sizeGroups_.size()-2].x()
132 + sizeGroups_.last().x()
142 sizeGroups_.last().x()
147 delta_.append(
new PtrList<dimensionedScalar>());
173 void Foam::diameterModels::populationBalanceModel::createPhasePairs()
175 forAll(velocityGroups_, i)
177 const phaseModel&
phasei = velocityGroups_[i].phase();
179 forAll(velocityGroups_, j)
181 const phaseModel& phasej = velocityGroups_[j].phase();
185 const phasePairKey
key 192 if (!phasePairs_.found(
key))
213 void Foam::diameterModels::populationBalanceModel::correct()
217 forAll(velocityGroups_, v)
219 velocityGroups_[v].preSolve();
222 forAll(coalescence_, model)
224 coalescence_[model].correct();
229 breakup_[model].correct();
231 breakup_[model].dsdPtr()().
correct();
234 forAll(binaryBreakup_, model)
236 binaryBreakup_[model].correct();
241 drift_[model].correct();
244 forAll(nucleation_, model)
246 nucleation_[model].correct();
251 void Foam::diameterModels::populationBalanceModel::
258 const sizeGroup& fj = sizeGroups_[j];
259 const sizeGroup& fk = sizeGroups_[
k];
264 for (label i = j; i < sizeGroups_.size(); i++)
269 if (Gamma.value() == 0)
continue;
271 const sizeGroup& fi = sizeGroups_[i];
277 0.5*fi.x()*coalescenceRate_()*Gamma
278 *fj*fj.phase()/fj.x()
279 *fk*fk.phase()/fk.x();
284 fi.x()*coalescenceRate_()*Gamma
285 *fj*fj.phase()/fj.x()
286 *fk*fk.phase()/fk.x();
291 const phasePairKey pairij
297 if (pDmdt_.found(pairij))
299 const scalar dmdtSign
304 pDmdt_[pairij]->ref() += dmdtSign*fj.x()/v*Sui_*fi.phase().rho();
307 const phasePairKey pairik
313 if (pDmdt_.found(pairik))
315 const scalar dmdtSign
320 pDmdt_[pairik]->ref() += dmdtSign*fk.x()/v*Sui_*fi.phase().rho();
326 void Foam::diameterModels::populationBalanceModel::
333 const sizeGroup& fi = sizeGroups_[i];
334 const sizeGroup& fj = sizeGroups_[j];
336 SuSp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x();
340 SuSp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x();
345 void Foam::diameterModels::populationBalanceModel::
352 const sizeGroup& fk = sizeGroups_[
k];
354 for (label i = 0; i <=
k; i++)
356 const sizeGroup& fi = sizeGroups_[i];
359 fi.x()*breakupRate_()*breakup_[model].dsdPtr()().nik(i,
k)
360 *fk*fk.phase()/fk.x();
364 const phasePairKey pair
370 if (pDmdt_.found(pair))
372 const scalar dmdtSign
377 pDmdt_[pair]->ref() += dmdtSign*Sui_*fi.phase().rho();
383 void Foam::diameterModels::populationBalanceModel::deathByBreakup(
const label i)
385 const sizeGroup& fi = sizeGroups_[i];
387 SuSp_[i] += breakupRate_()*fi.phase();
391 void Foam::diameterModels::populationBalanceModel::calcDeltas()
395 if (delta_[i].empty())
397 for (label j = 0; j <= sizeGroups_.size() - 1; j++)
405 v_[i+1].value() - v_[i].value()
411 v_[i].value() < 0.5*sizeGroups_[j].
x().value()
413 0.5*sizeGroups_[j].
x().value() < v_[i+1].value()
416 delta_[i][j] =
mag(0.5*sizeGroups_[j].
x() - v_[i]);
418 else if (0.5*sizeGroups_[j].
x().value() <= v_[i].value())
420 delta_[i][j].value() = 0;
428 void Foam::diameterModels::populationBalanceModel::
435 const sizeGroup& fj = sizeGroups_[j];
436 const sizeGroup& fi = sizeGroups_[i];
438 Sui_ = fi.x()*binaryBreakupRate_()*delta_[i][j]*fj*fj.phase()/fj.x();
442 const phasePairKey pairij
448 if (pDmdt_.found(pairij))
450 const scalar dmdtSign
455 pDmdt_[pairij]->ref() += dmdtSign*Sui_*fi.phase().rho();
461 for (label
k = 0;
k <= j;
k++)
466 if (Gamma.value() == 0)
continue;
468 const sizeGroup& fk = sizeGroups_[
k];
473 sizeGroups_[
k].x()*binaryBreakupRate_()*delta_[i][j]*Gamma
474 *fj*fj.phase()/fj.x();
478 const phasePairKey pairkj
484 if (pDmdt_.found(pairkj))
486 const scalar dmdtSign
490 pDmdt_.find(pairkj).key(),
495 pDmdt_[pairkj]->ref() += dmdtSign*Suk*fi.phase().rho();
501 void Foam::diameterModels::populationBalanceModel::
510 SuSp_[i] += alphai*binaryBreakupRate_()*delta_[j][i];
514 void Foam::diameterModels::populationBalanceModel::drift(
const label i)
516 const sizeGroup& fp = sizeGroups_[i];
520 rx_() =
pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x();
522 else if (i == sizeGroups_.size() - 1)
524 rx_() =
neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
529 pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x()
530 +
neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
534 (
neg(1 - rx_()) +
neg(rx_() - rx_()/(1 - rx_())))*driftRate_()
535 *fp.phase()/((rx_() - 1)*fp.x());
540 if (i == sizeGroups_.size() - 2)
542 rx_() =
pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x();
546 *(sizeGroups_[i+1].x() - sizeGroups_[i].x())
547 /(sizeGroups_[i].
x() - sizeGroups_[i-1].x());
549 else if (i < sizeGroups_.size() - 2)
551 rx_() =
pos(driftRate_())*sizeGroups_[i+2].x()/sizeGroups_[i+1].x();
555 *(sizeGroups_[i+2].x() - sizeGroups_[i+1].x())
556 /(sizeGroups_[i+1].
x() - sizeGroups_[i].x());
561 rx_() +=
neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x();
565 *(sizeGroups_[i].x() - sizeGroups_[i-1].x())
566 /(sizeGroups_[i+1].
x() - sizeGroups_[i].x());
570 rx_() +=
neg(driftRate_())*sizeGroups_[i-2].x()/sizeGroups_[i-1].x();
574 *(sizeGroups_[i-1].x() - sizeGroups_[i-2].x())
575 /(sizeGroups_[i].
x() - sizeGroups_[i-1].x());
578 if (i != sizeGroups_.size() - 1)
580 const sizeGroup& fe = sizeGroups_[i+1];
584 pos(driftRate_())*driftRate_()*rdx_()
585 *fp*fp.phase()/fp.x()
590 const phasePairKey pairij
596 if (pDmdt_.found(pairij))
598 const scalar dmdtSign
603 pDmdt_[pairij]->ref() -= dmdtSign*Sue*fp.phase().rho();
609 const sizeGroup& fw = sizeGroups_[i-1];
613 neg(driftRate_())*driftRate_()*rdx_()
614 *fp*fp.phase()/fp.x()
619 const phasePairKey pairih
625 if (pDmdt_.found(pairih))
627 const scalar dmdtSign
632 pDmdt_[pairih]->ref() -= dmdtSign*Suw*fp.phase().rho();
638 void Foam::diameterModels::populationBalanceModel::nucleation(
const label i)
640 Su_[i] += sizeGroups_[i].x()*nucleationRate_();
644 void Foam::diameterModels::populationBalanceModel::sources()
659 pDmdt_(phasePairIter())->ref() *= 0.0;
667 const sizeGroup& fi = sizeGroups_[i];
669 if (coalescence_.size() != 0)
671 for (label j = 0; j <= i; j++)
673 const sizeGroup& fj = sizeGroups_[j];
675 if (fi.x() + fj.x() > sizeGroups_.last().x())
break;
677 coalescenceRate_() *= 0.0;
679 forAll(coalescence_, model)
681 coalescence_[model].addToCoalescenceRate
689 birthByCoalescence(i, j);
691 deathByCoalescence(i, j);
695 if (breakup_.size() != 0)
699 breakup_[model].setBreakupRate(breakupRate_(), i);
701 birthByBreakup(i, model);
707 if (binaryBreakup_.size() != 0)
711 while (delta_[j][i].value() != 0)
713 binaryBreakupRate_() *= 0.0;
715 forAll(binaryBreakup_, model)
717 binaryBreakup_[model].addToBinaryBreakupRate
719 binaryBreakupRate_(),
725 birthByBinaryBreakup(j, i);
727 deathByBinaryBreakup(j, i);
733 if (drift_.size() != 0)
739 drift_[model].addToDriftRate(driftRate_(), i);
745 if (nucleation_.size() != 0)
747 nucleationRate_() *= 0.0;
749 forAll(nucleation_, model)
751 nucleation_[model].addToNucleationRate(nucleationRate_(), i);
760 void Foam::diameterModels::populationBalanceModel::dmdt()
762 forAll(velocityGroups_, v)
764 velocityGroup& velGroup = velocityGroups_[v];
766 velGroup.dmdtRef() *= 0.0;
770 if (&sizeGroups_[i].phase() == &velGroup.phase())
772 sizeGroup& fi = sizeGroups_[i];
774 velGroup.dmdtRef() += fi.phase().rho()*(Su_[i] - SuSp_[i]*fi);
781 void Foam::diameterModels::populationBalanceModel::calcAlphas()
785 forAll(velocityGroups_, v)
787 const phaseModel& phase = velocityGroups_[v].phase();
789 alphas_() +=
max(phase, phase.residualAlpha());
795 Foam::diameterModels::populationBalanceModel::calcDsm()
797 tmp<volScalarField> tInvDsm
809 forAll(velocityGroups_, v)
811 const phaseModel& phase = velocityGroups_[v].phase();
813 invDsm +=
max(phase, phase.residualAlpha())/(phase.d()*alphas_());
820 void Foam::diameterModels::populationBalanceModel::calcVelocity()
824 forAll(velocityGroups_, v)
826 const phaseModel& phase = velocityGroups_[v].phase();
828 U_() +=
max(phase, phase.residualAlpha())*phase.U()/alphas_();
832 bool Foam::diameterModels::populationBalanceModel::updateSources()
834 const bool result = sourceUpdateCounter_ % sourceUpdateInterval() == 0;
836 ++ sourceUpdateCounter_;
862 mesh_(fluid_.
mesh()),
866 fluid_.subDict(
"populationBalanceCoeffs").subDict(name_)
868 pimple_(mesh_.lookupObject<
pimpleControl>(
"solutionControl")),
895 dict_.
lookup(
"coalescenceModels"),
901 dict_.
lookup(
"breakupModels"),
907 dict_.
lookup(
"binaryBreakupModels"),
910 binaryBreakupRate_(),
913 dict_.
lookup(
"driftModels"),
921 dict_.
lookup(
"nucleationModels"),
933 this->registerVelocityGroups();
935 this->createPhasePairs();
937 if (sizeGroups_.size() < 3)
940 <<
"The populationBalance " << name_
941 <<
" requires a minimum number of three sizeGroups to be specified." 946 if (coalescence_.size() != 0)
948 coalescenceRate_.reset
959 if (breakup_.size() != 0)
972 if (binaryBreakup_.size() != 0)
974 binaryBreakupRate_.reset
985 if (drift_.size() != 0)
1018 if (nucleation_.size() != 0)
1020 nucleationRate_.reset
1031 if (velocityGroups_.size() > 1)
1128 lowerBoundary = sizeGroups_[i-1].x();
1131 if (i == sizeGroups_.size() - 1)
1137 upperBoundary = sizeGroups_[i+1].x();
1140 if (v < lowerBoundary || v > upperBoundary)
1146 return (v - lowerBoundary)/(xi - lowerBoundary);
1150 return (upperBoundary - v)/(upperBoundary - xi);
1164 phasePair(dispersedPhase, continuousPhase_)
1178 continuousPhase_.name()
1186 const dictionary& solutionControls = mesh_.solverDict(name_);
1187 const bool solveOnFinalIterOnly =
1188 solutionControls.
getOrDefault(
"solveOnFinalIterOnly",
false);
1190 if (!solveOnFinalIterOnly || pimple_.finalIter())
1193 const scalar tolerance = solutionControls.
get<scalar>(
"tolerance");
1201 scalar maxInitialResidual = 1;
1203 while (++iCorr <= nCorr && maxInitialResidual > tolerance)
1205 Info<<
"populationBalance " 1211 if (updateSources())
1218 maxInitialResidual = 0;
1222 sizeGroup& fi = sizeGroups_[i];
1223 const phaseModel& phase = fi.phase();
1226 const tmp<volScalarField>
trho(phase.thermo().rho());
1232 + fi.VelocityGroup().mvConvection()->fvmDiv
1234 phase.alphaRhoPhi(),
1240 - fi.VelocityGroup().dmdt(),
1250 sizeGroupEqn.relax();
1252 maxInitialResidual =
max 1254 sizeGroupEqn.solve().initialResidual(),
1262 forAll(velocityGroups_, i)
1264 velocityGroups_[i].postSolve();
1268 if (velocityGroups_.size() > 1)
1277 sizeGroups_.first()*sizeGroups_.first().phase()
1282 sizeGroups_.last()*sizeGroups_.last().phase()
1285 Info<< this->
name() <<
" sizeGroup phase fraction first, last = " 1286 << fAlpha0.weightedAverage(this->
mesh().V()).value()
1287 <<
' ' << fAlphaN.weightedAverage(this->
mesh().V()).value()
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
const Type & value() const noexcept
Return const reference to value.
fvMatrix< scalar > fvScalarMatrix
Base class for drift models.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
populationBalanceModel(const phaseSystem &fluid, const word &name, HashPtrTable< volScalarField, phasePairKey, phasePairKey::hash > &pDmdt)
const word & name() const noexcept
Return the object name.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< volScalarField > trho
const tmp< volScalarField > sigmaWithContinuousPhase(const phaseModel &dispersedPhase) const
Return the surface tension coefficient between a given dispersed.
label k
Boltzmann constant.
const dimensionSet dimless
Dimensionless.
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T. FatalIOError if not found, or if the number of tokens is incorrect.
bool writeData(Ostream &) const
Dummy write for regIOobject.
Lookup type of boundary radiation properties.
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object.
GeometricField< vector, fvPatchField, volMesh > volVectorField
dimensionedScalar sigma("sigma", dimMass/sqr(dimTime), transportProperties)
dimensionedScalar neg(const dimensionedScalar &ds)
autoPtr< populationBalanceModel > clone() const
Return clone.
Base class for coalescence models.
#define forAll(list, i)
Loop across all elements in list.
constexpr const char *const group
Group name for atomic constants.
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Base class for breakup models which give a total breakup rate and a separate daughter size distributi...
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Calculate the first temporal derivative.
Base class for binary breakup models which give the breakup rate between a sizeGroup pair directly...
const phaseModelList & phases() const
Return the phase models.
dimensionedScalar pos(const dimensionedScalar &ds)
const dimensionSet dimVolume(pow3(dimLength))
static word groupName(StringType base, const word &group)
Create dot-delimited name.group string.
Base class for nucleation models.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
static int compare(const Pair< word > &a, const Pair< word > &b)
Compare Pairs.
Class to represent a system of phases and model interfacial transfers between them.
static const word propertiesName
Default name of the turbulence properties dictionary.
A class for handling words, derived from Foam::string.
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Calculate the matrix for the first temporal derivative.
zeroField SuSp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
Calculate the field for explicit evaluation of implicit and explicit sources.
Constant dispersed-phase particle diameter model.
zeroField Sp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
const Time & time() const noexcept
Return Time associated with the objectRegistry.
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, IOobjectOption::registerOption regOpt, const Mesh &mesh, const dimensionSet &dims, const word &patchFieldType=fvPatchField< scalar >::calculatedType())
Return tmp field (NO_READ, NO_WRITE) from name, mesh, dimensions and patch type. [Takes current timeN...
Calculate the divergence of the given field.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
static word timeName(const scalar t, const int precision=precision_)
Return a time name for the given scalar time value formatted with the given precision.
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/chef2/andy/OpenFOAM/release/v2406/OpenFOAM-v2406/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=clamp((rho - rholSat)/(rhovSat - rholSat), zero_one{});alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
OBJstream os(runTime.globalPath()/outputName)
friend class velocityGroup
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
PIMPLE control class to supply convergence information/checks for the PIMPLE loop.
Automatically write from objectRegistry::writeObject()
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
regIOobject is an abstract class derived from IOobject to handle automatic object registration with t...
messageStream Info
Information stream (stdout output on master, null elsewhere)
const phaseCompressibleTurbulenceModel & continuousTurbulence() const
Return reference to turbulence model of the continuous phase.
Internal & ref(const bool updateAccessTime=true)
Same as internalFieldRef()
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T, or return the given default value. FatalIOError if it is found and the number of...
IOobject newIOobject(const word &name, IOobjectOption ioOpt) const
Create an IOobject at the current time instance (timeName) with the specified options.
A class for managing temporary objects.
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
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)
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Defines the attributes of an object for which implicit objectRegistry management is supported...
const dimensionedScalar gamma(const label i, const dimensionedScalar &v) const
Return allocation coefficient.
Request registration (bool: true)
Single incompressible phase derived from the phase-fraction. Used as part of the multiPhaseMixture fo...
Description for mass transfer between a pair of phases. The direction of the mass transfer is from th...
Do not request registration (bool: false)
Calculate the finiteVolume matrix for implicit and explicit sources.
virtual ~populationBalanceModel()
Destructor.
static constexpr const zero Zero
Global zero (0)
const dimensionSet dimVelocity