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)

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

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.C

@@ -49,74 +49,33 @@ namespace diameterModels
 
 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField>
-Foam::diameterModels::velocityGroup::secondMoment() const
-{
-    tmp<volScalarField> 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::volScalarField>
-Foam::diameterModels::velocityGroup::thirdMoment() const
-{
-    tmp<volScalarField> 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::volScalarField> Foam::diameterModels::velocityGroup::dsm() const
 {
-    return
+    tmp<volScalarField> tInvDsm
-        max(min(phase_/m2_, sizeGroups_.last().d()), sizeGroups_.first().d());
+    (
+        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 = "

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<volScalarField> secondMoment() const;
-
-        tmp<volScalarField> thirdMoment() const;
-
         tmp<volScalarField> dsm() const;
 
         tmp<volScalarField> fSum() const;
@@ -180,12 +170,6 @@ public:
         //- Return sizeGroups belonging to this velocityGroup
         inline const PtrList<sizeGroup>& 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<fv::convectionScheme<scalar>>& mvConvection() const;
 

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::fv::convectionScheme<Foam::scalar>>&
 Foam::diameterModels::velocityGroup::mvConvection() const
 {

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_i         |  Volume of droplet i [m3]
-        v_j         |  Volume of bubble j [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]

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)));
             }
         }
     }

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;
 };
 
 

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::
+Foam::diameterModels::daughterSizeDistributionModels::uniformBinary::
-uniformBinaryDsd
+uniformBinary
 (
     const breakupModel& breakup,
     const dictionary& dict
@@ -62,81 +62,35 @@ uniformBinaryDsd
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::diameterModels::daughterSizeDistributionModels::uniformBinaryDsd::
+Foam::diameterModels::daughterSizeDistributionModels::uniformBinary::
-~uniformBinaryDsd()
+~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 dimensionedScalar& xi = breakup_.popBal().sizeGroups()[i]->x();
-    const sizeGroup& fk = *breakup_.popBal().sizeGroups()[k];
+    const dimensionedScalar& xk = breakup_.popBal().sizeGroups()[k]->x();
-
+    const List<sizeGroup*>& sizeGroups = breakup_.popBal().sizeGroups();
-    const dimensionedScalar& xi = fi.x();
-    const dimensionedScalar& xk = fk.x();
-    dimensionedScalar x("x", dimVolume, Zero);
-    dimensionedScalar xii("xii", dimVolume, Zero);
 
     if (i == 0)
     {
-       x = xi;
+        return (sizeGroups[i+1]->x() - xi)/xk;
-    }
-    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;
     }
     else if (i == k)
     {
-        return
+        return (xi - sizeGroups[i-1]->x())/xk;
-            1.0/(xi - x)
-           *(
-                (0.5*sqr(xi) - xi*x)
-              - (0.5*sqr(x) - x*x)
-            )
-           *2.0/xk;
     }
     else
     {
-        return
+        return (sizeGroups[i+1]->x() - xi)/xk + (xi - sizeGroups[i-1]->x())/xk;
-            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;
     }
 }
 

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 <http://www.gnu.org/licenses/>.
 
 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
+#ifndef uniformBinary_H
-#define uniformBinaryDsd_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;
 };
 
 

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
+            Sui_ =
-               /fk.x()*Gamma;
+                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()
+            Sui_ =
-               *Gamma;
+                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<word>::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<word>::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] +=
+    SuSp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x();
-        coalescenceRate_()
-       *alphai
-       *fj*alphaj/fj.x();
 
     if (i != j)
     {
-        SuSp_[j] +=
+        SuSp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x();
-            coalescenceRate_()
-           *alphaj
-           *fi*alphai/fi.x();
     }
 }
 
@@ -372,13 +364,12 @@ birthByBreakup
     {
         const sizeGroup& fi = *sizeGroups_[i];
 
-        Sui_ = fi.x()*breakupRate_()*breakup_[model].dsdPtr()().nik(i, k)*fk
+        Sui_ =
-           *fk.phase()/fk.x();
+            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<word>::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<word>::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
+        Suk =
-           /fj.x()*Gamma;
+            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()
+        rx_() += pos(driftRate_())*sizeGroups_[i+1]->x()/sizeGroups_[i]->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<word>::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<word>::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<velocityGroup>, 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<velocityGroup>, velocityGroups_, iter)
+    forAll(velocityGroups_, v)
     {
-        U_ += iter().phase().U()
+        const phaseModel& phase = velocityGroups_[v]->phase();
-           *max(iter().phase(), iter().phase().residualAlpha())
+
-           /max(alphas_, iter().phase().residualAlpha());
+        U_ += phase*phase.U()/max(alphas_, phase.residualAlpha());
     }
 }
 
 
-Foam::tmp<Foam::volScalarField>
-Foam::diameterModels::populationBalanceModel::dsm() const
-{
-    tmp<volScalarField> 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.time().timeName(),
-            fluid_.mesh()
+            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.time().timeName(),
-                    fluid_.mesh()
+                    fluid.mesh()
                 ),
-                fluid_.mesh(),
+                fluid.mesh(),
                 dimensionedScalar("r", dimless, Zero)
             )
         );
@@ -1140,10 +1047,10 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
                 IOobject
                 (
                     "r",
-                    fluid_.time().timeName(),
+                    fluid.time().timeName(),
-                    fluid_.mesh()
+                    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()

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H

@@ -275,9 +275,6 @@ class populationBalanceModel
         //- Average velocity
         volVectorField U_;
 
-        //- Mean Sauter diameter
-        autoPtr<volScalarField> 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<volScalarField> dsm() const;
-
 
 public: