From 376b51b58b217acba3fe49b082e18f8b1ecf6bb0 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 29 Apr 2022 16:18:03 +0100
Subject: [PATCH] multiphaseEulerFoam::populationBalanceModel: improved
 dilatation treatment

The population balance model considers dilatation originating from density
change and mass transfer via source terms describing nucleation as well as
"drift" of the size distribution to smaller or larger sizes. Numerically, the
treatment does not necessarily equal the total dilatation, hence a correction is
introduced to ensure boundedness of the size group fractions.

Patch contributed by Institute of Fluid Dynamics,
Helmholtz-Zentrum Dresden - Rossendorf (HZDR)
and VTT Technical Research Centre of Finland Ltd.
---
 .../PopulationBalancePhaseSystem.C            |  63 ++++++----
 .../sizeGroup/shapeModels/fractal/fractal.C   |  14 +--
 .../velocityGroup/sizeGroup/sizeGroup.C       |   6 +-
 .../velocityGroup/sizeGroup/sizeGroup.H       |   4 +-
 .../driftModels/constantDrift/constantDrift.C |   4 +-
 .../driftModels/densityChange/densityChange.C |   4 +-
 .../driftModels/phaseChange/phaseChange.C     |  23 ++--
 .../populationBalanceModel.C                  | 115 +++++++++++++-----
 .../populationBalanceModel.H                  |  20 +--
 .../populationBalanceModelI.H                 |  15 +--
 .../populationBalance/drift/constant/fvModels |  49 +-------
 .../drift/constant/physicalProperties.air     |   6 +-
 .../populationBalance/drift/system/fvSolution |   2 +-
 .../isothermalGrowth/validation/createGraphs  |   1 -
 .../negativeDrift/constant/fvModels           |  51 +-------
 .../constant/physicalProperties.air           |   6 +-
 .../negativeDrift/system/fvSolution           |   2 +-
 .../wallBoilingPolydisperse/system/fvSolution |   2 +-
 .../system/fvSolution                         |   2 +-
 19 files changed, 183 insertions(+), 206 deletions(-)

diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
index 593ed919ab..86499a7827 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
@@ -46,20 +46,23 @@ PopulationBalancePhaseSystem
         const diameterModels::populationBalanceModel& popBal =
             populationBalances_[popBali];
 
-        forAll(popBal.velocityGroups(), velGrp1i)
+        forAllConstIter
+        (
+            HashTable<const diameterModels::velocityGroup*>,
+            popBal.velocityGroupPtrs(),
+            iter1
+        )
         {
-            const diameterModels::velocityGroup& velGrp1 =
-                popBal.velocityGroups()[velGrp1i];
+            const diameterModels::velocityGroup& velGrp1 = *iter1();
 
-            for
+            forAllConstIter
             (
-                label velGrp2i = velGrp1i + 1;
-                velGrp2i < popBal.velocityGroups().size();
-                ++ velGrp2i
+                HashTable<const diameterModels::velocityGroup*>,
+                popBal.velocityGroupPtrs(),
+                iter2
             )
             {
-                const diameterModels::velocityGroup& velGrp2 =
-                    popBal.velocityGroups()[velGrp2i];
+                const diameterModels::velocityGroup& velGrp2 = *iter2();
 
                 const phaseInterface interface
                 (
@@ -67,28 +70,36 @@ PopulationBalancePhaseSystem
                     velGrp2.phase()
                 );
 
-                this->template validateMassTransfer
-                    <diameterModels::populationBalanceModel>(interface);
-
-                dmdtfs_.insert
+                if
                 (
-                    interface,
-                    new volScalarField
+                    &velGrp1 != &velGrp2
+                    &&
+                    !dmdtfs_.found(interface)
+                )
+                {
+                    this->template validateMassTransfer
+                        <diameterModels::populationBalanceModel>(interface);
+
+                    dmdtfs_.insert
                     (
-                        IOobject
+                        interface,
+                        new volScalarField
                         (
-                            IOobject::groupName
+                            IOobject
                             (
-                                "populationBalance:dmdtf",
-                                interface.name()
+                                IOobject::groupName
+                                (
+                                    "populationBalance:dmdtf",
+                                    interface.name()
+                                ),
+                                this->mesh().time().timeName(),
+                                this->mesh()
                             ),
-                            this->mesh().time().timeName(),
-                            this->mesh()
-                        ),
-                        this->mesh(),
-                        dimensionedScalar(dimDensity/dimTime, 0)
-                    )
-                );
+                            this->mesh(),
+                            dimensionedScalar(dimDensity/dimTime, 0)
+                        )
+                    );
+                }
             }
         }
     }
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 e56019fe73..b10e1f363c 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
@@ -26,9 +26,9 @@ License
 #include "fractal.H"
 #include "addToRunTimeSelectionTable.H"
 #include "sinteringModel.H"
