reactingEulerFoam: Corrected LaakkonenAlopaeusAittamaa daughter size distribution model

This fix also required a generalization of the corresponding base class,
which allows the user to specify the number of daughter particles per
breakup event separately.

Patch contributed by Institute of Fluid Dynamics,
Helmholtz-Zentrum Dresden - Rossendorf (HZDR)

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/Make/files

@@ -37,7 +37,7 @@ populationBalanceModel/breakupModels/LaakkonenAlopaeusAittamaa/LaakkonenAlopaeus
 populationBalanceModel/breakupModels/powerLaw/powerLaw.C
 
 populationBalanceModel/daughterSizeDistributionModels/daughterSizeDistributionModel/daughterSizeDistributionModel.C
-populationBalanceModel/daughterSizeDistributionModels/uniformBinary/uniformBinary.C
+populationBalanceModel/daughterSizeDistributionModels/uniform/uniform.C
 populationBalanceModel/daughterSizeDistributionModels/LaakkonenAlopaeusAittamaaDsd/LaakkonenAlopaeusAittamaaDsd.C
 
 populationBalanceModel/driftModels/driftModel/driftModel.C

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,8 @@ LaakkonenAlopaeusAittamaaDsd::LaakkonenAlopaeusAittamaaDsd
     const dictionary& dict
 )
 :
-    daughterSizeDistributionModel(breakup, dict)
+    daughterSizeDistributionModel(breakup, dict),
+    C4_(dimensionedScalar::lookupOrDefault("C4", dict, dimless, 4.3))
 {}
 
 
@@ -69,6 +70,26 @@ LaakkonenAlopaeusAittamaaDsd::~LaakkonenAlopaeusAittamaaDsd()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::dimensionedScalar
+Foam::diameterModels::daughterSizeDistributionModels::
+LaakkonenAlopaeusAittamaaDsd::antiderivative
+(
+    const dimensionedScalar& xk,
+    const dimensionedScalar& v,
+    const dimensionedScalar& bndr,
+    const dimensionedScalar range
+) const
+{
+    return
+        (
+            pow(xk, -C4_ - 3)*pow(xk - v, C4_)*(v - xk)*((C4_ + 1)*(C4_ + 2)
+           *(C4_ + 3)*pow3(v) - (C4_ + 1)*(C4_ + 2)*(bndr*(C4_ + 4) - 3.0*xk)
+           *sqr(v) - 2.0*xk*(C4_ + 1)*(bndr*(C4_ + 4) - 3.0*xk)*v - 2.0*bndr
+           *sqr(xk)*C4_ + 6.0*pow3(xk) - 8.0*bndr*sqr(xk))
+        )/(2.0*(range)*(C4_ + 4));
+}
+
+
 Foam::dimensionedScalar
 Foam::diameterModels::daughterSizeDistributionModels::
 LaakkonenAlopaeusAittamaaDsd::calcNik
