From 1a3ae93fb372e1d49aac51e7bd0b9cbfc4ceb0b6 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Thu, 3 Jun 2021 16:57:11 +0100 Subject: [PATCH 1/4] multiphaseEulerFoam: Ensure mass transfer rates stay registered --- ...terfaceCompositionPhaseChangePhaseSystem.C | 34 +++++-------------- ...terfaceCompositionPhaseChangePhaseSystem.H | 4 +-- 2 files changed, 10 insertions(+), 28 deletions(-) diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C index 59f742d3db..65ae173efa 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C @@ -32,16 +32,9 @@ License // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // template -Foam::autoPtr -Foam::InterfaceCompositionPhaseChangePhaseSystem:: -totalDmdtfs() const +void Foam::InterfaceCompositionPhaseChangePhaseSystem:: +correctDmdtfs() { - autoPtr totalDmdtfsPtr - ( - new phaseSystem::dmdtfTable - ); - phaseSystem::dmdtfTable& totalDmdtfs = totalDmdtfsPtr(); - forAllConstIter ( interfaceCompositionModelTable, @@ -52,6 +45,8 @@ totalDmdtfs() const const phasePair& pair = this->phasePairs_[interfaceCompositionModelIter.key()]; + *dmdtfs_[pair] = Zero; + forAllConstIter(phasePair, pair, pairIter) { const autoPtr& compositionModelPtr = @@ -73,28 +68,15 @@ totalDmdtfs() const { const word& specie = *specieIter; - tmp dmidtf - ( + *dmdtfs_[pair] += (pairIter.index() == 0 ? +1 : -1) *( *(*dmidtfSus_[pair])[specie] + *(*dmidtfSps_[pair])[specie]*phase.Y(specie) - ) - ); - - if (totalDmdtfs.found(pair)) - { - *totalDmdtfs[pair] += dmidtf; - } - else - { - totalDmdtfs.insert(pair, dmidtf.ptr()); - } + ); } } } - - return totalDmdtfsPtr; } @@ -610,7 +592,7 @@ correct() } } - dmdtfs_ = totalDmdtfs(); + correctDmdtfs(); } @@ -620,7 +602,7 @@ correctSpecies() { BasePhaseSystem::correctSpecies(); - dmdtfs_ = totalDmdtfs(); + correctDmdtfs(); } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H index 8ba98dc2e8..33f69aeb67 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H @@ -108,8 +108,8 @@ class InterfaceCompositionPhaseChangePhaseSystem // Private Member Functions - //- Return mass transfers across each interface - autoPtr totalDmdtfs() const; + //- Correct mass transfers across each interface + void correctDmdtfs(); //- Return species mass transfers across each interface autoPtr totalDmidtfs() const; From 8fec877d5bacfe0ae9985740b2a9ad60f94908c0 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Fri, 4 Jun 2021 10:51:38 +0100 Subject: [PATCH 2/4] multiphaseEulerFoam: populationBalance: Simplifed drift term implementation Patch contributed by Institute of Fluid Dynamics, Helmholtz-Zentrum Dresden - Rossendorf (HZDR) --- .../sizeGroup/shapeModels/fractal/fractal.C | 4 +- .../populationBalanceModel.C | 144 ++++-------------- .../populationBalanceModel.H | 14 +- .../populationBalanceModelI.H | 6 +- 4 files changed, 41 insertions(+), 127 deletions(-) diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C index 0229f7b35e..afd0c95eb8 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2019-2020 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2019-2021 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -216,7 +216,7 @@ void Foam::diameterModels::shapeModels::fractal::correct() == - sinteringModel_->R() + Su_ - - fvm::SuSp(popBal.SuSp(fi.i()())*fi, kappa_) + - fvm::Sp(popBal.Sp(fi.i()())*fi, kappa_) + fvc::ddt(fi.phase().residualAlpha(), kappa_) - fvm::ddt(fi.phase().residualAlpha(), kappa_) ); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C index 3240198eea..66751ccc56 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C @@ -167,13 +167,13 @@ void Foam::diameterModels::populationBalanceModel::registerSizeGroups ) ); - SuSp_.append + Sp_.append ( new volScalarField ( IOobject ( - "SuSp", + "Sp", fluid_.time().timeName(), mesh_ ), @@ -344,11 +344,11 @@ deathByCoalescence const sizeGroup& fi = sizeGroups_[i]; const sizeGroup& fj = sizeGroups_[j]; - SuSp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x(); + Sp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x(); if (i != j) { - SuSp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x(); + Sp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x(); } } @@ -397,7 +397,7 @@ void Foam::diameterModels::populationBalanceModel::deathByBreakup(const label i) { const sizeGroup& fi = sizeGroups_[i]; - SuSp_[i] += breakupRate_()*fi.phase(); + Sp_[i] += breakupRate_()*fi.phase(); } @@ -524,7 +524,7 @@ deathByBinaryBreakup { const volScalarField& alphai = sizeGroups_[i].phase(); - SuSp_[i] += alphai*binaryBreakupRate_()*delta_[j][i]; + Sp_[i] += alphai*binaryBreakupRate_()*delta_[j][i]; } @@ -536,77 +536,18 @@ void Foam::diameterModels::populationBalanceModel::drift { model.addToDriftRate(driftRate_(), i); - const sizeGroup& fp = sizeGroups_[i]; + const sizeGroup& fp = sizeGroups()[i]; - if (i == 0) + if (i < sizeGroups().size() - 1) { - rx_() = pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x(); - } - else if (i == sizeGroups_.size() - 1) - { - rx_() = neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x(); - } - else - { - rx_() = - pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x() - + neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x(); - } - - SuSp_[i] += - (neg(1 - rx_()) + neg(rx_() - rx_()/(1 - rx_())))*driftRate_() - *fp.phase()/((rx_() - 1)*fp.x()); - - rx_() = Zero; - rdx_() = Zero; - - if (i == sizeGroups_.size() - 2) - { - rx_() = pos(driftRate_())*sizeGroups_[i+1].x()/sizeGroups_[i].x(); - - rdx_() = - pos(driftRate_()) - *(sizeGroups_[i+1].x() - sizeGroups_[i].x()) - /(sizeGroups_[i].x() - sizeGroups_[i-1].x()); - } - else if (i < sizeGroups_.size() - 2) - { - rx_() = pos(driftRate_())*sizeGroups_[i+2].x()/sizeGroups_[i+1].x(); - - rdx_() = - pos(driftRate_()) - *(sizeGroups_[i+2].x() - sizeGroups_[i+1].x()) - /(sizeGroups_[i+1].x() - sizeGroups_[i].x()); - } - - if (i == 1) - { - rx_() += neg(driftRate_())*sizeGroups_[i-1].x()/sizeGroups_[i].x(); - - rdx_() += - neg(driftRate_()) - *(sizeGroups_[i].x() - sizeGroups_[i-1].x()) - /(sizeGroups_[i+1].x() - sizeGroups_[i].x()); - } - else if (i > 1) - { - rx_() += neg(driftRate_())*sizeGroups_[i-2].x()/sizeGroups_[i-1].x(); - - rdx_() += - neg(driftRate_()) - *(sizeGroups_[i-1].x() - sizeGroups_[i-2].x()) - /(sizeGroups_[i].x() - sizeGroups_[i-1].x()); - } - - if (i != sizeGroups_.size() - 1) - { - const sizeGroup& fe = sizeGroups_[i+1]; + const sizeGroup& fe = sizeGroups()[i+1]; volScalarField& Sue = Sui_; + Sp_[i] += 1/(fe.x() - fp.x())*pos(driftRate_())*driftRate_()*fp.phase(); + Sue = - pos(driftRate_())*driftRate_()*rdx_() - *fp*fp.phase()/fp.x() - /(rx_() - 1); + fe.x()/(fp.x()*(fe.x() - fp.x())) + *pos(driftRate_())*driftRate_()*fp*fp.phase(); Su_[i+1] += Sue; @@ -629,15 +570,21 @@ void Foam::diameterModels::populationBalanceModel::drift sizeGroups_[i+1].shapeModelPtr()->addDrift(Sue, fp, model); } - if (i != 0) + if (i == sizeGroups().size() - 1) { - const sizeGroup& fw = sizeGroups_[i-1]; + Sp_[i] -= pos(driftRate_())*driftRate_()*fp.phase()/fp.x(); + } + + if (i > 0) + { + const sizeGroup& fw = sizeGroups()[i-1]; volScalarField& Suw = Sui_; + Sp_[i] += 1/(fw.x() - fp.x())*neg(driftRate_())*driftRate_()*fp.phase(); + Suw = - neg(driftRate_())*driftRate_()*rdx_() - *fp*fp.phase()/fp.x() - /(rx_() - 1); + fw.x()/(fp.x()*(fw.x() - fp.x())) + *neg(driftRate_())*driftRate_()*fp*fp.phase(); Su_[i-1] += Suw; @@ -659,6 +606,11 @@ void Foam::diameterModels::populationBalanceModel::drift sizeGroups_[i-1].shapeModelPtr()->addDrift(Suw, fp, model); } + + if (i == 0) + { + Sp_[i] -= neg(driftRate_())*driftRate_()*fp.phase()/fp.x(); + } } @@ -687,7 +639,7 @@ void Foam::diameterModels::populationBalanceModel::sources() { sizeGroups_[i].shapeModelPtr()->reset(); Su_[i] = Zero; - SuSp_[i] = Zero; + Sp_[i] = Zero; } forAllConstIter @@ -883,7 +835,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel v_(), delta_(), Su_(), - SuSp_(), + Sp_(), Sui_ ( IOobject @@ -919,8 +871,6 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel driftModel::iNew(*this) ), driftRate_(), - rx_(), - rdx_(), nucleation_ ( dict_.lookup("nucleationModels"), @@ -1020,36 +970,6 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel dimensionedScalar(dimVolume/dimTime, Zero) ) ); - - rx_.set - ( - new volScalarField - ( - IOobject - ( - "r", - fluid_.time().timeName(), - mesh_ - ), - mesh_, - dimensionedScalar(dimless, Zero) - ) - ); - - rdx_.set - ( - new volScalarField - ( - IOobject - ( - "r", - fluid_.time().timeName(), - mesh_ - ), - mesh_, - dimensionedScalar(dimless, Zero) - ) - ); } if (nucleation_.size() != 0) @@ -1269,7 +1189,7 @@ void Foam::diameterModels::populationBalanceModel::solve() + fvm::div(phase.alphaPhi(), fi) == Su_[i] - - fvm::SuSp(SuSp_[i], fi) + - fvm::Sp(Sp_[i], fi) + fluid_.fvModels().source(alpha, rho, fi)/rho + fvc::ddt(residualAlpha, fi) - fvm::ddt(residualAlpha, fi) diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H index da83e1019e..87bc0f548d 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H @@ -233,8 +233,8 @@ private: //- Explicitly treated sources PtrList Su_; - //- Sources treated implicitly or explicitly depending on sign - PtrList SuSp_; + //- Implicitly treated sources + PtrList Sp_; //- Field for caching sources volScalarField Sui_; @@ -263,12 +263,6 @@ private: //- Drift rate autoPtr driftRate_; - //- Ratio between successive representative volumes - autoPtr rx_; - - //- Ratio between successive class widths - autoPtr rdx_; - //- Zeroeth order models PtrList nucleation_; @@ -416,8 +410,8 @@ public: //- Return the sizeGroups belonging to this populationBalance inline const UPtrList& sizeGroups() const; - //- Return semi-implicit source terms - inline const volScalarField& SuSp(const label i) const; + //- Return implicit source terms + inline const volScalarField& Sp(const label i) const; //- Return list of unordered phasePairs in this populationBalance inline const phasePairTable& phasePairs() const; diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H index e7666b2296..f554a6d1a7 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2017-2020 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2021 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -85,9 +85,9 @@ Foam::diameterModels::populationBalanceModel::sizeGroups() const inline const Foam::volScalarField& -Foam::diameterModels::populationBalanceModel::SuSp(const label i) const +Foam::diameterModels::populationBalanceModel::Sp(const label i) const { - return SuSp_[i]; + return Sp_[i]; } From 3f9435fb18412cfa5314bbc07fac5e0ef8d00b7a Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 8 Jun 2021 08:55:17 +0100 Subject: [PATCH 3/4] multiphaseEulerFoam: populationBalance: Revised algebraic operations Field algebra has been optimised by careful ordering to minimise the number of expensive operations; e.g., changing a/b/c to a/(b*c) in order to minimise the number of divisions. Some minor consistency improvements have also been made throughout population balance. Patch contributed by Institute of Fluid Dynamics, Helmholtz-Zentrum Dresden - Rossendorf (HZDR) --- .../sizeGroup/shapeModels/fractal/fractal.C | 20 +- .../KochFriedlander/KochFriedlander.C | 6 +- .../shapeModels/spherical/spherical.C | 4 +- .../velocityGroup/sizeGroup/sizeGroup.C | 2 +- .../velocityGroup/velocityGroup.C | 2 +- .../LehrMilliesMewes/LehrMilliesMewes.C | 4 +- .../LuoSvendsen/LuoSvendsen.C | 6 +- .../powerLawUniformBinary.C | 4 +- .../LehrMilliesMewesCoalescence.C | 4 +- .../coalescenceModels/Luo/Luo.C | 10 +- .../LaakkonenDaughterSizeDistribution.C | 2 +- .../uniformBinary/uniformBinary.C | 2 +- .../reactionDriven/reactionDriven.C | 4 +- .../wallBoiling/wallBoiling.C | 2 +- .../populationBalanceModel.C | 184 ++++++++---------- .../populationBalanceModel.H | 13 +- .../populationBalanceModelI.H | 7 - 17 files changed, 125 insertions(+), 151 deletions(-) diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C index afd0c95eb8..ff3f8310a4 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C @@ -184,8 +184,8 @@ Foam::diameterModels::shapeModels::fractal::dColl() const volScalarField& dColl = tDColl.ref(); dColl = - 6.0/kappa_ - *pow(sizeGroup_.x()*pow3(kappa_)/(36.0*pi*alphaC_), 1.0/Df_); + 6/kappa_ + *pow(sizeGroup_.x()*pow3(kappa_)/(36*pi*alphaC_), 1/Df_); return tDColl; } @@ -230,8 +230,8 @@ void Foam::diameterModels::shapeModels::fractal::correct() kappa_ = min ( - max(kappa_, 6.0/sizeGroup_.dSph()), - 6.0/popBal.sizeGroups().first().dSph() + max(kappa_, 6/sizeGroup_.dSph()), + 6/popBal.sizeGroups().first().dSph() ); kappa_.correctBoundaryConditions(); @@ -283,7 +283,7 @@ void Foam::diameterModels::shapeModels::fractal::addDrift fu.shapeModelPtr()() ); - volScalarField dp(6.0/sourceKappa); + volScalarField dp(6/sourceKappa); const volScalarField a(sourceKappa*fu.x()); const dimensionedScalar dv(sizeGroup_.x() - fu.x()); @@ -292,18 +292,18 @@ void Foam::diameterModels::shapeModels::fractal::addDrift (2.0/3.0)*dv *( sourceKappa - + sourceShape.Df_*(1.0/sourceShape.d() - 1.0/dp) + + sourceShape.Df_*(1/sourceShape.d() - 1/dp) ) ); - dp += 6.0*(dv*a - fu.x()*da1)/sqr(a); + dp += 6*(dv*a - fu.x()*da1)/sqr(a); - const volScalarField np(6.0*sizeGroup_.x()/pi/pow3(dp)); - const volScalarField dc(dp*pow(np/alphaC_, 1.0/Df_)); + const volScalarField np(6*sizeGroup_.x()/pi/pow3(dp)); + const volScalarField dc(dp*pow(np/alphaC_, 1/Df_)); const volScalarField da2 ( - dv*(4.0/dp + 2.0*Df_/3.0*(1.0/dc - 1.0/dp)) + dv*(4/dp + 2*Df_/3*(1/dc - 1/dp)) ); Su_ += (a + 0.5*da1 + 0.5*da2)/sizeGroup_.x()*Su; diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/sinteringModels/KochFriedlander/KochFriedlander.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/sinteringModels/KochFriedlander/KochFriedlander.C index 302bb7c33e..d8c881bb88 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/sinteringModels/KochFriedlander/KochFriedlander.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/sinteringModels/KochFriedlander/KochFriedlander.C @@ -82,7 +82,7 @@ Foam::diameterModels::shapeModels::sinteringModels::KochFriedlander::tau() const ( "tau", fractal_.SizeGroup().mesh(), - dimensionedScalar(dimTime, 0.0) + dimensionedScalar(dimTime, Zero) ) ); @@ -120,7 +120,7 @@ Foam::diameterModels::shapeModels::sinteringModels::KochFriedlander::R() const fi.mesh() ), fi.mesh(), - dimensionedScalar(inv(dimTime), 0) + dimensionedScalar(inv(dimTime), Zero) ); volScalarField::Internal tau(this->tau()); @@ -130,7 +130,7 @@ Foam::diameterModels::shapeModels::sinteringModels::KochFriedlander::R() const R[celli] = fi[celli]*alpha[celli]/tau[celli]; } - return fvm::Sp(R, kappai) - 6.0/fi.dSph()*R; + return fvm::Sp(R, kappai) - 6/fi.dSph()*R; } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/spherical/spherical.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/spherical/spherical.C index 8dffdcde16..ca9d87ee2a 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/spherical/spherical.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/spherical/spherical.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2019-2020 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2019-2021 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -79,7 +79,7 @@ Foam::diameterModels::shapeModels::spherical::a() const ( "a", sizeGroup_.mesh(), - 6.0/sizeGroup_.dSph()*sizeGroup_.x() + 6/sizeGroup_.dSph()*sizeGroup_.x() ) ); } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.C index bb860f3678..df016b0228 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.C @@ -67,7 +67,7 @@ Foam::diameterModels::sizeGroup::sizeGroup phase_(phase), velocityGroup_(velocityGroup), dSph_("dSph", dimLength, dict), - x_("x", pi/6.0*pow3(dSph_)), + x_("x", pi/6*pow3(dSph_)), value_(dict.lookup("value")) { // Adjust refValue at mixedFvPatchField boundaries diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.C index 23e20f27ae..3eba5bfa3b 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.C @@ -64,7 +64,7 @@ Foam::tmp Foam::diameterModels::velocityGroup::dsm() const invDsm += fi.a()*fi/fi.x(); } - return 6.0/tInvDsm; + return 6/tInvDsm; } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/LehrMilliesMewes/LehrMilliesMewes.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/LehrMilliesMewes/LehrMilliesMewes.C index d3c0abf887..2e5f15b410 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/LehrMilliesMewes/LehrMilliesMewes.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/LehrMilliesMewes/LehrMilliesMewes.C @@ -102,9 +102,9 @@ addToBinaryBreakupRate binaryBreakupRate += 0.5*pow(fj.dSph()/L, 5.0/3.0) *exp(-sqrt(2.0)/pow3(fj.dSph()/L)) - *6.0/pow(pi, 1.5)/pow3(fi.dSph()/L) + *6/pow(pi, 1.5)/pow3(fi.dSph()/L) *exp(-9.0/4.0*sqr(log(pow(2.0, 0.4)*fi.dSph()/L))) - /max(1.0 + erf(1.5*log(pow(2.0, 1.0/15.0)*fj.dSph()/L)), small) + /max(1 + erf(1.5*log(pow(2.0, 1.0/15.0)*fj.dSph()/L)), small) /(T*pow3(L)); } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/LuoSvendsen/LuoSvendsen.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/LuoSvendsen/LuoSvendsen.C index 7f0c6a492c..5e5e2f16be 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/LuoSvendsen/LuoSvendsen.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/LuoSvendsen/LuoSvendsen.C @@ -180,7 +180,7 @@ Foam::diameterModels::binaryBreakupModels::LuoSvendsen::addToBinaryBreakupRate const volScalarField b ( - 12.0*cf*popBal_.sigmaWithContinuousPhase(fi.phase()) + 12*cf*popBal_.sigmaWithContinuousPhase(fi.phase()) /( beta_*continuousPhase.rho()*pow(fj.dSph(), 5.0/3.0) *pow(popBal_.continuousTurbulence().epsilon(), 2.0/3.0) @@ -191,12 +191,12 @@ Foam::diameterModels::binaryBreakupModels::LuoSvendsen::addToBinaryBreakupRate const volScalarField tMin(b/pow(xiMin, 11.0/3.0)); - volScalarField integral(3.0/(11.0*pow(b, 8.0/11.0))); + volScalarField integral(3/(11*pow(b, 8.0/11.0))); forAll(integral, celli) { integral[celli] *= - 2.0*pow(b[celli], 3.0/11.0)*tgamma(5.0/11.0) + 2*pow(b[celli], 3.0/11.0)*tgamma(5.0/11.0) *( gammaUpperReg5by11_->value(b[celli]) - gammaUpperReg5by11_->value(tMin[celli]) diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/powerLawUniformBinary/powerLawUniformBinary.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/powerLawUniformBinary/powerLawUniformBinary.C index 07869d9615..a1ab27e14f 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/powerLawUniformBinary/powerLawUniformBinary.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/binaryBreakupModels/powerLawUniformBinary/powerLawUniformBinary.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2017-2020 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2021 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -74,7 +74,7 @@ addToBinaryBreakupRate const sizeGroup& fj = popBal_.sizeGroups()[j]; binaryBreakupRate.primitiveFieldRef() += - pow(fj.x().value(), power_)*2.0/fj.x().value(); + pow(fj.x().value(), power_)*2/fj.x().value(); } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/LehrMilliesMewesCoalescence/LehrMilliesMewesCoalescence.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/LehrMilliesMewesCoalescence/LehrMilliesMewesCoalescence.C index 076a076fd0..f31057536a 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/LehrMilliesMewesCoalescence/LehrMilliesMewesCoalescence.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/LehrMilliesMewesCoalescence/LehrMilliesMewesCoalescence.C @@ -95,11 +95,11 @@ addToCoalescenceRate ); coalescenceRate += - pi/4.0*sqr(fi.dSph() + fj.dSph())*min(uChar, uCrit_) + pi/4*sqr(fi.dSph() + fj.dSph())*min(uChar, uCrit_) *exp ( - sqr(cbrt(alphaMax_) - /cbrt(max(popBal_.alphas(), fi.phase().residualAlpha())) - 1.0) + /cbrt(max(popBal_.alphas(), fi.phase().residualAlpha())) - 1) ); } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/Luo/Luo.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/Luo/Luo.C index aadb3114d3..0b6067795a 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/Luo/Luo.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/Luo/Luo.C @@ -62,7 +62,7 @@ Luo : coalescenceModel(popBal, dict), beta_(dimensionedScalar::lookupOrDefault("beta", dict, dimless, 2.05)), - C1_(dimensionedScalar::lookupOrDefault("C1", dict, dimless, 1.0)) + C1_(dimensionedScalar::lookupOrDefault("C1", dict, dimless, 1)) {} @@ -102,18 +102,18 @@ addToCoalescenceRate ( sqrt(beta_) *cbrt(popBal_.continuousTurbulence().epsilon()*fi.dSph()) - *sqrt(1.0 + pow(xi, -2.0/3.0)) + *sqrt(1 + pow(xi, -2.0/3.0)) ); coalescenceRate += - pi/4.0*sqr(fi.dSph() + fj.dSph())*uij + pi/4*sqr(fi.dSph() + fj.dSph())*uij *exp ( - C1_ - *sqrt(0.75*(1.0 + sqr(xi))*(1.0 + pow3(xi))) + *sqrt(0.75*(1 + sqr(xi))*(1 + pow3(xi))) /( sqrt(fi.phase().rho()/continuousPhase.rho() - + vm.Cvm())*pow3(1.0 + xi) + + vm.Cvm())*pow3(1 + xi) ) *sqrt ( diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/LaakkonenDaughterSizeDistribution/LaakkonenDaughterSizeDistribution.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/LaakkonenDaughterSizeDistribution/LaakkonenDaughterSizeDistribution.C index a276e2c281..cb13192ace 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/LaakkonenDaughterSizeDistribution/LaakkonenDaughterSizeDistribution.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/LaakkonenDaughterSizeDistribution/LaakkonenDaughterSizeDistribution.C @@ -115,7 +115,7 @@ LaakkonenDaughterSizeDistribution::calcNik if (k == 0) { - return 1.0; + return 1; } return diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.C index ea317d7a45..76dfacd7b0 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.C @@ -85,7 +85,7 @@ Foam::diameterModels::daughterSizeDistributionModels::uniformBinary::calcNik { if (k == 0) { - return 1.0; + return 1; } return (sizeGroups[i+1].x() - xi)/(xk - x0); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/reactionDriven/reactionDriven.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/reactionDriven/reactionDriven.C index 20d439d0a2..db77555dd6 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/reactionDriven/reactionDriven.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/reactionDriven/reactionDriven.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2019-2020 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2019-2021 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -139,7 +139,7 @@ Foam::diameterModels::nucleationModels::reactionDriven::addToNucleationRate velGroup_.phase().name() == pair_.first() ? +1 : -1; nucleationRate += - popBal_.eta(i, pi/6.0*pow3(dNuc_))*dmidtfSign*dmidtf/rho/fi.x(); + popBal_.eta(i, pi/6*pow3(dNuc_))*dmidtfSign*dmidtf/rho/fi.x(); } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C index d0dca1d537..4ec4b5a505 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C @@ -174,7 +174,7 @@ Foam::diameterModels::nucleationModels::wallBoiling::addToNucleationRate const labelList& faceCells = alphatw.patch().faceCells(); - dimensionedScalar unitLength("unitLength", dimLength, 1.0); + dimensionedScalar unitLength("unitLength", dimLength, 1); forAll(alphatw, facei) { diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C index 66751ccc56..e3a575bf66 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C @@ -88,9 +88,9 @@ void Foam::diameterModels::populationBalanceModel::registerSizeGroups { if ( - sizeGroups_.size() != 0 + sizeGroups().size() != 0 && - group.x().value() <= sizeGroups_.last().x().value() + group.x().value() <= sizeGroups().last().x().value() ) { FatalErrorInFunction @@ -99,53 +99,44 @@ void Foam::diameterModels::populationBalanceModel::registerSizeGroups << exit(FatalError); } - sizeGroups_.resize(sizeGroups_.size() + 1); - sizeGroups_.set(sizeGroups_.size() - 1, &group); + sizeGroups_.resize(sizeGroups().size() + 1); + sizeGroups_.set(sizeGroups().size() - 1, &group); - // Grid generation over property space - if (sizeGroups_.size() == 1) + if (sizeGroups().size() == 1) { - // Set the first sizeGroup boundary to the representative value of - // the first sizeGroup v_.append ( new dimensionedScalar ( "v", - sizeGroups_.last().x() + sizeGroups().last().x() ) ); - // Set the last sizeGroup boundary to the representative size of the - // last sizeGroup v_.append ( new dimensionedScalar ( "v", - sizeGroups_.last().x() + sizeGroups().last().x() ) ); } else { - // Overwrite the next-to-last sizeGroup boundary to lie halfway between - // the last two sizeGroups v_.last() = 0.5 *( - sizeGroups_[sizeGroups_.size()-2].x() - + sizeGroups_.last().x() + sizeGroups()[sizeGroups().size()-2].x() + + sizeGroups().last().x() ); - // Set the last sizeGroup boundary to the representative size of the - // last sizeGroup v_.append ( new dimensionedScalar ( "v", - sizeGroups_.last().x() + sizeGroups().last().x() ) ); } @@ -228,21 +219,21 @@ void Foam::diameterModels::populationBalanceModel::precompute() { calcDeltas(); - forAll(coalescence_, model) + forAll(coalescenceModels_, model) { - coalescence_[model].precompute(); + coalescenceModels_[model].precompute(); } - forAll(breakup_, model) + forAll(breakupModels_, model) { - breakup_[model].precompute(); + breakupModels_[model].precompute(); - breakup_[model].dsdPtr()->precompute(); + breakupModels_[model].dsdPtr()->precompute(); } - forAll(binaryBreakup_, model) + forAll(binaryBreakupModels_, model) { - binaryBreakup_[model].precompute(); + binaryBreakupModels_[model].precompute(); } forAll(drift_, model) @@ -264,35 +255,31 @@ birthByCoalescence const label k ) { - const sizeGroup& fj = sizeGroups_[j]; - const sizeGroup& fk = sizeGroups_[k]; + const sizeGroup& fj = sizeGroups()[j]; + const sizeGroup& fk = sizeGroups()[k]; dimensionedScalar Eta; dimensionedScalar v = fj.x() + fk.x(); - for (label i = j; i < sizeGroups_.size(); i++) + for (label i = j; i < sizeGroups().size(); i++) { - // Calculate fraction for intra-interval events Eta = eta(i, v); if (Eta.value() == 0) continue; - const sizeGroup& fi = sizeGroups_[i]; + const sizeGroup& fi = sizeGroups()[i]; - // Avoid double counting of events if (j == k) { Sui_ = - 0.5*fi.x()*coalescenceRate_()*Eta - *fj*fj.phase()/fj.x() - *fk*fk.phase()/fk.x(); + 0.5*fi.x()/(fj.x()*fk.x())*Eta + *coalescenceRate_()*fj*fj.phase()*fk*fk.phase(); } else { Sui_ = - fi.x()*coalescenceRate_()*Eta - *fj*fj.phase()/fj.x() - *fk*fk.phase()/fk.x(); + fi.x()/(fj.x()*fk.x())*Eta + *coalescenceRate_()*fj*fj.phase()*fk*fk.phase(); } Su_[i] += Sui_; @@ -310,7 +297,7 @@ birthByCoalescence Pair::compare(pDmdt_.find(pairij).key(), pairij) ); - *pDmdt_[pairij] += dmdtSign*fj.x()/v*Sui_*fi.phase().rho(); + *pDmdt_[pairij] += dmdtSign*fj.x()/v*Sui_*fj.phase().rho(); } const phasePairKey pairik @@ -326,7 +313,7 @@ birthByCoalescence Pair::compare(pDmdt_.find(pairik).key(), pairik) ); - *pDmdt_[pairik] += dmdtSign*fk.x()/v*Sui_*fi.phase().rho(); + *pDmdt_[pairik] += dmdtSign*fk.x()/v*Sui_*fk.phase().rho(); } sizeGroups_[i].shapeModelPtr()->addCoalescence(Sui_, fj, fk); @@ -341,8 +328,8 @@ deathByCoalescence const label j ) { - const sizeGroup& fi = sizeGroups_[i]; - const sizeGroup& fj = sizeGroups_[j]; + const sizeGroup& fi = sizeGroups()[i]; + const sizeGroup& fj = sizeGroups()[j]; Sp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x(); @@ -360,15 +347,15 @@ birthByBreakup const label model ) { - const sizeGroup& fk = sizeGroups_[k]; + const sizeGroup& fk = sizeGroups()[k]; for (label i = 0; i <= k; i++) { - const sizeGroup& fi = sizeGroups_[i]; + const sizeGroup& fi = sizeGroups()[i]; Sui_ = - fi.x()*breakupRate_()*breakup_[model].dsdPtr()().nik(i, k) - *fk*fk.phase()/fk.x(); + fi.x()*breakupModels_[model].dsdPtr()().nik(i, k)/fk.x() + *breakupRate_()*fk*fk.phase(); Su_[i] += Sui_; @@ -385,7 +372,7 @@ birthByBreakup Pair::compare(pDmdt_.find(pair).key(), pair) ); - *pDmdt_[pair] += dmdtSign*Sui_*fi.phase().rho(); + *pDmdt_[pair] += dmdtSign*Sui_*fk.phase().rho(); } sizeGroups_[i].shapeModelPtr()->addBreakup(Sui_, fk); @@ -395,20 +382,20 @@ birthByBreakup void Foam::diameterModels::populationBalanceModel::deathByBreakup(const label i) { - const sizeGroup& fi = sizeGroups_[i]; - - Sp_[i] += breakupRate_()*fi.phase(); + Sp_[i] += breakupRate_()*sizeGroups()[i].phase(); } void Foam::diameterModels::populationBalanceModel::calcDeltas() { - forAll(sizeGroups_, i) + forAll(sizeGroups(), i) { if (delta_[i].empty()) { - for (label j = 0; j <= sizeGroups_.size() - 1; j++) + for (label j = 0; j <= sizeGroups().size() - 1; j++) { + const sizeGroup& fj = sizeGroups()[j]; + delta_[i].append ( new dimensionedScalar @@ -421,14 +408,14 @@ void Foam::diameterModels::populationBalanceModel::calcDeltas() if ( - v_[i].value() < 0.5*sizeGroups_[j].x().value() + v_[i].value() < 0.5*fj.x().value() && - 0.5*sizeGroups_[j].x().value() < v_[i+1].value() + 0.5*fj.x().value() < v_[i+1].value() ) { - delta_[i][j] = mag(0.5*sizeGroups_[j].x() - v_[i]); + delta_[i][j] = mag(0.5*fj.x() - v_[i]); } - else if (0.5*sizeGroups_[j].x().value() <= v_[i].value()) + else if (0.5*fj.x().value() <= v_[i].value()) { delta_[i][j].value() = 0; } @@ -445,10 +432,12 @@ birthByBinaryBreakup const label j ) { - const sizeGroup& fj = sizeGroups_[j]; - const sizeGroup& fi = sizeGroups_[i]; + const sizeGroup& fi = sizeGroups()[i]; + const sizeGroup& fj = sizeGroups()[j]; - Sui_ = fi.x()*binaryBreakupRate_()*delta_[i][j]*fj*fj.phase()/fj.x(); + const volScalarField Su(binaryBreakupRate_()*fj*fj.phase()); + + Sui_ = fi.x()*delta_[i][j]/fj.x()*Su; Su_[i] += Sui_; @@ -467,7 +456,7 @@ birthByBinaryBreakup Pair::compare(pDmdt_.find(pairij).key(), pairij) ); - *pDmdt_[pairij] += dmdtSign*Sui_*fi.phase().rho(); + *pDmdt_[pairij] += dmdtSign*Sui_*fj.phase().rho(); } dimensionedScalar Eta; @@ -475,18 +464,15 @@ birthByBinaryBreakup for (label k = 0; k <= j; k++) { - // Calculate fraction for intra-interval events Eta = eta(k, v); if (Eta.value() == 0) continue; - const sizeGroup& fk = sizeGroups_[k]; + const sizeGroup& fk = sizeGroups()[k]; volScalarField& Suk = Sui_; - Suk = - sizeGroups_[k].x()*binaryBreakupRate_()*delta_[i][j]*Eta - *fj*fj.phase()/fj.x(); + Suk = fk.x()*delta_[i][j]*Eta/fj.x()*Su; Su_[k] += Suk; @@ -507,10 +493,10 @@ birthByBinaryBreakup ) ); - *pDmdt_[pairkj] += dmdtSign*Suk*fi.phase().rho(); + *pDmdt_[pairkj] += dmdtSign*Suk*fj.phase().rho(); } - sizeGroups_[i].shapeModelPtr()->addBreakup(Suk, fj); + sizeGroups_[k].shapeModelPtr()->addBreakup(Suk, fj); } } @@ -522,9 +508,7 @@ deathByBinaryBreakup const label i ) { - const volScalarField& alphai = sizeGroups_[i].phase(); - - Sp_[i] += alphai*binaryBreakupRate_()*delta_[j][i]; + Sp_[i] += sizeGroups()[i].phase()*binaryBreakupRate_()*delta_[j][i]; } @@ -621,11 +605,11 @@ nucleation nucleationModel& model ) { - const sizeGroup& fi = sizeGroups_[i]; + const sizeGroup& fi = sizeGroups()[i]; model.addToNucleationRate(nucleationRate_(), i); - Sui_ = sizeGroups_[i].x()*nucleationRate_(); + Sui_ = fi.x()*nucleationRate_(); Su_[i] += Sui_; @@ -635,7 +619,7 @@ nucleation void Foam::diameterModels::populationBalanceModel::sources() { - forAll(sizeGroups_, i) + forAll(sizeGroups(), i) { sizeGroups_[i].shapeModelPtr()->reset(); Su_[i] = Zero; @@ -654,15 +638,15 @@ void Foam::diameterModels::populationBalanceModel::sources() forAll(sizeGroups_, i) { - if (coalescence_.size() != 0) + if (coalescenceModels_.size() != 0) { for (label j = 0; j <= i; j++) { coalescenceRate_() = Zero; - forAll(coalescence_, model) + forAll(coalescenceModels_, model) { - coalescence_[model].addToCoalescenceRate + coalescenceModels_[model].addToCoalescenceRate ( coalescenceRate_(), i, @@ -676,11 +660,11 @@ void Foam::diameterModels::populationBalanceModel::sources() } } - if (breakup_.size() != 0) + if (breakupModels_.size() != 0) { - forAll(breakup_, model) + forAll(breakupModels_, model) { - breakup_[model].setBreakupRate(breakupRate_(), i); + breakupModels_[model].setBreakupRate(breakupRate_(), i); birthByBreakup(i, model); @@ -688,7 +672,7 @@ void Foam::diameterModels::populationBalanceModel::sources() } } - if (binaryBreakup_.size() != 0) + if (binaryBreakupModels_.size() != 0) { label j = 0; @@ -696,9 +680,9 @@ void Foam::diameterModels::populationBalanceModel::sources() { binaryBreakupRate_() = Zero; - forAll(binaryBreakup_, model) + forAll(binaryBreakupModels_, model) { - binaryBreakup_[model].addToBinaryBreakupRate + binaryBreakupModels_[model].addToBinaryBreakupRate ( binaryBreakupRate_(), j, @@ -770,7 +754,7 @@ Foam::diameterModels::populationBalanceModel::calcDsm() invDsm += max(phase, phase.residualAlpha())/(phase.d()*alphas_()); } - return 1.0/tInvDsm; + return 1/tInvDsm; } @@ -847,19 +831,19 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel mesh_, dimensionedScalar(inv(dimTime), Zero) ), - coalescence_ + coalescenceModels_ ( dict_.lookup("coalescenceModels"), coalescenceModel::iNew(*this) ), coalescenceRate_(), - breakup_ + breakupModels_ ( dict_.lookup("breakupModels"), breakupModel::iNew(*this) ), breakupRate_(), - binaryBreakup_ + binaryBreakupModels_ ( dict_.lookup("binaryBreakupModels"), binaryBreakupModel::iNew(*this) @@ -886,7 +870,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel this->createPhasePairs(); - if (sizeGroups_.size() < 3) + if (sizeGroups().size() < 3) { FatalErrorInFunction << "The populationBalance " << name_ @@ -895,7 +879,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel } - if (coalescence_.size() != 0) + if (coalescenceModels_.size() != 0) { coalescenceRate_.set ( @@ -913,7 +897,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel ); } - if (breakup_.size() != 0) + if (breakupModels_.size() != 0) { breakupRate_.set ( @@ -931,7 +915,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel ); } - if (binaryBreakup_.size() != 0) + if (binaryBreakupModels_.size() != 0) { binaryBreakupRate_.set ( @@ -1078,17 +1062,17 @@ Foam::diameterModels::populationBalanceModel::eta const dimensionedScalar& v ) const { - const dimensionedScalar& x0 = sizeGroups_[0].x(); - const dimensionedScalar& xi = sizeGroups_[i].x(); - const dimensionedScalar& xm = sizeGroups_.last().x(); + const dimensionedScalar& x0 = sizeGroups()[0].x(); + const dimensionedScalar& xi = sizeGroups()[i].x(); + const dimensionedScalar& xm = sizeGroups().last().x(); dimensionedScalar lowerBoundary(x0); dimensionedScalar upperBoundary(xm); - if (i != 0) lowerBoundary = sizeGroups_[i-1].x(); + if (i != 0) lowerBoundary = sizeGroups()[i-1].x(); - if (i != sizeGroups_.size() - 1) upperBoundary = sizeGroups_[i+1].x(); + if (i != sizeGroups().size() - 1) upperBoundary = sizeGroups()[i+1].x(); - if ((i == 0 && v < x0) || (i == sizeGroups_.size() - 1 && v > xm)) + if ((i == 0 && v < x0) || (i == sizeGroups().size() - 1 && v > xm)) { return v/xi; } @@ -1098,7 +1082,7 @@ Foam::diameterModels::populationBalanceModel::eta } else if (v.value() == xi.value()) { - return 1.0; + return 1; } else if (v > xi) { @@ -1175,7 +1159,7 @@ void Foam::diameterModels::populationBalanceModel::solve() maxInitialResidual = 0; - forAll(sizeGroups_, i) + forAll(sizeGroups(), i) { sizeGroup& fi = sizeGroups_[i]; const phaseModel& phase = fi.phase(); @@ -1210,12 +1194,12 @@ void Foam::diameterModels::populationBalanceModel::solve() volScalarField fAlpha0 ( - sizeGroups_.first()*sizeGroups_.first().phase() + sizeGroups().first()*sizeGroups().first().phase() ); volScalarField fAlphaN ( - sizeGroups_.last()*sizeGroups_.last().phase() + sizeGroups().last()*sizeGroups().last().phase() ); Info<< this->name() << " sizeGroup phase fraction first, last = " diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H index 87bc0f548d..1f7b4e34fa 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H @@ -53,7 +53,7 @@ Description \verbatim Coalescence and breakup term formulation: Kumar, S., & Ramkrishna, D. (1996). - On the solution of population balance equations by discretisation-I. A + On the solution of population balance equations by discretization-I. A fixed pivot technique. Chemical Engineering Science, 51(8), 1311-1332. \endverbatim @@ -240,19 +240,19 @@ private: volScalarField Sui_; //- Coalescence models - PtrList coalescence_; + PtrList coalescenceModels_; //- Coalescence rate autoPtr coalescenceRate_; - //- BreakupModels - PtrList breakup_; + //- Breakup models + PtrList breakupModels_; //- Breakup rate autoPtr breakupRate_; //- Binary breakup models - PtrList binaryBreakup_; + PtrList binaryBreakupModels_; //- Binary breakup rate autoPtr binaryBreakupRate_; @@ -420,9 +420,6 @@ public: // belong to this populationBalance inline bool isVelocityGroupPair(const phasePair& pair) const; - //- Return the sizeGroup boundaries - inline const PtrList& v() const; - //- Return total void of phases belonging to this populationBalance inline const volScalarField& alphas() const; diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H index f554a6d1a7..24bdc9a72a 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H @@ -123,13 +123,6 @@ inline bool Foam::diameterModels::populationBalanceModel::isVelocityGroupPair } -inline const Foam::PtrList& -Foam::diameterModels::populationBalanceModel::v() const -{ - return v_; -} - - inline const Foam::volScalarField& Foam::diameterModels::populationBalanceModel::alphas() const { From bc844c383c96de737bfe0c44bc982b2958c9de36 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 8 Jun 2021 09:03:33 +0100 Subject: [PATCH 4/4] multiphaseEulerFoam: populationBalance: Store coalescence and breakup pair indexing Patch contributed by Institute of Fluid Dynamics, Helmholtz-Zentrum Dresden - Rossendorf (HZDR) --- .../populationBalanceModel.C | 156 +++++++++--------- .../populationBalanceModel.H | 25 ++- .../populationBalanceModelI.H | 64 ++++--- 3 files changed, 138 insertions(+), 107 deletions(-) diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C index e3a575bf66..9a20b1652f 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C @@ -217,8 +217,6 @@ void Foam::diameterModels::populationBalanceModel::createPhasePairs() void Foam::diameterModels::populationBalanceModel::precompute() { - calcDeltas(); - forAll(coalescenceModels_, model) { coalescenceModels_[model].precompute(); @@ -248,8 +246,7 @@ void Foam::diameterModels::populationBalanceModel::precompute() } -void Foam::diameterModels::populationBalanceModel:: -birthByCoalescence +void Foam::diameterModels::populationBalanceModel::birthByCoalescence ( const label j, const label k @@ -321,8 +318,7 @@ birthByCoalescence } -void Foam::diameterModels::populationBalanceModel:: -deathByCoalescence +void Foam::diameterModels::populationBalanceModel::deathByCoalescence ( const label i, const label j @@ -340,8 +336,7 @@ deathByCoalescence } -void Foam::diameterModels::populationBalanceModel:: -birthByBreakup +void Foam::diameterModels::populationBalanceModel::birthByBreakup ( const label k, const label model @@ -425,8 +420,7 @@ void Foam::diameterModels::populationBalanceModel::calcDeltas() } -void Foam::diameterModels::populationBalanceModel:: -birthByBinaryBreakup +void Foam::diameterModels::populationBalanceModel::birthByBinaryBreakup ( const label i, const label j @@ -501,8 +495,7 @@ birthByBinaryBreakup } -void Foam::diameterModels::populationBalanceModel:: -deathByBinaryBreakup +void Foam::diameterModels::populationBalanceModel::deathByBinaryBreakup ( const label j, const label i @@ -598,8 +591,7 @@ void Foam::diameterModels::populationBalanceModel::drift } -void Foam::diameterModels::populationBalanceModel:: -nucleation +void Foam::diameterModels::populationBalanceModel::nucleation ( const label i, nucleationModel& model @@ -636,84 +628,73 @@ void Foam::diameterModels::populationBalanceModel::sources() *pDmdt_(phasePairIter()) = Zero; } - forAll(sizeGroups_, i) + forAll(coalescencePairs_, coalescencePairi) { - if (coalescenceModels_.size() != 0) + label i = coalescencePairs_[coalescencePairi].first(); + label j = coalescencePairs_[coalescencePairi].second(); + + coalescenceRate_() = Zero; + + forAll(coalescenceModels_, model) { - for (label j = 0; j <= i; j++) - { - coalescenceRate_() = Zero; - - forAll(coalescenceModels_, model) - { - coalescenceModels_[model].addToCoalescenceRate - ( - coalescenceRate_(), - i, - j - ); - } - - birthByCoalescence(i, j); - - deathByCoalescence(i, j); - } + coalescenceModels_[model].addToCoalescenceRate + ( + coalescenceRate_(), + i, + j + ); } - if (breakupModels_.size() != 0) + birthByCoalescence(i, j); + + deathByCoalescence(i, j); + } + + forAll(binaryBreakupPairs_, binaryBreakupPairi) + { + label i = binaryBreakupPairs_[binaryBreakupPairi].first(); + label j = binaryBreakupPairs_[binaryBreakupPairi].second(); + + binaryBreakupRate_() = Zero; + + forAll(binaryBreakupModels_, model) { - forAll(breakupModels_, model) - { - breakupModels_[model].setBreakupRate(breakupRate_(), i); - - birthByBreakup(i, model); - - deathByBreakup(i); - } + binaryBreakupModels_[model].addToBinaryBreakupRate + ( + binaryBreakupRate_(), + j, + i + ); } - if (binaryBreakupModels_.size() != 0) + birthByBinaryBreakup(j, i); + + deathByBinaryBreakup(j, i); + } + + forAll(sizeGroups(), i) + { + forAll(breakupModels_, model) { - label j = 0; + breakupModels_[model].setBreakupRate(breakupRate_(), i); - while (delta_[j][i].value() != 0) - { - binaryBreakupRate_() = Zero; + birthByBreakup(i, model); - forAll(binaryBreakupModels_, model) - { - binaryBreakupModels_[model].addToBinaryBreakupRate - ( - binaryBreakupRate_(), - j, - i - ); - } - - birthByBinaryBreakup(j, i); - - deathByBinaryBreakup(j, i); - - j++; - } + deathByBreakup(i); } - if (drift_.size() != 0) + forAll(drift_, model) { - forAll(drift_, model) - { - driftRate_() = Zero; - drift(i, drift_[model]); - } + driftRate_() = Zero; + + drift(i, drift_[model]); } - if (nucleation_.size() != 0) + forAll(nucleation_, model) { - forAll(nucleation_, model) - { - nucleationRate_() = Zero; - nucleation(i, nucleation_[model]); - } + nucleationRate_() = Zero; + + nucleation(i, nucleation_[model]); } } } @@ -837,6 +818,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel coalescenceModel::iNew(*this) ), coalescenceRate_(), + coalescencePairs_(), breakupModels_ ( dict_.lookup("breakupModels"), @@ -849,6 +831,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel binaryBreakupModel::iNew(*this) ), binaryBreakupRate_(), + binaryBreakupPairs_(), drift_ ( dict_.lookup("driftModels"), @@ -895,6 +878,14 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel dimensionedScalar(dimVolume/dimTime, Zero) ) ); + + forAll(sizeGroups(), i) + { + for (label j = 0; j <= i; j++) + { + coalescencePairs_.append(labelPair(i, j)); + } + } } if (breakupModels_.size() != 0) @@ -936,6 +927,19 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel ) ) ); + + calcDeltas(); + + forAll(sizeGroups(), i) + { + label j = 0; + + while (delta_[j][i].value() != 0) + { + binaryBreakupPairs_.append(labelPair(i, j)); + j++; + } + } } if (drift_.size() != 0) diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H index 1f7b4e34fa..f7588a785c 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H @@ -154,6 +154,7 @@ SourceFiles #include "pimpleControl.H" #include "phaseDynamicMomentumTransportModel.H" #include "HashPtrTable.H" +#include "Pair.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -245,6 +246,9 @@ private: //- Coalescence rate autoPtr coalescenceRate_; + //- Coalescence relevant size group pairs + List coalescencePairs_; + //- Breakup models PtrList breakupModels_; @@ -257,6 +261,9 @@ private: //- Binary breakup rate autoPtr binaryBreakupRate_; + //- Binary breakup relevant size group pairs + List binaryBreakupPairs_; + //- Drift models PtrList drift_; @@ -410,15 +417,17 @@ public: //- Return the sizeGroups belonging to this populationBalance inline const UPtrList& sizeGroups() const; - //- Return implicit source terms - inline const volScalarField& Sp(const label i) const; - //- Return list of unordered phasePairs in this populationBalance inline const phasePairTable& phasePairs() const; - //- Returns true if both phases are velocity groups and - // belong to this populationBalance - inline bool isVelocityGroupPair(const phasePair& pair) const; + //- Return implicit source terms + inline const volScalarField& Sp(const label i) const; + + //- Return coalescence relevant size group pairs + inline const List& coalescencePairs() const; + + //- Return binary breakup relevant size group pairs + inline const List& binaryBreakupPairs() const; //- Return total void of phases belonging to this populationBalance inline const volScalarField& alphas() const; @@ -426,6 +435,10 @@ public: //- Return average velocity inline const volVectorField& U() const; + //- Returns true if both phases are velocity groups and + // belong to this populationBalance + inline bool isVelocityGroupPair(const phasePair& pair) const; + //- Return allocation coefficient const dimensionedScalar eta ( diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H index 24bdc9a72a..232321bcba 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H @@ -84,13 +84,6 @@ Foam::diameterModels::populationBalanceModel::sizeGroups() const } -inline const Foam::volScalarField& -Foam::diameterModels::populationBalanceModel::Sp(const label i) const -{ - return Sp_[i]; -} - - inline const Foam::diameterModels::populationBalanceModel::phasePairTable& Foam::diameterModels::populationBalanceModel::phasePairs() const { @@ -98,28 +91,24 @@ Foam::diameterModels::populationBalanceModel::phasePairs() const } -inline bool Foam::diameterModels::populationBalanceModel::isVelocityGroupPair -( - const phasePair& pair -) const +inline const Foam::volScalarField& +Foam::diameterModels::populationBalanceModel::Sp(const label i) const { - if - ( - isA(pair.phase1().dPtr()()) - && - isA(pair.phase2().dPtr()()) - ) - { - const velocityGroup& velGroup1 = - refCast(pair.phase1().dPtr()()); + return Sp_[i]; +} - const velocityGroup& velGroup2 = - refCast(pair.phase2().dPtr()()); - return velGroup1.popBalName() == velGroup2.popBalName(); - } +inline const Foam::List>& +Foam::diameterModels::populationBalanceModel::coalescencePairs() const +{ + return coalescencePairs_; +} - return false; + +inline const Foam::List>& +Foam::diameterModels::populationBalanceModel::binaryBreakupPairs() const +{ + return binaryBreakupPairs_; } @@ -151,4 +140,29 @@ Foam::diameterModels::populationBalanceModel::U() const } +inline bool Foam::diameterModels::populationBalanceModel::isVelocityGroupPair +( + const phasePair& pair +) const +{ + if + ( + isA(pair.phase1().dPtr()()) + && + isA(pair.phase2().dPtr()()) + ) + { + const velocityGroup& velGroup1 = + refCast(pair.phase1().dPtr()()); + + const velocityGroup& velGroup2 = + refCast(pair.phase2().dPtr()()); + + return velGroup1.popBalName() == velGroup2.popBalName(); + } + + return false; +} + + // ************************************************************************* //