reactingEulerFoam: improvements to population balance modeling

Removed possibility for the user to specify a driftRate in the constantDrift
model which is independent of a fvOptions mass source. The driftRate must be
calculated from/be consistent with the mass source in order to yield a particle
number conserving result.

Made calculation of the over-all Sauter mean diameter of an entire population
balance conditional on more than one velocityGroup being present. This diameter
field is for post-processing purposes only and would be redundant in case of one
velocityGroup being used.

Solution control is extended to allow for solution of the population balance
equation at the last PIMPLE loop only, using an optional switch. This can be
beneficial in terms of simulation time as well as coupling between the
population balance based diameter calculation and the rest of the equation
system.

Patch contributed by Institute of Fluid Dynamics, Helmholtz-Zentrum Dresden - Rossendorf
(HZDR) and VTT Technical Research Centre of Finland Ltd.

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.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
@@ -51,8 +51,6 @@ Foam::diameterModels::driftModels::constantDrift::constantDrift
 )
 :
     driftModel(popBal, dict),
-    calculateFromFvOptions_(dict.lookup("calculateFromFvOptions")),
-    rate_("rate", dimVolume/dimTime, dict.lookupOrDefault<scalar>("rate", 1.0)),
     N_
     (
         IOobject
@@ -64,12 +62,7 @@ Foam::diameterModels::driftModels::constantDrift::constantDrift
         popBal.mesh(),
         dimensionedScalar("Sui", inv(dimVolume), Zero)
     )
-{
-    if (calculateFromFvOptions_ == false)
-    {
-        rate_.read(dict);
-    }
-}
+{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -97,15 +90,7 @@ void Foam::diameterModels::driftModels::constantDrift::driftRate
     phaseModel& phase = const_cast<phaseModel&>(fi.phase());
     volScalarField& rho = phase.thermo().rho();
 
-    if (calculateFromFvOptions_ == true)
-    {
-        driftRate +=(popBal_.fluid().fvOptions()(phase, rho)&rho)
-           /(N_*rho);
-    }
-    else
-    {
-        driftRate += rate_;
-    }
+    driftRate += (popBal_.fluid().fvOptions()(phase, rho)&rho)/(N_*rho);
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.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
@@ -25,10 +25,10 @@ Class
     Foam::diameterModels::driftModels::constant
 
 Description
-    Constant growth rate. Used for verification and validation of the drift
-    formulation implemented in the populationBalanceModel class. Rate is either
-    given in the corresponding dictionary or calculated from fvOptions mass
-    source, depending on Switch calculateFromFvOptions.
+    Constant drift rate within all classes. Used for verification and
+    validation of the drift formulation implemented in the
+    populationBalanceModel class. Rate is calculated from fvOptions mass source.
+
 
 SourceFiles
     constant.C
@@ -59,12 +59,6 @@ class constantDrift
 {
     // Private data
 
-        //- If true, calculate drift rate from fvOptions mass source
-        Switch calculateFromFvOptions_;
-
-        //- Optional driftRate
-        dimensionedScalar rate_;
-
         //- Total number concentration
         volScalarField N_;
 

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

@@ -969,6 +969,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
     (
         fluid.subDict("populationBalanceCoeffs").subDict(name_)
     ),
+    pimple_(mesh_.lookupObject<pimpleControl>("solutionControl")),
     continuousPhase_
     (
         mesh_.lookupObject<phaseModel>
@@ -1044,24 +1045,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
             Zero
         )
     ),
-    d_
-    (
-        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
-        )
-    )
+    d_()
 {
     this->registerVelocityGroups();
 
@@ -1083,6 +1067,30 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
                 dimensionedScalar("Sui", inv(dimTime), Zero)
             )
         );
+
+        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)
@@ -1201,13 +1209,27 @@ bool Foam::diameterModels::populationBalanceModel::writeData(Ostream& os) const
 
 void Foam::diameterModels::populationBalanceModel::solve()
 {
-    const dictionary& populationBalanceControls = mesh_.solverDict(name_);
-    label nCorr(readLabel(populationBalanceControls.lookup("nCorr")));
-    scalar tolerance(readScalar(populationBalanceControls.lookup("tolerance")));
+    const dictionary& solutionControls = mesh_.solverDict(name_);
+    bool solveOnFinalIterOnly
+        (
+            solutionControls.lookupOrDefault<bool>
+            (
+                "solveOnFinalIterOnly",
+                false
+            )
+        );
 
+    if (!solveOnFinalIterOnly || pimple_.finalIter())
+    {
         calcAlphas();
         calcVelocity();
 
+        label nCorr(readLabel(solutionControls.lookup("nCorr")));
+        scalar tolerance
+            (
+                readScalar(solutionControls.lookup("tolerance"))
+            );
+
         if (nCorr > 0)
         {
             preSolve();
@@ -1286,7 +1308,11 @@ void Foam::diameterModels::populationBalanceModel::solve()
             }
         }
 
-    d_ = dsm();
+        if (velocityGroups_.size() > 1)
+        {
+            d_() = dsm();
+        }
+    }
 }
 
 // ************************************************************************* //

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

@@ -136,6 +136,7 @@ SourceFiles
 
 #include "sizeGroup.H"
 #include "phasePair.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -184,6 +185,9 @@ class populationBalanceModel
         //- Dictionary
         dictionary dict_;
 
+        //- Reference to pimpleControl
+        const pimpleControl& pimple_;
+
         //- Continuous phase
         const phaseModel& continuousPhase_;
 
@@ -245,7 +249,7 @@ class populationBalanceModel
         volVectorField U_;
 
         //- Mean Sauter diameter
-        volScalarField d_;
+        autoPtr<volScalarField> d_;
 
 
     // Private member functions

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/constant/phaseProperties

@@ -122,9 +122,6 @@ populationBalanceCoeffs
         (
             densityChange{}
         );
-
-        zeroethOrderModels
-        ();
     }
 }
 

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnPolydisperse/constant/phaseProperties

@@ -104,9 +104,6 @@ populationBalanceCoeffs
         (
             densityChange{}
         );
-
-        zeroethOrderModels
-        ();
     }
 }
 

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnPolydisperse/system/fvSolution

@@ -28,6 +28,7 @@ solvers
         nCorr                1;
         tolerance            1e-4;
         renormalize          false;
+        solveOnFinalIterOnly true;
     }
 
     p_rgh