@@ -77,52 +98,42 @@ LaakkonenAlopaeusAittamaaDsd::calcNik
     const label k
 ) const
 {
-    const dimensionedScalar& xi = breakup_.popBal().sizeGroups()[i].x();
+    const dimensionedScalar& x0 = breakup_.popBal().sizeGroups()[0].x();
-    const dimensionedScalar& xk = breakup_.popBal().sizeGroups()[k].x();
+    dimensionedScalar xi = breakup_.popBal().sizeGroups()[i].x() - x0;
+    dimensionedScalar xk = breakup_.popBal().sizeGroups()[k].x() - x0;
     const UPtrList<sizeGroup>& sizeGroups = breakup_.popBal().sizeGroups();
 
     if (i == 0)
     {
-        const dimensionedScalar& xii = sizeGroups[i+1].x();
+        dimensionedScalar xii = sizeGroups[i+1].x() - x0;
+
+        if (k == 0)
+        {
+            return 1.0;
+        }
 
         return
-            (
+            antiderivative(xk, xi, xii, (xii-xi))
-                5.0*pow4(xii)*sqr(xk) - 12.0*pow5(xi)*xii + 10.0*pow6(xi)
+          - antiderivative(xk, xii, xii, (xii-xi));
-              - (20.0*pow3(xi)*xii - 15.0*pow4(xi))*sqr(xk) - 6.0*pow5(xii)*xk
-              - (24*pow5(xi) - 30*pow4(xi)*xii)*xk + 2*pow6(xii)
-            )
-           /((xii - xi)*pow5(xk));
     }
     else if (i == k)
     {
-        const dimensionedScalar& x = sizeGroups[i-1].x();
+        dimensionedScalar x = sizeGroups[i-1].x() - x0;
 
         return
-            (
+            antiderivative(xk, xi, x, (xi-x))
-                (15.0*pow4(xi) - 20.0*x*pow3(xi))*sqr(xk)
+          - antiderivative(xk, x, x, (xi-x));
-              + 5.0*pow4(x)*sqr(xk) + (30.0*x*pow4(xi) - 24.0*pow5(xi))*xk
-              - 6.0*pow5(x)*xk + 10.0*pow6(xi) - 12.0*x*pow5(xi) + 2.0*pow6(x)
-            )
-           /((xi - x)*pow5(xk));
     }
     else
     {
-        const dimensionedScalar& x = sizeGroups[i-1].x();
+        dimensionedScalar x = sizeGroups[i-1].x() - x0;
-        const dimensionedScalar& xii = sizeGroups[i+1].x();
+        dimensionedScalar xii = sizeGroups[i+1].x() - x0;
 
         return
-            (
+            antiderivative(xk, xi, xii, (xii-xi))
-                5.0*pow4(xii)*sqr(xk) - 12.0*pow5(xi)*xii + 10.0*pow6(xi)
+          - antiderivative(xk, xii, xii, (xii-xi))
-              - (20.0*pow3(xi)*xii - 15.0*pow4(xi))*sqr(xk) - 6.0*pow5(xii)*xk
+          + antiderivative(xk, xi, x, (xi-x))
-              - (24*pow5(xi) - 30*pow4(xi)*xii)*xk + 2*pow6(xii)
+          - antiderivative(xk, x, x, (xi-x));
-            )
-           /((xii - xi)*pow5(xk))
-          + (
-                (15.0*pow4(xi) - 20.0*x*pow3(xi))*sqr(xk)
-              + 5.0*pow4(x)*sqr(xk) + (30.0*x*pow4(xi) - 24.0*pow5(xi))*xk
-              - 6.0*pow5(x)*xk + 10.0*pow6(xi) - 12.0*x*pow5(xi) + 2.0*pow6(x)
-            )
-           /((xi - x)*pow5(xk));
     }
 }
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,14 +27,13 @@ Class
 
 Description
     Daughter size distribution model used by Laakkonen et al. (2006). Note that
-    it must be multiplied by the number of daughter particles (2). Also, the
+    the diameters in the original expression was substituted by bubble volumes
-    coefficient C4 needs to be set to the value of 2 as to ensure number and
+    giving
-    mass conservation following Kumar and Ramkrishna (1996). Also the internal
-    coordinate was changed to particle volume, which gives
 
     \f[
-        \frac{60}{v_j} \left(\frac{v_i}{v_j}\right)^{2}
+        \left(9 + \frac{33}{2}C_4 + 9C_4^2 + \frac{3}{2}C_4^3\right)
-        \left(1 - \frac{v_i}{v_j}\right)^{2}
+        \times \left(\frac{1}{3}\right) \left(\frac{1}{v_j}\right)
+        \left(\frac{v_i}{v_j}\right)^{2} \left(1 - \frac{v_i}{v_j}\right)^{C_4}
     \f]
 
     where
@@ -53,20 +52,10 @@ Description
         Eq. 3, p. 220.
     \endverbatim
 
-    \verbatim
-        Kumar, S., & Ramkrishna, D. (1996).
-        On the solution of population balance equations by discretization-I.
-        A fixed pivot technique.
-        Chemical Engineering Science, 51(8), 1311-1332.
-        Eq. 33, p. 1318.
-    \endverbatim
-
 Usage
     \table
         Property     | Description             | Required    | Default value
-        C1           | Coefficient C1          | no          | 6.0
+        C4           | Coefficient C4          | no          | 4.3
-        C2           | Coefficient C2          | no          | 0.04
-        C3           | Coefficient C3          | no          | 0.01
     \endtable
 
 SourceFiles