-#include "fvm.H"
-#include "fvcDdt.H"
-#include "fvcDiv.H"
+#include "fvmDdt.H"
+#include "fvmDiv.H"
+#include "fvmSup.H"
 #include "mixedFvPatchField.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -122,7 +122,7 @@ Foam::diameterModels::shapeModels::fractal::fractal
 
         if
         (
-            isA<mixedFvPatchScalarField>(kappa_.boundaryFieldRef()[patchi])
+            isA<const mixedFvPatchScalarField>(kappa_.boundaryField()[patchi])
         )
         {
             mixedFvPatchScalarField& kappa =
@@ -205,12 +205,12 @@ void Foam::diameterModels::shapeModels::fractal::correct()
 
     fvScalarMatrix kappaEqn
     (
-        fvm::ddt(alpha, fi, kappa_) + fvm::div(fAlphaPhi, kappa_)
+        fvm::ddt(alpha, fi, kappa_)
+      + fvm::div(fAlphaPhi, kappa_)
       ==
       - sinteringModel_->R()
       + Su_
-      - fvm::Sp(popBal.Sp(fi.i()())*fi, kappa_)
-
+      - fvm::Sp(popBal.Sp(fi.i())*fi, kappa_)
       - correction
         (
             fvm::Sp
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 7db66a5957..0f5a8b3477 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.C
@@ -75,7 +75,7 @@ Foam::diameterModels::sizeGroup::sizeGroup
 
         if
         (
-            isA<mixedFvPatchScalarField>(this->boundaryFieldRef()[patchi])
+            isA<const mixedFvPatchScalarField>(this->boundaryField()[patchi])
         )
         {
             mixedFvPatchScalarField& f =
@@ -108,7 +108,7 @@ Foam::diameterModels::sizeGroup::clone() const
 }
 
 
-const Foam::autoPtr<Foam::label>& Foam::diameterModels::sizeGroup::i() const
+const Foam::label& Foam::diameterModels::sizeGroup::i() const
 {
     if (!i_.valid())
     {
@@ -127,7 +127,7 @@ const Foam::autoPtr<Foam::label>& Foam::diameterModels::sizeGroup::i() const
         }
     }
 
-    return i_;
+    return i_();
 }
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.H
index 1fd317d8e4..43121086eb 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.H
@@ -192,8 +192,8 @@ public:
         //- Return const-reference to diameterModel of the phase
         inline const autoPtr<shapeModel>& shapeModelPtr() const;
 
-        //- Return label of the sizeGroup within the population balance
-        const autoPtr<label>& i() const;
+        //- Return index of the size group within the population balance
+        const label& i() const;
 
         //- Return representative surface area of the sizeGroup
         const tmp<volScalarField> a() const;
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C
index de960a56b0..562dfe4f31 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.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-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,7 +63,7 @@ void Foam::diameterModels::driftModels::constantDrift::addToDriftRate
     const label i
 )
 {
-    driftRate += rate_;
+    driftRate += popBal_.sizeGroups()[i].phase()*rate_;
 }
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/densityChange/densityChange.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/densityChange/densityChange.C
index e3da300864..d01d7344f0 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/densityChange/densityChange.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/densityChange/densityChange.C
@@ -70,11 +70,11 @@ void Foam::diameterModels::driftModels::densityChangeDrift::addToDriftRate
     const volScalarField& rho = phase.thermo().rho();
 
     driftRate -=
-        fi.x()/(rho*max(alpha, phase.residualAlpha()))
+        fi.x()
        *(
             fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
           - fvc::Sp(fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), rho)
-        );
+        )/rho;
 }
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/phaseChange/phaseChange.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/phaseChange/phaseChange.C
index 940096ce33..6e73374ba7 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/phaseChange/phaseChange.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/phaseChange/phaseChange.C
@@ -97,15 +97,20 @@ void Foam::diameterModels::driftModels::phaseChange::precompute()
 
     forAll(interfaces_, k)
     {
-        forAll(popBal_.velocityGroups(), j)
+        forAllConstIter
+        (
+            HashTable<const velocityGroup*>,
+            popBal_.velocityGroupPtrs(),
+            iter
+        )
         {
-            const velocityGroup& vgj = popBal_.velocityGroups()[j];
+            const velocityGroup& velGrp = *iter();
 
-            if (interfaces_[k].contains(vgj.phase()))
+            if (interfaces_[k].contains(velGrp.phase()))
             {
-                forAll(vgj.sizeGroups(), i)
+                forAll(velGrp.sizeGroups(), i)
                 {
-                    const sizeGroup& fi = vgj.sizeGroups()[i];
+                    const sizeGroup& fi = velGrp.sizeGroups()[i];
 
                     if (numberWeighted_)
                     {
@@ -128,11 +133,11 @@ void Foam::diameterModels::driftModels::phaseChange::addToDriftRate
     const label i
 )
 {
-    const velocityGroup& vg = popBal_.sizeGroups()[i].VelocityGroup();
+    const velocityGroup& velGrp = popBal_.sizeGroups()[i].VelocityGroup();
 
     forAll(interfaces_, k)
     {
-        if (interfaces_[k].contains(vg.phase()))
+        if (interfaces_[k].contains(velGrp.phase()))
         {
             const volScalarField& dmidtf =
                 popBal_.mesh().lookupObject<volScalarField>
@@ -149,7 +154,7 @@ void Foam::diameterModels::driftModels::phaseChange::addToDriftRate
                 );
 
             const scalar dmidtfSign =
-                interfaces_[k].index(vg.phase()) == 0 ? +1 : -1;
+                interfaces_[k].index(velGrp.phase()) == 0 ? +1 : -1;
 
             const sizeGroup& fi = popBal_.sizeGroups()[i];
 
@@ -163,7 +168,7 @@ void Foam::diameterModels::driftModels::phaseChange::addToDriftRate
                 dDriftRate.ref() *= fi.a();
             }
 
-            driftRate += dDriftRate;
+            driftRate += velGrp.phase()*dDriftRate;
         }
     }
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
index e442a878eb..e4f7116dbe 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
@@ -31,9 +31,12 @@ License
 #include "driftModel.H"
 #include "nucleationModel.H"
 #include "surfaceTensionModel.H"
-#include "fvm.H"
+#include "fvmDdt.H"
+#include "fvmDiv.H"
+#include "fvmSup.H"
 #include "fvcDdt.H"
 #include "fvcDiv.H"
+#include "fvcSup.H"
 #include "shapeModel.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -60,12 +63,26 @@ void Foam::diameterModels::populationBalanceModel::registerVelocityGroups()
 
             if (velGroup.popBalName() == this->name())
             {
-                velocityGroups_.resize(velocityGroups_.size() + 1);
+                velocityGroupPtrs_.insert(velGroup.phase().name(), &velGroup);
 
-                velocityGroups_.set
+                dilatationErrors_.insert
                 (
-                    velocityGroups_.size() - 1,
-                    &const_cast<velocityGroup&>(velGroup)
+                    velGroup.phase().name(),
+                    volScalarField
+                    (
+                        IOobject
+                        (
+                            IOobject::groupName
+                            (
+                                "dilatationError",
+                                velGroup.phase().name()
+                            ),
+                            fluid_.time().timeName(),
+                            mesh_
+                        ),
+                        mesh_,
+                        dimensionedScalar(inv(dimTime), 0)
+                    )
                 );
 
                 forAll(velGroup.sizeGroups(), i)
@@ -446,22 +463,21 @@ void Foam::diameterModels::populationBalanceModel::drift
         const sizeGroup& fe = sizeGroups()[i+1];
         volScalarField& Sue = Sui_;
 
-        Sp_[i] += 1/(fe.x() - fp.x())*pos(driftRate_())*driftRate_()*fp.phase();
+        Sp_[i] += 1/(fe.x() - fp.x())*pos(driftRate_())*driftRate_();
 
         Sue =
-            fe.x()/(fp.x()*(fe.x() - fp.x()))
-           *pos(driftRate_())*driftRate_()*fp*fp.phase();
+            fe.x()/(fp.x()*(fe.x() - fp.x()))*pos(driftRate_())*driftRate_()*fp;
 
         Su_[i+1] += Sue;
 
-        const phaseInterface interfaceij(fp.phase(), fe.phase());
+        const phaseInterface interfacepe(fp.phase(), fe.phase());
 
-        if (pDmdt_.found(interfaceij))
+        if (pDmdt_.found(interfacepe))
         {
             const scalar dmdtSign =
-                interfaceij.index(fp.phase()) == 0 ? +1 : -1;
+                interfacepe.index(fp.phase()) == 0 ? +1 : -1;
 
-            *pDmdt_[interfaceij] -= dmdtSign*Sue*fp.phase().rho();
+            *pDmdt_[interfacepe] -= dmdtSign*Sue*fp.phase().rho();
         }
 
         sizeGroups_[i+1].shapeModelPtr()->addDrift(Sue, fp, model);
@@ -469,7 +485,7 @@ void Foam::diameterModels::populationBalanceModel::drift
 
     if (i == sizeGroups().size() - 1)
     {
-        Sp_[i] -= pos(driftRate_())*driftRate_()*fp.phase()/fp.x();
+        Sp_[i] -= pos(driftRate_())*driftRate_()/fp.x();
     }
 
     if (i > 0)
@@ -477,22 +493,21 @@ void Foam::diameterModels::populationBalanceModel::drift
         const sizeGroup& fw = sizeGroups()[i-1];
         volScalarField& Suw = Sui_;
 
-        Sp_[i] += 1/(fw.x() - fp.x())*neg(driftRate_())*driftRate_()*fp.phase();
+        Sp_[i] += 1/(fw.x() - fp.x())*neg(driftRate_())*driftRate_();
 
         Suw =
-            fw.x()/(fp.x()*(fw.x() - fp.x()))
-           *neg(driftRate_())*driftRate_()*fp*fp.phase();
+            fw.x()/(fp.x()*(fw.x() - fp.x()))*neg(driftRate_())*driftRate_()*fp;
 
         Su_[i-1] += Suw;
 
-        const phaseInterface interfaceih(fp.phase(), fw.phase());
+        const phaseInterface interfacepw(fp.phase(), fw.phase());
 
-        if (pDmdt_.found(interfaceih))
+        if (pDmdt_.found(interfacepw))
         {
             const scalar dmdtSign =
-                interfaceih.index(fp.phase()) == 0 ? +1 : -1;
+                interfacepw.index(fp.phase()) == 0 ? +1 : -1;
 
-            *pDmdt_[interfaceih] -= dmdtSign*Suw*fp.phase().rho();
+            *pDmdt_[interfacepw] -= dmdtSign*Suw*fp.phase().rho();
         }
 
         sizeGroups_[i-1].shapeModelPtr()->addDrift(Suw, fp, model);
@@ -500,7 +515,7 @@ void Foam::diameterModels::populationBalanceModel::drift
 
     if (i == 0)
     {
-        Sp_[i] -= neg(driftRate_())*driftRate_()*fp.phase()/fp.x();
+        Sp_[i] -= neg(driftRate_())*driftRate_()/fp.x();
     }
 }
 
@@ -609,13 +624,43 @@ void Foam::diameterModels::populationBalanceModel::sources()
 }
 
 
+void Foam::diameterModels::populationBalanceModel::correctDilatationError()
+{
+    forAllIter
+    (
+        HashTable<volScalarField>,
+        dilatationErrors_,
+        iter
+    )
+    {
+        volScalarField& dilatationError = iter();
+        const word& phaseName = iter.key();
+        const phaseModel& phase = fluid_.phases()[phaseName];
+        const velocityGroup& velGroup = *velocityGroupPtrs_[phaseName];
+        const volScalarField& alpha = phase;
+        const volScalarField& rho = phase.thermo().rho();
+
+        dilatationError =
+            fvc::ddt(alpha) + fvc::div(phase.alphaPhi())
+          - (fluid_.fvModels().source(alpha, rho) & rho)/rho;
+
+        forAll(velGroup.sizeGroups(), i)
+        {
+            const sizeGroup& fi = velGroup.sizeGroups()[i];
+
+            dilatationError -= Su_[fi.i()] - fvc::Sp(Sp_[fi.i()], fi);
+        }
+    }
+}
+
+
 void Foam::diameterModels::populationBalanceModel::calcAlphas()
 {
     alphas_() = Zero;
 
-    forAll(velocityGroups_, v)
+    forAllConstIter(HashTable<const velocityGroup*>, velocityGroupPtrs_, iter)
     {
-        const phaseModel& phase = velocityGroups_[v].phase();
+        const phaseModel& phase = iter()->phase();
 
         alphas_() += max(phase, phase.residualAlpha());
     }
@@ -637,9 +682,9 @@ Foam::diameterModels::populationBalanceModel::calcDsm()
 
     volScalarField& invDsm = tInvDsm.ref();
 
-    forAll(velocityGroups_, v)
+    forAllConstIter(HashTable<const velocityGroup*>, velocityGroupPtrs_, iter)
     {
-        const phaseModel& phase = velocityGroups_[v].phase();
+        const phaseModel& phase = iter()->phase();
 
         invDsm += max(phase, phase.residualAlpha())/(phase.d()*alphas_());
     }
@@ -652,9 +697,9 @@ void Foam::diameterModels::populationBalanceModel::calcVelocity()
 {
     U_() = Zero;
 
-    forAll(velocityGroups_, v)
+    forAllConstIter(HashTable<const velocityGroup*>, velocityGroupPtrs_, iter)
     {
-        const phaseModel& phase = velocityGroups_[v].phase();
+        const phaseModel& phase = iter()->phase();
 
         U_() += max(phase, phase.residualAlpha())*phase.U()/alphas_();
     }
@@ -704,7 +749,6 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
             IOobject::groupName("alpha", dict_.lookup("continuousPhase"))
         )
     ),
-    velocityGroups_(),
     sizeGroups_(),
     v_(),
     delta_(),
@@ -890,7 +934,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
         );
     }
 
