reactingEulerFoam: Added population balance modeling capability

This patch enables the reactingEulerFoam solvers to simulate polydisperse flow
situations, i.e. flows where the disperse phase is subject to a size
distribution.

The newly added populationBalanceModel class solves the integro-partial
differential population balance equation (PBE) by means of a class method, also
called discrete or sectional method. This approach is based on discretizing the
PBE over its internal coordinate, the particle volume. This yields a set of
transport equations for the number concentration of particles in classes with a
different representative size. These are coupled through their source-terms and
solved in a segregated manner. The implementation is done in a way, that the
total particle number and mass is preserved for coalescence, breakup and drift
(i.e. isothermal growth or phase change) processes, irrespective of the chosen
discretization over the internal coordinate.

A population balance can be split over multiple velocity (temperature) fields,
using the capability of reactingMultiphaseEulerFoam to solve for n momentum
(energy) equations. To a certain degree, this takes into account the dependency
of heat- and momentum transfer on the disperse phase diameter. It is also possible
to define multiple population balances, e.g. bubbles and droplets simultaneously.

The functionality can be switched on by choosing the appropriate phaseSystem
type, e.g. populationBalanceMultiphaseSystem and the newly added diameterModel
class called velocityGroup. To illustrate the use of the functionality, a
bubbleColumnPolydisperse tutorial was added for reactingTwoPhaseEulerFoam and
reactingMultiphaseEulerFoam.

Furthermore, a reactingEulerFoam-specific functionObject called sizeDistribution
was added to allow post-Processing of the size distribution, e.g. to obtain the
number density function in a specific region.

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

applications/solvers/multiphase/reactingEulerFoam/Allwclean

@@ -7,5 +7,6 @@ wclean libso interfacialCompositionModels
 wclean libso derivedFvPatchFields
 reactingTwoPhaseEulerFoam/Allwclean
 reactingMultiphaseEulerFoam/Allwclean
+wclean libso functionObjects
 
 #------------------------------------------------------------------------------

applications/solvers/multiphase/reactingEulerFoam/Allwmake

@@ -14,5 +14,6 @@ wmakeLnInclude reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels
 wmakeLnInclude reactingMultiphaseEulerFoam/multiphaseCompressibleTurbulenceModels
 reactingTwoPhaseEulerFoam/Allwmake $targetType $*
 reactingMultiphaseEulerFoam/Allwmake $targetType $*
+wmake $targetType functionObjects
 
 #------------------------------------------------------------------------------

applications/solvers/multiphase/reactingEulerFoam/functionObjects/Make/files

@@ -0,0 +1,3 @@
+sizeDistribution/sizeDistribution.C
+
+LIB = $(FOAM_LIBBIN)/libreactingEulerFoamFunctionObjects

applications/solvers/multiphase/reactingEulerFoam/functionObjects/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/functionObjects/field/lnInclude \
+    -I../phaseSystems/lnInclude
+
+LIB_LIBS = \
+    -lfieldFunctionObjects \
+    -lfiniteVolume

applications/solvers/multiphase/reactingEulerFoam/functionObjects/sizeDistribution/sizeDistribution.C