@@ -96,6 +85,12 @@ class LaakkonenAlopaeusAittamaaDsd
 :
     public daughterSizeDistributionModel
 {
+    // Private Data
+
+        //- Optional coefficient C4, defaults to 4.3
+        dimensionedScalar C4_;
+
+
 public:
 
     //- Runtime type information
@@ -117,8 +112,16 @@ public:
 
     // Member Functions
 
-        //- Return total number of particles assigned to class i when a particle
+        //- Calculate antiderivative
-        //  of class k breaks
+        dimensionedScalar antiderivative
+        (
+            const dimensionedScalar& xk,
+            const dimensionedScalar& v,
+            const dimensionedScalar& bndr,
+            const dimensionedScalar range
+        ) const;
+
+        //- Calculate contribution to sizeGroup i due to breakup in sizeGroup k
         virtual dimensionedScalar calcNik(const label i, const label k) const;
 };
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,7 +79,18 @@ daughterSizeDistributionModel
     const dictionary& dict
 )
 :
-    breakup_(breakup)
+    breakup_(breakup),
+    particleNumber_
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "particleNumber",
+            dict,
+            dimless,
+            2.0
+        )
+    ),
+    nik_()
 {}
 
 
@@ -93,8 +104,7 @@ Foam::diameterModels::daughterSizeDistributionModel::
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const Foam::dimensionedScalar&
-Foam::diameterModels::daughterSizeDistributionModel::
+Foam::diameterModels::daughterSizeDistributionModel::nik
-nik
 (
     const label i,
     const label k
@@ -114,7 +124,12 @@ void Foam::diameterModels::daughterSizeDistributionModel::correct()
 
             for (label i = 0; i <= k; i++)
             {
-                nik_[k].append(new dimensionedScalar (this->calcNik(i, k)));
+                nik_[k].append
+                (
+                    new dimensionedScalar
+                    (
+                        particleNumber_*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     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,14 @@ Class
     Foam::diameterModels::daughterSizeDistributionModel
 
 Description
-    Base class for daughter size distribution models.
+    Base class for daughter size distribution models. Calculates and stores
+    the contribution to a sizeGroup i due to breakup in a larger sizeGroup k.
+
+Usage
+    \table
+        Property     | Description              | Required    | Default value
+        p            | daughter particle number | no          | 2
+    \endtable
 
 SourceFiles
     daughterSizeDistributionModel.C
@@ -60,7 +67,10 @@ protected:
         //- Reference to the breakupModel
         const breakupModel& breakup_;
 
-        //- nik
+        //- Number of daughter particles, defaults to binary breakup
+        dimensionedScalar particleNumber_;
+
+        //- Contribution to sizeGroup i due to breakup in sizeGroup k
         PtrList<PtrList<dimensionedScalar>> nik_;
 
 
@@ -109,15 +119,13 @@ public:
 
     // Member Functions
 
-        //- Return total number of particles assigned to class i when a particle
+        //- Return contribution to sizeGroup i due to breakup in sizeGroup k
-        //  of class k breaks
         const dimensionedScalar& nik(const label i, const label k) const;
 
         //- Correct field independent expressions
         void correct();
 
-        //- Calculate and return total number of particles assigned to class i
+        //- Calculate contribution to sizeGroup i due to breakup in sizeGroup k
-        //  when a particle of class k breaks
         virtual dimensionedScalar calcNik
         (
             const label i,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,7 +23,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "uniformBinary.H"
+#include "uniform.H"
 #include "addToRunTimeSelectionTable.H"
 #include "breakupModel.H"
 
@@ -35,11 +35,11 @@ namespace diameterModels
 {
 namespace daughterSizeDistributionModels
 {
-    defineTypeNameAndDebug(uniformBinary, 0);
+    defineTypeNameAndDebug(uniform, 0);
     addToRunTimeSelectionTable
     (
         daughterSizeDistributionModel,
-        uniformBinary,
+        uniform,
         dictionary
     );
 }
@@ -49,8 +49,8 @@ namespace daughterSizeDistributionModels
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::diameterModels::daughterSizeDistributionModels::uniformBinary::
+Foam::diameterModels::daughterSizeDistributionModels::uniform::
-uniformBinary
+uniform
 (
     const breakupModel& breakup,
     const dictionary& dict
@@ -62,35 +62,41 @@ uniformBinary
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::diameterModels::daughterSizeDistributionModels::uniformBinary::
+Foam::diameterModels::daughterSizeDistributionModels::uniform::
-~uniformBinary()
+~uniform()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::dimensionedScalar
-Foam::diameterModels::daughterSizeDistributionModels::uniformBinary::calcNik
+Foam::diameterModels::daughterSizeDistributionModels::uniform::calcNik
 (
     const label i,
     const label k
 ) const
 {
+    const dimensionedScalar& x0 = breakup_.popBal().sizeGroups()[0].x();
     const dimensionedScalar& xi = breakup_.popBal().sizeGroups()[i].x();
     const dimensionedScalar& xk = breakup_.popBal().sizeGroups()[k].x();
     const UPtrList<sizeGroup>& sizeGroups = breakup_.popBal().sizeGroups();
 
     if (i == 0)
     {
-        return (sizeGroups[i+1].x() - xi)/xk;
+        if (k == 0)
+        {
+            return 1.0;
+        }
+
+        return (sizeGroups[i+1].x() - xi)/(2.0*(xk - x0));
     }
     else if (i == k)
     {
-        return (xi - sizeGroups[i-1].x())/xk;
+        return (xi - sizeGroups[i-1].x())/(2.0*(xk - x0));
     }
     else
     {
-        return (sizeGroups[i+1].x() - xi)/xk + (xi - sizeGroups[i-1].x())/xk;
+        return (sizeGroups[i+1].x() - sizeGroups[i-1].x())/(2.0*(xk - x0));
     }
 }
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,20 +22,18 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::diameterModels::daughterSizeDistributionModels::uniformBinary
+    Foam::diameterModels::daughterSizeDistributionModels::uniform
 
 Description
-    Uniform binary daughter size distribution. Used for verification and
+    Uniform daughter size distribution.
-    validation of the breakup formulation implemented in the
-    populationBalanceModel class.
 
 SourceFiles
-    uniformBinary.C
+    uniform.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef uniformBinary_H
+#ifndef uniform_H
-#define uniformBinary_H
+#define uniform_H
 
 #include "daughterSizeDistributionModel.H"
 
@@ -49,22 +47,22 @@ namespace daughterSizeDistributionModels
 {
 
 /*---------------------------------------------------------------------------*\
-                      Class uniformBinary Declaration
+                           Class uniform Declaration
 \*---------------------------------------------------------------------------*/
 
-class uniformBinary
+class uniform
 :
     public daughterSizeDistributionModel
 {
 public:
 
     //- Runtime type information
-    TypeName("uniformBinary");
+    TypeName("uniform");
 
 
     // Constructor
 
-        uniformBinary
+        uniform
         (
             const breakupModel& breakup,
             const dictionary& dict
@@ -72,13 +70,12 @@ public:
 
 
     //- Destructor
-    virtual ~uniformBinary();
+    virtual ~uniform();
 
 
     // Member Functions
 
-        //- Return total number of particles assigned to class i when a particle
+        //- Calculate contribution to sizeGroup i due to breakup in sizeGroup k
-        //  of class k breaks
         virtual dimensionedScalar calcNik(const label i, const label k) const;
 };
 

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

@@ -122,7 +122,7 @@ Usage
                 {
                     C 0.5;
                     exponent 0.01;
-                    daughterSizeDistributionModel uniformBinary;
+                    daughterSizeDistributionModel uniform;
                 }
             );
 

test/multiphase/reactingMultiphaseEulerFoam/populationBalanceModeling/breakup/constant/phaseProperties

@@ -133,7 +133,7 @@ populationBalanceCoeffs
             powerLaw
             {
                 power   2.0;
-                daughterSizeDistributionModel uniformBinary;
+                daughterSizeDistributionModel uniform;
             }
         );
 

test/multiphase/reactingMultiphaseEulerFoam/populationBalanceModeling/simultaneousCoalescenceAndBreakup/constant/phaseProperties

@@ -138,7 +138,7 @@ populationBalanceCoeffs
             {
                 C 0.1;
                 exponent 0.01;
-                daughterSizeDistributionModel uniformBinary;
+                daughterSizeDistributionModel uniform;
             }
         );
 

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

@@ -114,7 +114,7 @@ populationBalanceCoeffs
             {
                 C 0.5;
                 exponent 0.01;
-                daughterSizeDistributionModel uniformBinary;
+                daughterSizeDistributionModel uniform;
             }
         );
 

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

@@ -93,7 +93,7 @@ populationBalanceCoeffs
             {
                 C 10.0;
                 exponent 1.0;
-                daughterSizeDistributionModel uniformBinary;
+                daughterSizeDistributionModel uniform;
             }
         );
 

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoilingPolyDisperse/constant/phaseProperties

@@ -91,7 +91,7 @@ populationBalanceCoeffs
             {
                 C 0.5;
                 exponent 0.01;
-                daughterSizeDistributionModel uniformBinary;
+                daughterSizeDistributionModel uniform;
             }
         );