From ae459d60e0e841ee02028761e2e265097d9e7f30 Mon Sep 17 00:00:00 2001
From: "Schlegel, Dr. Fabian (FWDC) - 107795" <f.schlegel@hzdr.de>
Date: Mon, 6 Apr 2020 18:33:59 +0200
Subject: [PATCH] Adding multiphase solver cipsaMultiphaseEulerFoam

    Solver, which makes use of the compact momentum interpolation method,
    according to:

    Cubero, A. and Fueyo, N. (2007).
    A compact momentum interpolation procedure for unsteady flows and
    relaxation.
    Numerical Heat Transfer, Part B: Fundamentals, 52(6):507-529.

    * A partial elimination sum formulation is added to
      approximately resolve phase coupling due to the large drag coefficient for
      an arbitrary number of phases
    * The virtual mass force is handled consistently
    * Necessary infrastructure in phaseSystems is added
    * New drag models are added: StrubeljTiselj, IshiiZuberSparse
    * A new blending method is added: definitePhases
---
 .../multiphase/reactingEulerFoam/Allwclean    |   1 +
 .../multiphase/reactingEulerFoam/Allwmake     |   1 +
 .../cipsaMultiphaseEulerFoam/CourantNo.H      |  60 +++++
 .../cipsaMultiphaseEulerFoam/Make/files       |   3 +
 .../cipsaMultiphaseEulerFoam/Make/options     |  25 ++
 .../cipsaMultiphaseEulerFoam.C                | 125 ++++++++++
 .../createFieldRefs.H                         |  46 ++++
 .../cipsaMultiphaseEulerFoam/createFields.H   |  79 +++++++
 .../initZeroPtrList.H                         |  82 +++++++
 .../cipsaMultiphaseEulerFoam/pU/HUEqn.H       |  79 +++++++
 .../cipsaMultiphaseEulerFoam/pU/UEqns.H       | 205 ++++++++++++++++
 .../cipsaMultiphaseEulerFoam/pU/ddtUEqn.H     |  47 ++++
 .../cipsaMultiphaseEulerFoam/pU/dragDiag.H    |  65 ++++++
 .../cipsaMultiphaseEulerFoam/pU/dragSource.H  |  68 ++++++
 .../cipsaMultiphaseEulerFoam/pU/dragSourcef.H |  61 +++++
 .../pU/forcesSource.H                         |  36 +++
 .../pU/forcesSourcef.H                        |  36 +++
 .../cipsaMultiphaseEulerFoam/pU/pEqn.H        | 218 ++++++++++++++++++
 .../pU/stabilizeHUEqn.H                       |  43 ++++
 .../cipsaMultiphaseEulerFoam/pU/vmDiag.H      |  16 ++
 .../cipsaMultiphaseEulerFoam/pU/vmSource.H    |  22 ++
 .../cipsaMultiphaseEulerFoam/pU/vmSourcef.H   |  23 ++
 .../interfacialModels/Make/files              |   2 +
 .../IshiiZuberSparse/IshiiZuberSparse.C       |  83 +++++++
 .../IshiiZuberSparse/IshiiZuberSparse.H       | 105 +++++++++
 .../StrubeljTiselj/StrubeljTiselj.C           |  94 ++++++++
 .../StrubeljTiselj/StrubeljTiselj.H           | 121 ++++++++++
 .../BlendedInterfacialModel.C                 |   4 +-
 .../definitePhases/definitePhases.C           | 140 +++++++++++
 .../definitePhases/definitePhases.H           | 152 ++++++++++++
 .../reactingEulerFoam/phaseSystems/Make/files |   1 +
 .../MomentumTransferPhaseSystem.C             | 203 +++++++++++++++-
 .../MomentumTransferPhaseSystem.H             | 139 ++++++++++-
 .../MovingPhaseModel/MovingPhaseModel.C       |  12 +-
 .../MovingPhaseModel/MovingPhaseModel.H       |   7 +-
 .../StationaryPhaseModel.C                    |  19 +-
 .../StationaryPhaseModel.H                    |   7 +-
 .../phaseModel/phaseModel/phaseModel.C        |  15 +-
 .../phaseModel/phaseModel/phaseModel.H        |  20 +-
 .../phasePair/phasePair/phasePair.C           |  23 +-
 .../phaseSystems/phaseSystem/phaseSystem.C    |  89 ++++++-
 .../phaseSystems/phaseSystem/phaseSystem.H    |  30 ++-
 .../phaseSystems/phaseSystem/phaseSystemI.H   |  11 +-
 .../multiphaseSystem/multiphaseSystem.H       | 132 ++++++++++-
 44 files changed, 2727 insertions(+), 23 deletions(-)
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/CourantNo.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/Make/files
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/Make/options
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/cipsaMultiphaseEulerFoam.C
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/createFieldRefs.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/createFields.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/initZeroPtrList.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/HUEqn.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/UEqns.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/ddtUEqn.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/dragDiag.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSource.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSourcef.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSource.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSourcef.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/pEqn.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/stabilizeHUEqn.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/vmDiag.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSource.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSourcef.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.C
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.C
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.H
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.C
 create mode 100644 applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.H

diff --git a/applications/solvers/multiphase/reactingEulerFoam/Allwclean b/applications/solvers/multiphase/reactingEulerFoam/Allwclean
index 40a8acade..a5f64a484 100755
--- a/applications/solvers/multiphase/reactingEulerFoam/Allwclean
+++ b/applications/solvers/multiphase/reactingEulerFoam/Allwclean
@@ -7,6 +7,7 @@ wclean libso interfacialCompositionModels
 wclean libso derivedFvPatchFields
 reactingTwoPhaseEulerFoam/Allwclean
 reactingMultiphaseEulerFoam/Allwclean
+wclean libso cipsaMultiphaseEulerFoam
 wclean libso functionObjects
 
 #------------------------------------------------------------------------------
diff --git a/applications/solvers/multiphase/reactingEulerFoam/Allwmake b/applications/solvers/multiphase/reactingEulerFoam/Allwmake
index 77f61dd78..2d77123ee 100755
--- a/applications/solvers/multiphase/reactingEulerFoam/Allwmake
+++ b/applications/solvers/multiphase/reactingEulerFoam/Allwmake
@@ -10,6 +10,7 @@ wmake $targetType interfacialCompositionModels
 wmake $targetType derivedFvPatchFields
 reactingTwoPhaseEulerFoam/Allwmake $targetType $*
 reactingMultiphaseEulerFoam/Allwmake $targetType $*