@@ -0,0 +1,602 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sizeDistribution.H"
+#include "sizeGroup.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(sizeDistribution, 0);
+    addToRunTimeSelectionTable(functionObject, sizeDistribution, dictionary);
+}
+}
+
+template<>
+const char*
+Foam::NamedEnum
+<
+    Foam::functionObjects::sizeDistribution::regionTypes,
+    2
+>::names[] = {"cellZone", "all"};
+
+template<>
+const char*
+Foam::NamedEnum
+<
+    Foam::functionObjects::sizeDistribution::functionTypes,
+    4
+>::names[] =
+    {
+        "numberDensity",
+        "volumeDensity",
+        "numberConcentration",
+        "moments"
+    };
+
+template<>
+const char*
+Foam::NamedEnum
+<
+    Foam::functionObjects::sizeDistribution::abszissaTypes,
+    2
+>::names[] = {"diameter", "volume"};
+
+const Foam::NamedEnum
+<
+    Foam::functionObjects::sizeDistribution::regionTypes,
+    2
+> Foam::functionObjects::sizeDistribution::regionTypeNames_;
+
+const Foam::NamedEnum
+<
+    Foam::functionObjects::sizeDistribution::functionTypes,
+    4
+> Foam::functionObjects::sizeDistribution::functionTypeNames_;
+
+const Foam::NamedEnum
+<
+    Foam::functionObjects::sizeDistribution::abszissaTypes,
+    2
+> Foam::functionObjects::sizeDistribution::abszissaTypeNames_;
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::sizeDistribution::initialise
+(
+    const dictionary& dict
+)
+{
+    switch (functionType_)
+    {
+        case ftNdf:
+        {
+            break;
+        }
+
+        case ftVdf:
+        {
+            break;
+        }
+
+        case ftNc:
+        {
+            break;
+        }
+
+        case ftMom:
+        {
+            break;
+        }
+
+        default:
+        {
+            FatalErrorInFunction
+               << "Unknown function type. Valid function types are:"
+                << functionTypeNames_ << nl << exit(FatalError);
+        }
+    }
+
+    switch (abszissaType_)
+    {
+        case atDiameter:
+        {
+            break;
+        }
+
+        case atVolume:
+        {
+            break;
+        }
+
+        default:
+        {
+            FatalErrorInFunction
+               << "Unknown abszissa type. Valid abszissa types are:"
+                << abszissaTypeNames_ << nl << exit(FatalError);
+        }
+    }
+
+    setCellZoneCells();
+
+    if (nCells_ == 0)
+    {
+        FatalErrorInFunction
+            << type() << " " << name() << ": "
+            << regionTypeNames_[regionType_] << "(" << regionName_ << "):" << nl
+            << "    Region has no cells" << exit(FatalError);
+    }
+
+    volume_ = volume();
+
+    Info<< type() << " " << name() << ":"
+        << regionTypeNames_[regionType_] << "(" << regionName_ << "):" << nl
+        << "    total cells  = " << nCells_ << nl
+        << "    total volume = " << volume_
+        << nl << endl;
+
+    Info<< nl << endl;
+}
+
+
+void Foam::functionObjects::sizeDistribution::setCellZoneCells()
+{
+    switch (regionType_)
+    {
+        case rtCellZone:
+        {
+            dict().lookup("name") >> regionName_;
+
+            label zoneId = mesh().cellZones().findZoneID(regionName_);
+
+            if (zoneId < 0)
+            {
+                FatalErrorInFunction
+                    << "Unknown cell zone name: " << regionName_
+                    << ". Valid cell zones are: " << mesh().cellZones().names()
+                    << nl << exit(FatalError);
+            }
+
+            cellId_ = mesh().cellZones()[zoneId];
+            nCells_ = returnReduce(cellId_.size(), sumOp<label>());
+            break;
+        }
+
+        case rtAll:
+        {
+            cellId_ = identity(mesh().nCells());
+            nCells_ = returnReduce(cellId_.size(), sumOp<label>());
+            break;
+        }
+
+        default:
+        {
+            FatalErrorInFunction
+               << "Unknown region type. Valid region types are:"
+                << regionTypeNames_ << nl << exit(FatalError);
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "Selected region size = " << cellId_.size() << endl;
+    }
+}
+
+
+Foam::scalar Foam::functionObjects::sizeDistribution::volume() const
+{
+    return gSum(filterField(mesh().V()));
+}
+
+
+void Foam::functionObjects::sizeDistribution::combineFields(scalarField& field)
+{
+    List<scalarField> allValues(Pstream::nProcs());
+
+    allValues[Pstream::myProcNo()] = field;
+
+    Pstream::gatherList(allValues);
+
+    if (Pstream::master())
+    {
+        field =
+            ListListOps::combine<scalarField>
+            (
+                allValues,
+                accessOp<scalarField>()
+            );
+    }
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::functionObjects::sizeDistribution::filterField
+(
+    const scalarField& field
+) const
+{
+    return tmp<scalarField>(new scalarField(field, cellId_));
+}
+
+
+void Foam::functionObjects::sizeDistribution::writeFileHeader
+(
+    const label i
+)
+{
+    OFstream& file = this->file();
+
+    switch (functionType_)
+    {
+        case ftNdf:
+        {
+            writeHeader(file, "Number density function");
+            break;
+        }
+
+        case ftVdf:
+        {
+            writeHeader(file, "Volume density function");
+            break;
+        }
+
+        case ftNc:
+        {
+            writeHeader(file, "Number concentration");
+            break;
+        }
+
+        case ftMom:
+        {
+            writeHeader(file, "Moments");
+            break;
+        }
+    }
+
+    switch (abszissaType_)
+    {
+        case atVolume:
+        {
+            writeCommented(file, "Time/volume");
+            break;
+        }
+
+        case atDiameter:
+        {
+            writeCommented(file, "Time/diameter");
+            break;
+        }
+    }
+
+    switch (functionType_)
+    {
+        case ftMom:
+        {
+            for (label i = 0; i <= momentOrder_; i++)
+            {
+                file() << tab << i;
+            }
+
+            break;
+        }
+        default:
+        {
+            forAll(popBal_.sizeGroups(), sizeGroupi)
+            {
+                const diameterModels::sizeGroup& fi =
+                    *popBal_.sizeGroups()[sizeGroupi];
+
+                switch (abszissaType_)
+                {
+                    case atDiameter:
+                    {
+                        file() << tab  << fi.d().value();
+
+                        break;
+                    }
+
+                    case atVolume:
+                    {
+                        file() << tab  << fi.x().value();
+
+                        break;
+                    }
+                }
+            }
+
+            break;
+        }
+    }
+
+    file << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::sizeDistribution::sizeDistribution
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    logFiles(obr_, name),
+    dict_(dict),
+    regionType_(regionTypeNames_.read(dict.lookup("regionType"))),
+    regionName_(word::null),
+    functionType_(functionTypeNames_.read(dict.lookup("functionType"))),
+    abszissaType_(abszissaTypeNames_.read(dict.lookup("abszissaType"))),
+    nCells_(0),
+    cellId_(),
+    volume_(0.0),
+    writeVolume_(dict.lookupOrDefault("writeVolume", false)),
+    popBal_
+    (
+        obr_.lookupObject<Foam::diameterModels::populationBalanceModel>
+        (
+            dict.lookup("populationBalance")
+        )
+    ),
+    N_(popBal_.sizeGroups().size()),
+    momentOrder_(dict.lookupOrDefault<label>("momentOrder", 0)),
+    normalize_(dict.lookupOrDefault("normalize", false)),
+    sumN_(0.0),
+    sumV_(0.0)
+{
+    read(dict);
+    resetName(name);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::sizeDistribution::~sizeDistribution()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::sizeDistribution::read(const dictionary& dict)
+{
+    if (dict != dict_)
+    {
+        dict_ = dict;
+    }
+
+    fvMeshFunctionObject::read(dict);
+
+    initialise(dict);
+
+    return true;
+}
+
+
+bool Foam::functionObjects::sizeDistribution::execute()
+{
+    return true;
+}
+
+
+bool Foam::functionObjects::sizeDistribution::write()
+{
+    logFiles::write();
+
+    if (Pstream::master())
+    {
+        writeTime(file());
+    }
+
+    Log << type() << " " << name() << " write" << nl;
+
+    scalarField V(filterField(mesh().V()));
+    combineFields(V);
+
+    sumN_ *= 0.0;
+    sumV_ *= 0.0;
+
+    forAll(N_, i)
+    {
+        const Foam::diameterModels::sizeGroup& fi =
+            *popBal_.sizeGroups()[i];
+
+        const volScalarField& alpha = fi.VelocityGroup().phase();
+
+        scalarField Ni(fi*alpha/fi.x());
+        scalarField values(filterField(Ni));
+        scalarField V(filterField(mesh().V()));
+
+        // Combine onto master
+        combineFields(values);
+        combineFields(V);
+
+        if (Pstream::master())
+        {
+            // Calculate volume-averaged number concentration
+            N_[i] = sum(V*values)/sum(V);
+        }
+
+        sumN_ += N_[i];
+
+        sumV_ += N_[i]*fi.x().value();
+    }
+
+    if (Pstream::master())
+    {
+        switch (functionType_)
+        {
+            case ftMom:
+            {
+                for (label m = 0; m <= momentOrder_; m++)
+                {
+                    scalar result(0.0);
+
+                    forAll(N_, i)
+                    {
+                        const Foam::diameterModels::sizeGroup& fi =
+                            *popBal_.sizeGroups()[i];
+
+                        switch (abszissaType_)
+                        {
+                            case atVolume:
+                            {
+                                result += pow(fi.x().value(), m)*N_[i];
+
+                                break;
+                            }
+
+                            case atDiameter:
+                            {
+                                result += pow(fi.d().value(), m)*N_[i];
+
+                                break;
+                            }
+                        }
+                    }
+
+                    file() << tab << result;
+                }
+
+                break;
+            }
+
+            default:
+            {
+                forAll(popBal_.sizeGroups(), i)
+                {
+                    const Foam::diameterModels::sizeGroup& fi =
+                        *popBal_.sizeGroups()[i];
+
+                    scalar result(0.0);
+                    scalar delta(0.0);
+
+                    switch (abszissaType_)
+                    {
+                        case atVolume:
+                        {
+                            delta = popBal_.v()[i+1].value()
+                              - popBal_.v()[i].value();
+
+                            break;
+                        }
+
+                        case atDiameter:
+                        {
+                            const scalar& formFactor =
+                                fi.VelocityGroup().formFactor().value();
+
+                            delta =
+                                pow
+                                (
+                                    popBal_.v()[i+1].value()
+                                   /formFactor,
+                                    1.0/3.0
+                                )
+                              - pow
+                                (
+                                    popBal_.v()[i].value()
+                                   /formFactor,
+                                    1.0/3.0
+                                );
+
+                            break;
+                        }
+                    }
+
+                    switch (functionType_)
+                    {
+                        case ftNdf:
+                        {
+                            if (normalize_ == true)
+                            {
+                                result = N_[i]/delta/sumN_;
+                            }
+                            else
+                            {
+                                result = N_[i]/delta;
+                            }
+
+                            break;
+                        }
+
+                        case ftVdf:
+                        {
+                            if (normalize_ == true)
+                            {
+                                result = N_[i]*fi.x().value()/delta/sumV_;
+                            }
+                            else
+                            {
+                                result = N_[i]*fi.x().value()/delta;
+                            }
+
+                            break;
+                        }
+
+                        case ftNc:
+                        {
+                            if (normalize_ == true)
+                            {
+                                result = N_[i]/sumN_;
+                            }
+                            else
+                            {
+                                result = N_[i];
+                            }
+
+                            break;
+                        }
+
+                        default:
+                        {
+                            break;
+                        }
+                    }
+
+                    file()<< tab << result;
+                }
+            }
+        }
+    }
+
+    if (Pstream::master())
+    {
+        file()<< endl;
+    }
+
+    Log << endl;
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/functionObjects/sizeDistribution/sizeDistribution.H

@@ -0,0 +1,280 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::sizeDistribution
+
+Group
+    grpreactingEulerFoamFunctionObjects
+
+Description
+    This function object calculates and outputs information on the
+    sizeDistribution as calculated by the populationBalanceModel class
+    implemented in the reactionEulerFoam solver suite. It is build up very
+    similar to the fieldValue/volRegion functionObjects, since it must be
+    possible to obtain information about the populationBalanceModel within a
+    defined region.
+
+    Example of function object specification:
+    \verbatim
+    box.all.numberDensity.volume.bubbles
+    {
+        type sizeDistribution;
+        libs ("libreactingEulerFoamFunctionObjects.so");
+        writeControl    outputTime;
+        writeInterval   1;
+        log             true;
+        ...
+        functionType        numberDensity;
+        abszissaType        volume;
+        regionType          all;
+        populationBalanceModel   bubbles;
+        normalize           true;
+    }
+    \endverbatim
+
+Usage
+    \table
+        Property        | Description           | Required    | Default value
+        type            | type name: sizeDistribution | yes |
+        functionType    | numberDensity, volumeDensity, numberConcentration,
+                          moments | yes |
+        abszissaType    | volume, diameter | yes |
+        momentOrder     | Write moment up to given order | no | 0
+        regionType      | Evaluate for cellZone or entire mesh | yes |
+        cellZoneName    | Required if regionType is cellZone |   |
+        populationBalanceModel    | Respective populationBalanceModel | yes |
+        normalize | Normalization | no |
+    \endtable
+
+See also
+    Foam::diameterModels::populationBalanceModel
+    Foam::functionObjects::fvMeshFunctionObject
+    Foam::functionObjects::logFiles
+    Foam::functionObject
+
+SourceFiles
+    sizeDistribution.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_sizeDistribution_H
+#define functionObjects_sizeDistribution_H
+
+#include "fvMeshFunctionObject.H"
+#include "logFiles.H"
+#include "populationBalanceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class fvMesh;
+
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class sizeDistribution Declaration
+\*---------------------------------------------------------------------------*/
+
+class sizeDistribution
+:
+    public fvMeshFunctionObject,
+    public logFiles
+{
+
+public:
+
+    // Public data types
+
+        //- Region type enumeration
+        enum regionTypes
+        {
+            rtCellZone,
+            rtAll
+        };
+
+        //- Region type names
+        static const NamedEnum<regionTypes, 2> regionTypeNames_;
+
+
+        //- Function type enumeration
+        enum functionTypes
+        {
+            ftNdf,
+            ftVdf,
+            ftNc,
+            ftMom
+        };
+
+        //- Function type names
+        static const NamedEnum<functionTypes, 4> functionTypeNames_;
+
+
+        //- abszissa type enumeration
+        enum abszissaTypes
+        {
+            atDiameter,
+            atVolume,
+        };
+
+        //- Abszissa type names
+        static const NamedEnum<abszissaTypes, 2> abszissaTypeNames_;
+
+
+protected:
+
+    // Protected data
+
+        //- Construction dictionary
+        dictionary dict_;
+
+        //- Region type
+        regionTypes regionType_;
+
+        //- Name of region
+        word regionName_;
+
+        //- Function type
+        functionTypes functionType_;
+
+        //- Abszissa type
+        abszissaTypes abszissaType_;
+
+        //- Global number of cells
+        label nCells_;
+
+        //- Local list of cell IDs
+        labelList cellId_;
+
+        //- Total volume of the evaluated region
+        scalar volume_;
+
+        //- Optionally write the volume of the sizeDistribution
+        bool writeVolume_;
+
+        //- PopulationBalance
+        const Foam::diameterModels::populationBalanceModel& popBal_;
+
+        //- Number concentrations
+        List<scalar> N_;
+
+        //- Write moments up to spezified order with respect to abszissaType
+        label momentOrder_;
+
+        //- Normalization switch
+        const Switch normalize_;
+
+        //- Sum of number concentrations
+        scalar sumN_;
+
+        //- Volumertic sum
+        scalar sumV_;
+
+
+    // Protected Member Functions
+
+        //- Initialise, e.g. cell addressing
+        void initialise(const dictionary& dict);
+
+        //- Set cells to evaluate based on a cell zone
+        void setCellZoneCells();
+
+        //- Calculate and return volume of the evaluated cell zone
+        scalar volume() const;
+
+        //- Combine fields from all processor domains into single field
+        void combineFields(scalarField& field);
+
+        //- Filter field according to cellIds
+        tmp<scalarField> filterField(const scalarField& field) const;
+
+        //- Output file header information
+        virtual void writeFileHeader(const label i);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("sizeDistribution");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        sizeDistribution
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~sizeDistribution();
+
+
+    // Member Functions
+
+        //- Return the reference to the construction dictionary
+        const dictionary& dict() const
+        {
+            return dict_;
+        }
+
+        //- Return the local list of cell IDs
+        const labelList& cellId() const
+        {
+            return cellId_;
+        }
+
+        //- Helper function to return the reference to the mesh
+        const fvMesh& mesh() const
+        {
+            return refCast<const fvMesh>(obr_);
+        }
+
+        //- Read from dictionary
+        virtual bool read(const dictionary& dict);
+
+        //- Execute
+        virtual bool execute();
+
+        //- Write
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -12,6 +12,28 @@ diameterModels/diameterModel/diameterModel.C
 diameterModels/diameterModel/newDiameterModel.C
 diameterModels/constantDiameter/constantDiameter.C
 diameterModels/isothermalDiameter/isothermalDiameter.C
+diameterModels/velocityGroup/velocityGroup.C
+
+populationBalanceModel/populationBalanceModel/populationBalanceModel.C
+diameterModels/velocityGroup/sizeGroup/sizeGroup.C
+
+populationBalanceModel/coalescenceModels/coalescenceModel/coalescenceModel.C
+populationBalanceModel/coalescenceModels/constantCoalescence/constantCoalescence.C
+populationBalanceModel/coalescenceModels/hydrodynamic/hydrodynamic.C
+
+populationBalanceModel/binaryBreakupModels/binaryBreakupModel/binaryBreakupModel.C
+populationBalanceModel/binaryBreakupModels/powerLawUniformBinary/powerLawUniformBinary.C
+
+populationBalanceModel/breakupModels/breakupModel/breakupModel.C
+populationBalanceModel/breakupModels/exponential/exponential.C
+populationBalanceModel/breakupModels/powerLaw/powerLaw.C
+
+populationBalanceModel/daughterSizeDistributionModels/daughterSizeDistributionModel/daughterSizeDistributionModel.C
+populationBalanceModel/daughterSizeDistributionModels/uniformBinaryDsd/uniformBinaryDsd.C
+
+populationBalanceModel/driftModels/driftModel/driftModel.C
+populationBalanceModel/driftModels/constantDrift/constantDrift.C
+populationBalanceModel/driftModels/densityChange/densityChange.C
 
 BlendedInterfacialModel/blendingMethods/blendingMethod/blendingMethod.C
 BlendedInterfacialModel/blendingMethods/blendingMethod/newBlendingMethod.C

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C

@@ -0,0 +1,378 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PopulationBalancePhaseSystem.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasePhaseSystem>
+Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::
+PopulationBalancePhaseSystem
+(
+    const fvMesh& mesh
+)
+:
+    BasePhaseSystem(mesh),
+
+    populationBalances_
+    (
+        this->lookup("populationBalances"),
+        diameterModels::populationBalanceModel::iNew(*this, pDmdt_)
+    )
+{
+    forAll(populationBalances_, i)
+    {
+        const Foam::diameterModels::populationBalanceModel& popBal =
+            populationBalances_[i];
+
+        forAllConstIter(phaseSystem::phasePairTable, popBal.phasePairs(), iter)
+        {
+            const phasePairKey& key = iter.key();
+
+            if (!this->phasePairs_.found(key))
+            {
+                this->phasePairs_.insert
+                (
+                    key,
+                    autoPtr<phasePair>
+                    (
+                        new phasePair
+                        (
+                            this->phaseModels_[key.first()],
+                            this->phaseModels_[key.second()]
+                        )
+                    )
+                );
+            }
+        }
+    }
+
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        // Initially assume no mass transfer
+        pDmdt_.insert
+        (
+            pair,
+            new volScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("pDmdt", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar("zero", dimDensity/dimTime, 0)
+            )
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class BasePhaseSystem>
+Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::
+~PopulationBalancePhaseSystem()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class BasePhaseSystem>
+Foam::tmp<Foam::volScalarField>
+Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::pDmdt
+(
+    const phasePairKey& key
+) const
+{
+    if (pDmdt_.found(key))
+    {
+        const scalar pDmdtSign
+        (
+            Pair<word>::compare(pDmdt_.find(key).key(), key)
+        );
+
+        return pDmdtSign**pDmdt_[key];
+    }
+
+    tmp<volScalarField> tpDmdt
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "dmdt",
+                this->mesh_.time().timeName(),
+                this->mesh_
+            ),
+            this->mesh_,
+            dimensionedScalar("zero", dimDensity/dimTime, 0)
+        )
+    );
+
+    return tpDmdt;
+}
+
+
+template<class BasePhaseSystem>
+Foam::tmp<Foam::volScalarField>
+Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::dmdt
+(
+    const phaseModel& phase
+) const
+{
+    tmp<volScalarField> tDmdtFromKey
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("dmdtFromKey", phase.name()),
+                this->mesh_.time().timeName(),
+                this->mesh_
+            ),
+            this->mesh_,
+            dimensionedScalar("zero", dimDensity/dimTime, 0)
+        )
+    );
+
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        const phaseModel* phase1 = &pair.phase1();
+        const phaseModel* phase2 = &pair.phase2();
+
+        forAllConstIter(phasePair, pair, iter)
+        {
+            if (phase1 == &phase)
+            {
+                tDmdtFromKey.ref() += this->dmdt
+                (
+                    phasePairKey(phase1->name(), phase2->name(),false)
+                );
+            }
+
+            Swap(phase1, phase2);
+        }
+    }
+
+    return tDmdtFromKey;
+}
+
+template<class BasePhaseSystem>
+Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
+Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::momentumTransfer() const
+{
+    autoPtr<phaseSystem::momentumTransferTable>
+        eqnsPtr(BasePhaseSystem::momentumTransfer());
+
+    phaseSystem::momentumTransferTable& eqns = eqnsPtr();
+
+    // Source term due to mass trasfer
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        const volVectorField& U1(pair.phase1().U());
+        const volVectorField& U2(pair.phase2().U());
+
+        const volScalarField dmdt(this->pDmdt(pair));
+        const volScalarField dmdt21(posPart(dmdt));
+        const volScalarField dmdt12(negPart(dmdt));
+
+        *eqns[pair.phase1().name()] += dmdt21*U2 - fvm::Sp(dmdt21, U1);
+        *eqns[pair.phase2().name()] -= dmdt12*U1 - fvm::Sp(dmdt12, U2);
+    }
+
+    return eqnsPtr;
+}
+
+
+template<class BasePhaseSystem>
+Foam::autoPtr<Foam::phaseSystem::massTransferTable>
+Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::heatTransfer() const
+{
+    autoPtr<phaseSystem::heatTransferTable> eqnsPtr =
+        BasePhaseSystem::heatTransfer();
+
+    phaseSystem::heatTransferTable& eqns = eqnsPtr();
+
+    // Source term due to mass trasfer
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        const volScalarField& he1(pair.phase1().thermo().he());
+        const volScalarField& he2(pair.phase2().thermo().he());
+
+        const volScalarField dmdt(this->pDmdt(pair));
+        const volScalarField dmdt21(posPart(dmdt));
+        const volScalarField dmdt12(negPart(dmdt));
+
+        *eqns[pair.phase1().name()] += dmdt21*he2 - fvm::Sp(dmdt21, he1);
+        *eqns[pair.phase2().name()] -= dmdt12*he1 - fvm::Sp(dmdt12, he2);
+    }
+
+    return eqnsPtr;
+}
+
+
+template<class BasePhaseSystem>
+Foam::autoPtr<Foam::phaseSystem::massTransferTable>
+Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::massTransfer() const
+{
+    autoPtr<phaseSystem::massTransferTable> eqnsPtr =
+        BasePhaseSystem::massTransfer();
+
+    phaseSystem::massTransferTable& eqns = eqnsPtr();
+
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        const phaseModel& phase = pair.phase1();
+        const phaseModel& otherPhase = pair.phase2();
+
+        const volScalarField dmdt(this->pDmdt(pair));
+        const volScalarField dmdt12(posPart(dmdt));
+        const volScalarField dmdt21(negPart(dmdt));
+
+        const PtrList<volScalarField>& Yi = phase.Y();
+
+        forAll(Yi, i)
+        {
+            const word name
+            (
+                IOobject::groupName(Yi[i].member(), phase.name())
+            );
+
+            const word otherName
+            (
+                IOobject::groupName(Yi[i].member(), otherPhase.name())
+            );
+
+            *eqns[name] +=
+                dmdt21*eqns[otherName]->psi()
+              - fvm::Sp(dmdt21, eqns[name]->psi());
+
+            *eqns[otherName] -=
+                dmdt12*eqns[name]->psi()
+              - fvm::Sp(dmdt12, eqns[otherName]->psi());
+        }
+    }
+
+    return eqnsPtr;
+}
+
+
+template<class BasePhaseSystem>
+bool Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::read()
+{
+    if (BasePhaseSystem::read())
+    {
+        bool readOK = true;
+
+        // Models ...
+
+        return readOK;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasePhaseSystem>
+void Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::solve()
+{
+    BasePhaseSystem::solve();
+
+    BasePhaseSystem::correct();
+
+    forAll(populationBalances_, i)
+    {
+        populationBalances_[i].solve();
+    }
+}
+
+
+// ************************************************************************* //1

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H

@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PopulationBalancePhaseSystem
+
+Description
+    Class which provides population balance functionality.
+
+SourceFiles
+    PopulationBalancePhaseSystem.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PopulationBalancePhaseSystem_H
+#define PopulationBalancePhaseSystem_H
+
+#include "HeatAndMassTransferPhaseSystem.H"
+#include "saturationModel.H"
+#include "Switch.H"
+
+#include "populationBalanceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class PopulationBalancePhaseSystem Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasePhaseSystem>
+class PopulationBalancePhaseSystem
+:
+    public BasePhaseSystem
+{
+
+protected:
+
+    // Protected data
+
+        //- populationBalanceModels
+        PtrList<diameterModels::populationBalanceModel> populationBalances_;
+
+        //- Interfacial Mass transfer rate
+        HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash> pDmdt_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh
+        PopulationBalancePhaseSystem(const fvMesh&);
+
+
+    //- Destructor
+    virtual ~PopulationBalancePhaseSystem();
+
+
+    // Member Functions
+
+        //- Return the interfacial mass flow rate
+        virtual tmp<volScalarField> pDmdt(const phasePairKey& key) const;
+
+        //- Return the interfacial mass flow rate
+        virtual tmp<volScalarField> dmdt(const phasePairKey& key) const
+        {
+            return BasePhaseSystem::dmdt(key) + this->pDmdt(key);
+        };
+
+        //- Return the total interfacial mass transfer rate for phase
+        virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;
+
+        //- Return the momentum transfer matrices
+        virtual autoPtr<phaseSystem::momentumTransferTable>
+            momentumTransfer() const;
+
+        //- Return the heat transfer matrices
+        virtual autoPtr<phaseSystem::heatTransferTable> heatTransfer() const;
+
+        //- Return the mass transfer matrices
+        virtual autoPtr<phaseSystem::massTransferTable> massTransfer() const;
+
+        //- Read base phaseProperties dictionary
+        virtual bool read();
+
+        //- Solve all population balance equations
+        virtual void solve();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "PopulationBalancePhaseSystem.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -54,7 +54,20 @@ Foam::diameterModels::isothermal::isothermal
 :
     diameterModel(diameterProperties, phase),
     d0_("d0", dimLength, diameterProperties_),
-    p0_("p0", dimPressure, diameterProperties_)
+    p0_("p0", dimPressure, diameterProperties_),
+    d_
+    (
+        IOobject
+        (
+            IOobject::groupName("d", phase.name()),
+            phase_.time().timeName(),
+            phase_.mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        phase_.mesh(),
+        dimensionedScalar("d", dimLength, 0.0)
+    )
 {}
 
 
@@ -77,6 +90,12 @@ Foam::tmp<Foam::volScalarField> Foam::diameterModels::isothermal::d() const
 }
 
 
+void Foam::diameterModels::isothermal::correct()
+{
+    d_ = d();
+}
+
+
 bool Foam::diameterModels::isothermal::read(const dictionary& phaseProperties)
 {
     diameterModel::read(phaseProperties);

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/isothermalDiameter/isothermalDiameter.H

@@ -60,6 +60,9 @@ class isothermal
         //- Reference pressure for the isothermal expansion
         dimensionedScalar p0_;
 
+        //- Actual diameter field
+        volScalarField d_;
+
 
 public:
 
@@ -86,6 +89,9 @@ public:
         //- Return the diameter field
         virtual tmp<volScalarField> d() const;
 
+        //- Correct the diameter field
+        virtual void correct();
+
         //- Read phaseProperties dictionary
         virtual bool read(const dictionary& phaseProperties);
 };

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

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sizeGroup.H"
+#include "populationBalanceModel.H"
+#include "mixedFvPatchField.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::sizeGroup::sizeGroup
+(
+    const word& name,
+    const dictionary& dict,
+    const phaseModel& phase,
+    const velocityGroup& velocityGroup,
+    const fvMesh& mesh
+)
+:
+    volScalarField
+    (
+        IOobject
+        (
+            IOobject::groupName
+            (
+                name,
+                IOobject::groupName
+                (
+                    velocityGroup.phase().name(),
+                    velocityGroup.popBalName()
+                )
+            ),
+            mesh.time().timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar(name, dimless, readScalar(dict.lookup("value"))),
+        velocityGroup.f().boundaryField().types()
+    ),
+    phase_(phase),
+    velocityGroup_(velocityGroup),
+    d_("d", dimLength, dict),
+    x_("x", velocityGroup.formFactor()*pow3(d_)),
+    value_(readScalar(dict.lookup("value")))
+{
+    // Adjust refValue at mixedFvPatchField boundaries
+    forAll(this->boundaryField(), patchi)
+    {
+        typedef mixedFvPatchField<scalar> mixedFvPatchScalarField;
+
+        if
+        (
+            isA<mixedFvPatchScalarField>(this->boundaryFieldRef()[patchi])
+        )
+        {
+            mixedFvPatchScalarField& f =
+                refCast<mixedFvPatchScalarField>
+                (
+                    this->boundaryFieldRef()[patchi]
+                );
+
+            f.refValue() = value_;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::diameterModels::sizeGroup::~sizeGroup()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::diameterModels::sizeGroup>
+Foam::diameterModels::sizeGroup::clone() const
+{
+    notImplemented("sizeGroup::clone() const");
+    return autoPtr<sizeGroup>(nullptr);
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.H

@@ -0,0 +1,185 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::sizeGroup
+
+Description
+    This class represents a single sizeGroup belonging to a velocityGroup.
+    The main property of a sizeGroup is its representative diameter. The
+    corresponding volScalarField f<number>.<phaseName>.<populationBalanceName>
+    gives the volume fraction of the sizeGroup such that all sizeGroup fractions
+    over a velocityGroup must sum to unity. The field is either read from the
+    startTime directory if present or constructed from a reference field called
+    f.<phaseName>.<populationBalanceName> where the boundary condition types
+    must be specified. All field and boundary condition values are reset to
+    match the "value" given in the sizeGroup dictionary.
+
+Usage
+    \table
+        Property     | Description             | Required    | Default value
+        d            | Representative diameter | yes         |
+        value        | Field and BC value      | yes         |
+    \endtable
+    Example
+    \verbatim
+    f1
+    {
+        d            3e-3;
+        value        1.0;
+    }
+    \endverbatim
+
+SourceFiles
+    sizeGroup.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sizeGroup_H
+#define sizeGroup_H
+
+#include "dictionaryEntry.H"
+#include "velocityGroup.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class sizeGroup Declaration
+\*---------------------------------------------------------------------------*/
+
+class sizeGroup
+:
+    public volScalarField
+{
+    // Private data
+
+        //- Phase this  sizeGroup belongs to
+        const phaseModel& phase_;
+
+        //- VelocityGroup this sizeGroup belongs to
+        const velocityGroup& velocityGroup_;
+
+        //- Representative diameter of the sizeGroup
+        const dimensionedScalar d_;
+
+        //- Representative volume of the sizeGroup
+        const dimensionedScalar x_;
+
+        //- Initial value and value at boundaries
+        const scalar value_;
+
+
+public:
+
+    // Constructors
+
+        sizeGroup
+        (
+            const word& name,
+            const dictionary& dict,
+            const phaseModel& phase,
+            const velocityGroup& velocityGroup,
+            const fvMesh& mesh
+        );
+
+        //- Return clone
+        autoPtr<sizeGroup> clone() const;
+
+        //- Return a pointer to a new sizeGroup created on freestore
+        //  from Istream
+        class iNew
+        {
+            const phaseModel& phase_;
+            const velocityGroup& velocityGroup_;
+
+        public:
+
+            iNew
+            (
+                const phaseModel& phase,
+                const velocityGroup& velocityGroup
+            )
+            :
+                phase_(phase),
+                velocityGroup_(velocityGroup)
+            {}
+
+            autoPtr<sizeGroup> operator()(Istream& is) const
+            {
+                dictionaryEntry ent(dictionary::null, is);
+                return autoPtr<sizeGroup>
+                (
+                    new sizeGroup
+                    (
+                        ent.keyword(),
+                        ent,
+                        phase_,
+                        velocityGroup_,
+                        phase_.mesh()
+                    )
+                );
+            }
+        };
+
+
+    //- Destructor
+    virtual ~sizeGroup();
+
+
+    // Member Functions
+
+        inline const word& keyword() const;
+
+        //- Return const-reference to the phase
+        inline const phaseModel&  phase() const;
+
+        //- Return const-reference to the velocityGroup
+        inline const velocityGroup& VelocityGroup() const;
+
+        //- Return representative diameter of the sizeGroup
+        inline const dimensionedScalar& d() const;
+
+        //- Return representative volume of the sizeGroup
+        inline const dimensionedScalar& x() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "sizeGroupI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroupI.H

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::word&
+Foam::diameterModels::sizeGroup::keyword() const
+{
+    return name();
+}
+
+
+inline const Foam::phaseModel&
+Foam::diameterModels::sizeGroup::phase() const
+{
+    return phase_;
+}
+
+
+inline const Foam::diameterModels::velocityGroup&
+Foam::diameterModels::sizeGroup::VelocityGroup() const
+{
+    return velocityGroup_;
+}
+
+
+inline const Foam::dimensionedScalar&
+Foam::diameterModels::sizeGroup::d() const
+{
+    return d_;
+}
+
+
+inline const Foam::dimensionedScalar&
+Foam::diameterModels::sizeGroup::x() const
+{
+    return x_;
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,401 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "velocityGroup.H"
+#include "sizeGroup.H"
+#include "populationBalanceModel.H"
+#include "addToRunTimeSelectionTable.H"
+#include "zeroGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+    defineTypeNameAndDebug(velocityGroup, 0);
+
+    addToRunTimeSelectionTable
+    (
+        diameterModel,
+        velocityGroup,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * 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(), fi.phase().residualAlpha())/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(), fi.phase().residualAlpha())/fi.x();
+    }
+
+    return tm3;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::velocityGroup::dsm() const
+{
+    return
+        max(min(phase_/m2_, sizeGroups_.last().d()), sizeGroups_.first().d());
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::diameterModels::velocityGroup::fSum() const
+{
+    tmp<volScalarField> tsumSizeGroups
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "sumSizeGroups",
+                phase_.time().timeName(),
+                phase_.mesh()
+            ),
+            phase_.mesh(),
+            dimensionedScalar("sumSizeGroups", dimless, 0)
+        )
+    );
+
+    volScalarField& sumSizeGroups = tsumSizeGroups.ref();
+
+    forAll(sizeGroups_,i)
+    {
+        sumSizeGroups += sizeGroups_[i];
+    }
+
+    Info<< phase_.name() << " sizeGroups-sum volume fraction, min, max = "
+        << sumSizeGroups.weightedAverage(phase_.mesh().V()).value()
+        << ' ' << min(sumSizeGroups).value()
+        << ' ' << max(sumSizeGroups).value()
+        << endl;
+
+    return tsumSizeGroups;
+}
+
+
+void Foam::diameterModels::velocityGroup::renormalize()
+{
+        Info<< phase_.name()
+            << " renormalizing sizeGroups"
+            << endl;
+
+        // Set negative values to zero
+        forAll(sizeGroups_, i)
+        {
+            sizeGroups_[i] *= pos(sizeGroups_[i]);
+        };
+
+        forAll(sizeGroups_, i)
+        {
+            sizeGroups_[i] /= fSum_;
+        };
+}
+
+
+Foam::tmp<Foam::fv::convectionScheme<Foam::scalar>>
+Foam::diameterModels::velocityGroup::mvconvection() const
+{
+    tmp<fv::convectionScheme<Foam::scalar> > mvConvection
+    (
+        fv::convectionScheme<Foam::scalar>::New
+        (
+            phase_.mesh(),
+            fields_,
+            phase_.alphaRhoPhi(),
+            phase_.mesh().divScheme
+            (
+                "div(" + phase_.alphaRhoPhi()().name() + ",f)"
+            )
+        )
+    );
+
+    return mvConvection;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::velocityGroup::velocityGroup
+(
+    const dictionary& diameterProperties,
+    const phaseModel& phase
+)
+:
+    diameterModel(diameterProperties, phase),
+    popBalName_(diameterProperties.lookup("populationBalance")),
+    f_
+    (
+        IOobject
+        (
+            IOobject::groupName
+            (
+                "f",
+                IOobject::groupName
+                (
+                    phase.name(),
+                    popBalName_
+                )
+            ),
+            phase.time().timeName(),
+            phase.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        phase.mesh()
+    ),
+    formFactor_("formFactor", dimless, diameterProperties.lookup("formFactor")),
+    sizeGroups_
+    (
+        diameterProperties.lookup("sizeGroups"),
+        sizeGroup::iNew(phase, *this)
+    ),
+    fSum_
+    (
+        IOobject
+        (
+            IOobject::groupName
+            (
+                "fsum",
+                IOobject::groupName
+                (
+                    phase.name(),
+                    popBalName_
+                )
+            ),
+            phase.time().timeName(),
+            phase.mesh()
+        ),
+        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
+        (
+            IOobject::groupName("d", phase.name()),
+            phase.time().timeName(),
+            phase.mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        phase.mesh(),
+        dimensionedScalar("d", dimLength, Zero)
+    ),
+    dmdt_
+    (
+        IOobject
+        (
+            IOobject::groupName("source", phase.name()),
+            phase.time().timeName(),
+            phase.mesh()
+        ),
+        phase.mesh(),
+        dimensionedScalar("dmdt", dimDensity/dimTime, Zero)
+    )
+{
+    if
+    (
+        mag
+        (
+            1.0
+          - fSum_.weightedAverage(fSum_.mesh().V()).value()
+        )
+      >= 1e-5
+     || mag(1.0 - max(fSum_).value()) >= 1e-5
+     || mag(1.0 - min(fSum_).value()) >= 1e-5
+    )
+    {
+        FatalErrorInFunction
+            << "Initial values of the sizeGroups belonging to velocityGroup "
+            << this->phase().name()
+            << " must add to unity." << endl
+            << exit(FatalError);
+    }
+
+    forAll(sizeGroups_, i)
+    {
+        fields_.add(sizeGroups_[i]);
+    }
+
+    m2_ = secondMoment();
+
+    m3_ = thirdMoment();
+
+    d_ = dsm();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::diameterModels::velocityGroup::~velocityGroup()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+void Foam::diameterModels::velocityGroup::preSolve()
+{
+    mvConvection_ = mvconvection();
+}
+
+
+void Foam::diameterModels::velocityGroup::postSolve()
+{
+    m2_ = secondMoment();
+
+    m3_ = thirdMoment();
+
+    d_ = dsm();
+
+    fSum_ = fSum();
+
+    if
+    (
+        phase_.mesh().solverDict(popBalName_).lookupOrDefault<Switch>
+        (
+            "renormalize",
+            false
+        )
+    )
+    {
+        renormalize();
+    }
+}
+
+
+bool Foam::diameterModels::velocityGroup::
+read(const dictionary& phaseProperties)
+{
+    diameterModel::read(phaseProperties);
+
+    return true;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::diameterModels::velocityGroup::d() const
+{
+    return d_;
+}
+
+// ************************************************************************* //

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

@@ -0,0 +1,221 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::velocityGroup
+
+Description
+    This diameterModel is intended for use with a populationBalanceModel in
+    order to simulate polydispersed bubbly or particulate flows. Can hold any
+    number of sizeGroups from which the Sauter mean diameter is calculated. It
+    can also be used as a diameterModel without a populationBalance and would
+    then behave like a constantDiameter model. In this case, some arbitrary name
+    must be entered for the populationBalance keyword.
+
+Usage
+    \table
+        Property          | Description
+        populationBalance | Name of the corresponding populationBalance
+        formFactor        | Form factor for converting diameter into volume
+        sizeGroups        | List of sizeGroups
+    \endtable
+    Example
+    \verbatim
+    diameterModel   velocityGroup;
+    velocityGroupCoeffs
+    {
+        populationBalance    bubbles;
+
+        formFactor      0.5235987756;
+
+        sizeGroups
+        (
+            f0{d  1.00e-3; value 0;}
+            f1{d  1.08e-3; value 0;}
+            f2{d  1.16e-3; value 0.25;}
+            f3{d  1.25e-3; value 0.5;}
+            f4{d  1.36e-3; value 0.25;}
+            f5{d  1.46e-3; value 0;}
+            ...
+        );
+    }
+    \endverbatim
+
+SourceFiles
+    velocityGroup.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef velocityGroup_H
+#define velocityGroup_H
+
+#include "diameterModel.H"
+#include "multivariateScheme.H"
+#include "convectionScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+
+class populationBalanceModelModel;
+class sizeGroup;
+
+/*---------------------------------------------------------------------------*\
+                        Class velocityGroup Declaration
+\*---------------------------------------------------------------------------*/
+
+class velocityGroup
+:
+    public diameterModel
+{
+    // Private data
+
+        //- Name of the populationBalance this velocityGroup belongs to
+        word popBalName_;
+
+        //- Reference field from which the sizeGroup fields are derived
+        volScalarField f_;
+
+        //- Form factor relating diameter and volume
+        dimensionedScalar formFactor_;
+
+        //- sizeGroups belonging to this velocityGroup
+        PtrList<sizeGroup> sizeGroups_;
+
+        //- 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_;
+
+        //- Multivariate convection scheme
+        tmp<fv::convectionScheme<scalar>> mvConvection_;
+
+        //- Table of fields for multivariate convection
+        multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields_;
+
+        //- Mass transfer rate
+        volScalarField dmdt_;
+
+
+    // Private member functions
+
+        tmp<volScalarField> secondMoment() const;
+
+        tmp<volScalarField> thirdMoment() const;
+
+        tmp<volScalarField> dsm() const;
+
+        tmp<volScalarField> fSum() const;
+
+        void renormalize();
+
+        tmp<Foam::fv::convectionScheme<Foam::scalar>> mvconvection() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("velocityGroup");
+
+
+    // Constructors
+
+        //- Construct from components
+        velocityGroup
+        (
+            const dictionary& diameterProperties,
+            const phaseModel& phase
+        );
+
+
+    //- Destructor
+    virtual ~velocityGroup();
+
+
+    // Member Functions
+
+        //- Return name of populationBalance this velocityGroup belongs to
+        inline const word& popBalName() const;
+
+        //- Return reference field for sizeGroup's
+        inline const volScalarField& f() const;
+
+        //- Return the form factor
+        inline const dimensionedScalar& formFactor() const;
+
+        //- 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;
+
+        //- Return const-reference to mass transfer rate
+        inline const volScalarField& dmdt() const;
+
+        //- Return reference to mass transfer rate
+        inline volScalarField& dmdt();
+
+        //- Corrections before populationBalanceModel::solve()
+        void preSolve();
+
+        //- Corrections after populationBalanceModel::solve()
+        void postSolve();
+
+        //- Read diameterProperties dictionary
+        virtual bool read(const dictionary& diameterProperties);
+
+        //- Return the Sauter-mean diameter
+        virtual tmp<volScalarField> d() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "velocityGroupI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroupI.H

@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::word&
+Foam::diameterModels::velocityGroup::popBalName() const
+{
+    return popBalName_;
+}
+
+
+inline const Foam::volScalarField&
+Foam::diameterModels::velocityGroup::f() const
+{
+    return f_;
+}
+
+
+inline const Foam::dimensionedScalar&
+Foam::diameterModels::velocityGroup::formFactor() const
+{
+    return formFactor_;
+}
+
+
+inline const Foam::PtrList<Foam::diameterModels::sizeGroup>&
+Foam::diameterModels::velocityGroup::sizeGroups() const
+{
+    return sizeGroups_;
+}
+
+
+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
+{
+    return mvConvection_;
+}
+
+
+inline const Foam::volScalarField& Foam::diameterModels::velocityGroup::
+dmdt() const
+{
+    return dmdt_;
+}
+
+
+inline Foam::volScalarField& Foam::diameterModels::velocityGroup::dmdt()
+{
+    return dmdt_;
+}
+
+
+// ************************************************************************* //

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

@@ -131,6 +131,12 @@ Foam::tmp<Foam::volScalarField> Foam::phaseModel::d() const
 }
 
 
+const Foam::autoPtr<Foam::diameterModel>& Foam::phaseModel::dPtr() const
+{
+    return diameterModel_;
+}
+
+
 void Foam::phaseModel::correct()
 {
     diameterModel_->correct();

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

@@ -180,6 +180,9 @@ public:
         //- Return the Sauter-mean diameter
         tmp<volScalarField> d() const;
 
+        //- Return const-reference to diameterModel of the phase
+        const autoPtr<diameterModel>& dPtr() const;
+
         //- Correct the phase properties
         virtual void correct();
 

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

@@ -283,6 +283,30 @@ void Foam::phaseSystem::solve()
 {}
 
 
+Foam::tmp<Foam::volScalarField> Foam::phaseSystem::dmdt
+(
+const Foam::phaseModel& phase
+) const
+{
+    tmp<volScalarField> tDmdt
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("dmdt", phase.name()),
+                this->mesh_.time().timeName(),
+                this->mesh_
+            ),
+            this->mesh_,
+            dimensionedScalar("zero", dimDensity/dimTime, 0)
+        )
+    );
+
+    return tDmdt;
+}
+
+
 void Foam::phaseSystem::correct()
 {
     forAll(phaseModels_, phasei)

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

@@ -352,6 +352,9 @@ public:
         //- Solve for the phase fractions
         virtual void solve();
 
+        //- Return the total interfacial mass transfer rate for a phase
+        virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;
+
         //- Correct the fluid properties other than the thermo and turbulence
         virtual void correct();
 

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

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "binaryBreakupModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+    defineTypeNameAndDebug(binaryBreakupModel, 0);
+    defineRunTimeSelectionTable(binaryBreakupModel, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Selector  * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::diameterModels::binaryBreakupModel>
+Foam::diameterModels::binaryBreakupModel::New
+(
+    const word& type,
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+{
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(type);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown binary breakup model type "
+            << type << nl << nl
+            << "Valid binary breakup model types : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<binaryBreakupModel>(cstrIter()(popBal, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::binaryBreakupModel::binaryBreakupModel
+(
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+:
+    popBal_(popBal),
+    C_("C", dimless, dict.lookupOrDefault("C", 1.0))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::diameterModels::populationBalanceModel&
+Foam::diameterModels::binaryBreakupModel::popBal() const
+{
+    return popBal_;
+}
+
+
+void Foam::diameterModels::binaryBreakupModel::correct()
+{}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::binaryBreakupModel
+
+Description
+    Base class for binary breakup models.
+
+SourceFiles
+    binaryBreakupModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef binaryBreakupModel_H
+#define binaryBreakupModel_H
+
+#include "populationBalanceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class binaryBreakupModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class binaryBreakupModel
+{
+protected:
+
+    // Protected data
+
+        //- Reference to the populationBalance
+        const populationBalanceModel& popBal_;
+
+        //- Optional coefficient
+        const dimensionedScalar C_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("binaryBreakupModel");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            binaryBreakupModel,
+            dictionary,
+            (
+                const populationBalanceModel& popBal,
+                const dictionary& dict
+            ),
+            (popBal, dict)
+        );
+
+
+    //- Class used for the read-construction of
+    //  PtrLists of binary breakup models
+    class iNew
+    {
+        const populationBalanceModel& popBal_;
+
+    public:
+
+        iNew(const populationBalanceModel& popBal)
+        :
+            popBal_(popBal)
+        {}
+
+        autoPtr<binaryBreakupModel> operator()(Istream& is) const
+        {
+            word type(is);
+            dictionary dict(is);
+            return binaryBreakupModel::New(type, popBal_, dict);
+        }
+    };
+
+
+    // Constructor
+
+        binaryBreakupModel
+        (
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+        autoPtr<binaryBreakupModel> clone() const
+        {
+            NotImplemented;
+            return autoPtr<binaryBreakupModel>(nullptr);
+        }
+
+
+    // Selector
+
+        static autoPtr<binaryBreakupModel> New
+        (
+            const word& type,
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~binaryBreakupModel()
+    {}
+
+
+    // Member Functions
+
+        //- Return reference to populationBalance
+        const populationBalanceModel& popBal() const;
+
+        //- Correct diameter independent expressions
+        virtual void correct();
+
+        //- Add to binary breakupRate
+        virtual void binaryBreakupRate
+        (
+            volScalarField& binaryBreakupRate,
+            const label i,
+            const label j
+        ) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "powerLawUniformBinary.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace binaryBreakupModels
+{
+    defineTypeNameAndDebug(powerLawUniformBinary, 0);
+    addToRunTimeSelectionTable
+    (
+        binaryBreakupModel,
+        powerLawUniformBinary,
+        dictionary
+    );
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::binaryBreakupModels::powerLawUniformBinary::
+powerLawUniformBinary
+(
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+:
+    binaryBreakupModel(popBal, dict),
+    power_(readScalar(dict.lookup("power")))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::diameterModels::binaryBreakupModels::powerLawUniformBinary::
+binaryBreakupRate
+(
+    volScalarField& binaryBreakupRate,
+    const label i,
+    const label j
+)
+{
+    const sizeGroup& fj = *popBal_.sizeGroups()[j];
+
+    binaryBreakupRate.primitiveFieldRef() +=
+        C_.value()
+       *pow(fj.x().value(), power_)
+       *2.0/fj.x().value();
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::binaryBreakupModels::powerLawUniformBinary
+
+Description
+    Powerlaw kernel with a uniform daughter size distribution. Used for
+    verification and validation of the binary breakup formulation implemented
+    in the populationBalanceModel class.
+
+SourceFiles
+    powerLawUniformBinary.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef powerLawUniformBinary_H
+#define powerLawUniformBinary_H
+
+#include "binaryBreakupModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace binaryBreakupModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class powerLawUniformBinary Declaration
+\*---------------------------------------------------------------------------*/
+
+class powerLawUniformBinary
+:
+    public binaryBreakupModel
+{
+    // Private data
+
+        //- Power
+        scalar power_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("powerLawUniformBinary");
+
+    // Constructor
+
+        powerLawUniformBinary
+        (
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~powerLawUniformBinary()
+    {}
+
+
+    // Member Functions
+
+        //- Add to binary breakupRate
+        virtual void binaryBreakupRate
+        (
+            volScalarField& binaryBreakupRate,
+            const label i,
+            const label j
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace binaryBreakupModels
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "breakupModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+    defineTypeNameAndDebug(breakupModel, 0);
+    defineRunTimeSelectionTable(breakupModel, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Selector  * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::diameterModels::breakupModel>
+Foam::diameterModels::breakupModel::New
+(
+    const word& type,
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+{
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(type);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown breakup model type "
+            << type << nl << nl
+            << "Valid breakup model types : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<breakupModel>(cstrIter()(popBal, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::breakupModel::breakupModel
+(
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+:
+    popBal_(popBal),
+    dict_(dict),
+    C_("C", dimless, dict.lookupOrDefault("C", 1.0))
+{
+    dsd_ = daughterSizeDistributionModel::New(*this, dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::diameterModels::populationBalanceModel&
+Foam::diameterModels::breakupModel::popBal() const
+{
+    return popBal_;
+}
+
+
+void Foam::diameterModels::breakupModel::correct()
+{}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,182 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::breakupModel
+
+Description
+    Base class for breakup models.
+
+SourceFiles
+    breakupModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef breakupModel_H
+#define breakupModel_H
+
+#include "populationBalanceModel.H"
+#include "daughterSizeDistributionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class breakupModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class breakupModel
+{
+protected:
+
+    // Protected data
+
+        //- Reference to the populationBalance
+        const populationBalanceModel& popBal_;
+
+        //- Dictionary
+        dictionary dict_;
+
+        //- Optional coefficient
+        const dimensionedScalar C_;
+
+        //- Daughter size distribution model
+        autoPtr<daughterSizeDistributionModel> dsd_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("breakupModel");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            breakupModel,
+            dictionary,
+            (
+                const populationBalanceModel& popBal,
+                const dictionary& dict
+            ),
+            (popBal, dict)
+        );
+
+
+    //- Class used for the read-construction of
+    //  PtrLists of breakup models
+    class iNew
+    {
+        const populationBalanceModel& popBal_;
+
+    public:
+
+        iNew(const populationBalanceModel& popBal)
+        :
+            popBal_(popBal)
+        {}
+
+        autoPtr<breakupModel> operator()(Istream& is) const
+        {
+            word type(is);
+            dictionary dict(is);
+            return breakupModel::New(type, popBal_, dict);
+        }
+    };
+
+
+    // Constructor
+
+        breakupModel
+        (
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+        autoPtr<breakupModel> clone() const
+        {
+            NotImplemented;
+            return autoPtr<breakupModel>(nullptr);
+        }
+
+
+    // Selector
+
+        static autoPtr<breakupModel> New
+        (
+            const word& type,
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~breakupModel()
+    {}
+
+
+    // Member Functions
+
+        //- Return reference to the populationBalance
+        const populationBalanceModel& popBal() const;
+
+        //- Return const reference to daughter size distribution pointer
+        const autoPtr<daughterSizeDistributionModel>& dsdPtr() const
+        {
+            return dsd_;
+        }
+
+        //- Return reference to daughter size distribution pointer
+        autoPtr<daughterSizeDistributionModel>& dsdPtr()
+        {
+            return dsd_;
+        }
+
+        //- Correct diameter independent expressions
+        virtual void correct();
+
+        //- Add to breakupRate
+        virtual void breakupRate
+        (
+            volScalarField& breakupRate,
+            const label i
+        ) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,71 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "exponential.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace breakupModels
+{
+    defineTypeNameAndDebug(exponential, 0);
+    addToRunTimeSelectionTable(breakupModel, exponential, dictionary);
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::breakupModels::exponential::exponential
+(
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+:
+    breakupModel(popBal, dict),
+    exponent_(readScalar(dict.lookup("exponent")))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::diameterModels::breakupModels::exponential::breakupRate
+(
+    volScalarField& breakupRate,
+    const label i
+)
+{
+    const sizeGroup& fi = *popBal_.sizeGroups()[i];
+
+    breakupRate.primitiveFieldRef() = C_.value()*exp(exponent_*fi.x().value());
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::breakupModels::exponential
+
+Description
+    Exponential kernel. Used for verification and validation of the breakup
+    formulation implemented in the populationBalanceModel class.
+
+SourceFiles
+    exponential.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef exponential_H
+#define exponential_H
+
+#include "breakupModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace breakupModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class exponential Declaration
+\*---------------------------------------------------------------------------*/
+
+class exponential
+:
+    public breakupModel
+{
+    // Private data
+
+        //- Power
+        scalar exponent_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("exponential");
+
+    // Constructor
+
+        exponential
+        (
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~exponential()
+    {}
+
+
+    // Member Functions
+
+        //- Add to total breakupRate
+        virtual void breakupRate
+        (
+            volScalarField& breakupRate,
+            const label i
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace breakupModels
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,71 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "powerLaw.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace breakupModels
+{
+    defineTypeNameAndDebug(powerLaw, 0);
+    addToRunTimeSelectionTable(breakupModel, powerLaw, dictionary);
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::breakupModels::powerLaw::powerLaw
+(
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+:
+    breakupModel(popBal, dict),
+    power_(readScalar(dict.lookup("power")))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::diameterModels::breakupModels::powerLaw::breakupRate
+(
+    volScalarField& breakupRate,
+    const label i
+)
+{
+    const sizeGroup& fi = *popBal_.sizeGroups()[i];
+
+    breakupRate.primitiveFieldRef() = C_.value()*pow(fi.x().value(), power_);
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::breakupModels::powerLaw
+
+Description
+    Powerlaw kernel. Used for verification and validation of the breakup
+    formulation implemented in the populationBalanceModel class.
+
+SourceFiles
+    powerLaw.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef powerLaw_H
+#define powerLaw_H
+
+#include "breakupModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace breakupModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class powerLaw Declaration
+\*---------------------------------------------------------------------------*/
+
+class powerLaw
+:
+    public breakupModel
+{
+    // Private data
+
+        //- Power
+        scalar power_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("powerLaw");
+
+    // Constructor
+
+        powerLaw
+        (
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~powerLaw()
+    {}
+
+
+    // Member Functions
+
+        //- Add to total breakupRate
+        virtual void breakupRate
+        (
+            volScalarField& breakupRate,
+            const label i
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace breakupModels
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "coalescenceModel.H"
+#include "phaseSystem.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+    defineTypeNameAndDebug(coalescenceModel, 0);
+    defineRunTimeSelectionTable(coalescenceModel, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Selector  * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::diameterModels::coalescenceModel>
+Foam::diameterModels::coalescenceModel::New
+(
+    const word& type,
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+{
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(type);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown coalescence model type "
+            << type << nl << nl
+            << "Valid coalescence model types : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<coalescenceModel>(cstrIter()(popBal, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::coalescenceModel::coalescenceModel
+(
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+:
+    popBal_(popBal),
+    C_("C", dimless, dict.lookupOrDefault("C", 1.0))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::diameterModels::populationBalanceModel&
+Foam::diameterModels::coalescenceModel::popBal() const
+{
+    return popBal_;
+}
+
+
+void Foam::diameterModels::coalescenceModel::correct()
+{}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::coalescenceModel
+
+Description
+    Base class for coalescence models.
+
+SourceFiles
+    coalescenceModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef coalescenceModel_H
+#define coalescenceModel_H
+
+#include "populationBalanceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class coalescenceModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class coalescenceModel
+{
+protected:
+
+    // Protected data
+
+        //- Reference to the populationBalance
+        const populationBalanceModel& popBal_;
+
+        //- Optional coefficient
+        const dimensionedScalar C_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("coalescenceModel");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            coalescenceModel,
+            dictionary,
+            (
+                const populationBalanceModel& popBal,
+                const dictionary& dict
+            ),
+            (popBal, dict)
+        );
+
+
+    //- Class used for the read-construction of
+    //  PtrLists of coalescence models
+    class iNew
+    {
+        const populationBalanceModel& popBal_;
+
+    public:
+
+        iNew(const populationBalanceModel& popBal)
+        :
+            popBal_(popBal)
+        {}
+
+        autoPtr<coalescenceModel> operator()(Istream& is) const
+        {
+            word type(is);
+            dictionary dict(is);
+            return coalescenceModel::New(type, popBal_, dict);
+        }
+    };
+
+
+    // Constructor
+
+        coalescenceModel
+        (
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+        autoPtr<coalescenceModel> clone() const
+        {
+            NotImplemented;
+            return autoPtr<coalescenceModel>(nullptr);
+        }
+
+
+    // Selector
+
+        static autoPtr<coalescenceModel> New
+        (
+            const word& type,
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~coalescenceModel()
+    {}
+
+
+    // Member Functions
+
+        //- Return reference to populationBalance
+        const populationBalanceModel& popBal() const;
+
+        //- Correct diameter independent expressions
+        virtual void correct();
+
+        //- Add to coalescenceRate
+        virtual void coalescenceRate
+        (
+            volScalarField& coalescenceRate,
+            const label i,
+            const label j
+        ) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "constantCoalescence.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace coalescenceModels
+{
+    defineTypeNameAndDebug(constantCoalescence, 0);
+    addToRunTimeSelectionTable
+    (
+        coalescenceModel,
+        constantCoalescence,
+        dictionary
+    );
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::coalescenceModels::constantCoalescence::
+constantCoalescence
+(
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+:
+    coalescenceModel(popBal, dict),
+    rate_("rate", dimVolume/dimTime, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::diameterModels::coalescenceModels::constantCoalescence::
+coalescenceRate
+(
+    volScalarField& coalescenceRate,
+    const label i,
+    const label j
+)
+{
+    coalescenceRate += C_.value()*rate_;
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::coalescenceModels::constantCoalescence
+
+Description
+    Constant coalescence kernel. Used for verification and validation of the
+    coalescence formulation implemented in the populationBalanceModel class.
+
+SourceFiles
+    constantCoalescence.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantCoalescence_H
+#define constantCoalescence_H
+
+#include "coalescenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace coalescenceModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class constantCoalescence Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantCoalescence
+:
+    public coalescenceModel
+{
+    // Private data
+
+        //- constant coalescenceRate
+        const dimensionedScalar rate_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constant");
+
+    // Constructor
+
+        constantCoalescence
+        (
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~constantCoalescence()
+    {}
+
+
+    // Member Functions
+
+        //- Add to coalescenceRate
+        virtual void coalescenceRate
+        (
+            volScalarField& coalescenceRate,
+            const label i,
+            const label j
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace coalescenceModels
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "hydrodynamic.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace coalescenceModels
+{
+    defineTypeNameAndDebug(hydrodynamic, 0);
+    addToRunTimeSelectionTable(coalescenceModel, hydrodynamic, dictionary);
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::coalescenceModels::hydrodynamic::hydrodynamic
+(
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+:
+    coalescenceModel(popBal, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::diameterModels::coalescenceModels::hydrodynamic::coalescenceRate
+(
+    volScalarField& coalescenceRate,
+    const label i,
+    const label j
+)
+{
+    const sizeGroup& fi = *popBal_.sizeGroups()[i];
+    const sizeGroup& fj = *popBal_.sizeGroups()[j];
+
+    coalescenceRate.primitiveFieldRef() +=
+        C_.value()*pow3((fi.d().value() + fj.d().value()));
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::coalescenceModels::hydrodynamic
+
+Description
+    Hydrodynamic kernel. Used for verification and validation of the coalescence
+    formulation implemented in the populationBalanceModel class.
+
+SourceFiles
+    hydrodynamic.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef hydrodynamic_H
+#define hydrodynamic_H
+
+#include "coalescenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace coalescenceModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class hydrodynamic Declaration
+\*---------------------------------------------------------------------------*/
+
+class hydrodynamic
+:
+    public coalescenceModel
+{
+public:
+
+    //- Runtime type information
+    TypeName("hydrodynamic");
+
+    // Constructor
+
+        hydrodynamic
+        (
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~hydrodynamic()
+    {}
+
+
+    // Member Functions
+
+        //- Add to coalescenceRate
+        virtual void coalescenceRate
+        (
+            volScalarField& coalescenceRate,
+            const label i,
+            const label j
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace coalescenceModels
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "daughterSizeDistributionModel.H"
+#include "breakupModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+    defineTypeNameAndDebug(daughterSizeDistributionModel, 0);
+    defineRunTimeSelectionTable(daughterSizeDistributionModel, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Selector  * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::diameterModels::daughterSizeDistributionModel>
+Foam::diameterModels::daughterSizeDistributionModel::New
+(
+    const breakupModel& breakup,
+    const dictionary& dict
+)
+{
+    word daughterSizeDistributionModelType
+    (
+        dict.lookup("daughterSizeDistributionModel")
+    );
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(daughterSizeDistributionModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown daughter size distribution model type "
+            << daughterSizeDistributionModelType << endl << endl
+            << "Valid daughter size distribution model types are : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return cstrIter()(breakup, dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+Foam::diameterModels::daughterSizeDistributionModel::
+daughterSizeDistributionModel
+(
+    const breakupModel& breakup,
+    const dictionary& dict
+)
+:
+    breakup_(breakup)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::diameterModels::daughterSizeDistributionModel::
+~daughterSizeDistributionModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::dimensionedScalar&
+Foam::diameterModels::daughterSizeDistributionModel::
+nik
+(
+    const label i,
+    const label k
+) const
+{
+    return nik_[k][i];
+}
+
+
+void Foam::diameterModels::daughterSizeDistributionModel::correct()
+{
+    if (nik_.size() == 0)
+    {
+        forAll(breakup_.popBal().sizeGroups(), k)
+        {
+            nik_.append(new PtrList<dimensionedScalar>());
+
+            for (label i = 0; i <= k; i++)
+            {
+                nik_[k].append(new dimensionedScalar (this->n(i, k)));
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::daughterSizeDistributionModel
+
+Description
+    Base class for daughter size distribution models.
+
+SourceFiles
+    daughterSizeDistributionModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef daughterSizeDistributionModel_H
+#define daughterSizeDistributionModel_H
+
+#include "dictionary.H"
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+
+class breakupModel;
+
+/*---------------------------------------------------------------------------*\
+                Class daughterSizeDistributionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class daughterSizeDistributionModel
+{
+protected:
+
+    // Protected data
+
+        //- Reference to the breakupModel
+        const breakupModel& breakup_;
+
+        //- nik
+        PtrList<PtrList<dimensionedScalar>> nik_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("daughterSizeDistributionModel");
+
+
+    // Declare runtime construction
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            daughterSizeDistributionModel,
+            dictionary,
+            (
+                const breakupModel& breakup,
+                const dictionary& dict
+            ),
+            (breakup, dict)
+        );
+
+
+    // Constructor
+
+        daughterSizeDistributionModel
+        (
+            const breakupModel& breakup,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~daughterSizeDistributionModel();
+
+
+    // Selectors
+
+        static autoPtr<daughterSizeDistributionModel> New
+        (
+            const breakupModel& breakup,
+            const dictionary& dict
+        );
+
+
+    // Member Functions
+
+        //- Return total number of particles assigned to class i when a particle
+        //  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
+        //  when a particle of class k breaks
+        virtual dimensionedScalar n(const label i, const label k) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "uniformBinaryDsd.H"
+#include "addToRunTimeSelectionTable.H"
+#include "breakupModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace daughterSizeDistributionModels
+{
+    defineTypeNameAndDebug(uniformBinaryDsd, 0);
+    addToRunTimeSelectionTable
+    (
+        daughterSizeDistributionModel,
+        uniformBinaryDsd,
+        dictionary
+    );
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::daughterSizeDistributionModels::uniformBinaryDsd::
+uniformBinaryDsd
+(
+    const breakupModel& breakup,
+    const dictionary& dict
+)
+:
+    daughterSizeDistributionModel(breakup, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::diameterModels::daughterSizeDistributionModels::uniformBinaryDsd::
+~uniformBinaryDsd()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::dimensionedScalar
+Foam::diameterModels::daughterSizeDistributionModels::uniformBinaryDsd::n
+(
+    const label i,
+    const label k
+) const
+{
+    const sizeGroup& fi = *breakup_.popBal().sizeGroups()[i];
+    const sizeGroup& fk = *breakup_.popBal().sizeGroups()[k];
+
+    const dimensionedScalar& xi = fi.x();
+    const dimensionedScalar& xk = fk.x();
+    dimensionedScalar x("x", dimVolume, Zero);
+    dimensionedScalar xii("xii", dimVolume, Zero);
+
+    if (i == 0)
+    {
+       x = xi;
+    }
+    else
+    {
+       x = breakup_.popBal().sizeGroups()[i-1]->x();
+    }
+
+    if (i == breakup_.popBal().sizeGroups().size() - 1)
+    {
+        xii = fi.x();
+    }
+    else
+    {
+        xii = breakup_.popBal().sizeGroups()[i+1]->x();
+    }
+
+    if (i == 0)
+    {
+        return
+            1.0/(xii - xi)
+           *(
+                (xii*xii - 0.5*sqr(xii))
+              - (xi*xii - 0.5*sqr(xi))
+            )
+           *2.0/xk;
+    }
+    else if (i == k)
+    {
+        return
+            1.0/(xi - x)
+           *(
+                (0.5*sqr(xi) - xi*x)
+              - (0.5*sqr(x) - x*x)
+            )
+           *2.0/xk;
+    }
+    else
+    {
+        return
+            1.0/(xii - xi)
+           *(
+                (xii*xii - 0.5*sqr(xii))
+              - (xi*xii - 0.5*sqr(xi))
+            )
+           *2.0/xk
+          + 1.0/(xi - x)
+           *(
+                (0.5*sqr(xi) - xi*x)
+              - (0.5*sqr(x) - x*x)
+            )
+           *2.0/xk;
+    }
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::daughterSizeDistributionModels::uniformBinaryDsd
+
+Description
+    Uniform binary daughter size distribution. Used for verification and
+    validation of the breakup formulation implemented in the
+    populationBalanceModel class.
+
+SourceFiles
+    uniformBinaryDsd.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef uniformBinaryDsd_H
+#define uniformBinaryDsd_H
+
+#include "daughterSizeDistributionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace daughterSizeDistributionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class uniformBinaryDsd Declaration
+\*---------------------------------------------------------------------------*/
+
+class uniformBinaryDsd
+:
+    public daughterSizeDistributionModel
+{
+public:
+
+    //- Runtime type information
+    TypeName("uniformBinary");
+
+
+    // Constructor
+
+        uniformBinaryDsd
+        (
+            const breakupModel& breakup,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~uniformBinaryDsd();
+
+
+    // 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;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace daughterSizeDistributionModels
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "constantDrift.H"
+#include "phaseSystem.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace driftModels
+{
+    defineTypeNameAndDebug(constantDrift, 0);
+    addToRunTimeSelectionTable(driftModel, constantDrift, dictionary);
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::driftModels::constantDrift::constantDrift
+(
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+:
+    driftModel(popBal, dict),
+    calculateFromFvOptions_(dict.lookup("calculateFromFvOptions")),
+    rate_("rate", dimVolume/dimTime, dict.lookupOrDefault<scalar>("rate", 1.0)),
+    N_
+    (
+        IOobject
+        (
+            "N",
+            popBal.mesh().time().timeName(),
+            popBal.mesh()
+        ),
+        popBal.mesh(),
+        dimensionedScalar("Sui", inv(dimVolume), Zero)
+    )
+{
+    if (calculateFromFvOptions_ == false)
+    {
+        rate_.read(dict);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::diameterModels::driftModels::constantDrift::correct()
+{
+    N_ *= 0.0;
+
+    forAll(popBal_.sizeGroups(), i)
+    {
+        const sizeGroup& fi = *popBal_.sizeGroups()[i];
+
+        N_ += fi*fi.phase()/fi.x();
+    }
+}
+
+
+void Foam::diameterModels::driftModels::constantDrift::driftRate
+(
+    volScalarField& driftRate,
+    const label i
+)
+{
+    const sizeGroup& fi = *popBal_.sizeGroups()[i];
+    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_;
+    }
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+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.
+
+SourceFiles
+    constant.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantDrift_H
+#define constantDrift_H
+
+#include "driftModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace driftModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class constantDrift Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantDrift
+:
+    public driftModel
+{
+    // Private data
+
+        //- If true, calculate drift rate from fvOptions mass source
+        Switch calculateFromFvOptions_;
+
+        //- Optional driftRate
+        dimensionedScalar rate_;
+
+        //- Total number concentration
+        volScalarField N_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constant");
+
+    // Constructor
+
+        constantDrift
+        (
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~constantDrift()
+    {}
+
+
+    // Member Functions
+
+        //- Correct diameter independent expressions
+        virtual void correct();
+
+        //- Add to driftRate
+        virtual void driftRate
+        (
+            volScalarField& driftRate,
+            const label i
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace driftModels
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "densityChange.H"
+#include "addToRunTimeSelectionTable.H"
+#include "phaseSystem.H"
+#include "fvcDdt.H"
+#include "fvcGrad.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace driftModels
+{
+    defineTypeNameAndDebug(densityChangeDrift, 0);
+    addToRunTimeSelectionTable(driftModel, densityChangeDrift, dictionary);
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::driftModels::densityChangeDrift::densityChangeDrift
+(
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+:
+    driftModel(popBal, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+void Foam::diameterModels::driftModels::densityChangeDrift::driftRate
+(
+    volScalarField& driftRate,
+    const label i
+)
+{
+    const sizeGroup& fi = *popBal_.sizeGroups()[i];
+    volScalarField& rho = const_cast<volScalarField&>(fi.phase().rho()());
+
+    driftRate -= (fvc::ddt(rho) + (fvc::grad(rho)&popBal_.U()))
+       *popBal_.sizeGroups()[i]->x()/rho;
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::driftModels::densityChange
+
+Description
+    Drift rate induced by changes in density.
+
+SourceFiles
+    densityChange.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef densityChange_H
+#define densityChange_H
+
+#include "driftModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace driftModels
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class densityChangeDrift Declaration
+\*---------------------------------------------------------------------------*/
+
+class densityChangeDrift
+:
+    public driftModel
+{
+public:
+
+    //- Runtime type information
+    TypeName("densityChange");
+
+    // Constructor
+
+        densityChangeDrift
+        (
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~densityChangeDrift()
+    {}
+
+
+    // Member Functions
+
+        //- Add to driftRate
+        virtual void driftRate
+        (
+            volScalarField& driftRate,
+            const label i
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace driftModels
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "driftModel.H"
+#include "phaseSystem.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+    defineTypeNameAndDebug(driftModel, 0);
+    defineRunTimeSelectionTable(driftModel, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Selector  * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::diameterModels::driftModel>
+Foam::diameterModels::driftModel::New
+(
+    const word& type,
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+{
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(type);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown drift model type "
+            << type << nl << nl
+            << "Valid drift model types : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<driftModel>(cstrIter()(popBal, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::driftModel::driftModel
+(
+    const populationBalanceModel& popBal,
+    const dictionary& dict
+)
+:
+    popBal_(popBal)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::diameterModels::populationBalanceModel&
+Foam::diameterModels::driftModel::popBal() const
+{
+    return popBal_;
+}
+
+
+void Foam::diameterModels::driftModel::correct()
+{}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::driftModel
+
+Description
+    Base class for drift models.
+
+SourceFiles
+    driftModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef driftModel_H
+#define driftModel_H
+
+#include "populationBalanceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class driftModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class driftModel
+{
+protected:
+
+    // Protected data
+
+        //- Reference to the populationBalanceModel
+        const populationBalanceModel& popBal_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("driftModel");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            driftModel,
+            dictionary,
+            (
+                const populationBalanceModel& popBal,
+                const dictionary& dict
+            ),
+            (popBal, dict)
+        );
+
+
+    //- Class used for the read-construction of
+    //  PtrLists of drift models
+    class iNew
+    {
+        const populationBalanceModel& popBal_;
+
+    public:
+
+        iNew(const populationBalanceModel& popBal)
+        :
+            popBal_(popBal)
+        {}
+
+        autoPtr<driftModel> operator()(Istream& is) const
+        {
+            word type(is);
+            dictionary dict(is);
+            return driftModel::New(type, popBal_, dict);
+        }
+    };
+
+
+    // Constructor
+
+        driftModel
+        (
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+        autoPtr<driftModel> clone() const
+        {
+            NotImplemented;
+            return autoPtr<driftModel>(nullptr);
+        }
+
+
+    // Selector
+
+        static autoPtr<driftModel> New
+        (
+            const word& type,
+            const populationBalanceModel& popBal,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~driftModel()
+    {}
+
+
+    // Member Functions
+
+        //- Return reference to populationBalanceModel
+        const populationBalanceModel& popBal() const;
+
+        //- Correct diameter independent expressions
+        virtual void correct();
+
+        //- Add to driftRate
+        virtual void driftRate
+        (
+            volScalarField& driftRate,
+            const label i
+        ) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,403 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::populationBalance
+
+Description
+    Class that solves the univariate population balance equation by means of
+    a class method (also called sectional or discrete method). The internal
+    coordinate is set to the particle volume, so the equation is based on
+    a transport equation of the volume-based number density function. The
+    discretization is done using the fixed pivot technique of Kumar and
+    Ramkrishna (1996). The source terms for coalescence and breakup are written
+    in a way that particle number and mass are preserved. An additional breakup
+    formulation that is restricted to binary breakup is also included which
+    allows a more efficient handling of breakup models with an integrated
+    daughter size distribution. It is given by Liao et al. (2017). The
+    implementation allows to split the population balance over multiple velocity
+    fields using the capability of reactingMultiphaseEulerFoam to solve for n
+    momentum equations. It is also possible to define multiple population
+    balances, e.g. bubbles and droplets simultaneously.
+
+    References:
+    \verbatim
+        Coalescence and breakup term formulation:
+            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.
+    \endverbatim
+
+    \verbatim
+        Binary breakup term formulation:
+            Liao, Y.; Oertel, R.; Kriebitzsch, S.; Schlegel, F.; Lucas, D.
+            (2017). A Discrete Population Balance Equation for Binary Breakage.
+            International Journal for Numerical Methods in Fluids, In review.
+    \endverbatim
+
+Usage
+    Example excerpt from a phaseProperties dictionary.
+    \verbatim
+    type populationBalanceTwoPhaseSystem;
+
+    phases (air water);
+
+    populationBalances (bubbles);
+
+    air
+    {
+        type            purePhaseModel;
+        diameterModel   velocityGroup;
+        velocityGroupCoeffs
+        {
+            populationBalance    bubbles;
+
+            formFactor      0.5235987756;
+
+            sizeGroups
+            (
+                f0{d  1.00e-3; value 0;}
+                f1{d  1.08e-3; value 0;}
+                f2{d  1.16e-3; value 0.25;}
+                f3{d  1.25e-3; value 0.5;}
+                f4{d  1.36e-3; value 0.25;}
+                f5{d  1.46e-3; value 0;}
+                ...
+            );
+        }
+
+        residualAlpha   1e-6;
+    }
+
+    populationBalanceCoeffs
+    {
+        bubbles
+        {
+            continuousPhase water;
+
+            coalescenceModels
+            (
+                hydrodynamic
+                {
+                    C 0.25;
+                }
+            );
+
+            binaryBreakupModels
+            ();
+
+            breakupModels
+            (
+                exponential
+                {
+                    C 0.5;
+                    exponent 0.01;
+                    daughterSizeDistributionModel uniformBinary;
+                }
+            );
+
+            driftModels
+            (
+                densityChange{}
+            );
+        }
+    }
+    \endverbatim
+
+SourceFiles
+    populationBalanceModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef populationBalanceModel_H
+#define populationBalanceModel_H
+
+#include "sizeGroup.H"
+#include "phasePair.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class phaseSystem;
+
+namespace diameterModels
+{
+
+class coalescenceModel;
+class breakupModel;
+class binaryBreakupModel;
+class driftModel;
+
+/*---------------------------------------------------------------------------*\
+                   Class populationBalanceModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class populationBalanceModel
+:
+    public regIOobject
+{
+    // Private typedefs
+
+        typedef
+            HashTable<autoPtr<phasePair>, phasePairKey, phasePairKey::hash>
+            phasePairTable;
+
+
+    // Private data
+
+        //- Reference to the phaseSystem
+        const phaseSystem& fluid_;
+
+        //- Interfacial Mass transfer rate between velocityGroups
+        HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>& pDmdt_;
+
+        //- Reference to the mesh
+        const fvMesh& mesh_;
+
+        //- Name of the populationBalance
+        word name_;
+
+        //- Dictionary
+        dictionary dict_;
+
+        //- Continuous phase
+        const phaseModel& continuousPhase_;
+
+        //- velocityGroups belonging to this populationBalance
+        List<velocityGroup*> velocityGroups_;
+
+        //- sizeGroups belonging to this populationBalance
+        List<sizeGroup*> sizeGroups_;
+
+        //- List of unordered phasePairs in this populationBalance
+        phasePairTable phasePairs_;
+
+        //- sizeGroup boundaries
+        PtrList<dimensionedScalar> v_;
+
+        //- Section width required for binary breakup formulation
+        PtrList<PtrList<dimensionedScalar>> delta_;
+
+        //- Explicitly treated sources
+        PtrList<volScalarField> Su_;
+
+        //- Sources treated implicitly or explicitly depending on sign
+        PtrList<volScalarField> SuSp_;
+
+        //- Field for caching pDmdt sources
+        autoPtr<volScalarField> Sui_;
+
+        //- Coalescence models
+        PtrList<coalescenceModel> coalescence_;
+
+        //- Coalescence rate
+        autoPtr<volScalarField> coalescenceRate_;
+
+        //- BreakupModels
+        PtrList<breakupModel> breakup_;
+
+        //- Breakup rate
+        autoPtr<volScalarField> breakupRate_;
+
+        //- Binary breakup models
+        PtrList<binaryBreakupModel> binaryBreakup_;
+
+        //- Binary breakup rate
+        autoPtr<volScalarField> binaryBreakupRate_;
+
+        //- Drift models
+        PtrList<driftModel> drift_;
+
+        //- Drift rate
+        autoPtr<volScalarField> driftRate_;
+
+        //- Ratio between sucessive reprasentative volumes
+        autoPtr<volScalarField> r_;
+
+        //- Total void fraction of phases belonging to this populationBalance
+        volScalarField alphas_;
+
+        //- Average velocity
+        volVectorField U_;
+
+        //- Mean Sauter diameter
+        volScalarField d_;
+
+
+    // Private member functions
+
+        void registerVelocityGroups();
+
+        void add(sizeGroup* group);
+
+        void createPhasePairs();
+
+        void preSolve();
+
+        dimensionedScalar gamma
+        (
+            const label i,
+            const dimensionedScalar& v
+        ) const;
+
+        void birthByCoalescence(const label j, const label k);
+
+        void deathByCoalescence(const label i, const label j);
+
+        void birthByBreakup(const label k, const label model);
+
+        void deathByBreakup(const label i);
+
+        void calcDeltas();
+
+        void birthByBinaryBreakup(const label i,const label j);
+
+        void deathByBinaryBreakup(const label j, const label i);
+
+        void drift(const label i);
+
+        void sources();
+
+        void dmdt();
+
+        void calcAlphas();
+
+        void calcVelocity();
+
+        tmp<volScalarField> dsm() const;
+
+
+public:
+
+    friend class sizeGroup;
+    friend class velocityGroup;
+
+    // Constructor
+
+        populationBalanceModel
+        (
+            const phaseSystem& fluid,
+            const word& name,
+            HashPtrTable
+            <
+                volScalarField,
+                phasePairKey,
+                phasePairKey::hash
+            >& pDmdt
+        );
+
+        //- Return clone
+        autoPtr<populationBalanceModel> clone() const;
+
+        //- Return a pointer to a new populationBalanceModel object created on
+        //  freestore from Istream
+        class iNew
+        {
+            const phaseSystem& fluid_;
+
+            HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>&
+                pDmdt_;
+
+        public:
+
+            iNew
+            (
+                const phaseSystem& fluid,
+                HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>&
+                    pDmdt
+            )
+            :
+                fluid_(fluid),
+                pDmdt_(pDmdt)
+            {}
+
+            autoPtr<populationBalanceModel> operator()(Istream& is) const
+            {
+                return autoPtr<populationBalanceModel>
+                (
+                    new populationBalanceModel(fluid_, word(is), pDmdt_)
+                );
+            }
+        };
+
+
+    //- Destructor
+    virtual ~populationBalanceModel();
+
+    // Member Functions
+
+        //- Dummy write for regIOobject
+        bool writeData(Ostream&) const;
+
+        //- Return reference to the phaseSystem
+        inline const phaseSystem& fluid() const;
+
+        //- Return reference to the mesh
+        inline const fvMesh& mesh() const;
+
+        //- Return populationBalanceCoeffs dictionary
+        inline const dictionary& dict() const;
+
+        //- Return continuous phase
+        inline const phaseModel& continuousPhase() const;
+
+        //- Return the velocityGroups belonging to this populationBalance
+        inline const List<velocityGroup*>& velocityGroups() const;
+
+        //- Return the sizeGroups belonging to this populationBalance
+        inline const List<sizeGroup*>& sizeGroups() const;
+
+        //- Return list of unordered phasePairs in this populationBalance
+        inline const phasePairTable& phasePairs() const;
+
+        //- Return the sizeGroup boundaries
+        inline const PtrList<dimensionedScalar>& v() const;
+
+        //- Return total void of phases belonging to this populationBalance
+        inline const volScalarField& alphas() const;
+
+        //- Return average velocity
+        inline const volVectorField& U() const;
+
+        //- Solve the population balance equation
+        void solve();
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "populationBalanceModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::phaseSystem&
+Foam::diameterModels::populationBalanceModel::fluid() const
+{
+    return fluid_;
+}
+
+
+inline const Foam::fvMesh&
+Foam::diameterModels::populationBalanceModel::mesh() const
+{
+    return mesh_;
+}
+
+
+inline const Foam::dictionary&
+Foam::diameterModels::populationBalanceModel::dict() const
+{
+    return dict_;
+}
+
+
+inline const Foam::phaseModel&
+Foam::diameterModels::populationBalanceModel::continuousPhase() const
+{
+    return continuousPhase_;
+}
+
+
+inline const Foam::List<Foam::diameterModels::velocityGroup*>&
+Foam::diameterModels::populationBalanceModel::velocityGroups() const
+{
+    return velocityGroups_;
+}
+
+
+inline const Foam::List<Foam::diameterModels::sizeGroup*>&
+Foam::diameterModels::populationBalanceModel::sizeGroups() const
+{
+    return sizeGroups_;
+}
+
+
+inline const Foam::diameterModels::populationBalanceModel::phasePairTable&
+Foam::diameterModels::populationBalanceModel::phasePairs() const
+{
+    return phasePairs_;
+}
+
+
+inline const Foam::PtrList<Foam::dimensionedScalar>&
+Foam::diameterModels::populationBalanceModel::v() const
+{
+    return v_;
+}
+
+
+inline const Foam::volScalarField&
+Foam::diameterModels::populationBalanceModel::alphas() const
+{
+    return alphas_;
+}
+
+
+inline const Foam::volVectorField&
+Foam::diameterModels::populationBalanceModel::U() const
+{
+    return U_;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystems.C

@@ -30,6 +30,7 @@ License
 #include "MomentumTransferPhaseSystem.H"
 #include "HeatTransferPhaseSystem.H"
 #include "HeatAndMassTransferPhaseSystem.H"
+#include "PopulationBalancePhaseSystem.H"
 #include "InterfaceCompositionPhaseChangePhaseSystem.H"
 #include "ThermalPhaseChangePhaseSystem.H"
 
@@ -102,6 +103,29 @@ namespace Foam
         dictionary,
         thermalPhaseChangeMultiphaseSystem
     );
+
+    typedef
+
+        PopulationBalancePhaseSystem
+        <
+            HeatAndMassTransferPhaseSystem
+            <
+                MomentumTransferPhaseSystem<multiphaseSystem>
+            >
+        >
+        populationBalanceMultiphaseSystem;
+
+    addNamedToRunTimeSelectionTable
+    (
+        multiphaseSystem,
+        populationBalanceMultiphaseSystem,
+        dictionary,
+        populationBalanceMultiphaseSystem
+    );
+
+
+
+
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/phaseChange/phaseChange.C

@@ -56,7 +56,7 @@ Foam::diameterModels::IATEsources::phaseChange::phaseChange
 :
     IATEsource(iate),
     pairName_(dict.lookup("pairName")),
-    iDmdtPtr_(nullPtr)
+    iDmdtPtr_(nullptr)
 {}
 
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystems.C

@@ -30,6 +30,7 @@ License
 #include "MomentumTransferPhaseSystem.H"
 #include "HeatTransferPhaseSystem.H"
 #include "HeatAndMassTransferPhaseSystem.H"
+#include "PopulationBalancePhaseSystem.H"
 #include "InterfaceCompositionPhaseChangePhaseSystem.H"
 #include "ThermalPhaseChangePhaseSystem.H"
 
@@ -88,10 +89,10 @@ namespace Foam
     typedef
         ThermalPhaseChangePhaseSystem
         <
-            HeatAndMassTransferPhaseSystem
+        HeatAndMassTransferPhaseSystem
-            <
+        <
-                MomentumTransferPhaseSystem<twoPhaseSystem>
+            MomentumTransferPhaseSystem<twoPhaseSystem>
-            >
+        >
         >
         thermalPhaseChangeTwoPhaseSystem;
 
@@ -102,6 +103,24 @@ namespace Foam
         dictionary,
         thermalPhaseChangeTwoPhaseSystem
     );
+
+    typedef
+        PopulationBalancePhaseSystem
+        <
+            HeatAndMassTransferPhaseSystem
+            <
+                MomentumTransferPhaseSystem<twoPhaseSystem>
+            >
+        >
+        populationBalanceTwoPhaseSystem;
+
+    addNamedToRunTimeSelectionTable
+    (
+        twoPhaseSystem,
+        populationBalanceTwoPhaseSystem,
+        dictionary,
+        populationBalanceTwoPhaseSystem
+    );
 }
 
 

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/T.airI

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T.airI;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 300;
+
+boundaryField
+{
+    walls
+    {
+        type               zeroGradient;
+    }
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.airI;
+        inletValue         $internalField;
+        value              $internalField;
+    }
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/T.airII

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T.airII;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 300;
+
+boundaryField
+{
+    walls
+    {
+        type               zeroGradient;
+    }
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.airII;
+        inletValue         $internalField;
+        value              $internalField;
+    }
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/T.water

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 350;
+
+boundaryField
+{
+    walls
+    {
+        type               zeroGradient;
+    }
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.water;
+        inletValue         $internalField;
+        value              $internalField;
+    }
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/Theta

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Theta;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 1.0e-7;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1.0e-7;
+        value           uniform 1.0e-7;
+    }
+
+    walls
+    {
+        type            zeroGradient;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/U.airI

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.airI;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0.1 0);
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               pressureInletOutletVelocity;
+        phi                phi.airI;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedValue;
+        value              uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/U.airII

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.airII;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0.1 0);
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               pressureInletOutletVelocity;
+        phi                phi.airII;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedValue;
+        value              uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/U.water

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               pressureInletOutletVelocity;
+        phi                phi.water;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedValue;
+        value              uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/alpha.airI

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.airI;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.25;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.airI;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/alpha.airI.orig

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.airI;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.25;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.airI;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/alpha.airII.orig

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.airII;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.25;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.airII;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/alpha.water.orig

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.5;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.water;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/p

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [1 -1 -2 0 0 0 0];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type               calculated;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               calculated;
+        value              $internalField;
+    }
+    walls
+    {
+        type               calculated;
+        value              $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/0/p_rgh

@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [1 -1 -2 0 0 0 0];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            prghTotalPressure;
+        p0              $internalField;
+        U               U.airII;
+        phi             phi.airII;
+        rho             thermo:rho.airII;
+        value           $internalField;
+    }
+    walls
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/Allrun

@@ -0,0 +1,14 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# Set application name
+application=`getApplication`
+
+runApplication blockMesh
+runApplication setFields
+runApplication $application
+
+# ----------------------------------------------------------------- end-of-file

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/constant/g

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 -9.81 0);
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/constant/phaseProperties

@@ -0,0 +1,307 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      phaseProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+type heatAndMomentumTransferMultiphaseSystem;
+
+phases (airI airII water);
+
+airI
+{
+    type            purePhaseModel;
+    diameterModel   isothermal;
+    isothermalCoeffs
+    {
+        d0              3e-3;
+        p0              1e5;
+    }
+
+    residualAlpha   1e-6;
+}
+
+airII
+{
+    type            purePhaseModel;
+    diameterModel   isothermal;
+    isothermalCoeffs
+    {
+        d0              5e-3;
+        p0              1e5;
+    }
+
+    residualAlpha   1e-6;
+}
+
+water
+{
+    type            purePhaseModel;
+    diameterModel   constant;
+    constantCoeffs
+    {
+        d               1e-4;
+    }
+
+    residualAlpha   1e-6;
+}
+
+blending
+{
+    default
+    {
+        type            linear;
+        minFullyContinuousAlpha.airI 0.7;
+        minPartlyContinuousAlpha.airI 0.5;
+        minFullyContinuousAlpha.airII 0.7;
+        minPartlyContinuousAlpha.airII 0.5;
+        minFullyContinuousAlpha.water 0.7;
+        minPartlyContinuousAlpha.water 0.5;
+    }
+}
+
+surfaceTension
+(
+    (airI and water)
+    {
+        type            constant;
+        sigma           0.07;
+    }
+
+    (airII and water)
+    {
+        type            constant;
+        sigma           0.07;
+    }
+);
+
+interfaceCompression
+();
+
+aspectRatio
+(
+    (airI in water)
+    {
+        type            constant;
+        E0              1.0;
+    }
+
+    (water in airI)
+    {
+        type            constant;
+        E0              1.0;
+    }
+
+    (airII in water)
+    {
+        type            constant;
+        E0              1.0;
+    }
+
+    (water in airII)
+    {
+        type            constant;
+        E0              1.0;
+    }
+);
+
+drag
+(
+    (airI in water)
+    {
+        type            SchillerNaumann;
+        residualRe      1e-3;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+
+    (water in airI)
+    {
+        type            SchillerNaumann;
+        residualRe      1e-3;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+
+    (airI and water)
+    {
+        type            segregated;
+        m               0.5;
+        n               8;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+
+    (airII in water)
+    {
+        type            SchillerNaumann;
+        residualRe      1e-3;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+
+    (water in airII)
+    {
+        type            SchillerNaumann;
+        residualRe      1e-3;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+
+    (airII and water)
+    {
+        type            segregated;
+        m               0.5;
+        n               8;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+
+    (airI in airII)
+    {
+        type            SchillerNaumann;
+        residualRe      1e-3;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+
+    (airII in airI)
+    {
+        type            SchillerNaumann;
+
+        residualRe      1e-3;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+
+    (airI and airII)
+    {
+        type            segregated;
+        m               0.5;
+        n               8;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+);
+
+virtualMass
+(
+    (airI in water)
+    {
+        type            constantCoefficient;
+        Cvm             0.5;
+    }
+
+    (water in airI)
+    {
+        type            constantCoefficient;
+        Cvm             0.5;
+    }
+
+    (airII in water)
+    {
+        type            constantCoefficient;
+        Cvm             0.5;
+    }
+
+    (water in airII)
+    {
+        type            constantCoefficient;
+        Cvm             0.5;
+    }
+
+    (airI in airII)
+    {
+        type            constantCoefficient;
+        Cvm             0.5;
+    }
+
+    (airI in airII)
+    {
+        type            constantCoefficient;
+        Cvm             0.5;
+    }
+);
+
+heatTransfer
+(
+    (airI in water)
+    {
+        type            RanzMarshall;
+        residualAlpha   1e-4;
+    }
+
+    (water in airI)
+    {
+        type            RanzMarshall;
+        residualAlpha   1e-4;
+    }
+
+    (airII in water)
+    {
+        type            RanzMarshall;
+        residualAlpha   1e-4;
+    }
+
+    (water in airII)
+    {
+        type            RanzMarshall;
+        residualAlpha   1e-4;
+    }
+
+    (airI in airII)
+    {
+        type            RanzMarshall;
+        residualAlpha   1e-4;
+    }
+
+    (airII in airI)
+    {
+        type            RanzMarshall;
+        residualAlpha   1e-4;
+    }
+);
+
+lift
+();
+
+wallLubrication
+();
+
+turbulentDispersion
+();
+
+// Minimum allowable pressure
+pMin            10000;
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/constant/thermophysicalProperties.airI

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.airI;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        molWeight   28.9;
+    }
+    thermodynamics
+    {
+        Cp          1007;
+        Hf          0;
+    }
+    transport
+    {
+        mu          1.84e-05;
+        Pr          0.7;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/constant/thermophysicalProperties.airII

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.airII;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        molWeight   28.9;
+    }
+    thermodynamics
+    {
+        Cp          1007;
+        Hf          0;
+    }
+    transport
+    {
+        mu          1.84e-05;
+        Pr          0.7;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/constant/thermophysicalProperties.water

@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectFluid;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        molWeight   18;
+    }
+    equationOfState
+    {
+        R           3000;
+        rho0        1027;
+    }
+    thermodynamics
+    {
+        Cp          4195;
+        Hf          0;
+    }
+    transport
+    {
+        mu          3.645e-4;
+        Pr          2.289;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/constant/turbulenceProperties.airI

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties.airI;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/constant/turbulenceProperties.airII

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties.airII;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/constant/turbulenceProperties.water

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/system/blockMeshDict

@@ -0,0 +1,61 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    (0 0 0)
+    (0.15 0 0)
+    (0.15 1 0)
+    (0 1 0)
+    (0 0 0.1)
+    (0.15 0 0.1)
+    (0.15 1 0.1)
+    (0 1 0.1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (25 75 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+patches
+(
+    patch inlet
+    (
+        (1 5 4 0)
+    )
+    patch outlet
+    (
+        (3 7 6 2)
+    )
+    wall walls
+    (
+        (0 4 7 3)
+        (2 6 5 1)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/system/controlDict

@@ -0,0 +1,109 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     reactingMultiphaseEulerFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         100;
+
+deltaT          0.005;
+
+writeControl    runTime;
+
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  no;
+
+maxCo           0.5;
+
+maxDeltaT       1;
+
+functions0
+{
+    fieldAverage1
+    {
+        type            fieldAverage;
+        libs            ("libfieldFunctionObjects.so");
+        writeControl    writeTime;
+        fields
+        (
+            U.airI
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+
+            U.airII
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+
+            U.water
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+
+            alpha.airI
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+
+            alpha.airII
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+
+            p
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+        );
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/system/fvSchemes

@@ -0,0 +1,62 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default                         none;
+
+    "div\(phi,alpha.*)"             Gauss vanLeer;
+    "div\(phir,alpha.*,alpha.*\)"   Gauss vanLeer;
+    "div\(phir,alpha.*\)"           Gauss vanLeer;
+
+    "div\(alphaRhoPhi.*,U.*\)"      Gauss limitedLinearV 1;
+    "div\(phi.*,U.*\)"              Gauss limitedLinearV 1;
+
+    "div\(alphaRhoPhi.*,(h|e).*\)"  Gauss limitedLinear 1;
+    "div\(alphaRhoPhi.*,K.*\)"      Gauss limitedLinear 1;
+    "div\(alphaPhi.*,p\)"           Gauss limitedLinear 1;
+
+    "div\(\(\(\(alpha.*\*thermo:rho.*\)*nuEff.*\)\*dev2\(T\(grad\(U.*\)\)\)\)\)" Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/system/fvSolution

@@ -0,0 +1,76 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "alpha.*"
+    {
+        nAlphaCorr      1;
+        nAlphaSubCycles 2;
+    }
+
+    p_rgh
+    {
+        solver          GAMG;
+        smoother        DIC;
+        tolerance       1e-8;
+        relTol          0;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        relTol          0;
+    }
+
+    "U.*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-5;
+        relTol          0;
+        minIter         1;
+    }
+
+    "e.*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-8;
+        relTol          0;
+        minIter         1;
+    }
+}
+
+PIMPLE
+{
+    nOuterCorrectors 3;
+    nCorrectors      1;
+    nNonOrthogonalCorrectors 0;
+    faceMomentum     yes;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        ".*"            1;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnFixedPolydisperse/system/setFieldsDict

@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha.airI 0
+    volScalarFieldValue alpha.airII 1
+    volScalarFieldValue alpha.water 0
+);
+
+regions
+(
+    boxToCell
+    {
+        box (0 0 -0.1) (0.15 0.701 0.1);
+        fieldValues
+        (
+            volScalarFieldValue alpha.airII 0
+            volScalarFieldValue alpha.water 1
+        );
+    }
+);
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/T.airI

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T.airI;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 300;
+
+boundaryField
+{
+    walls
+    {
+        type               zeroGradient;
+    }
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.airI;
+        inletValue         $internalField;
+        value              $internalField;
+    }
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/T.airII

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T.airII;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 300;
+
+boundaryField
+{
+    walls
+    {
+        type               zeroGradient;
+    }
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.airII;
+        inletValue         $internalField;
+        value              $internalField;
+    }
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/T.water

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 350;
+
+boundaryField
+{
+    walls
+    {
+        type               zeroGradient;
+    }
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.water;
+        inletValue         $internalField;
+        value              $internalField;
+    }
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/Theta

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Theta;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 1.0e-7;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1.0e-7;
+        value           uniform 1.0e-7;
+    }
+
+    walls
+    {
+        type            zeroGradient;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/U.airI

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.airI;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0.1 0);
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               pressureInletOutletVelocity;
+        phi                phi.airI;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedValue;
+        value              uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/U.airII

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.airII;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0.1 0);
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               pressureInletOutletVelocity;
+        phi                phi.airII;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedValue;
+        value              uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/U.water

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               pressureInletOutletVelocity;
+        phi                phi.water;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedValue;
+        value              uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/alpha.airI

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.airI;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.25;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.airI;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/alpha.airI.orig

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.airI;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.25;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.airI;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/alpha.airII.orig

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.airII;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.25;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.airII;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/alpha.water.orig

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.5;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.water;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/f.airI.bubbles

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      f.airI.bubbles;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.airI;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/f.airII.bubbles

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                    |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      f.airII.bubbles;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.airII;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/p

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [1 -1 -2 0 0 0 0];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type               calculated;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               calculated;
+        value              $internalField;
+    }
+    walls
+    {
+        type               calculated;
+        value              $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/0/p_rgh

@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [1 -1 -2 0 0 0 0];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            prghTotalPressure;
+        p0              $internalField;
+        U               U.airII;
+        phi             phi.airII;
+        rho             thermo:rho.airII;
+        value           $internalField;
+    }
+    walls
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/Allrun

@@ -0,0 +1,14 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# Set application name
+application=`getApplication`
+
+runApplication blockMesh
+runApplication setFields
+runApplication $application
+
+# ----------------------------------------------------------------- end-of-file