-    if (velocityGroups_.size() > 1)
+    if (velocityGroupPtrs_.size() > 1)
     {
         alphas_.set
         (
@@ -1060,6 +1104,8 @@ void Foam::diameterModels::populationBalanceModel::solve()
                 sources();
             }
 
+            correctDilatationError();
+
             maxInitialResidual = 0;
 
             forAll(sizeGroups(), i)
@@ -1068,15 +1114,18 @@ void Foam::diameterModels::populationBalanceModel::solve()
                 const phaseModel& phase = fi.phase();
                 const volScalarField& alpha = phase;
                 const volScalarField& rho = phase.thermo().rho();
+                const volScalarField& dilatationError =
+                    dilatationErrors_[phase.name()];
 
                 fvScalarMatrix sizeGroupEqn
                 (
-                    fvm::ddt(alpha, fi) + fvm::div(phase.alphaPhi(), fi)
+                    fvm::ddt(alpha, fi)
+                  + fvm::div(phase.alphaPhi(), fi)
+                  - fvm::Sp(dilatationError, fi)
                 ==
-                    Su_[i]
+                    fvc::Su(Su_[i], fi)
                   - fvm::Sp(Sp_[i], fi)
                   + fluid_.fvModels().source(alpha, rho, fi)/rho
-
                   - correction
                     (
                         fvm::Sp
@@ -1121,7 +1170,7 @@ void Foam::diameterModels::populationBalanceModel::solve()
 
 void Foam::diameterModels::populationBalanceModel::correct()
 {
-    if (velocityGroups_.size() > 1)
+    if (velocityGroupPtrs_.size() > 1)
     {
         calcAlphas();
         dsm_() = calcDsm();
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H
index 17bca4f326..62ea1f1b83 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H
@@ -257,10 +257,10 @@ private:
         //- Continuous phase
         const phaseModel& continuousPhase_;
 
-        //- velocityGroups belonging to this populationBalance
-        UPtrList<velocityGroup> velocityGroups_;
+        //- Velocity groups belonging to this populationBalance
+        HashTable<const velocityGroup*> velocityGroupPtrs_;
 
-        //- sizeGroups belonging to this populationBalance
+        //- Size groups belonging to this populationBalance
         UPtrList<sizeGroup> sizeGroups_;
 
         //- sizeGroup boundaries
@@ -275,6 +275,9 @@ private:
         //- Implicitly treated sources
         PtrList<volScalarField> Sp_;
 
+        //- Dilatation errors
+        HashTable<volScalarField> dilatationErrors_;
+
         //- Field for caching sources
         volScalarField Sui_;
 
@@ -357,7 +360,7 @@ private:
 
         void sources();
 
-        void dmdt();
+        void correctDilatationError();
 
         void calcAlphas();
 
@@ -454,15 +457,18 @@ public:
         //- Return continuous phase
         inline const phaseModel& continuousPhase() const;
 
-        //- Return the velocityGroups belonging to this populationBalance
-        inline const UPtrList<velocityGroup>& velocityGroups() const;
+        //- Return the velocity groups belonging to this populationBalance
+        inline const HashTable<const velocityGroup*>& velocityGroupPtrs() const;
 
-        //- Return the sizeGroups belonging to this populationBalance
+        //- Return the size groups belonging to this populationBalance
         inline const UPtrList<sizeGroup>& sizeGroups() const;
 
         //- Return implicit source terms
         inline const volScalarField& Sp(const label i) const;
 
+        //- Return dilatation obtained from source terms
+        inline const HashTable<volScalarField>& sourceDilatation() const;
+
         //- Return coalescence relevant size group pairs
         inline const List<labelPair>& coalescencePairs() const;
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
index f38dd139e8..48b5c8be46 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
@@ -82,10 +82,10 @@ Foam::diameterModels::populationBalanceModel::continuousPhase() const
 }
 
 
-inline const Foam::UPtrList<Foam::diameterModels::velocityGroup>&
-Foam::diameterModels::populationBalanceModel::velocityGroups() const
+inline const Foam::HashTable<const Foam::diameterModels::velocityGroup*>&
+Foam::diameterModels::populationBalanceModel::velocityGroupPtrs() const
 {
-    return velocityGroups_;
+    return velocityGroupPtrs_;
 }
 
 
@@ -120,13 +120,13 @@ Foam::diameterModels::populationBalanceModel::binaryBreakupPairs() const
 inline const Foam::volScalarField&
 Foam::diameterModels::populationBalanceModel::alphas() const
 {
-    if (velocityGroups_.size() > 1)
+    if (velocityGroupPtrs_.size() > 1)
     {
         return alphas_();
     }
     else
     {
-        return velocityGroups_.first().phase();
+        return velocityGroupPtrs_[velocityGroupPtrs_.begin().key()]->phase();
     }
 }
 
@@ -134,13 +134,14 @@ Foam::diameterModels::populationBalanceModel::alphas() const
 inline const Foam::volVectorField&
 Foam::diameterModels::populationBalanceModel::U() const
 {
-    if (velocityGroups_.size() > 1)
+    if (velocityGroupPtrs_.size() > 1)
     {
         return U_();
     }
     else
     {
-        return velocityGroups_.first().phase().U();
+        return
+            velocityGroupPtrs_[velocityGroupPtrs_.begin().key()]->phase().U();
     }
 }
 
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/fvModels b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/fvModels
index 941923760a..cc8a2e57e6 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/fvModels
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/fvModels
@@ -14,51 +14,6 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-airSource
-{
-    type            massSource;
-
-    selectionMode   all;
-
-    massFlowRate    0.02186682;
-
-    phase           air;
-    rho             thermo:rho.air;
-
-    fieldValues
-    {
-        f0.air      0;
-        f1.air      0;
-        f2.air      0;
-        f3.air      0;
-        f4.air      0;
-        f5.air      0;
-        f6.air      0;
-        f7.air      0;
-        f8.air      0;
-        f9.air      0;
-        f10.air     0;
-        f11.air     0;
-        f12.air     0;
-        f13.air     0;
-        f14.air     0;
-        f15.air     0;
-        f16.air     0;
-        f17.air     0;
-        f18.air     0;
-        f19.air     0;
-        f20.air     0;
-        f21.air     0;
-        f22.air     0;
-        f23.air     0;
-        f24.air     0;
-        f25.air     0;
-        f26.air     0;
-        f27.air     0;
-        f28.air     0;
-    }
-}
-
 waterSink
 {
     type            massSource;
@@ -71,7 +26,9 @@ waterSink
     rho             thermo:rho.water;
 
     fieldValues
-    {}
+    {
+        U.water     (0 0 0);
+    }
 }
 
 // ************************************************************************* //
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/physicalProperties.air b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/physicalProperties.air
index 2c9bc0082d..ec67d27ce9 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/physicalProperties.air
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/physicalProperties.air
@@ -20,7 +20,7 @@ thermoType
     mixture         pureMixture;
     transport       const;
     thermo          hConst;
-    equationOfState rhoConst;
+    equationOfState perfectGas;
     specie          specie;
     energy          sensibleInternalEnergy;
 }
@@ -31,10 +31,6 @@ mixture
     {
         molWeight   28.9;
     }
-    equationOfState
-    {
-        rho         1;
-    }
     thermodynamics
     {
         Cp          1007;
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/fvSolution b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/fvSolution
index 5f62647fd7..bab15a15b5 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/fvSolution
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/fvSolution
@@ -60,7 +60,7 @@ PIMPLE
     pRefCell        0;
     pRefValue       1e5;
 
-    flow            no;
+    flow            yes;
     thermophysics   yes;
     models          yes;
 }
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/isothermalGrowth/validation/createGraphs b/test/multiphase/multiphaseEulerFoam/populationBalance/isothermalGrowth/validation/createGraphs
index 7f97180f8d..843ccf6778 100755
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/isothermalGrowth/validation/createGraphs
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/isothermalGrowth/validation/createGraphs
@@ -12,7 +12,6 @@ gnuplot<<EOF
 
     set xlabel 'd (mm)'
     set ylabel 'y (m)'
-    set yrange [0.05:]
     set key bottom center
 
     latestTime = system("foamListTimes -case .. -latestTime")
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/fvModels b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/fvModels
index c165b0f115..81d711b964 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/fvModels
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/fvModels
@@ -14,52 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-airSource
-{
-    type            massSource;
-
-    selectionMode   all;
-
-    massFlowRate    -0.030785654;
-
-    phase           air;
-    rho             thermo:rho.air;
-
-    fieldValues
-    {
-        f0.air      0;
-        f1.air      0;
-        f2.air      0;
-        f3.air      0;
-        f4.air      0;
-        f5.air      0;
-        f6.air      0;
-        f7.air      0;
-        f8.air      0;
-        f9.air      0;
-        f10.air     0;
-        f11.air     0;
-        f12.air     0;
-        f13.air     0;
-        f14.air     0;
-        f15.air     0;
-        f16.air     0;
-        f17.air     0;
-        f18.air     0;
-        f19.air     0;
-        f20.air     0;
-        f21.air     0;
-        f22.air     0;
-        f23.air     0;
-        f24.air     0;
-        f25.air     0;
-        f26.air     0;
-        f27.air     0;
-        f28.air     0;
-    }
-}
-
-waterSink
+waterSource
 {
     type            massSource;
 
@@ -71,7 +26,9 @@ waterSink
     rho             thermo:rho.water;
 
     fieldValues
-    {}
+    {
+        U.water     (0 0 0);
+    }
 }
 
 // ************************************************************************* //
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/physicalProperties.air b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/physicalProperties.air
index 2c9bc0082d..ec67d27ce9 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/physicalProperties.air
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/physicalProperties.air
@@ -20,7 +20,7 @@ thermoType
     mixture         pureMixture;
     transport       const;
     thermo          hConst;
-    equationOfState rhoConst;
+    equationOfState perfectGas;
     specie          specie;
     energy          sensibleInternalEnergy;
 }
@@ -31,10 +31,6 @@ mixture
     {
         molWeight   28.9;
     }
-    equationOfState
-    {
-        rho         1;
-    }
     thermodynamics
     {
         Cp          1007;
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/fvSolution b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/fvSolution
index 5f62647fd7..bab15a15b5 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/fvSolution
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/fvSolution
@@ -60,7 +60,7 @@ PIMPLE
     pRefCell        0;
     pRefValue       1e5;
 
-    flow            no;
+    flow            yes;
     thermophysics   yes;
     models          yes;
 }
diff --git a/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperse/system/fvSolution b/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperse/system/fvSolution
index 772b122868..abc735062b 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperse/system/fvSolution
+++ b/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperse/system/fvSolution
@@ -28,7 +28,7 @@ solvers
         tolerance               1e-4;
         scale                   true;
         solveOnFinalIterOnly    true;
-        sourceUpdateInterval    1;
+        sourceUpdateInterval    20;
     }
 
     p_rgh
diff --git a/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperseTwoGroups/system/fvSolution b/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperseTwoGroups/system/fvSolution
index 772b122868..abc735062b 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperseTwoGroups/system/fvSolution
+++ b/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperseTwoGroups/system/fvSolution
@@ -28,7 +28,7 @@ solvers
         tolerance               1e-4;
         scale                   true;
         solveOnFinalIterOnly    true;
-        sourceUpdateInterval    1;
+        sourceUpdateInterval    20;
     }
 
     p_rgh