+wmake $targetType cipsaMultiphaseEulerFoam
 wmake $targetType functionObjects
 
 #------------------------------------------------------------------------------
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/CourantNo.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/CourantNo.H
new file mode 100644
index 000000000..d1d28ab5c
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/CourantNo.H
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2020 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/>.
+
+Global
+    CourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar CoNum = 0.0;
+scalar meanCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    scalarField sumPhi
+    (
+        fvc::surfaceSum(mag(phi))().primitiveField()
+    );
+
+    forAll(phases, phasei)
+    {
+        sumPhi = max
+        (
+            sumPhi,
+            fvc::surfaceSum(mag(phases[phasei].phi()))().primitiveField()
+        );
+    }
+
+    CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Courant Number mean: " << meanCoNum
+    << " max: " << CoNum << endl;
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/Make/files b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/Make/files
new file mode 100644
index 000000000..3299c2137
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/Make/files
@@ -0,0 +1,3 @@
+cipsaMultiphaseEulerFoam.C
+
+EXE = $(FOAM_APPBIN)/cipsaMultiphaseEulerFoam
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/Make/options b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/Make/options
new file mode 100644
index 000000000..f8cc5e8d8
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/Make/options
@@ -0,0 +1,25 @@
+EXE_INC = \
+    -I../reactingMultiphaseEulerFoam/multiphaseSystem/lnInclude \
+    -I../phaseSystems/lnInclude \
+    -I../interfacialModels/lnInclude \
+    -I../interfacialCompositionModels/lnInclude \
+    -I../reactingMultiphaseEulerFoam/multiphaseCompressibleTurbulenceModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -lreactingPhaseSystem \
+    -lreactingMultiphaseSystem \
+    -lreactingEulerianInterfacialModels \
+    -lreactingEulerianInterfacialCompositionModels \
+    -lmultiphaseReactingTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/cipsaMultiphaseEulerFoam.C b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/cipsaMultiphaseEulerFoam.C
new file mode 100644
index 000000000..fc465e032
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/cipsaMultiphaseEulerFoam.C
@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 2018-2020 OpenFOAM Foundation
+  ==  == ====== ====   ==  ==  | and Helmholtz-Zentrum Dresden-Rossendorf
+-------------------------------------------------------------------------------
+License
+    This file is a derivative work 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/>.
+
+Application
+    cipsaMultiphaseEulerFoam
+
+Description
+    Solver for a system of any number of incompressible fluid phases with a
+    common pressure, but otherwise separate properties. The solver uses the
+    consistent momentum interpolation method of Cubero et al. (2014) and a
+    partial-elimination algorithm.
+
+    Reference:
+    \verbatim
+        Cubero, A., Sánchez-Insa, A., & Fueyo, N. (2014).
+        A consistent momentum interpolation method for steady and unsteady
+        multiphase flows.
+        Computers & Chemical Engineering, 62, 96-107.
+    \endverbatim
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "multiphaseSystem.H"
+#include "phaseCompressibleTurbulenceModel.H"
+#include "pimpleControl.H"
+#include "initZeroPtrList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    Switch faceMomentum
+    (
+        pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
+    );
+
+    Switch partialElimination
+    (
+        pimple.dict().lookupOrDefault<Switch>("partialElimination", true)
+    );
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (pimple.run(runTime))
+    {
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity corrector loop (PIMPLE)
+        while (pimple.loop())
+        {
+            fluid.solve(rAUs, rHAUfs);
+            fluid.correct();
+
+            #include "pU/UEqns.H"
+
+            // --- Pressure corrector loop (PISO)
+            while (pimple.correct())
+            {
+                #include "pU/pEqn.H"
+            }
+
+            // Update material derivitives
+            fluid.correctKinematics();
+
+            if (pimple.turbCorr())
+            {
+                fluid.correctTurbulence();
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/createFieldRefs.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/createFieldRefs.H
new file mode 100644
index 000000000..249867420
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/createFieldRefs.H
@@ -0,0 +1,46 @@
+surfaceScalarField& phi = fluid.phi();
+
+List<scalar> relaxU(phases.size());
+PtrList<surfaceScalarField> alphafs(phases.size());
+
+// Equations
+PtrList<fvVectorMatrix> HUEqns(phases.size());
+PtrList<fvVectorMatrix> ddtUEqns(phases.size());
+PtrList<surfaceScalarField> ddtPhiCoeffs(phases.size());
+
+// Cell-based coefficients
+PtrList<volScalarField> rHUs(phases.size());
+PtrList<volScalarField> rAs(phases.size());
+PtrList<volScalarField> rAUs; // Dummy argument for solveAlphas()
+
+PtrList<volScalarField> dragDiag(phases.size());
+PtrList<volVectorField> dragSource(phases.size());
+
+PtrList<volScalarField> vmDiag(phases.size());
+PtrList<volVectorField> vmSource(phases.size());
+
+PtrList<volVectorField> forcesSource(phases.size());
+
+// Face-based coefficients
+PtrList<surfaceScalarField> rHUfs(phases.size());
+PtrList<surfaceScalarField> rAfs(phases.size());
+PtrList<surfaceScalarField> rAUfs(phases.size());
+PtrList<surfaceScalarField> rHAUfs(phases.size()); // Argument for solveAlphas()
+PtrList<surfaceScalarField> alpharAUfs(phases.size());
+
+PtrList<surfaceScalarField> dragDiagf(phases.size());
+PtrList<surfaceScalarField> dragSourcef(phases.size());
+PtrList<surfaceScalarField> dragPCoefff(phases.size());
+
+PtrList<surfaceScalarField> vmDiagf(phases.size());
+PtrList<surfaceScalarField> vmSourcef(phases.size());
+
+PtrList<surfaceScalarField> forcesSourcef(phases.size());
+
+// Cell velocities
+PtrList<volVectorField> HbyAs(phases.size());
+
+// Face fluxes
+PtrList<surfaceScalarField> HbyASumf(phases.size());
+PtrList<surfaceScalarField> phigs(phases.size());
+PtrList<surfaceScalarField> phiHbyAs(phases.size());
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/createFields.H
new file mode 100644
index 000000000..3f659cd22
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/createFields.H
@@ -0,0 +1,79 @@
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+
+Info<< "\nCreating phaseSystem\n" << endl;
+
+autoPtr<multiphaseSystem> fluidPtr
+(
+    multiphaseSystem::New(mesh)
+);
+multiphaseSystem& fluid = fluidPtr();
+multiphaseSystem::phaseModelList& phases = fluid.phases();
+
+// Create mixture density and velocity
+volScalarField rho("rho", fluid.rho());
+volVectorField U("U", fluid.U());
+
+dimensionedScalar pMin
+(
+    "pMin",
+    dimPressure,
+    fluid
+);
+
+#include "gh.H"
+
+volScalarField& p = phases[0].thermoRef().p();
+
+Info<< "Reading field p_rgh\n" << endl;
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+if (fluid.incompressible())
+{
+    p = max(p_rgh + fluid.rho()*gh, pMin);
+
+    if (p_rgh.needReference())
+    {
+        setRefCell
+        (
+            p,
+            p_rgh,
+            pimple.dict(),
+            pRefCell,
+            pRefValue
+        );
+
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - fluid.rho()*gh;
+    }
+}
+
+mesh.setFluxRequired(p_rgh.name());
+
+// Do some checks and abort if necessary
+if (Switch(pimple.dict().lookup("correctFixedFluxBCs")[0].wordToken()))
+{
+    FatalErrorInFunction
+        << "Option correctFixedFluxBCs "
+        << "- required for BlendedInterfacialModel - "
+        << "has to be turned off for this solver."
+        << exit(FatalError);
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/initZeroPtrList.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/initZeroPtrList.H
new file mode 100644
index 000000000..c7816d470
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/initZeroPtrList.H
@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 2018-2020 OpenFOAM Foundation
+  ==  == ====== ====   ==  ==  | and Helmholtz-Zentrum Dresden-Rossendorf
+-------------------------------------------------------------------------------
+License
+    This file is a derivative work 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/>.
+
+InNamespace
+    Foam
+
+Description
+    Initialize the fields of a PtrList with zero
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef initZeroPtrList_H
+#define initZeroPtrList_H
+
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type,  template< class > class PatchField, class GeoMesh>
+void initZeroPtrList
+(
+    PtrList<GeometricField<Type, PatchField, GeoMesh>>& a,
+    const fvMesh& mesh,
+    const dimensionSet& ds,
+    const word& name
+)
+{
+    using T = Foam::GeometricField<Type, PatchField, GeoMesh>;
+
+    forAll(a, i)
+    {
+        if (a.set(i))
+        {
+            a[i] == Zero;
+        }
+        else
+        {
+            a.set
+            (
+                i,
+                T::New(name, mesh, dimensioned<Type>(ds, Zero))
+            );
+        }
+    }
+
+    return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/HUEqn.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/HUEqn.H
new file mode 100644
index 000000000..07419a8b5
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/HUEqn.H
@@ -0,0 +1,79 @@
+forAll(phases, phasei)
+{
+    phaseModel& phase = phases[phasei];
+    volVectorField& UAlpha = phase.URef();
+
+    HUEqns.set
+    (
+        phasei,
+        (
+            fvm::div(phase.alphaRhoPhi(), UAlpha)
+          + phase.divDevRhoReff(UAlpha)
+          + fvm::SuSp(- fvc::div(phase.alphaRhoPhi()), UAlpha)
+        )
+    );
+}
+
+// Create U & grad(U) fields
+PtrList<fvVectorMatrix> UgradUs(phases.size());
+forAll(phases, phasei)
+{
+    const phaseModel& phase = phases[phasei];
+    const volVectorField& UAlpha = phase.U();
+    const surfaceScalarField& phiAlpha = phase.phi();
+
+    UgradUs.set
+    (
+        phasei,
+        new fvVectorMatrix
+        (
+            fvm::div(phiAlpha, UAlpha)
+          + fvm::SuSp(- fvc::div(phiAlpha), UAlpha)
+        )
+    );
+}
+
+// Add convective part of the virtual mass force
+forAllConstIter(multiphaseSystem::VmTable, fluid.Vms(), VmIter)
+{
+    const volScalarField& Vm(*VmIter());
+    const phasePair& pair(fluid.phasePairs()[VmIter.key()]);
+
+    forAllConstIter(phasePair, pair, iter)
+    {
+        const phaseModel& phase = iter();
+        const phaseModel& otherPhase = iter.otherPhase();
+
+        HUEqns[phase.index()] +=
+            Vm
+           *(
+                UgradUs[phase.index()]
+              - (UgradUs[otherPhase.index()] & otherPhase.U())
+            );
+    }
+}
+
+// Relaxation factor
+forAll(phases, phasei)
+{
+    word name = HUEqns[phasei].psi().select(pimple.finalPimpleIter());
+
+    relaxU[phasei] = 1;
+
+    if (HUEqns[phasei].psi().mesh().relaxField(name))
+    {
+        relaxU[phasei] =
+            HUEqns[phasei].psi().mesh().fieldRelaxationFactor(name);
+
+        if
+        (
+            pimple.finalPimpleIter()
+         && HUEqns[phasei].psi().mesh().fieldRelaxationFactor(name) != 1
+        )
+        {
+            Warning << "Field relaxation factor of " << relaxU[phasei]
+                    << " for solution variable '" << name
+                    << "' is not recommended." << endl << endl;
+        }
+    }
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/UEqns.H
new file mode 100644
index 000000000..9583f9915
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/UEqns.H
@@ -0,0 +1,205 @@
+// Update mixture properties
+rho = fluid.rho();
+U = fluid.U();
+
+// Store velocity and flux of previous iteration
+forAll(phases, phasei)
+{
+    phases[phasei].URef().storePrevIter();
+    phases[phasei].phiRef().storePrevIter();
+}
+
+// Mixture density gradient times gravitational acceleration
+const surfaceScalarField ghSnGradRho
+(
+    "ghSnGradRho",
+    ghf*fvc::snGrad(rho)*mesh.magSf()
+);
+
+// Update drag and virtual mass coefficient table and interfacial forces tables
+fluid.momentumTransfer();
+
+// Time derivative equation
+#include "pU/ddtUEqn.H"
+
+// Convection, diffusion equation
+#include "pU/HUEqn.H"
+#include "pU/stabilizeHUEqn.H"
+
+// alphafs, rHUs, rHUfs, rAs, rAfs
+forAll(phases, phasei)
+{
+    phaseModel& phase = phases[phasei];
+    const volScalarField& alpha = phase;
+
+    // Face volume fractions
+    alphafs.set(phasei, fvc::interpolate(alpha));
+
+    rHUs.set(phasei, 1.0/HUEqns[phasei].A());
+    rHUfs.set(phasei, fvc::interpolate(rHUs[phasei]));
+
+    // Diagonal coefficients
+    rAs.set
+    (
+        phasei,
+        1.0/(HUEqns[phasei].A() + ddtUEqns[phasei].A())
+    );
+    rAfs.set
+    (
+        phasei,
+        1.0
+       /(
+            fvc::interpolate(HUEqns[phasei].A())
+          + fvc::interpolate(ddtUEqns[phasei].A())
+        )
+    );
+}
+
+// Implicit and pressure part of the drag
+#include "pU/dragDiag.H"
+
+// Implicit part of the time derivative of the virtual mass
+#include "pU/vmDiag.H"
+
+// Momentum predictor
+if (pimple.momentumPredictor())
+{
+    // Explicit part of the drag terms
+    #include "pU/dragSource.H"
+
+    // Explicit part of the time derivative of the virtual mass
+    #include "pU/vmSource.H"
+
+    // Add lift, wall lubrication and turbulent dispersion forces
+    #include "pU/forcesSource.H"
+
+    // rAUfs, alpharAUfs, phigs WITHOUT field relaxation of UAlpha
+    forAll(phases, phasei)
+    {
+        phaseModel& phase = phases[phasei];
+
+        rAUfs.set
+        (
+            phasei,
+            1.0
+           /(
+                1.0
+              + rHUfs[phasei]*fvc::interpolate(ddtUEqns[phasei].A())
+              + rHUfs[phasei]*dragDiagf[phasei]
+              + rHUfs[phasei]*vmDiagf[phasei]
+            )
+        );
+
+        // Phase diagonal coefficients
+        alpharAUfs.set
+        (
+            phasei,
+            rHUfs[phasei]*rAUfs[phasei]
+           *(alphafs[phasei] + dragPCoefff[phasei])
+        );
+
+        // Combined buoyancy and force fluxes
+        phigs.set
+        (
+            phasei,
+            alpharAUfs[phasei]
+           *(
+              - ghSnGradRho
+              + fluid.surfaceTension(phase)*mesh.magSf()
+            )
+        );
+    }
+
+    // Momentum predictor equation
+    forAll(phases, phasei)
+    {
+        phaseModel& phase = phases[phasei];
+        const volScalarField& alpha = phase;
+        volVectorField& UAlpha = phase.URef();
+
+        fvVectorMatrix UEqns
+        (
+            ddtUEqns[phasei]
+          + HUEqns[phasei]
+          + fvm::Sp(dragDiag[phasei], UAlpha)
+          + fvm::Sp(vmDiag[phasei], UAlpha)
+            ==
+            dragSource[phasei]
+          + alpha*(phase.rho() - rho)*g
+          + (
+                ddtUEqns[phasei].A()
+              + HUEqns[phasei].A()
+              + dragDiag[phasei]
+            )
+           *fvc::reconstruct
+            (
+                phigs[phasei]
+              - alpharAUfs[phasei]*fvc::snGrad(p_rgh)*mesh.magSf()
+            )
+          + vmSource[phasei]
+          + forcesSource[phasei]
+        );
+
+        UEqns.relax();
+        UEqns.solve();
+    }
+}
+
+// rAUfs, alpharAUfs, phigs WITH field relaxation of UAlpha
+forAll(phases, phasei)
+{
+    phaseModel& phase = phases[phasei];
+
+    rAUfs.set
+    (
+        phasei,
+        1.0
+       /(
+            1.0/max(relaxU[phasei], small)
+          + rHUfs[phasei]*fvc::interpolate(ddtUEqns[phasei].A())
+          + rHUfs[phasei]*dragDiagf[phasei]
+          + rHUfs[phasei]*vmDiagf[phasei]
+        )
+    );
+
+    rHAUfs.set(phasei, rHUfs[phasei]*rAUfs[phasei]);
+
+    // Phase diagonal coefficients
+    alpharAUfs.set
+    (
+        phasei,
+        rHUfs[phasei]*rAUfs[phasei]
+       *(alphafs[phasei] + dragPCoefff[phasei])
+    );
+
+    // Combined buoyancy and force fluxes
+    phigs.set
+    (
+        phasei,
+        alpharAUfs[phasei]
+       *(
+          - ghSnGradRho
+          + fluid.surfaceTension(phase)*mesh.magSf()
+        )
+    );
+}
+
+// Construct pressure "diffusivity"
+surfaceScalarField rAUf
+(
+    IOobject
+    (
+        "rAUf",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar(dimensionSet(-1, 3, 1, 0, 0), Zero)
+);
+
+forAll(phases, phasei)
+{
+    rAUf += alphafs[phasei]*alpharAUfs[phasei];
+}
+
+rAUf = mag(rAUf);
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/ddtUEqn.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/ddtUEqn.H
new file mode 100644
index 000000000..e9cfd2abb
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/ddtUEqn.H
@@ -0,0 +1,47 @@
+forAll(phases, phasei)
+{
+    phaseModel& phase = phases[phasei];
+    const volScalarField& alpha = phase;
+    const volScalarField& rhoAlpha = phase.rho();
+    volVectorField& UAlpha = phase.URef();
+
+    ddtUEqns.set
+    (
+        phasei,
+        fvm::ddt(alpha, rhoAlpha, UAlpha)
+      - fvm::Sp(fvc::ddt(alpha, rhoAlpha), UAlpha)
+    );
+
+    word ddtScheme
+    (
+        mesh.ddtScheme
+        (
+            "ddt("
+          + alpha.name()
+          + ',' + rhoAlpha.name() + ','
+          + UAlpha.name() + ')'
+        ).last().wordToken()
+    );
+
+    if (ddtScheme == "Euler")
+    {
+        ddtPhiCoeffs.set
+        (
+            phasei,
+            (
+                byDt
+                (
+                    fvc::interpolate(alpha.oldTime())
+                   *fvc::interpolate(rhoAlpha.oldTime())
+                   *phase.phiRef().oldTime()
+                )
+            )
+        );
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Only Euler is an allowed time scheme "
+            << "for this solver." << exit(FatalError);
+    }
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/dragDiag.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/dragDiag.H
new file mode 100644
index 000000000..2f564e07c
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/dragDiag.H
@@ -0,0 +1,65 @@
+initZeroPtrList(dragDiag, mesh, dimensionSet(1, -3, -1, 0, 0), "dragDiag");
+initZeroPtrList(dragDiagf, mesh, dimensionSet(1, -3, -1, 0, 0), "dragDiagf");
+initZeroPtrList(dragPCoefff, mesh, dimensionSet(0, 0, 0, 0, 0), "dragPCoefff");
+
+forAllConstIter(multiphaseSystem::KdTable, fluid.Kds(), KdIter)
+{
+    const volScalarField& K(*KdIter());
+    const phasePair& pair(fluid.phasePairs()[KdIter.key()]);
+
+    forAllConstIter(phasePair, pair, iter)
+    {
+        const phaseModel& phase = iter();
+        const phaseModel& otherPhase = iter.otherPhase();
+
+        if (partialElimination)
+        {
+            dragDiag[phase.index()] +=
+                K
+               *(
+                    1
+                  + rAs[otherPhase.index()]
+                   *(
+                        HUEqns[phase.index()].A()
+                      + ddtUEqns[phase.index()].A()
+                    )
+                );
+        }
+        else
+        {
+            dragDiag[phase.index()] += K;
+        }
+    }
+}
+
+forAllConstIter(multiphaseSystem::KdfTable, fluid.Kdfs(), KdfIter)
+{
+    const surfaceScalarField& Kf(*KdfIter());
+    const phasePair& pair(fluid.phasePairs()[KdfIter.key()]);
+
+    forAllConstIter(phasePair, pair, iter)
+    {
+        const phaseModel& phase = iter();
+        const phaseModel& otherPhase = iter.otherPhase();
+
+        if (partialElimination)
+        {
+            dragDiagf[phase.index()] +=
+                Kf
+               *(
+                    1
+                  + rAfs[otherPhase.index()]
+                   *(
+                        fvc::interpolate(HUEqns[phase.index()].A())
+                      + fvc::interpolate(ddtUEqns[phase.index()].A())
+                    )
+                );
+
+            dragPCoefff[phase.index()] += Kf*rAfs[otherPhase.index()];
+        }
+        else
+        {
+            dragDiagf[phase.index()] += Kf;
+        }
+    }
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSource.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSource.H
new file mode 100644
index 000000000..4ffda234e
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSource.H
@@ -0,0 +1,68 @@
+// Summed RHS of cell-based momentum balance equation
+volVectorField HSum
+(
+    IOobject
+    (
+        "HSum",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedVector(dimensionSet(1, -2, -2, 0, 0), Zero)
+);
+
+if (partialElimination)
+{
+    forAll(phases, phasei)
+    {
+        HSum +=
+            HUEqns[phasei].H()
+          + ddtUEqns[phasei].H();
+    }
+}
+
+initZeroPtrList
+(
+    dragSource, mesh, dimensionSet(1, -2, -2, 0, 0), "dragSource"
+);
+
+
+forAllConstIter(multiphaseSystem::KdTable, fluid.Kds(), KdIter)
+{
+    const volScalarField& K(*KdIter());
+    const phasePair& pair(fluid.phasePairs()[KdIter.key()]);
+
+    forAllConstIter(phasePair, pair, iter)
+    {
+        const phaseModel& phase = iter();
+        const phaseModel& otherPhase = iter.otherPhase();
+
+        if (partialElimination)
+        {
+            dragSource[phase.index()] += K*rAs[otherPhase.index()]*HSum;
+
+            forAll(phases, phasei)
+            {
+                if
+                (
+                    (phasei == phase.index())
+                 || (phasei == otherPhase.index())
+                ) continue;
+
+                dragSource[phase.index()] -=
+                (
+                    K*rAs[otherPhase.index()]
+                   *(
+                        HUEqns[phasei].A()
+                      + ddtUEqns[phasei].A()
+                    )
+                   *phases[phasei].U()
+                );
+            }
+        }
+        else
+        {
+            dragSource[phase.index()] += K*phases[otherPhase.index()].U();
+        }
+    }
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSourcef.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSourcef.H
new file mode 100644
index 000000000..7691b08f4
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSourcef.H
@@ -0,0 +1,61 @@
+// Summed RHS of face-based momentum balance equation
+if (partialElimination)
+{
+    initZeroPtrList(HbyASumf, mesh, dimensionSet(0, 3, -1, 0, 0), "HbyASumf");
+    forAll(phases, phasei)
+    {
+        forAll(phases, phasej)
+        {
+            HbyASumf[phasei] +=
+                rHUfs[phasei]
+               *(
+                    fvc::flux(HUEqns[phasej].H())
+                  + ddtPhiCoeffs[phasej]
+                );
+        }
+    }
+}
+
+initZeroPtrList(dragSourcef, mesh, dimensionSet(0, 3, -1, 0, 0), "dragSourcef");
+
+forAllConstIter(multiphaseSystem::KdfTable, fluid.Kdfs(), KdfIter)
+{
+    const surfaceScalarField& Kf(*KdfIter());
+    const phasePair& pair(fluid.phasePairs()[KdfIter.key()]);
+
+    forAllConstIter(phasePair, pair, iter)
+    {
+        const phaseModel& phase = iter();
+        const phaseModel& otherPhase = iter.otherPhase();
+
+        if (partialElimination)
+        {
+            dragSourcef[phase.index()]
+                += Kf*rAfs[otherPhase.index()]*HbyASumf[phase.index()];
+
+            forAll(phases, phasei)
+            {
+                if
+                (
+                    (phasei == phase.index())
+                 || (phasei == otherPhase.index())
+                ) continue;
+
+                dragSourcef[phase.index()] -=
+                (
+                    Kf*rAfs[otherPhase.index()]
+                   *(
+                        fvc::interpolate(HUEqns[phasei].A())
+                      + fvc::interpolate(ddtUEqns[phasei].A())
+                    )
+                   *rHUfs[phase.index()]*phases[phasei].phi()
+                );
+            }
+        }
+        else
+        {
+            dragSourcef[phase.index()]
+                += Kf*rHUfs[phase.index()]*phases[otherPhase.index()].phi();
+        }
+    }
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSource.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSource.H
new file mode 100644
index 000000000..792f35611
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSource.H
@@ -0,0 +1,36 @@
+initZeroPtrList
+(
+    forcesSource, mesh, dimensionSet(1, -2, -2, 0, 0), "forcesSource"
+);
+
+// Add lift force
+// Inspired by MomentumTransferPhaseSystem::phiFfs()
+forAllConstIter(multiphaseSystem::FlTable, fluid.Fls(), FlIter)
+{
+    const volVectorField& F(*FlIter());
+    const phasePair& pair(fluid.phasePairs()[FlIter.key()]);
+
+    forcesSource[pair.phase1().index()] -= F;
+    forcesSource[pair.phase2().index()] += F;
+}
+
+// Add wall lubrication force
+// Inspired by MomentumTransferPhaseSystem::phiFfs()
+forAllConstIter(multiphaseSystem::FwTable, fluid.Fws(), FwIter)
+{
+    const volVectorField& F(*FwIter());
+    const phasePair& pair(fluid.phasePairs()[FwIter.key()]);
+
+    forcesSource[pair.phase1().index()] -= F;
+    forcesSource[pair.phase2().index()] += F;
+}
+
+// Add turbulent dispersion force
+forAllConstIter(multiphaseSystem::FtTable, fluid.Fts(), FtIter)
+{
+    const volVectorField& F(*FtIter());
+    const phasePair& pair(fluid.phasePairs()[FtIter.key()]);
+
+    forcesSource[pair.phase1().index()] -= F;
+    forcesSource[pair.phase2().index()] += F;
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSourcef.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSourcef.H
new file mode 100644
index 000000000..43059959a
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSourcef.H
@@ -0,0 +1,36 @@
+initZeroPtrList
+(
+    forcesSourcef, mesh, dimensionSet(1, 0, -2, 0, 0), "forcesSourcef"
+);
+
+// Add lift force
+// Inspired by MomentumTransferPhaseSystem::phiFfs()
+forAllConstIter(multiphaseSystem::FlfTable, fluid.Flfs(), FlfIter)
+{
+    const surfaceScalarField& Ff(*FlfIter());
+    const phasePair& pair(fluid.phasePairs()[FlfIter.key()]);
+
+    forcesSourcef[pair.phase1().index()] -= Ff;
+    forcesSourcef[pair.phase2().index()] += Ff;
+}
+
+// Add wall lubrication force
+// Inspired by MomentumTransferPhaseSystem::phiFfs()
+forAllConstIter(multiphaseSystem::FwfTable, fluid.Fwfs(), FwfIter)
+{
+    const surfaceScalarField& Ff(*FwfIter());
+    const phasePair& pair(fluid.phasePairs()[FwfIter.key()]);
+
+    forcesSourcef[pair.phase1().index()] -= Ff;
+    forcesSourcef[pair.phase2().index()] += Ff;
+}
+
+// Add turbulent dispersion force
+forAllConstIter(multiphaseSystem::FtfTable, fluid.Ftfs(), FtfIter)
+{
+    const surfaceScalarField& Ff(*FtfIter());
+    const phasePair& pair(fluid.phasePairs()[FtfIter.key()]);
+
+    forcesSourcef[pair.phase1().index()] -= Ff;
+    forcesSourcef[pair.phase2().index()] += Ff;
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/pEqn.H
new file mode 100644
index 000000000..a8de62007
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/pEqn.H
@@ -0,0 +1,218 @@
+// Correct fixed-flux BCs to be consistent with the velocity BCs
+// functionality of MRF.correctBoundaryFlux(U, phi)
+forAll(fluid.movingPhases(), movingPhasei)
+{
+    phaseModel& phase = fluid.movingPhases()[movingPhasei];
+
+    FieldField<fvsPatchField, scalar> Uf
+    (
+        mesh.Sf().boundaryField() & phase.U()().boundaryField()
+    );
+
+    surfaceScalarField::Boundary& phibf = phase.phiRef().boundaryFieldRef();
+
+    forAll(mesh.boundary(), patchi)
+    {
+        if
+        (
+            isA<fixedValueFvsPatchScalarField>(phibf[patchi])
+        )
+        {
+            phibf[patchi] == Uf[patchi];
+        }
+    }
+}
+
+// Update explicit drag terms for PISO loop
+#include "pU/dragSource.H"
+#include "pU/dragSourcef.H"
+
+// Explicit part of the time derivative of the virtual mass for PISO loop
+#include "pU/vmSource.H"
+#include "pU/vmSourcef.H"
+
+// Update lift, wall lubrication and turbulent dispersion forces for PISO loop
+#include "pU/forcesSource.H"
+#include "pU/forcesSourcef.H"
+
+// Pseudo convected velocities
+forAll(phases, phasei)
+{
+    const volVectorField& UAlpha = phases[phasei].U();
+
+    HbyAs.set
+    (
+        phasei,
+        constrainHbyA
+        (
+            rHUs[phasei]*HUEqns[phasei].H(),
+            UAlpha,
+            p_rgh
+        )
+    );
+}
+
+// Predicted fluxes for each phase
+forAll(phases, phasei)
+{
+    phiHbyAs.set
+    (
+        phasei,
+        rAUfs[phasei]
+       *(
+            fvc::flux(HbyAs[phasei])
+          + rHUfs[phasei]*ddtPhiCoeffs[phasei]
+          + dragSourcef[phasei]
+          + rHUfs[phasei]*alphafs[phasei]*(g & mesh.Sf())
+           *(
+                fvc::interpolate(phases[phasei].rho()) - fvc::interpolate(rho)
+            )
+          + (1 - relaxU[phasei])/max(relaxU[phasei], small)
+           *phases[phasei].phiRef().prevIter()
+          + vmSourcef[phasei]
+          + rHUfs[phasei]*forcesSourcef[phasei]
+        )
+    );
+}
+
+// Total predicted flux
+surfaceScalarField phiHbyA
+(
+    IOobject
+    (
+        "phiHbyA",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(dimArea*dimVelocity, 0)
+);
+
+forAll(phases, phasei)
+{
+    phiHbyA += alphafs[phasei]*(phiHbyAs[phasei] + phigs[phasei]);
+}
+
+// Update the fixedFluxPressure BCs to ensure flux consistency
+{
+    surfaceScalarField::Boundary phib(phi.boundaryField());
+    phib = 0;
+    forAll(phases, phasei)
+    {
+        phaseModel& phase = phases[phasei];
+        phib += alphafs[phasei].boundaryField()*phase.phiRef().boundaryField();
+    }
+
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryFieldRef(),
+        (
+            phiHbyA.boundaryField() - phib
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+}
+
+// Iterate over the pressure equation to correct for non-orthogonality
+while (pimple.correctNonOrthogonal())
+{
+    // Construct the transport part of the pressure equation
+    fvScalarMatrix p_rghEqnIncomp
+    (
+        fvc::div(phiHbyA)
+      - fvm::laplacian(rAUf, p_rgh)
+    );
+
+    // Solve
+    {
+        fvScalarMatrix p_rghEqn(p_rghEqnIncomp);
+
+        if (fluid.incompressible())
+        {
+            p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        }
+
+        p_rghEqn.solve();
+    }
+
+    // Correct fluxes and velocities on last non-orthogonal iteration
+    if (pimple.finalNonOrthogonalIter())
+    {
+        phi = phiHbyA + p_rghEqnIncomp.flux();
+
+        surfaceScalarField mSfGradp("mSfGradp", p_rghEqnIncomp.flux()/rAUf);
+
+        forAll(phases, phasei)
+        {
+            phaseModel& phase = phases[phasei];
+
+            phase.phiRef() =
+                phiHbyAs[phasei]
+              + phigs[phasei]
+              + alpharAUfs[phasei]*mSfGradp;
+        }
+
+        // Optionally relax pressure for velocity correction
+        p_rgh.relax();
+
+        mSfGradp = p_rghEqnIncomp.flux()/rAUf;
+
+        // Momentum corrector
+        forAll(phases, phasei)
+        {
+            phaseModel& phase = phases[phasei];
+            const volScalarField& alpha = phase;
+            volVectorField& UAlpha = phase.URef();
+
+            if (faceMomentum)
+            {
+                UAlpha =
+                    fvc::reconstruct(phase.phi());
+            }
+            else
+            {
+                UAlpha =
+                    1.0
+                   /(
+                        ddtUEqns[phasei].A()
+                      + HUEqns[phasei].A()/max(relaxU[phasei], small)
+                      + dragDiag[phasei]
+                      + vmDiag[phasei]
+                    )
+                   *(
+                        HUEqns[phasei].A()*HbyAs[phasei]
+                      + ddtUEqns[phasei].H()
+                      + dragSource[phasei]
+                      + alpha*(phase.rho() - rho)*g
+                      + (1 - relaxU[phasei])/max(relaxU[phasei], small)
+                       *HUEqns[phasei].A()*UAlpha.prevIter()
+                      + vmSource[phasei]
+                      + forcesSource[phasei]
+                    )
+                  + fvc::reconstruct
+                    (
+                        phigs[phasei]
+                      + alpharAUfs[phasei]*mSfGradp
+                    );
+            }
+
+            UAlpha.correctBoundaryConditions();
+        }
+    }
+}
+
+// Update the mixture properties and static pressure
+U = fluid.U();
+p = p_rgh + rho*gh;
+
+// Account for static pressure reference
+if (p_rgh.needReference() && fluid.incompressible())
+{
+    p += dimensionedScalar
+    (
+        "p",
+        p.dimensions(),
+        pRefValue - getRefCellValue(p, pRefCell)
+    );
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/stabilizeHUEqn.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/stabilizeHUEqn.H
new file mode 100644
index 000000000..93de1d4be
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/stabilizeHUEqn.H
@@ -0,0 +1,43 @@
+forAll(phases, phasei)
+{
+    phaseModel& phase = phases[phasei];
+    const volScalarField& alpha = phase;
+    const volScalarField& rhoAlpha = phase.rho();
+    volVectorField& UAlpha = phase.URef();
+
+    // Stabilize due to zero flux
+    const dimensionedScalar residualDivU(dimensionSet(1, -3, -1, 0, 0), vSmall);
+
+    HUEqns[phasei] +=
+        fvm::Sp
+        (
+            max
+            (
+                residualDivU - fvc::div(phase.alphaRhoPhi()),
+                dimensionedScalar(dimensionSet(1, -3, -1, 0, 0), Zero)
+            ),
+            UAlpha
+        )
+      - fvc::Sp
+        (
+            max
+            (
+                residualDivU - fvc::div(phase.alphaRhoPhi()),
+                dimensionedScalar(dimensionSet(1, -3, -1, 0, 0), Zero)
+            ),
+            UAlpha
+        );
+
+    // Stabilize due to zero phase fraction
+    HUEqns[phasei] +=
+        fvm::ddt
+        (
+            max(phase.residualAlpha() - alpha, scalar(0))*rhoAlpha,
+            UAlpha
+        )
+      - fvc::ddt
+        (
+            max(phase.residualAlpha() - alpha, scalar(0))*rhoAlpha,
+            UAlpha
+        );
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/vmDiag.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/vmDiag.H
new file mode 100644
index 000000000..b209496ed
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/vmDiag.H
@@ -0,0 +1,16 @@
+initZeroPtrList(vmDiag, mesh, dimensionSet(1, -3, -1, 0, 0), "vmDiag");
+initZeroPtrList(vmDiagf, mesh, dimensionSet(1, -3, -1, 0, 0), "vmDiagf");
+
+forAllConstIter(multiphaseSystem::VmTable, fluid.Vms(), VmIter)
+{
+    const volScalarField& Vm(*VmIter());
+    const phasePair& pair(fluid.phasePairs()[VmIter.key()]);
+
+    forAllConstIter(phasePair, pair, iter)
+    {
+        const phaseModel& phase = iter();
+
+        vmDiag[phase.index()] += byDt(Vm);
+        vmDiagf[phase.index()] += byDt(fvc::interpolate(Vm));
+    }
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSource.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSource.H
new file mode 100644
index 000000000..f670aee1a
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSource.H
@@ -0,0 +1,22 @@
+initZeroPtrList(vmSource, mesh, dimensionSet(1, -2, -2, 0, 0), "vmSource");
+
+forAllConstIter(multiphaseSystem::VmTable, fluid.Vms(), VmIter)
+{
+    const volScalarField& Vm(*VmIter());
+    const phasePair& pair(fluid.phasePairs()[VmIter.key()]);
+
+    forAllConstIter(phasePair, pair, iter)
+    {
+        const phaseModel& phase = iter();
+        const phaseModel& otherPhase = iter.otherPhase();
+
+        vmSource[phase.index()] +=
+            Vm
+           *byDt
+            (
+                phase.U()().oldTime()
+              + otherPhase.U()
+              - otherPhase.U()().oldTime()
+            );
+    }
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSourcef.H b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSourcef.H
new file mode 100644
index 000000000..b107f94a5
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSourcef.H
@@ -0,0 +1,23 @@
+initZeroPtrList(vmSourcef, mesh, dimensionSet(0, 3, -1, 0, 0), "vmSourcef");
+
+forAllConstIter(multiphaseSystem::VmTable, fluid.Vms(), VmIter)
+{
+    const volScalarField& Vm(*VmIter());
+    const phasePair& pair(fluid.phasePairs()[VmIter.key()]);
+
+    forAllConstIter(phasePair, pair, iter)
+    {
+        const phaseModel& phase = iter();
+        const phaseModel& otherPhase = iter.otherPhase();
+
+        vmSourcef[phase.index()] +=
+            rHUfs[phase.index()]
+           *fvc::interpolate(Vm)
+           *byDt
+            (
+                phase.phi()().oldTime()
+              + otherPhase.phi()
+              - otherPhase.phi()().oldTime()
+            );
+    }
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/Make/files b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/Make/files
index a40b522c9..ed5796496 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/Make/files
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/Make/files
@@ -18,6 +18,8 @@ dragModels/IshiiZuber/IshiiZuber.C
 dragModels/AttouFerschneider/AttouFerschneider.C
 dragModels/aerosolDrag/aerosolDrag.C
 dragModels/timeScaleFilteredDrag/timeScaleFilteredDrag.C
+dragModels/IshiiZuberSparse/IshiiZuberSparse.C
+dragModels/StrubeljTiselj/StrubeljTiselj.C
 
 swarmCorrections/swarmCorrection/swarmCorrection.C
 swarmCorrections/swarmCorrection/newSwarmCorrection.C
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.C
new file mode 100644
index 000000000..412f482d2
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.C
@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 2018-2020 OpenFOAM Foundation
+  ==  == ====== ====   ==  ==  | and Helmholtz-Zentrum Dresden-Rossendorf
+-------------------------------------------------------------------------------
+License
+    This file is a derivative work 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 "IshiiZuberSparse.H"
+#include "phasePair.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace dragModels
+{
+    defineTypeNameAndDebug(IshiiZuberSparse, 0);
+    addToRunTimeSelectionTable(dragModel, IshiiZuberSparse, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dragModels::IshiiZuberSparse::IshiiZuberSparse
+(
+    const dictionary& dict,
+    const phasePair& pair,
+    const bool registerObject
+)
+:
+    dragModel(dict, pair, registerObject)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::dragModels::IshiiZuberSparse::~IshiiZuberSparse()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::dragModels::IshiiZuberSparse::CdRe() const
+{
+    const volScalarField Re(pair_.Re());
+    const volScalarField Eo(pair_.Eo());
+
+    const volScalarField CdRe
+    (
+        24*(1 + 0.1*pow(Re, 0.75))
+    );
+
+    const volScalarField CdReEllipse(0.6666*sqrt(Eo)*Re);
+
+    return
+        pos0(CdReEllipse - CdRe)
+       *min(CdReEllipse, Re*2.66667)
+      + neg(CdReEllipse - CdRe)*CdRe;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.H b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.H
new file mode 100644
index 000000000..5bc44b029
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.H
@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 2018-2020 OpenFOAM Foundation
+  ==  == ====== ====   ==  ==  | and Helmholtz-Zentrum Dresden-Rossendorf
+-------------------------------------------------------------------------------
+License
+    This file is a derivative work 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::dragModels::IshiiZuberSparse
+
+Description
+    Ishii and Zuber (1979) drag model for sparse dispersed bubbly flows. Due to
+    the simplification of a sparse flow regime no correction for swarm and/or
+    viscous effects in the elliptical and cap shape regime is included. To
+    account for such effects add a swarmCorrection model or use IshiiZuber for
+    dense dispersed bubbly flows.
+
+    Reference:
+    \verbatim
+        Ishii, M., & Zuber, N. (1979).
+        Drag coefficient and relative velocity in bubbly, droplet or particulate
+        flows.
+        AIChE Journal, 25(5), 843-855.
+    \endverbatim
+
+SourceFiles
+    IshiiZuberSparse.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IshiiZuberSparse_H
+#define IshiiZuberSparse_H
+
+#include "dragModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace dragModels
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class IshiiZuberSparse Declaration
+\*---------------------------------------------------------------------------*/
+
+class IshiiZuberSparse
+:
+    public dragModel
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("IshiiZuberSparse");
+
+
+    // Constructors
+
+        //- Construct from a dictionary and a phase pair
+        IshiiZuberSparse
+        (
+            const dictionary& dict,
+            const phasePair& pair,
+            const bool registerObject
+        );
+
+
+    //- Destructor
+    virtual ~IshiiZuberSparse();
+
+
+    // Member Functions
+
+        //- Drag coefficient
+        virtual tmp<volScalarField> CdRe() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace dragModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.C
new file mode 100644
index 000000000..aa1c8f26b
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.C
@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 2018-2020 OpenFOAM Foundation
+  ==  == ====== ====   ==  ==  | and Helmholtz-Zentrum Dresden-Rossendorf
+-------------------------------------------------------------------------------
+License
+    This file is a derivative work 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 "StrubeljTiselj.H"
+#include "phasePair.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceInterpolate.H"
+#include "fvcGrad.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace dragModels
+{
+    defineTypeNameAndDebug(StrubeljTiselj, 0);
+    addToRunTimeSelectionTable(dragModel, StrubeljTiselj, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dragModels::StrubeljTiselj::StrubeljTiselj
+(
+    const dictionary& dict,
+    const phasePair& pair,
+    const bool registerObject
+)
+:
+    dragModel(dict, pair, registerObject),
+    tau_("Tau", dimless, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::dragModels::StrubeljTiselj::~StrubeljTiselj()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::dragModels::StrubeljTiselj::Ki() const
+{
+    return
+        max(pair_.phase1(), pair_.phase1().residualAlpha())
+       *max(pair_.phase2(), pair_.phase2().residualAlpha())
+       *pair_.rho()*tau_/pair_.phase1().mesh().time().deltaT();
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::dragModels::StrubeljTiselj::K() const
+{
+    return Ki();
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>Foam::dragModels::StrubeljTiselj::Kf() const
+{
+    return fvc::interpolate(Ki());
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::dragModels::StrubeljTiselj::CdRe() const
+{
+    notImplemented("dragModels::StrubeljTiselj::CdRe() const");
+
+    return tmp<volScalarField>(nullptr);
+}
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.H b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.H
new file mode 100644
index 000000000..6c5374cb8
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.H
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 2018-2020 OpenFOAM Foundation
+  ==  == ====== ====   ==  ==  | and Helmholtz-Zentrum Dresden-Rossendorf
+-------------------------------------------------------------------------------
+License
+    This file is a derivative work 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::dragModels::StrubeljTiselj
+
+Description
+    Numerical drag for handling of resolved interfaces in multi-fluid model
+
+    \verbatim
+        Strubelj, L. and Tiselj, I. 'Two-fluid model with interface sharpening'
+        International Journal for Numerical Methods in Engineering, Vol. 85,
+        pp. 575-590 (2011)
+    \endverbatim
+
+SourceFiles
+    StrubeljTiselj.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef StrubeljTiselj_H
+#define StrubeljTiselj_H
+
+#include "dragModel.H"
+#include "phaseModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace dragModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class StrubeljTiselj Declaration
+\*---------------------------------------------------------------------------*/
+
+class StrubeljTiselj
+:
+    public dragModel
+{
+    // Private Data
+
+        //- Ratio of physical time step to relaxation time
+        const dimensionedScalar tau_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("StrubeljTiselj");
+
+
+    // Constructors
+
+        //- Construct from a dictionary and a phase pair
+        StrubeljTiselj
+        (
+            const dictionary& dict,
+            const phasePair& pair,
+            const bool registerObject
+        );
+
+
+    //- Destructor
+    virtual ~StrubeljTiselj();
+
+
+    // Member Functions
+
+        //- Return the phase-intensive drag coefficient Ki
+        //  used in the momentum equations
+        //    ddt(alpha1*rho1*U1) + ... = ... alphad*K*(U1-U2)
+        //    ddt(alpha2*rho2*U2) + ... = ... alphad*K*(U2-U1)
+        virtual tmp<volScalarField> Ki() const;
+
+        //- Return the drag coefficient K as a scalar
+        //  used in the momentum equations
+        //    ddt(alpha1*rho1*U1) + ... = ... K*(U1-U2)
+        //    ddt(alpha2*rho2*U2) + ... = ... K*(U2-U1)
+        virtual tmp<volScalarField> K() const;
+
+        //- Return the drag coefficient Kf
+        //  used in the face-momentum equations
+        virtual tmp<surfaceScalarField> Kf() const;
+
+        //- Drag coefficient
+        virtual tmp<volScalarField> CdRe() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace dragModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C
index abc7e3a52..299bb2e27 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -264,7 +264,7 @@ Foam::BlendedInterfacialModel<ModelType>::evaluate
 
     if (model2In1_.valid())
     {
-        addToXs(f2, (model1In2_().*method)(args ...));
+        addToXs(f2, (model2In1_().*method)(args ...));
     }
 
     if
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.C
new file mode 100644
index 000000000..d3b7e5dc3
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.C
@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 2019-2020 OpenFOAM Foundation
+  ==  == ====== ====   ==  ==  | and Helmholtz-Zentrum Dresden-Rossendorf
+-------------------------------------------------------------------------------
+License
+    This file is a derivative work 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 "definitePhases.H"
+#include "phaseSystem.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace blendingMethods
+{
+    defineTypeName(definitePhases);
+
+    addToRunTimeSelectionTable
+    (
+        blendingMethod,
+        definitePhases,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::blendingMethods::definitePhases::definitePhases
+(
+    const dictionary& dict,
+    const wordList& phaseNames
+)
+:
+    blendingMethod(dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::blendingMethods::definitePhases::~definitePhases()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::blendingMethods::definitePhases::f1
+(
+    const phaseModel& phase1,
+    const phaseModel& phase2
+) const
+{
+    if(phase1.fluid().continuousPhases().empty())
+    {
+        FatalErrorInFunction
+            << "The definitePhases blending requires at least one continuous "
+            << "phase, specified by keyword 'continuous' in phaseProperties."
+            << exit(FatalError);
+    }
+
+    const word phases12(IOobject::groupName(phase1.name(), phase2.name()));
+
+    tmp<volScalarField> f1
+    (
+        volScalarField::New
+        (
+            IOobject::groupName("f1."+typeName, phases12),
+            (!(phase1.continuous()))
+           *phase2
+           /max(phase2.fluid().sumAlphaContinuous(), phase2.residualAlpha())
+        )
+    );
+
+    if (debug && phase1.mesh().time().outputTime())
+    {
+        f1().write();
+    }
+
+    return f1;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::blendingMethods::definitePhases::f2
+(
+    const phaseModel& phase1,
+    const phaseModel& phase2
+) const
+{
+    if(phase1.fluid().continuousPhases().empty())
+    {
+        FatalErrorInFunction
+        << "The definitePhases blending requires at least one continuous"
+        << "phase, specified by keyword 'continuous' in phaseProperties."
+            << exit(FatalError);
+    }
+
+    const word phases12(IOobject::groupName(phase1.name(), phase2.name()));
+
+    tmp<volScalarField> f2
+    (
+        volScalarField::New
+        (
+            IOobject::groupName("f2."+typeName, phases12),
+            (!(phase2.continuous()))
+           *phase1
+           /max(phase1.fluid().sumAlphaContinuous(), phase1.residualAlpha())
+        )
+    );
+
+    if (debug && phase1.mesh().time().outputTime())
+    {
+        f2().write();
+    }
+
+    return f2;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.H
new file mode 100644
index 000000000..681c90c63
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.H
@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 2019-2020 OpenFOAM Foundation
+  ==  == ====== ====   ==  ==  | and Helmholtz-Zentrum Dresden-Rossendorf
+-------------------------------------------------------------------------------
+License
+    This file is a derivative work 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::blendingMethods::definitePhases
+
+Description
+    This blending method is intended to work similar to the 'none' method, but
+    for an arbitrary number of disperse and of continuous phases.
+
+    For every phase, the function 'isDisperse' indicates, whether this phase is
+    disperse (1) or not (0). This information is taken from the list called
+    continuousPhases, which has to be provided by the user.
+
+    For a pair of continuous phases, the blending coefficient for any kind of
+    unordered phase pair model (e.g. '(water and oil)') is always equal to 1,
+    whereas f1 and f2 are always zero for such a pair.
+
+    For a pair of a continuous and a disperse phase, the blending coefficient
+    (f1 or f2) corresponding to the ordered phase pair model '(disperse in
+    continuous)' is one, wherever this specific continuous phase is present and
+    no other continuous phase is. This blending coefficient becomes zero, where
+    this specific continuous phase is not present at al. In regions, where
+    multiple continuous phases are present, the blending coefficient is
+    proportional to the ratio of the volume fraction of this specific continuous
+    phase and the sum of volume fractions of all continuous phases at this
+    point.
+
+    Blending factor f1 is corresponding to an ordered phase pair model
+    (phase1 in phase2). f2 is defined vice versa, accordingly
+    \f[
+        f_1 = \mathit{isDisperse}(\mathrm{phase1})
+            \frac{\alpha_2}{\alpha_{\Sigma c}}
+    \f]
+    \f[
+        f_2 = \mathit{isDisperse}(\mathrm{phase2})
+            \frac{\alpha_1}{\alpha_{\Sigma c}}
+    \f]
+
+    with the sum volume fraction of all phases declared as continuous
+    \f[
+        \alpha_{\Sigma c} = \Sigma_{i\in c} \alpha_i
+    \f]
+
+    Note, that the blending factors could be stored during runtime by enabling
+    the debug switch 'blendingMethod'. Calls to functions f1() and f2() may be
+    enforced with a debug switch for a specific model type,
+    e.g. 'BlendedInterfacialModel<dragModel> '.
+
+Usage
+    \table
+        Name              | Description                     | Example
+        continuous phases | declare phases to be continuous | (water oil)
+                          | under all circumstances         |
+    \endtable
+
+    Example for a phase system with an arbitrary number of disperse air phases
+    and both continuous phases of water and oil
+
+SourceFiles
+    definitePhases.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef definitePhases_H
+#define definitePhases_H
+
+#include "blendingMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace blendingMethods
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class definitePhases Declaration
+\*---------------------------------------------------------------------------*/
+
+class definitePhases
+:
+    public blendingMethod
+{
+public:
+
+    //- Runtime type information
+    TypeNameNoDebug("definitePhases");
+
+
+    // Constructors
+
+        //- Construct from a dictionary and two phases
+        definitePhases
+        (
+            const dictionary& dict,
+            const wordList& phaseNames
+        );
+
+
+    //- Destructor
+    ~definitePhases();
+
+
+    // Member Functions
+
+        //- Factor for primary phase
+        virtual tmp<volScalarField> f1
+        (
+            const phaseModel& phase1,
+            const phaseModel& phase2
+        ) const;
+
+        //- Factor for secondary phase
+        virtual tmp<volScalarField> f2
+        (
+            const phaseModel& phase1,
+            const phaseModel& phase2
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace blendingMethods
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/Make/files b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/Make/files
index 9e3385057..787240118 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/Make/files
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/Make/files
@@ -76,6 +76,7 @@ BlendedInterfacialModel/blendingMethods/blendingMethod/newBlendingMethod.C
 BlendedInterfacialModel/blendingMethods/noBlending/noBlending.C
 BlendedInterfacialModel/blendingMethods/linear/linear.C
 BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.C
+BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.C
 
 reactionThermo/hRefConstThermos.C
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index 9df985ec0..f637286e9 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -176,38 +176,185 @@ MomentumTransferPhaseSystem
     const fvMesh& mesh
 )
 :
-    BasePhaseSystem(mesh)
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+    BasePhaseSystem(mesh),
+    correctFixedFluxBCs_
+    (
+        mesh.solutionDict().subDict("PIMPLE").lookupOrDefault<Switch>
+        (
+            "correctFixedFluxBCs",
+            true
+        )
+    )
 {
     this->generatePairsAndSubModels
     (
         "drag",
-        dragModels_
+        dragModels_,
+        correctFixedFluxBCs_
     );
 
     this->generatePairsAndSubModels
     (
         "virtualMass",
-        virtualMassModels_
+        virtualMassModels_,
+        correctFixedFluxBCs_
     );
 
     this->generatePairsAndSubModels
     (
         "lift",
-        liftModels_
+        liftModels_,
+        correctFixedFluxBCs_
     );
 
     this->generatePairsAndSubModels
     (
         "wallLubrication",
-        wallLubricationModels_
+        wallLubricationModels_,
+        correctFixedFluxBCs_
     );
 
     this->generatePairsAndSubModels
     (
         "turbulentDispersion",
-        turbulentDispersionModels_
+        turbulentDispersionModels_,
+        correctFixedFluxBCs_
     );
 
+    forAllConstIter
+    (
+        liftModelTable,
+        liftModels_,
+        liftModelIter
+    )
+    {
+        const phasePair& pair(this->phasePairs_[liftModelIter.key()]);
+
+        Fls_.insert
+        (
+            pair,
+            new volVectorField
+            (
+                IOobject
+                (
+                    IOobject::groupName("Fl", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh()
+                ),
+                this->mesh(),
+                dimensionedVector(liftModel::dimF, Zero)
+            )
+        );
+
+        Flfs_.insert
+        (
+            pair,
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("Flf", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh()
+                ),
+                this->mesh(),
+                dimensionedScalar(liftModel::dimF*dimArea, 0)
+            )
+        );
+    }
+
+    forAllConstIter
+    (
+        wallLubricationModelTable,
+        wallLubricationModels_,
+        wallLubricationModelIter
+    )
+    {
+        const phasePair& pair
+        (
+            this->phasePairs_[wallLubricationModelIter.key()]
+        );
+
+        Fws_.insert
+        (
+            pair,
+            new volVectorField
+            (
+                IOobject
+                (
+                    IOobject::groupName("Fw", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh()
+                ),
+                this->mesh(),
+                dimensionedVector(wallLubricationModel::dimF, Zero)
+            )
+        );
+
+        Fwfs_.insert
+        (
+            pair,
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("Fwf", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh()
+                ),
+                this->mesh(),
+                dimensionedScalar(wallLubricationModel::dimF*dimArea, 0)
+            )
+        );
+    }
+
+    forAllConstIter
+    (
+        turbulentDispersionModelTable,
+        turbulentDispersionModels_,
+        turbulentDispersionModelIter
+    )
+    {
+        const phasePair& pair
+        (
+            this->phasePairs_[turbulentDispersionModelIter.key()]
+        );
+
+        Fts_.insert
+        (
+            pair,
+            new volVectorField
+            (
+                IOobject
+                (
+                    IOobject::groupName("Ft", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh()
+                ),
+                this->mesh(),
+                dimensionedVector(turbulentDispersionModel::dimF, Zero)
+            )
+        );
+
+        Ftfs_.insert
+        (
+            pair,
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("Ftf", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh()
+                ),
+                this->mesh(),
+                dimensionedScalar(turbulentDispersionModel::dimF*dimArea, 0)
+            )
+        );
+    }
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
     forAllConstIter
     (
         dragModelTable,
@@ -399,6 +546,48 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer()
         }
     }
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+    // Update the lift forces
+    forAllConstIter
+    (
+        liftModelTable,
+        liftModels_,
+        liftModelIter
+    )
+    {
+        *Fls_[liftModelIter.key()] = liftModelIter()->F<vector>();
+        *Flfs_[liftModelIter.key()] = liftModelIter()->Ff();
+    }
+
+    // Update the wall lubrication forces
+    forAllConstIter
+    (
+        wallLubricationModelTable,
+        wallLubricationModels_,
+        wallLubricationModelIter
+    )
+    {
+        *Fws_[wallLubricationModelIter.key()] =
+            wallLubricationModelIter()->F<vector>();
+        *Fwfs_[wallLubricationModelIter.key()] =
+            wallLubricationModelIter()->Ff();
+    }
+
+    // Update the turbulent dispersion coefficients
+    forAllConstIter
+    (
+        turbulentDispersionModelTable,
+        turbulentDispersionModels_,
+        turbulentDispersionModelIter
+    )
+    {
+        *Fts_[turbulentDispersionModelIter.key()] =
+            turbulentDispersionModelIter()->F<vector>();
+        *Ftfs_[turbulentDispersionModelIter.key()] =
+            turbulentDispersionModelIter()->Ff();
+    }
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
     return eqnsPtr;
 }
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
index be5d48ea1..8966e3a78 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -134,8 +134,83 @@ private:
         > turbulentDispersionModelTable;
 
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+protected:
+
+    // Protected typedefs
+
+        typedef HashPtrTable
+        <
+            volVectorField,
+            phasePairKey,
+            phasePairKey::hash
+        > FlTable;
+
+        typedef HashPtrTable
+        <
+            surfaceScalarField,
+            phasePairKey,
+            phasePairKey::hash
+        > FlfTable;
+
+        typedef HashPtrTable
+        <
+            volVectorField,
+            phasePairKey,
+            phasePairKey::hash
+        > FwTable;
+
+        typedef HashPtrTable
+        <
+            surfaceScalarField,
+            phasePairKey,
+            phasePairKey::hash
+        > FwfTable;
+
+        typedef HashPtrTable
+        <
+            volVectorField,
+            phasePairKey,
+            phasePairKey::hash
+        > FtTable;
+
+        typedef HashPtrTable
+        <
+            surfaceScalarField,
+            phasePairKey,
+            phasePairKey::hash
+        > FtfTable;
+
+
+private:
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
+
     // Private Data
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+        //- Switch for correctFixedFluxBCs in BlendedInterfacialModel
+        Switch correctFixedFluxBCs_;
+
+        //- Lift forces
+        FlTable Fls_;
+
+        //- Face lift forces
+        FlfTable Flfs_;
+
+        //- Wall lubrication forces
+        FwTable Fws_;
+
+        //- Face wall lubrication forces
+        FwfTable Fwfs_;
+
+        //- Turbulent dispersion forces
+        FtTable Fts_;
+
+        //- Face Turbulent dispersion forces
+        FtfTable Ftfs_;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
         //- Drag coefficients
         KdTable Kds_;
 
@@ -299,6 +374,68 @@ public:
             const PtrList<surfaceScalarField>& rAUfs
         ) const;
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+        //- Return the drag coefficient for all phase-pairs
+        virtual const KdTable& Kds() const
+        {
+            return Kds_;
+        };
+
+        //- Return the face drag coefficient for all phase-pairs
+        virtual const KdfTable& Kdfs() const
+        {
+            return Kdfs_;
+        };
+
+        //- Return the virtual mass force coefficient for all phase-pairs
+        virtual const VmTable& Vms() const
+        {
+            return Vms_;
+        };
+
+        //- Return the virtual face mass force coefficient for all phase-pairs
+        virtual const VmfTable& Vmfs() const
+        {
+            return Vmfs_;
+        };
+
+        //- Return the lift force for all phase-pairs
+        virtual const FlTable& Fls() const
+        {
+            return Fls_;
+        };
+
+        //- Return the face lift force for all phase-pairs
+        virtual const FlfTable& Flfs() const
+        {
+            return Flfs_;
+        };
+
+        //- Return the wall lubrication force for all phase-pairs
+        virtual const FwTable& Fws() const
+        {
+            return Fws_;
+        };
+
+        //- Return the face wall lubrication force for all phase-pairs
+        virtual const FwfTable& Fwfs() const
+        {
+            return Fwfs_;
+        };
+
+        //- Return the turbulent dispersion force for all phase-pairs
+        virtual const FtTable& Fts() const
+        {
+            return Fts_;
+        };
+
+        //- Return the face turbulent dispersion force for all phase-pairs
+        virtual const FtfTable& Ftfs() const
+        {
+            return Ftfs_;
+        };
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
         //- Read base phaseProperties dictionary
         virtual bool read();
 };
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
index 8b2069d67..2445f1ce1 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -559,4 +559,14 @@ Foam::MovingPhaseModel<BasePhaseModel>::pPrime() const
 }
 
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+template<class BasePhaseModel>
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::MovingPhaseModel<BasePhaseModel>::divDevRhoReff(volVectorField& U) const
+{
+    return turbulence_->divDevRhoReff(U);
+}
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
index 546b62da3..5697fbbbf 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -246,6 +246,11 @@ public:
             //- Return the phase-pressure'
             //  (derivative of phase-pressure w.r.t. phase-fraction)
             virtual tmp<volScalarField> pPrime() const;
+
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+            //- Return the turbulent term for the momentum equation
+            virtual tmp<fvVectorMatrix> divDevRhoReff(volVectorField& U) const;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
 };
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C
index 7c6e9c368..fda4b3407 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -377,4 +377,21 @@ Foam::StationaryPhaseModel<BasePhaseModel>::pPrime() const
 }
 
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+template<class BasePhaseModel>
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::StationaryPhaseModel<BasePhaseModel>::divDevRhoReff
+(
+    volVectorField& U
+) const
+{
+    FatalErrorInFunction
+        << "Cannot access the turbulent term of a stationary phase"
+        << exit(FatalError);
+
+    return fvm::Su(U, U);
+}
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H
index 6279e0cef..91a93c2b4 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -178,6 +178,11 @@ public:
             //- Return the phase-pressure'
             //  (derivative of phase-pressure w.r.t. phase-fraction)
             virtual tmp<volScalarField> pPrime() const;
+
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+            //- Return the turbulent term for the momentum equation
+            virtual tmp<fvVectorMatrix> divDevRhoReff(volVectorField& U) const;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
 };
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
index e6c030632..e478f6c5f 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,7 +68,10 @@ Foam::phaseModel::phaseModel
         dimless,
         fluid.subDict(phaseName).lookup("residualAlpha")
     ),
-    alphaMax_(fluid.subDict(phaseName).lookupOrDefault("alphaMax", 1.0))
+    alphaMax_(fluid.subDict(phaseName).lookupOrDefault("alphaMax", 1.0)),
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+    continuous_(fluid.subDict(phaseName).lookupOrDefault("continuous", false))
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
 {
     diameterModel_ = diameterModel::New(fluid.subDict(phaseName), *this);
 }
@@ -125,6 +128,14 @@ Foam::scalar Foam::phaseModel::alphaMax() const
 }
 
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+bool Foam::phaseModel::continuous() const
+{
+    return continuous_;
+}
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
+
 Foam::tmp<Foam::volScalarField> Foam::phaseModel::d() const
 {
     return diameterModel_().d();
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
index 77db9e7b8..60c1fad1f 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -75,6 +75,11 @@ class phaseModel
         //- Optional maximum phase-fraction (e.g. packing limit)
         scalar alphaMax_;
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+        //- Morphology of phase (default is false, i.e. dispersed)
+        const Switch continuous_;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
         //- Diameter model
         autoPtr<diameterModel> diameterModel_;
 
@@ -177,6 +182,11 @@ public:
         //- Return the maximum phase-fraction (e.g. packing limit)
         scalar alphaMax() const;
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+        //- Return true if phase is continuous
+        bool continuous() const;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
         //- Return the Sauter-mean diameter
         tmp<volScalarField> d() const;
 
@@ -417,6 +427,14 @@ public:
             //- Return the phase-pressure'
             //  (derivative of phase-pressure w.r.t. phase-fraction)
             virtual tmp<volScalarField> pPrime() const = 0;
+
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+            //- Return the turbulent term for the momentum equation
+            virtual tmp<fvVectorMatrix> divDevRhoReff
+            (
+                volVectorField& U
+            ) const = 0;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
 };
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePair.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePair.C
index 079f95714..e16ab2306 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePair.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePair.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,7 +104,26 @@ Foam::word Foam::phasePair::otherName() const
 
 Foam::tmp<Foam::volScalarField> Foam::phasePair::rho() const
 {
-    return phase1()*phase1().rho() + phase2()*phase2().rho();
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+    volScalarField rhoMin(min(phase1().rho(), phase2().rho()));
+    volScalarField rhoMax(max(phase1().rho(), phase2().rho()));
+
+    return
+        min
+        (
+            max
+            (
+                phase1()*phase1().rho() + phase2()*phase2().rho()
+               /max
+                (
+                    phase1() + phase2(),
+                    phase1().residualAlpha() + phase2().residualAlpha()
+                ),
+                rhoMin
+            ),
+            rhoMax
+        );
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
 }
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
index 81d474936..98dc8fe79 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -163,6 +163,9 @@ Foam::phaseSystem::phaseSystem
     label stationaryPhasei = 0;
     label anisothermalPhasei = 0;
     label multiComponentPhasei = 0;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+    label continuousPhasei = 0;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
     forAll(phaseModels_, phasei)
     {
         phaseModel& phase = phaseModels_[phasei];
@@ -170,16 +173,25 @@ Foam::phaseSystem::phaseSystem
         stationaryPhasei += phase.stationary();
         anisothermalPhasei += !phase.isothermal();
         multiComponentPhasei += !phase.pure();
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+        continuousPhasei += phase.continuous();
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
     }
     movingPhaseModels_.resize(movingPhasei);
     stationaryPhaseModels_.resize(stationaryPhasei);
     anisothermalPhaseModels_.resize(anisothermalPhasei);
     multiComponentPhaseModels_.resize(multiComponentPhasei);
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+    continuousPhaseModels_.resize(continuousPhasei);
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
 
     movingPhasei = 0;
     stationaryPhasei = 0;
     anisothermalPhasei = 0;
     multiComponentPhasei = 0;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+    continuousPhasei = 0;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
     forAll(phaseModels_, phasei)
     {
         phaseModel& phase = phaseModels_[phasei];
@@ -199,7 +211,36 @@ Foam::phaseSystem::phaseSystem
         {
             multiComponentPhaseModels_.set(multiComponentPhasei ++, &phase);
         }
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+        if (phase.continuous())
+        {
+            continuousPhaseModels_.set(continuousPhasei ++, &phase);
+        }
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+    }
+
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+    if (!continuousPhaseModels_.empty())
+    {
+        Info<< "Continuous phases:";
+        forAll(continuousPhaseModels_, phasei)
+        {
+            Info<< " " << continuousPhaseModels_[phasei].name();
+        }
+        Info<< endl;
+    }
+
+    // Detect incompressibility
+    incompressible_ = true;
+    forAll(phaseModels_, phasei)
+    {
+        if (!phaseModels_[phasei].incompressible())
+        {
+            incompressible_ = false;
+            break;
+        }
     }
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
 
     // Write phi
     phi_.writeOpt() = IOobject::AUTO_WRITE;
@@ -263,6 +304,29 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::rho() const
 }
 
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+Foam::tmp<Foam::volScalarField> Foam::phaseSystem::sumAlphaContinuous() const
+{
+    const label nPhases = movingPhaseModels_.size();
+
+    tmp<volScalarField> alphaCont
+    (
+         movingPhaseModels_[0].continuous()
+        *movingPhaseModels_[0]
+    );
+
+    for (label phasei = 1; phasei < nPhases; ++ phasei)
+    {
+        alphaCont.ref() +=
+            movingPhaseModels_[phasei].continuous()
+           *movingPhaseModels_[phasei];
+    }
+
+    return alphaCont;
+}
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
+
 Foam::tmp<Foam::volVectorField> Foam::phaseSystem::U() const
 {
     const label nMovingPhases = movingPhaseModels_.size();
@@ -371,6 +435,14 @@ Foam::PtrList<Foam::volScalarField> Foam::phaseSystem::dmdts() const
 }
 
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+bool Foam::phaseSystem::incompressible() const
+{
+    return incompressible_;
+}
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
+
 bool Foam::phaseSystem::implicitPhasePressure(const phaseModel& phase) const
 {
     return false;
@@ -494,4 +566,19 @@ Foam::tmp<Foam::surfaceScalarField> Foam::byDt(const surfaceScalarField& sf)
 }
 
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+Foam::tmp<Foam::volVectorField> Foam::byDt(const volVectorField& vf)
+{
+    if (fv::localEulerDdt::enabled(vf.mesh()))
+    {
+        return fv::localEulerDdt::localRDeltaT(vf.mesh())*vf;
+    }
+    else
+    {
+        return vf/vf.mesh().time().deltaT();
+    }
+}
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
index 776e4dc6c..504787b4c 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -153,6 +153,11 @@ protected:
         //- Multi-component phase models
         phaseModelPartialList multiComponentPhaseModels_;
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+        //- Continuous phase models
+        phaseModelPartialList continuousPhaseModels_;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
         //- Phase pairs
         phasePairTable phasePairs_;
 
@@ -168,6 +173,11 @@ protected:
         //- Blending methods
         blendingMethodTable blendingMethods_;
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+        //- Incompressibility
+        bool incompressible_;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
 
         // Sub Models
 
@@ -304,6 +314,11 @@ public:
             //- Access the models for phases that have multiple species
             inline phaseModelPartialList& multiComponentPhases();
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+            //- Return the models for phases that are continuous
+            inline const phaseModelPartialList& continuousPhases() const;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
             //- Return the phase pairs
             inline const phasePairTable& phasePairs() const;
 
@@ -399,6 +414,11 @@ public:
             //- Return the mixture density
             tmp<volScalarField> rho() const;
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+            //- Return the total volume fraction of continuous phases
+            tmp<volScalarField> sumAlphaContinuous() const;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
             //- Return the mixture velocity
             tmp<volVectorField> U() const;
 
@@ -421,6 +441,11 @@ public:
             //- Return the mass transfer rates for each phase
             virtual PtrList<volScalarField> dmdts() const;
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+            //- Return incompressibility
+            bool incompressible() const;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
 
         // Transfers
 
@@ -544,6 +569,9 @@ public:
 
 tmp<volScalarField> byDt(const volScalarField& vf);
 tmp<surfaceScalarField> byDt(const surfaceScalarField& sf);
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+tmp<volVectorField> byDt(const volVectorField& vf);
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
index 906d013cd..7c1606f03 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -101,6 +101,15 @@ Foam::phaseSystem::multiComponentPhases()
 }
 
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+inline const Foam::phaseSystem::phaseModelPartialList&
+Foam::phaseSystem::continuousPhases() const
+{
+    return continuousPhaseModels_;
+}
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
+
 inline const Foam::phaseSystem::phasePairTable&
 Foam::phaseSystem::phasePairs() const
 {
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
index ef121a10a..6074dc518 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,6 +53,103 @@ class multiphaseSystem
 :
     public phaseSystem
 {
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+public:
+
+    // Public Typedefs
+
+        typedef
+            HashPtrTable
+            <
+                volScalarField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            KdTable;
+
+        typedef
+            HashPtrTable
+            <
+                surfaceScalarField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            KdfTable;
+
+        typedef
+            HashPtrTable
+            <
+                volScalarField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            VmTable;
+
+        typedef
+            HashPtrTable
+            <
+                surfaceScalarField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            VmfTable;
+
+        typedef
+            HashPtrTable
+            <
+                volVectorField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            FlTable;
+
+        typedef
+            HashPtrTable
+            <
+                surfaceScalarField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            FlfTable;
+
+        typedef
+            HashPtrTable
+            <
+                volVectorField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            FwTable;
+
+        typedef
+            HashPtrTable
+            <
+                surfaceScalarField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            FwfTable;
+
+        typedef
+            HashPtrTable
+            <
+                volVectorField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            FtTable;
+
+        typedef
+            HashPtrTable
+            <
+                surfaceScalarField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            FtfTable;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
+
 private:
 
     // Private Typedefs
@@ -178,6 +275,39 @@ public:
         //  Field values are 1 near and 0 away for the interface.
         tmp<volScalarField> nearInterface() const;
 
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Begin HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+        //- Return the drag coefficient for all phase-pairs
+        virtual const KdTable& Kds() const = 0;
+
+        //- Return the face drag coefficient for all phase-pairs
+        virtual const KdfTable& Kdfs() const = 0;
+
+        //- Return the virtual mass force coefficient for all phase-pairs
+        virtual const VmTable& Vms() const = 0;
+
+        //- Return the face based virtual mass force coefficient
+        //  for all phase-pairs
+        virtual const VmfTable& Vmfs() const = 0;
+
+        //- Return the lift force for all phase-pairs
+        virtual const FlTable& Fls() const = 0;
+
+        //- Return the face lift force for all phase-pairs
+        virtual const FlfTable& Flfs() const = 0;
+
+        //- Return the wall lubrication force for all phase-pairs
+        virtual const FwTable& Fws() const = 0;
+
+        //- Return the face wall lubrication force for all phase-pairs
+        virtual const FwfTable& Fwfs() const = 0;
+
+        //- Return the turbulent dispersion force for all phase-pairs
+        virtual const FtTable& Fts() const = 0;
+
+        //- Return the face turbulent dispersion force for all phase-pairs
+        virtual const FtfTable& Ftfs() const = 0;
+//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End HZDR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
+
         //- Solve for the phase fractions
         virtual void solve
         (
-- 
2.25.0

