From ee6216b4e75675f4e77e7a8704c8db3b7f76cb70 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Mon, 23 Apr 2018 09:02:14 +0100 Subject: [PATCH] reactingEulerFoam: Improved Sauter mean diameter calculation The Sauter mean diameter calculation has been modified to be more stable in the limit of vanishing phase fraction. The calculation of the overall Sauter mean diameter for a populationBalance involving more than one velocityGroup has been removed. This calculation depends upon the phase fraction and it is not stable as the fractions tend to zero. The overall Sauter mean diameter is only used for post-processing and can still be recovered from the individual diameter fields of the involved velocityGroups. Some parts of the population balance modeling have also been renamed and refactored. Patch contributed by Institute of Fluid Dynamics, Helmholtz-Zentrum Dresden - Rossendorf (HZDR) --- .../reactingEulerFoam/phaseSystems/Make/files | 2 +- .../velocityGroup/velocityGroup.C | 137 ++-------- .../velocityGroup/velocityGroup.H | 16 -- .../velocityGroup/velocityGroupI.H | 16 +- .../CoulaloglouTavlaridesCoalescence.H | 4 +- .../daughterSizeDistributionModel.C | 4 +- .../daughterSizeDistributionModel.H | 8 +- .../uniformBinary.C} | 76 ++---- .../uniformBinary.H} | 20 +- .../populationBalanceModel.C | 240 ++++++------------ .../populationBalanceModel.H | 7 +- 11 files changed, 134 insertions(+), 396 deletions(-) rename applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/{uniformBinaryDsd/uniformBinaryDsd.C => uniformBinary/uniformBinary.C} (60%) rename applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/{uniformBinaryDsd/uniformBinaryDsd.H => uniformBinary/uniformBinary.H} (87%) diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/Make/files b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/Make/files index 00bf5cbb4..2e4652fdd 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/Make/files +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/Make/files @@ -30,7 +30,7 @@ populationBalanceModel/breakupModels/exponential/exponential.C populationBalanceModel/breakupModels/powerLaw/powerLaw.C populationBalanceModel/daughterSizeDistributionModels/daughterSizeDistributionModel/daughterSizeDistributionModel.C -populationBalanceModel/daughterSizeDistributionModels/uniformBinaryDsd/uniformBinaryDsd.C +populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.C populationBalanceModel/driftModels/driftModel/driftModel.C populationBalanceModel/driftModels/constantDrift/constantDrift.C diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.C index 6ac0b4e16..2b42100f8 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.C @@ -49,74 +49,33 @@ namespace diameterModels // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // -Foam::tmp -Foam::diameterModels::velocityGroup::secondMoment() const -{ - tmp tm2 - ( - new volScalarField - ( - IOobject - ( - "m2", - phase_.time().timeName(), - phase_.mesh() - ), - phase_.mesh(), - dimensionedScalar("m2", inv(dimLength), Zero) - ) - ); - - volScalarField& m2 = tm2.ref(); - - forAll(sizeGroups_, i) - { - const sizeGroup& fi = sizeGroups_[i]; - - m2 += sqr(fi.d())*formFactor()*fi - *max(fi.phase(), small)/fi.x(); - } - - return tm2; -} - - -Foam::tmp -Foam::diameterModels::velocityGroup::thirdMoment() const -{ - tmp tm3 - ( - new volScalarField - ( - IOobject - ( - "m3", - phase_.time().timeName(), - phase_.mesh() - ), - phase_.mesh(), - dimensionedScalar("m3", dimless, Zero) - ) - ); - - volScalarField& m3 = tm3.ref(); - - forAll(sizeGroups_, i) - { - const sizeGroup& fi = sizeGroups_[i]; - - m3 += pow3(fi.d())*formFactor()*fi - *max(fi.phase(), small)/fi.x(); - } - - return tm3; -} - - Foam::tmp Foam::diameterModels::velocityGroup::dsm() const { - return - max(min(phase_/m2_, sizeGroups_.last().d()), sizeGroups_.first().d()); + tmp tInvDsm + ( + new volScalarField + ( + IOobject + ( + "invDsm", + phase_.time().timeName(), + phase_.mesh() + ), + phase_.mesh(), + dimensionedScalar("invDsm", inv(dimLength), Zero) + ) + ); + + volScalarField& invDsm = tInvDsm.ref(); + + forAll(sizeGroups_, i) + { + const sizeGroup& fi = sizeGroups_[i]; + + invDsm += fi/fi.d(); + } + + return 1.0/tInvDsm; } @@ -243,44 +202,6 @@ Foam::diameterModels::velocityGroup::velocityGroup ), fSum() ), - m2_ - ( - IOobject - ( - IOobject::groupName - ( - "m2", - IOobject::groupName - ( - phase.name(), - popBalName_ - ) - ), - phase.time().timeName(), - phase.mesh() - ), - phase.mesh(), - dimensionedScalar("m2", inv(dimLength), Zero) - ), - m3_ - ( - IOobject - ( - IOobject::groupName - ( - "m3", - IOobject::groupName - ( - phase.name(), - popBalName_ - ) - ), - phase.time().timeName(), - phase.mesh() - ), - phase.mesh(), - dimensionedScalar("m3", dimless, Zero) - ), d_ ( IOobject @@ -354,10 +275,6 @@ Foam::diameterModels::velocityGroup::velocityGroup fields_.add(sizeGroups_[i]); } - m2_ = secondMoment(); - - m3_ = thirdMoment(); - d_ = dsm(); } @@ -379,10 +296,6 @@ void Foam::diameterModels::velocityGroup::preSolve() void Foam::diameterModels::velocityGroup::postSolve() { - m2_ = secondMoment(); - - m3_ = thirdMoment(); - d_ = dsm(); Info<< this->phase().name() << " Sauter mean diameter, min, max = " diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.H index 88e3c31f8..c9cbcbcc0 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.H @@ -112,12 +112,6 @@ class velocityGroup //- Sum of sizeGroup volume fractions volScalarField fSum_; - //- Second moment of the distribution - volScalarField m2_; - - //- Third moment of the distribution - volScalarField m3_; - //- Number-based Sauter-mean diameter of the phase volScalarField d_; @@ -133,10 +127,6 @@ class velocityGroup // Private member functions - tmp secondMoment() const; - - tmp thirdMoment() const; - tmp dsm() const; tmp fSum() const; @@ -180,12 +170,6 @@ public: //- Return sizeGroups belonging to this velocityGroup inline const PtrList& sizeGroups() const; - //- Return second moment of the distribution - inline const volScalarField& m2() const; - - //- Return third moment of the distribution - inline const volScalarField& m3() const; - //- Return const-reference to multivariate convectionScheme inline const tmp>& mvConvection() const; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroupI.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroupI.H index ccf918ad7..3a930c917 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroupI.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroupI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -53,20 +53,6 @@ Foam::diameterModels::velocityGroup::sizeGroups() const } -inline const Foam::volScalarField& -Foam::diameterModels::velocityGroup::m2() const -{ - return m2_; -} - - -inline const Foam::volScalarField& -Foam::diameterModels::velocityGroup::m3() const -{ - return m3_; -} - - inline const Foam::tmp>& Foam::diameterModels::velocityGroup::mvConvection() const { diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/CoulaloglouTavlaridesCoalescence/CoulaloglouTavlaridesCoalescence.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/CoulaloglouTavlaridesCoalescence/CoulaloglouTavlaridesCoalescence.H index 9e7a24668..59bb46e86 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/CoulaloglouTavlaridesCoalescence/CoulaloglouTavlaridesCoalescence.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/CoulaloglouTavlaridesCoalescence/CoulaloglouTavlaridesCoalescence.H @@ -45,8 +45,8 @@ Description \vartable \sigma | Surface tension [N/m] - v_i | Volume of bubble i [m3] - v_j | Volume of bubble j [m3] + v_i | Volume of droplet i [m3] + v_j | Volume of droplet j [m3] \epsilon_c | Turbulent dissipation rate of continuous phase [m2/s3] \alpha_d | Total void fraction of disperse phase [-] \mu_c | Molecular dynamic viscosity of liquid phase [Pa s] diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/daughterSizeDistributionModel/daughterSizeDistributionModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/daughterSizeDistributionModel/daughterSizeDistributionModel.C index 2641d388e..4370e474e 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/daughterSizeDistributionModel/daughterSizeDistributionModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/daughterSizeDistributionModel/daughterSizeDistributionModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -114,7 +114,7 @@ void Foam::diameterModels::daughterSizeDistributionModel::correct() for (label i = 0; i <= k; i++) { - nik_[k].append(new dimensionedScalar (this->n(i, k))); + nik_[k].append(new dimensionedScalar (this->calcNik(i, k))); } } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/daughterSizeDistributionModel/daughterSizeDistributionModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/daughterSizeDistributionModel/daughterSizeDistributionModel.H index 8bea3f088..25c91caec 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/daughterSizeDistributionModel/daughterSizeDistributionModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/daughterSizeDistributionModel/daughterSizeDistributionModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -118,7 +118,11 @@ public: //- Calculate and return total number of particles assigned to class i // when a particle of class k breaks - virtual dimensionedScalar n(const label i, const label k) const = 0; + virtual dimensionedScalar calcNik + ( + const label i, + const label k + ) const = 0; }; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinaryDsd/uniformBinaryDsd.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.C similarity index 60% rename from applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinaryDsd/uniformBinaryDsd.C rename to applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.C index d74ec2b44..3862b4ee2 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinaryDsd/uniformBinaryDsd.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -23,7 +23,7 @@ License \*---------------------------------------------------------------------------*/ -#include "uniformBinaryDsd.H" +#include "uniformBinary.H" #include "addToRunTimeSelectionTable.H" #include "breakupModel.H" @@ -35,11 +35,11 @@ namespace diameterModels { namespace daughterSizeDistributionModels { - defineTypeNameAndDebug(uniformBinaryDsd, 0); + defineTypeNameAndDebug(uniformBinary, 0); addToRunTimeSelectionTable ( daughterSizeDistributionModel, - uniformBinaryDsd, + uniformBinary, dictionary ); } @@ -49,8 +49,8 @@ namespace daughterSizeDistributionModels // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::diameterModels::daughterSizeDistributionModels::uniformBinaryDsd:: -uniformBinaryDsd +Foam::diameterModels::daughterSizeDistributionModels::uniformBinary:: +uniformBinary ( const breakupModel& breakup, const dictionary& dict @@ -62,81 +62,35 @@ uniformBinaryDsd // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // -Foam::diameterModels::daughterSizeDistributionModels::uniformBinaryDsd:: -~uniformBinaryDsd() +Foam::diameterModels::daughterSizeDistributionModels::uniformBinary:: +~uniformBinary() {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // Foam::dimensionedScalar -Foam::diameterModels::daughterSizeDistributionModels::uniformBinaryDsd::n +Foam::diameterModels::daughterSizeDistributionModels::uniformBinary::calcNik ( const label i, const label k ) const { - const sizeGroup& fi = *breakup_.popBal().sizeGroups()[i]; - const sizeGroup& fk = *breakup_.popBal().sizeGroups()[k]; - - const dimensionedScalar& xi = fi.x(); - const dimensionedScalar& xk = fk.x(); - dimensionedScalar x("x", dimVolume, Zero); - dimensionedScalar xii("xii", dimVolume, Zero); + const dimensionedScalar& xi = breakup_.popBal().sizeGroups()[i]->x(); + const dimensionedScalar& xk = breakup_.popBal().sizeGroups()[k]->x(); + const List& sizeGroups = breakup_.popBal().sizeGroups(); if (i == 0) { - x = xi; - } - else - { - x = breakup_.popBal().sizeGroups()[i-1]->x(); - } - - if (i == breakup_.popBal().sizeGroups().size() - 1) - { - xii = fi.x(); - } - else - { - xii = breakup_.popBal().sizeGroups()[i+1]->x(); - } - - if (i == 0) - { - return - 1.0/(xii - xi) - *( - (xii*xii - 0.5*sqr(xii)) - - (xi*xii - 0.5*sqr(xi)) - ) - *2.0/xk; + return (sizeGroups[i+1]->x() - xi)/xk; } else if (i == k) { - return - 1.0/(xi - x) - *( - (0.5*sqr(xi) - xi*x) - - (0.5*sqr(x) - x*x) - ) - *2.0/xk; + return (xi - sizeGroups[i-1]->x())/xk; } else { - return - 1.0/(xii - xi) - *( - (xii*xii - 0.5*sqr(xii)) - - (xi*xii - 0.5*sqr(xi)) - ) - *2.0/xk - + 1.0/(xi - x) - *( - (0.5*sqr(xi) - xi*x) - - (0.5*sqr(x) - x*x) - ) - *2.0/xk; + return (sizeGroups[i+1]->x() - xi)/xk + (xi - sizeGroups[i-1]->x())/xk; } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinaryDsd/uniformBinaryDsd.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.H similarity index 87% rename from applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinaryDsd/uniformBinaryDsd.H rename to applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.H index 246c5c11f..85e50f39a 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinaryDsd/uniformBinaryDsd.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -22,7 +22,7 @@ License along with OpenFOAM. If not, see . Class - Foam::diameterModels::daughterSizeDistributionModels::uniformBinaryDsd + Foam::diameterModels::daughterSizeDistributionModels::uniformBinary Description Uniform binary daughter size distribution. Used for verification and @@ -30,12 +30,12 @@ Description populationBalanceModel class. SourceFiles - uniformBinaryDsd.C + uniformBinary.C \*---------------------------------------------------------------------------*/ -#ifndef uniformBinaryDsd_H -#define uniformBinaryDsd_H +#ifndef uniformBinary_H +#define uniformBinary_H #include "daughterSizeDistributionModel.H" @@ -49,10 +49,10 @@ namespace daughterSizeDistributionModels { /*---------------------------------------------------------------------------*\ - Class uniformBinaryDsd Declaration + Class uniformBinary Declaration \*---------------------------------------------------------------------------*/ -class uniformBinaryDsd +class uniformBinary : public daughterSizeDistributionModel { @@ -64,7 +64,7 @@ public: // Constructor - uniformBinaryDsd + uniformBinary ( const breakupModel& breakup, const dictionary& dict @@ -72,14 +72,14 @@ public: //- Destructor - virtual ~uniformBinaryDsd(); + virtual ~uniformBinary(); // Member Functions //- Return total number of particles assigned to class i when a particle // of class k breaks - virtual dimensionedScalar n(const label i, const label k) const; + virtual dimensionedScalar calcNik(const label i, const label k) const; }; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C index 0250cc614..4406d39eb 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C @@ -143,9 +143,9 @@ void Foam::diameterModels::populationBalanceModel::add(sizeGroup* group) ( "Su", fluid_.time().timeName(), - fluid_.mesh() + mesh_ ), - fluid_.mesh(), + mesh_, dimensionedScalar ( "Su", @@ -163,9 +163,9 @@ void Foam::diameterModels::populationBalanceModel::add(sizeGroup* group) ( "SuSp", fluid_.time().timeName(), - fluid_.mesh() + mesh_ ), - fluid_.mesh(), + mesh_, dimensionedScalar ( "SuSp", @@ -217,7 +217,7 @@ void Foam::diameterModels::populationBalanceModel::createPhasePairs() } -void Foam::diameterModels::populationBalanceModel::preSolve() +void Foam::diameterModels::populationBalanceModel::correct() { calcDeltas(); @@ -264,8 +264,6 @@ birthByCoalescence { const sizeGroup& fj = *sizeGroups_[j]; const sizeGroup& fk = *sizeGroups_[k]; - const volScalarField& alphaj = fj.phase(); - const volScalarField& alphak = fk.phase(); dimensionedScalar Gamma; dimensionedScalar v = fj.x() + fk.x(); @@ -282,21 +280,23 @@ birthByCoalescence // Avoid double counting of events if (j == k) { - Sui_ = 0.5*fi.x()*coalescenceRate_()*fj*alphaj/fj.x()*fk*alphak - /fk.x()*Gamma; + Sui_ = + 0.5*fi.x()*coalescenceRate_()*Gamma + *fj*fj.phase()/fj.x() + *fk*fk.phase()/fk.x(); } else { - Sui_ = fi.x()*coalescenceRate_()*fj*alphaj/fj.x()*fk*alphak/fk.x() - *Gamma; + Sui_ = + fi.x()*coalescenceRate_()*Gamma + *fj*fj.phase()/fj.x() + *fk*fk.phase()/fk.x(); } Su_[i] += Sui_; dimensionedScalar ratio = fj.x()/fi.x(); - const volScalarField& rho = fi.phase().rho(); - const phasePairKey pairij ( fi.phase().name(), @@ -310,7 +310,7 @@ birthByCoalescence Pair::compare(pDmdt_.find(pairij).key(), pairij) ); - pDmdt_[pairij]->ref() += dmdtSign*ratio*Sui_*rho; + pDmdt_[pairij]->ref() += dmdtSign*ratio*Sui_*fi.phase().rho(); } const phasePairKey pairik @@ -326,7 +326,7 @@ birthByCoalescence Pair::compare(pDmdt_.find(pairik).key(), pairik) ); - pDmdt_[pairik]->ref() += dmdtSign*(1 - ratio)*Sui_*rho; + pDmdt_[pairik]->ref() += dmdtSign*(1 - ratio)*Sui_*fi.phase().rho(); } } } @@ -341,20 +341,12 @@ deathByCoalescence { const sizeGroup& fi = *sizeGroups_[i]; const sizeGroup& fj = *sizeGroups_[j]; - const volScalarField& alphai = fi.phase(); - const volScalarField& alphaj = fj.phase(); - SuSp_[i] += - coalescenceRate_() - *alphai - *fj*alphaj/fj.x(); + SuSp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x(); if (i != j) { - SuSp_[j] += - coalescenceRate_() - *alphaj - *fi*alphai/fi.x(); + SuSp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x(); } } @@ -372,13 +364,12 @@ birthByBreakup { const sizeGroup& fi = *sizeGroups_[i]; - Sui_ = fi.x()*breakupRate_()*breakup_[model].dsdPtr()().nik(i, k)*fk - *fk.phase()/fk.x(); + Sui_ = + fi.x()*breakupRate_()*breakup_[model].dsdPtr()().nik(i, k) + *fk*fk.phase()/fk.x(); Su_[i] += Sui_; - const volScalarField& rho = fi.phase().rho(); - const phasePairKey pair ( fi.phase().name(), @@ -392,7 +383,7 @@ birthByBreakup Pair::compare(pDmdt_.find(pair).key(), pair) ); - pDmdt_[pair]->ref() += dmdtSign*Sui_*rho; + pDmdt_[pair]->ref() += dmdtSign*Sui_*fi.phase().rho(); } } } @@ -452,10 +443,8 @@ birthByBinaryBreakup { const sizeGroup& fj = *sizeGroups_[j]; const sizeGroup& fi = *sizeGroups_[i]; - const volScalarField& alphaj = fj.phase(); - const volScalarField& rho = fj.phase().rho(); - Sui_ = fi.x()*binaryBreakupRate_()*delta_[i][j]*fj*alphaj/fj.x(); + Sui_ = fi.x()*binaryBreakupRate_()*delta_[i][j]*fj*fj.phase()/fj.x(); Su_[i] += Sui_; @@ -472,7 +461,7 @@ birthByBinaryBreakup Pair::compare(pDmdt_.find(pairij).key(), pairij) ); - pDmdt_[pairij]->ref() += dmdtSign*Sui_*rho; + pDmdt_[pairij]->ref() += dmdtSign*Sui_*fi.phase().rho(); } dimensionedScalar Gamma; @@ -489,8 +478,9 @@ birthByBinaryBreakup volScalarField& Suk = Sui_; - Suk = sizeGroups_[k]->x()*binaryBreakupRate_()*delta_[i][j]*fj*alphaj - /fj.x()*Gamma; + Suk = + sizeGroups_[k]->x()*binaryBreakupRate_()*delta_[i][j]*Gamma + *fj*fj.phase()/fj.x(); Su_[k] += Suk; @@ -511,7 +501,7 @@ birthByBinaryBreakup ) ); - pDmdt_[pairkj]->ref() += dmdtSign*Suk*rho; + pDmdt_[pairkj]->ref() += dmdtSign*Suk*fi.phase().rho(); } } } @@ -533,7 +523,6 @@ deathByBinaryBreakup void Foam::diameterModels::populationBalanceModel::drift(const label i) { const sizeGroup& fi = *sizeGroups_[i]; - const volScalarField& rho = fi.phase().rho(); if (i == 0) { @@ -545,11 +534,13 @@ void Foam::diameterModels::populationBalanceModel::drift(const label i) } else { - rx_() = pos(driftRate_())*sizeGroups_[i+1]->x()/sizeGroups_[i]->x() + 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(1 - rx_()/(1 - rx_())))*driftRate_() + SuSp_[i] += + (neg(1 - rx_()) + neg(1 - rx_()/(1 - rx_())))*driftRate_() *fi.phase()/((rx_() - 1)*sizeGroups_[i]->x()); rx_() *= 0.0; @@ -559,26 +550,29 @@ void Foam::diameterModels::populationBalanceModel::drift(const label i) { rx_() += pos(driftRate_())*sizeGroups_[i+2]->x()/sizeGroups_[i+1]->x(); - rdx_() += pos(driftRate_()) + rdx_() += + pos(driftRate_()) *(sizeGroups_[i+2]->x() - sizeGroups_[i+1]->x()) /(sizeGroups_[i+1]->x() - sizeGroups_[i]->x()); } else if (i == sizeGroups_.size() - 2) { - rx_() += pos(driftRate_())*sizeGroups_[i+1]->x() - /sizeGroups_[i]->x(); + rx_() += pos(driftRate_())*sizeGroups_[i+1]->x()/sizeGroups_[i]->x(); - rdx_() += pos(driftRate_()) + rdx_() += + pos(driftRate_()) *(sizeGroups_[i+1]->x() - sizeGroups_[i]->x()) /(sizeGroups_[i]->x() - sizeGroups_[i-1]->x()); } if (i == 1) { - rx_() += neg(driftRate_())*sizeGroups_[i-1]->x() + rx_() += + neg(driftRate_())*sizeGroups_[i-1]->x() /sizeGroups_[i]->x(); - rdx_() += neg(driftRate_()) + rdx_() += + neg(driftRate_()) *(sizeGroups_[i]->x() - sizeGroups_[i-1]->x()) /(sizeGroups_[i+1]->x() - sizeGroups_[i]->x()); } @@ -586,7 +580,8 @@ void Foam::diameterModels::populationBalanceModel::drift(const label i) { rx_() += neg(driftRate_())*sizeGroups_[i-2]->x()/sizeGroups_[i-1]->x(); - rdx_() += neg(driftRate_()) + rdx_() += + neg(driftRate_()) *(sizeGroups_[i-1]->x() - sizeGroups_[i-2]->x()) /(sizeGroups_[i]->x() - sizeGroups_[i-1]->x()); } @@ -596,7 +591,9 @@ void Foam::diameterModels::populationBalanceModel::drift(const label i) const sizeGroup& fj = *sizeGroups_[i+1]; volScalarField& Suj = Sui_; - Suj = pos(driftRate_())*driftRate_()*rdx_()*fi*fi.phase()/fi.x() + Suj = + pos(driftRate_())*driftRate_()*rdx_() + *fi*fi.phase()/fi.x() /(rx_() - 1); Su_[i+1] += Suj; @@ -614,7 +611,7 @@ void Foam::diameterModels::populationBalanceModel::drift(const label i) Pair::compare(pDmdt_.find(pairij).key(), pairij) ); - pDmdt_[pairij]->ref() -= dmdtSign*Suj*rho; + pDmdt_[pairij]->ref() -= dmdtSign*Suj*fi.phase().rho(); } } @@ -623,7 +620,9 @@ void Foam::diameterModels::populationBalanceModel::drift(const label i) const sizeGroup& fh = *sizeGroups_[i-1]; volScalarField& Suh = Sui_; - Suh = neg(driftRate_())*driftRate_()*rdx_()*fi*fi.phase()/fi.x() + Suh = + neg(driftRate_())*driftRate_()*rdx_() + *fi*fi.phase()/fi.x() /(rx_() - 1); Su_[i-1] += Suh; @@ -641,7 +640,7 @@ void Foam::diameterModels::populationBalanceModel::drift(const label i) Pair::compare(pDmdt_.find(pairih).key(), pairih) ); - pDmdt_[pairih]->ref() -= dmdtSign*Suh*rho; + pDmdt_[pairih]->ref() -= dmdtSign*Suh*fi.phase().rho(); } } } @@ -796,9 +795,9 @@ void Foam::diameterModels::populationBalanceModel::calcAlphas() { alphas_ *= 0.0; - forAllIter(PtrListDictionary, velocityGroups_, iter) + forAll(velocityGroups_, v) { - alphas_ += iter().phase(); + alphas_ += velocityGroups_[v]->phase(); } } @@ -807,80 +806,15 @@ void Foam::diameterModels::populationBalanceModel::calcVelocity() { U_ *= 0.0; - forAllIter(PtrListDictionary, velocityGroups_, iter) + forAll(velocityGroups_, v) { - U_ += iter().phase().U() - *max(iter().phase(), iter().phase().residualAlpha()) - /max(alphas_, iter().phase().residualAlpha()); + const phaseModel& phase = velocityGroups_[v]->phase(); + + U_ += phase*phase.U()/max(alphas_, phase.residualAlpha()); } } -Foam::tmp -Foam::diameterModels::populationBalanceModel::dsm() const -{ - tmp tDsm - ( - new volScalarField - ( - IOobject - ( - "dsm", - fluid_.time().timeName(), - fluid_.mesh() - ), - fluid_.mesh(), - dimensionedScalar("dsm", dimLength, Zero) - ) - ); - - volScalarField& dsm = tDsm.ref(); - - volScalarField m2 - ( - IOobject - ( - "m2", - fluid_.time().timeName(), - fluid_.mesh() - ), - fluid_.mesh(), - dimensionedScalar("m2", inv(dimLength), Zero) - ); - - volScalarField m3 - ( - IOobject - ( - "m3", - fluid_.time().timeName(), - fluid_.mesh() - ), - fluid_.mesh(), - dimensionedScalar("m3", dimless, Zero) - ); - - forAll(velocityGroups_, i) - { - const velocityGroup& velGroup = *velocityGroups_[i]; - - m2 += velGroup.m2(); - - m3 += velGroup.m3(); - } - - dsm = m3/m2; - - Info<< this->name() << " Sauter mean diameter, min, max = " - << dsm.weightedAverage(dsm.mesh().V()).value() - << ' ' << min(dsm).value() - << ' ' << max(dsm).value() - << endl; - - return tDsm; -} - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::diameterModels::populationBalanceModel::populationBalanceModel @@ -926,10 +860,10 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel IOobject ( "Sui", - fluid_.time().timeName(), - fluid_.mesh() + fluid.time().timeName(), + fluid.mesh() ), - fluid_.mesh(), + mesh_, dimensionedScalar("Sui", inv(dimTime), Zero) ), coalescence_ @@ -986,7 +920,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel ( IOobject ( - IOobject::groupName("alpha", this->name()), + IOobject::groupName("U", this->name()), fluid.time().timeName(), fluid.mesh(), IOobject::NO_READ, @@ -999,8 +933,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel dimVelocity, Zero ) - ), - d_() + ) { this->registerVelocityAndSizeGroups(); @@ -1015,32 +948,6 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel << exit(FatalError); } - if (velocityGroups_.size() > 1) - { - d_.reset - ( - new volScalarField - ( - IOobject - ( - IOobject::groupName("d", this->name()), - fluid.time().timeName(), - fluid.mesh(), - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - fluid.mesh(), - dimensionedScalar - ( - IOobject::groupName("d", this->name()), - dimLength, - Zero - ) - ) - ); - - - } if (coalescence_.size() != 0) { @@ -1125,10 +1032,10 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel IOobject ( "r", - fluid_.time().timeName(), - fluid_.mesh() + fluid.time().timeName(), + fluid.mesh() ), - fluid_.mesh(), + fluid.mesh(), dimensionedScalar("r", dimless, Zero) ) ); @@ -1140,10 +1047,10 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel IOobject ( "r", - fluid_.time().timeName(), - fluid_.mesh() + fluid.time().timeName(), + fluid.mesh() ), - fluid_.mesh(), + fluid.mesh(), dimensionedScalar("r", dimless, Zero) ) ); @@ -1250,11 +1157,11 @@ void Foam::diameterModels::populationBalanceModel::solve() ) ); + calcAlphas(); + calcVelocity(); + if (!solveOnFinalIterOnly || pimple_.finalIter()) { - calcAlphas(); - calcVelocity(); - label nCorr(readLabel(solutionControls.lookup("nCorr"))); scalar tolerance ( @@ -1263,7 +1170,7 @@ void Foam::diameterModels::populationBalanceModel::solve() if (nCorr > 0) { - preSolve(); + correct(); } int iCorr = 0; @@ -1339,11 +1246,6 @@ void Foam::diameterModels::populationBalanceModel::solve() } } - if (velocityGroups_.size() > 1) - { - d_() = dsm(); - } - volScalarField fAlpha0 ( *sizeGroups_.first()*sizeGroups_.first()->phase() diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H index e7a02f2ae..2010ef09e 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H @@ -275,9 +275,6 @@ class populationBalanceModel //- Average velocity volVectorField U_; - //- Mean Sauter diameter - autoPtr d_; - // Private member functions @@ -287,7 +284,7 @@ class populationBalanceModel void createPhasePairs(); - void preSolve(); + void correct(); void birthByCoalescence(const label j, const label k); @@ -315,8 +312,6 @@ class populationBalanceModel void calcVelocity(); - tmp dsm() const; - public: