From 5bdc07e9339078f0198979b657b8538fe4141164 Mon Sep 17 00:00:00 2001
From: "Schlegel, Dr. Fabian (FWDC) - 107795" <f.schlegel@hzdr.de>
Date: Tue, 8 Dec 2020 15:31:26 +0100
Subject: [PATCH] multiphaseEulerFoam: Add cipsaMultiphaseEulerFoam solver

    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
    * Interfacial momentum and mass exchanges are handled consistently
    * Necessary infrastructure in phaseSystems is added
    * New drag models are added: StrubeljTiselj, IshiiZuberSparse
    * A new blending method is added: definitePhases
    * ResolvedPhaseSystem added: calculation of surfaceTension and
        interfaceCompression
    * InterfaceCompression becomes a blendedInterfacialModel
    * InterfaceIndicatorMethods added for detection of resolved interfaces
    * Added SGSdeltas: copy of LESdeltas, which can be accessed by
        interfaceIndicatorMethods
    * Support of PID time step control for cipsaMultiphaseEulerFoam solver
    * Mixture density of phasePair: formulation for arbitrary number of phases
---
 .../multiphase/multiphaseEulerFoam/Allwclean  |   1 +
 .../multiphase/multiphaseEulerFoam/Allwmake   |   1 +
 .../cipsaMultiphaseEulerFoam/CourantNo.H      |  60 +++
 .../cipsaMultiphaseEulerFoam/EEqns.H          |  61 +++
 .../cipsaMultiphaseEulerFoam/Make/files       |   3 +
 .../cipsaMultiphaseEulerFoam/Make/options     |  25 ++
 .../cipsaMultiphaseEulerFoam.C                | 135 +++++++
 .../createFieldRefs.H                         |  49 +++
 .../cipsaMultiphaseEulerFoam/createFields.H   |  69 ++++
 .../initZeroPtrList.H                         |  82 +++++
 .../cipsaMultiphaseEulerFoam/pU/HUEqn.H       |  84 +++++
 .../cipsaMultiphaseEulerFoam/pU/UEqns.H       | 234 ++++++++++++
 .../cipsaMultiphaseEulerFoam/pU/ddtUEqn.H     |  47 +++
 .../cipsaMultiphaseEulerFoam/pU/dmdtUEqn.H    |  74 ++++
 .../cipsaMultiphaseEulerFoam/pU/dragDiag.H    |  71 ++++
 .../cipsaMultiphaseEulerFoam/pU/dragSource.H  |  75 ++++
 .../cipsaMultiphaseEulerFoam/pU/dragSourcef.H |  84 +++++
 .../pU/forcesSourcef.H                        |  36 ++
 .../pU/fvOptionsEqn.H                         |  18 +
 .../cipsaMultiphaseEulerFoam/pU/pEqn.H        | 347 ++++++++++++++++++
 .../pU/stabilizeHUEqn.H                       |  43 +++
 .../cipsaMultiphaseEulerFoam/pU/vmDiag.H      |  16 +
 .../cipsaMultiphaseEulerFoam/pU/vmSource.H    |  22 ++
 .../cipsaMultiphaseEulerFoam/pU/vmSourcef.H   |  23 ++
 .../interfacialCompositionModels/Make/files   |   5 +
 .../SGSdeltas/SGSdelta/SGSdelta.C             |  93 +++++
 .../SGSdeltas/SGSdelta/SGSdelta.H             | 143 ++++++++
 .../cubeRootVolDelta/cubeRootVolDelta.C       | 126 +++++++
 .../cubeRootVolDelta/cubeRootVolDelta.H       | 119 ++++++
 .../SGSdeltas/maxDeltaxyz/maxDeltaxyz.C       | 131 +++++++
 .../SGSdeltas/maxDeltaxyz/maxDeltaxyz.H       | 124 +++++++
 .../oneDimensionalDelta/oneDimensionalDelta.C | 125 +++++++
 .../oneDimensionalDelta/oneDimensionalDelta.H | 112 ++++++
 .../interfacialModels/Make/files              |   2 +
 .../IshiiZuberSparse/IshiiZuberSparse.C       |  83 +++++
 .../IshiiZuberSparse/IshiiZuberSparse.H       | 105 ++++++
 .../StrubeljTiselj/StrubeljTiselj.C           |  94 +++++
 .../StrubeljTiselj/StrubeljTiselj.H           | 121 ++++++
 .../multiphaseSystems/multiphaseSystems.C     |  26 +-
 .../definitePhases/definitePhases.C           | 140 +++++++
 .../definitePhases/definitePhases.H           | 142 +++++++
 .../phaseSystems/Make/files                   |  16 +
 .../MomentumTransferPhaseSystem.C             | 118 ++++++
 .../MomentumTransferPhaseSystem.H             |  80 ++++
 .../PhaseTransferPhaseSystem.H                |   5 +
 .../ResolvedPhaseSystem/ResolvedPhaseSystem.C | 220 +++++++++++
 .../ResolvedPhaseSystem/ResolvedPhaseSystem.H | 141 +++++++
 .../constantInterfaceCompression.C            |  91 +++++
 .../constantInterfaceCompression.H            | 100 +++++
 .../interfaceCompressionModel.C               | 102 +++++
 .../interfaceCompressionModel.H               | 138 +++++++
 .../interfaceCompressionModelNew.C            |  60 +++
 .../noInterfaceCompression.C                  | 116 ++++++
 .../noInterfaceCompression.H                  |  96 +++++
 .../gradAlpha/gradAlpha.C                     | 107 ++++++
 .../gradAlpha/gradAlpha.H                     | 142 +++++++
 .../interfaceIndicatorMethod.C                | 118 ++++++
 .../interfaceIndicatorMethod.H                | 137 +++++++
 .../interfaceIndicatorMethodNew.C             | 112 ++++++
 .../useBlendingMethod/useBlendingMethod.C     | 106 ++++++
 .../useBlendingMethod/useBlendingMethod.H     | 134 +++++++
 .../continuumSurfaceForce.C                   | 241 ++++++++++++
 .../continuumSurfaceForce.H                   | 119 ++++++
 .../noSurfaceTensionForce.C                   | 102 +++++
 .../noSurfaceTensionForce.H                   |  96 +++++
 .../surfaceTensionForceModel.C                |  82 +++++
 .../surfaceTensionForceModel.H                | 129 +++++++
 .../surfaceTensionForceModelNew.C             |  78 ++++
 .../MovingPhaseModel/MovingPhaseModel.C       |   8 +
 .../MovingPhaseModel/MovingPhaseModel.H       |   3 +
 .../StationaryPhaseModel.C                    |  15 +
 .../StationaryPhaseModel.H                    |   3 +
 .../phaseModel/phaseModel/phaseModel.C        |   9 +-
 .../phaseModel/phaseModel/phaseModel.H        |  12 +
 .../phasePair/phasePair/phasePair.C           |  19 +-
 .../phaseSystems/phaseSystem/phaseSystem.C    | 303 +++++----------
 .../phaseSystems/phaseSystem/phaseSystem.H    | 180 ++++++---
 .../phaseSystems/phaseSystem/phaseSystemI.H   |   7 +
 .../phaseSystem/phaseSystemSolve.C            |  90 +----
 .../general/include/createPIDTimeControls.H   |  60 +++
 .../general/include/readPIDTimeControls.H     |  91 +++++
 .../general/include/setInitialPIDDeltaT.H     |  53 +++
 .../cfdTools/general/include/setPIDDeltaT.H   |  82 +++++
 83 files changed, 6935 insertions(+), 337 deletions(-)
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/CourantNo.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/EEqns.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/files
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/options
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/cipsaMultiphaseEulerFoam.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFieldRefs.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFields.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/initZeroPtrList.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/HUEqn.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/UEqns.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/ddtUEqn.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dmdtUEqn.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragDiag.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSource.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSourcef.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSourcef.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/fvOptionsEqn.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/pEqn.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/stabilizeHUEqn.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmDiag.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSource.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSourcef.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/maxDeltaxyz/maxDeltaxyz.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/maxDeltaxyz/maxDeltaxyz.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/oneDimensionalDelta/oneDimensionalDelta.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/oneDimensionalDelta/oneDimensionalDelta.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/ResolvedPhaseSystem.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/ResolvedPhaseSystem.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/constantInterfaceCompression/constantInterfaceCompression.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/constantInterfaceCompression/constantInterfaceCompression.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModel.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModel.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModelNew.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/noInterfaceCompression/noInterfaceCompression.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/noInterfaceCompression/noInterfaceCompression.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/gradAlpha/gradAlpha.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/gradAlpha/gradAlpha.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethod.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethod.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethodNew.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/useBlendingMethod/useBlendingMethod.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/useBlendingMethod/useBlendingMethod.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/continuumSurfaceForce/continuumSurfaceForce.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/continuumSurfaceForce/continuumSurfaceForce.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/noSurfaceTensionForce/noSurfaceTensionForce.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/noSurfaceTensionForce/noSurfaceTensionForce.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModel.C
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModel.H
 create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModelNew.C
 create mode 100644 src/finiteVolume/cfdTools/general/include/createPIDTimeControls.H
 create mode 100644 src/finiteVolume/cfdTools/general/include/readPIDTimeControls.H
 create mode 100644 src/finiteVolume/cfdTools/general/include/setInitialPIDDeltaT.H
 create mode 100644 src/finiteVolume/cfdTools/general/include/setPIDDeltaT.H

diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/Allwclean b/applications/solvers/multiphase/multiphaseEulerFoam/Allwclean
index b8fe60678..858f85f76 100755
--- a/applications/solvers/multiphase/multiphaseEulerFoam/Allwclean
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/Allwclean
@@ -7,6 +7,7 @@ wclean libso interfacialCompositionModels
 wclean libso multiphaseCompressibleMomentumTransportModels
 wclean libso multiphaseThermophysicalTransportModels
 multiphaseEulerFoam/Allwclean
+wclean libso cipsaMultiphaseEulerFoam
 wclean libso functionObjects
 
 #------------------------------------------------------------------------------
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/Allwmake b/applications/solvers/multiphase/multiphaseEulerFoam/Allwmake
index db3b01550..9c4618de9 100755
--- a/applications/solvers/multiphase/multiphaseEulerFoam/Allwmake
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/Allwmake
@@ -10,6 +10,7 @@ wmake $targetType interfacialCompositionModels
 wmake $targetType multiphaseCompressibleMomentumTransportModels
 wmake $targetType multiphaseThermophysicalTransportModels
 multiphaseEulerFoam/Allwmake $targetType $*
+wmake $targetType cipsaMultiphaseEulerFoam
 wmake $targetType functionObjects
 
 #------------------------------------------------------------------------------
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/CourantNo.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/CourantNo.H
new file mode 100644
index 000000000..d1d28ab5c
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/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/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/EEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/EEqns.H
new file mode 100644
index 000000000..a6e4c5256
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/EEqns.H
@@ -0,0 +1,61 @@
+for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
+{
+    fluid.correctEnergyTransport();
+
+    autoPtr<phaseSystem::heatTransferTable>
+        heatTransferPtr(fluid.heatTransfer());
+
+    phaseSystem::heatTransferTable& heatTransfer = heatTransferPtr();
+
+    forAll(fluid.anisothermalPhases(), anisothermalPhasei)
+    {
+        phaseModel& phase = fluid.anisothermalPhases()[anisothermalPhasei];
+
+        const volScalarField& alpha = phase;
+        tmp<volScalarField> tRhoAlpha = phase.rho();
+        const volScalarField& rhoAlpha = tRhoAlpha();
+        tmp<volVectorField> tUAlpha = phase.U();
+        const volVectorField& UAlpha = tUAlpha();
+
+        fvScalarMatrix EEqn
+        (
+            phase.heEqn()
+         ==
+           *heatTransfer[phase.name()]
+          + alpha*rhoAlpha*(UAlpha & g)
+          + fvOptions(alpha, rhoAlpha, phase.thermoRef().he())
+        );
+
+        EEqn +=
+            fvm::ddt
+            (
+                max(phase.residualAlpha() - alpha, Zero)*rhoAlpha,
+                phase.thermoRef().he()
+            )
+          - fvc::ddt
+            (
+                max(phase.residualAlpha() - alpha, Zero)*rhoAlpha,
+                phase.thermoRef().he()
+            );
+
+        EEqn.relax();
+        fvOptions.constrain(EEqn);
+        EEqn.solve();
+        fvOptions.correct(phase.thermoRef().he());
+    }
+
+    fluid.correctThermo();
+    fluid.correctContinuityError();
+}
+
+
+forAll(phases, phasei)
+{
+    phaseModel& phase = phases[phasei];
+
+    Info<< phase.name() << " min/max T "
+        << min(phase.thermo().T()).value()
+        << " - "
+        << max(phase.thermo().T()).value()
+        << endl;
+}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/files b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/files
new file mode 100644
index 000000000..3299c2137
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/files
@@ -0,0 +1,3 @@
+cipsaMultiphaseEulerFoam.C
+
+EXE = $(FOAM_APPBIN)/cipsaMultiphaseEulerFoam
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/options b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/options
new file mode 100644
index 000000000..915392055
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/options
@@ -0,0 +1,25 @@
+EXE_INC = \
+    -I../multiphaseEulerFoam/multiphaseSystems \
+    -I../phaseSystems/lnInclude \
+    -I../interfacialModels/lnInclude \
+    -I../interfacialCompositionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
+    -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -lphaseSystem \
+    -lmultiphaseSystems \
+    -leulerianInterfacialModels \
+    -leulerianInterfacialCompositionModels \
+    -lmultiphaseMomentumTransportModels \
+    -lmultiphaseThermophysicalTransportModels \
+    -lthermophysicalTransportModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/cipsaMultiphaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/cipsaMultiphaseEulerFoam.C
new file mode 100644
index 000000000..2a4bb8d72
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/cipsaMultiphaseEulerFoam.C
@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | 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 "phaseSystem.H"
+#include "phaseCompressibleMomentumTransportModel.H"
+#include "pimpleControl.H"
+#include "initZeroPtrList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createPIDTimeControls.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "CourantNo.H"
+    #include "setInitialPIDDeltaT.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 "readPIDTimeControls.H"
+
+        int nEnergyCorrectors
+        (
+            pimple.dict().lookupOrDefault<int>("nEnergyCorrectors", 1)
+        );
+
+        #include "CourantNo.H"
+        #include "setPIDDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity corrector loop (PIMPLE)
+        while (pimple.loop())
+        {
+            fluid.solve(rAUs, rHAUfs);
+            fluid.correct();
+            // Update continuityError in MovingPhaseModel for heEqn in
+            // AnisothermalPhaseModel
+            fluid.correctContinuityError();
+
+            #include "pU/UEqns.H"
+            #include "EEqns.H"
+
+            // --- Pressure corrector loop (PISO)
+            while (pimple.correct())
+            {
+                #include "pU/pEqn.H"
+            }
+
+            // Update material derivatives
+            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/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFieldRefs.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFieldRefs.H
new file mode 100644
index 000000000..1803f8d75
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFieldRefs.H
@@ -0,0 +1,49 @@
+surfaceScalarField& phi = fluid.phi();
+
+const IOMRFZoneList& MRF = fluid.MRF();
+const fv::options& fvOptions = fluid.fvOptions();
+
+List<scalar> relaxU(phases.size());
+PtrList<surfaceScalarField> alphafs(phases.size());
+
+// Equations
+PtrList<fvVectorMatrix> HUEqns(phases.size());
+PtrList<fvVectorMatrix> ddtUEqns(phases.size());
+PtrList<fvVectorMatrix> dmdtUEqns(phases.size());
+PtrList<fvVectorMatrix> fvOptionsEqns(phases.size());
+PtrList<surfaceScalarField> ddtPhiCoeffs(phases.size());
+PtrList<surfaceScalarField> dmdtPhiCoeffs(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());
+
+// 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> phigs(phases.size());
+PtrList<surfaceScalarField> phiHbyAs(phases.size());
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFields.H
new file mode 100644
index 000000000..d24413a6f
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFields.H
@@ -0,0 +1,69 @@
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+
+Info<< "Creating phaseSystem\n" << endl;
+
+autoPtr<phaseSystem> fluidPtr
+(
+    phaseSystem::New(mesh)
+);
+phaseSystem& fluid = fluidPtr();
+phaseSystem::phaseModelList& phases = fluid.phases();
+phaseSystem::phaseModelPartialList& movingPhases = fluid.movingPhases();
+
+// 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());
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/initZeroPtrList.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/initZeroPtrList.H
new file mode 100644
index 000000000..c7816d470
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/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/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/HUEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/HUEqn.H
new file mode 100644
index 000000000..d7ad17416
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/HUEqn.H
@@ -0,0 +1,84 @@
+forAll(movingPhases, movingPhasei)
+{
+    phaseModel& phase = movingPhases[movingPhasei];
+    volVectorField& UAlpha = phase.URef();
+    const volScalarField& alpha = phase;
+    const volScalarField& rhoAlpha = phase.thermoRef().rho();
+
+    HUEqns.set
+    (
+        phase.index(),
+        (
+            fvm::div(phase.alphaRhoPhi(), UAlpha)
+          + phase.divDevTau(UAlpha)
+          + fvm::SuSp(- fvc::div(phase.alphaRhoPhi()), UAlpha)
+          + MRF.DDt(alpha*rhoAlpha, UAlpha)
+        )
+    );
+}
+
+// Create U & grad(U) fields
+PtrList<fvVectorMatrix> UgradUs(phases.size());
+forAll(movingPhases, movingPhasei)
+{
+    const phaseModel& phase = movingPhases[movingPhasei];
+    const volVectorField& UAlpha = phase.U();
+    const surfaceScalarField& phiAlpha = phase.phi();
+
+    UgradUs.set
+    (
+        phase.index(),
+        new fvVectorMatrix
+        (
+            fvm::div(phiAlpha, UAlpha)
+          + fvm::SuSp(- fvc::div(phiAlpha), UAlpha)
+        )
+    );
+}
+
+// Add convective part of the virtual mass force
+forAllConstIter(phaseSystem::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(movingPhases, movingPhasei)
+{
+    const phaseModel& phase = movingPhases[movingPhasei];
+
+    word name = HUEqns[phase.index()].psi().select(pimple.finalPimpleIter());
+
+    relaxU[phase.index()] = 1;
+
+    if (HUEqns[phase.index()].psi().mesh().relaxField(name))
+    {
+        relaxU[phase.index()] =
+            HUEqns[phase.index()].psi().mesh().fieldRelaxationFactor(name);
+
+        if
+        (
+            pimple.finalPimpleIter()
+         && HUEqns[phase.index()].psi().mesh().fieldRelaxationFactor(name) != 1
+        )
+        {
+            Warning << "Field relaxation factor of " << relaxU[phase.index()]
+                    << " for solution variable '" << name
+                    << "' is not recommended." << endl << endl;
+        }
+    }
+}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/UEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/UEqns.H
new file mode 100644
index 000000000..90a7d8cf5
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/UEqns.H
@@ -0,0 +1,234 @@
+// Update mixture properties
+rho = fluid.rho();
+U = fluid.U();
+
+// Store velocity and flux of previous iteration
+forAll(movingPhases, movingPhasei)
+{
+    movingPhases[movingPhasei].U()().storePrevIter();
+    movingPhases[movingPhasei].phi()().storePrevIter();
+}
+
+// Update drag and virtual mass coefficient table and interfacial forces tables
+fluid.momentumTransfer();
+
+// Time derivative equation
+#include "pU/ddtUEqn.H"
+
+// Momentum transfer by mass transfer equation
+#include "pU/dmdtUEqn.H"
+
+// Convection, diffusion equation
+#include "pU/HUEqn.H"
+#include "pU/stabilizeHUEqn.H"
+
+// fvOptions for momentum equation
+#include "pU/fvOptionsEqn.H"
+
+// Face volume fractions
+forAll(phases, phasei)
+{
+    const phaseModel& phase = phases[phasei];
+    const volScalarField& alpha = phase;
+
+    alphafs.set(phasei, fvc::interpolate(alpha));
+}
+
+// rHUs, rHUfs, rAs, rAfs
+forAll(movingPhases, movingPhasei)
+{
+    const phaseModel& phase = movingPhases[movingPhasei];
+
+    rHUs.set(phase.index(), 1.0/HUEqns[phase.index()].A());
+    rHUfs.set(phase.index(), fvc::interpolate(rHUs[phase.index()]));
+
+    // Diagonal coefficients for partial elimination
+    rAs.set
+    (
+        phase.index(),
+        1.0
+       /(
+            HUEqns[phase.index()].A()
+          + ddtUEqns[phase.index()].A()
+          + fvOptionsEqns[phase.index()].A()
+        )
+    );
+
+    // Face diagonal coefficients for partial elimination
+    rAfs.set
+    (
+        phase.index(),
+        1.0
+       /(
+            fvc::interpolate(HUEqns[phase.index()].A())
+          + fvc::interpolate(ddtUEqns[phase.index()].A())
+          + fvc::interpolate(fvOptionsEqns[phase.index()].A())
+        )
+    );
+}
+fluid.fillFields("rHUs", dimless, rHUs);
+fluid.fillFields("rHUfs", dimless, rHUfs);
+fluid.fillFields("rAs", dimTime/dimDensity, rAs);
+fluid.fillFields("rAfs", dimTime/dimDensity, rAfs);
+
+// 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())
+{
+    // Correct p_rgh for consistency with p
+    p_rgh = p - rho*gh;
+
+    // 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/forcesSourcef.H"
+
+    // Mixture density gradient times gravitational acceleration
+    const surfaceScalarField ghSnGradRho
+    (
+        "ghSnGradRho",
+        ghf*fvc::snGrad(rho)*mesh.magSf()
+    );
+
+    // rAUfs, alpharAUfs, phigs WITHOUT field relaxation of UAlpha
+    forAll(movingPhases, movingPhasei)
+    {
+        phaseModel& phase = movingPhases[movingPhasei];
+        const volScalarField& alpha = phase;
+        const volScalarField& rhoAlpha = phase.rho();
+        volVectorField& UAlpha = phase.URef();
+
+        rAUfs.set
+        (
+            phase.index(),
+            1.0
+           /(
+                1.0
+              + rHUfs[phase.index()]
+               *(
+                    fvc::interpolate(ddtUEqns[phase.index()].A())
+                  + fvc::interpolate(dmdtUEqns[phase.index()].A())
+                  + fvc::interpolate(fvOptionsEqns[phase.index()].A())
+                  + dragDiagf[phase.index()]
+                  + vmDiagf[phase.index()]
+                )
+            )
+        );
+
+        // Phase diagonal coefficients
+        alpharAUfs.set
+        (
+            phase.index(),
+            rHUfs[phase.index()]*rAUfs[phase.index()]
+           *(alphafs[phase.index()] + dragPCoefff[phase.index()])
+        );
+
+        // Combined buoyancy and force fluxes
+        phigs.set
+        (
+            phase.index(),
+            alpharAUfs[phase.index()]
+           *(
+              - ghSnGradRho
+              + fluid.surfaceTension(phase)*mesh.magSf()
+            )
+        );
+
+        // Momentum predictor equation
+        fvVectorMatrix UEqns
+        (
+            ddtUEqns[phase.index()]
+          + dmdtUEqns[phase.index()]
+          + HUEqns[phase.index()]
+          + fvm::Sp(dragDiag[phase.index()], UAlpha)
+          + fvm::Sp(vmDiag[phase.index()], UAlpha)
+          + fvOptionsEqns[phase.index()]
+            ==
+            dragSource[phase.index()]
+          + alpha*(rhoAlpha - rho)*g
+          + (
+                ddtUEqns[phase.index()].A()
+              + dmdtUEqns[phase.index()].A()
+              + HUEqns[phase.index()].A()
+              + dragDiag[phase.index()]
+              + fvOptionsEqns[phase.index()].A()
+            )
+           *fvc::reconstruct
+            (
+                phigs[phase.index()]
+              - alpharAUfs[phase.index()]*fvc::snGrad(p_rgh)*mesh.magSf()
+              + rHUfs[phase.index()]*rAUfs[phase.index()]
+               *forcesSourcef[phase.index()]
+            )
+          + vmSource[phase.index()]
+        );
+
+        UEqns.relax();
+        fvOptions.constrain(UEqns);
+        UEqns.solve();
+        fvOptions.correct(UAlpha);
+    }
+}
+
+// Construct pressure "diffusivity"
+surfaceScalarField rAUf
+(
+    IOobject
+    (
+        "rAUf",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar(dimensionSet(-1, 3, 1, 0, 0), Zero)
+);
+
+// rAUfs, rHAUfs, alpharAUfs, rAUf WITH field relaxation of UAlpha
+initZeroPtrList(rAUfs, mesh, dimless, "rAUfs");
+initZeroPtrList(rHAUfs, mesh, dimTime/dimDensity, "rHAUfs");
+initZeroPtrList(alpharAUfs, mesh, dimTime/dimDensity, "alpharAUfs");
+forAll(movingPhases, movingPhasei)
+{
+    const phaseModel& phase = movingPhases[movingPhasei];
+
+    rAUfs.set
+    (
+        phase.index(),
+        1.0
+       /(
+            1.0/max(relaxU[phase.index()], small)
+          + rHUfs[phase.index()]
+           *(
+                fvc::interpolate(ddtUEqns[phase.index()].A())
+              + fvc::interpolate(dmdtUEqns[phase.index()].A())
+              + fvc::interpolate(fvOptionsEqns[phase.index()].A())
+              + dragDiagf[phase.index()]
+              + vmDiagf[phase.index()]
+            )
+        )
+    );
+
+    rHAUfs.set(phase.index(), rHUfs[phase.index()]*rAUfs[phase.index()]);
+
+    // Phase diagonal coefficients
+    alpharAUfs.set
+    (
+        phase.index(),
+        rHUfs[phase.index()]*rAUfs[phase.index()]
+       *(alphafs[phase.index()] + dragPCoefff[phase.index()])
+    );
+
+    // Pressure "diffusivity"
+    rAUf += alphafs[phase.index()]*alpharAUfs[phase.index()];
+}
+
+rAUf = mag(rAUf);
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/ddtUEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/ddtUEqn.H
new file mode 100644
index 000000000..fb81d8e0d
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/ddtUEqn.H
@@ -0,0 +1,47 @@
+forAll(movingPhases, movingPhasei)
+{
+    const phaseModel& phase = movingPhases[movingPhasei];
+    const volScalarField& alpha = phase;
+    const volScalarField& rhoAlpha = phase.rho();
+    const volVectorField& UAlpha = phase.U();
+
+    ddtUEqns.set
+    (
+        phase.index(),
+        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
+        (
+            phase.index(),
+            (
+                byDt
+                (
+                    fvc::interpolate(alpha.oldTime())
+                   *fvc::interpolate(rhoAlpha.oldTime())
+                   *MRF.absolute(phase.phi()().oldTime())
+                )
+            )
+        );
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Only Euler is an allowed time scheme "
+            << "for this solver." << exit(FatalError);
+    }
+}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dmdtUEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dmdtUEqn.H
new file mode 100644
index 000000000..b5b1d6bac
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dmdtUEqn.H
@@ -0,0 +1,74 @@
+forAll(movingPhases, movingPhasei)
+{
+    const phaseModel& phase = movingPhases[movingPhasei];
+    const volVectorField& UAlpha = phase.U();
+    const PtrList<volScalarField> dmdts = fluid.dmdts();
+
+    const dimensionedScalar ZERO(dimensionSet(1, -3, -1, 0, 0), Zero);
+
+    dmdtUEqns.set
+    (
+        phase.index(),
+        (
+            fvm::Sp(ZERO, UAlpha)
+          + fvm::Su(Zero, UAlpha)
+        )
+    );
+
+    volScalarField source
+    (
+        volScalarField::New
+        (
+            IOobject::groupName("source", phase.name()),
+            mesh,
+            dimensionedScalar(dimDensity/dimTime, 0)
+        )
+    );
+
+    if (fvOptions.appliesToField(rho.name()))
+    {
+        source += fvOptions(phase, rho)&rho;
+    }
+
+    if (dmdts.set(phase.index()))
+    {
+        source += dmdts[phase.index()];
+    }
+
+    dmdtUEqns[phase.index()] += fvm::SuSp(source, UAlpha);
+}
+
+initZeroPtrList
+(
+    dmdtPhiCoeffs,
+    mesh,
+    dimensionSet(1, 0, -2, 0, 0),
+    "dmdtPhiCoeff"
+);
+
+// Add mass sources
+// Inspired by MomentumTransferPhaseSystem::addDmdtUfs() but with fvm::SuSp()
+// instead of fvm::Sp()
+forAllConstIter(phaseSystem::dmdtfTable, fluid.dmdtfsPhaseTransfer(), dmdtfIter)
+{
+    const phasePairKey& key = dmdtfIter.key();
+    const phasePair& pair(fluid.phasePairs()[key]);
+
+    const volScalarField dmdtf(Pair<word>::compare(pair, key)**dmdtfIter());
+
+    phaseModel& phase1 = phases[pair.phase1().name()];
+    phaseModel& phase2 = phases[pair.phase2().name()];
+
+    const volScalarField dmdtf21(posPart(dmdtf));
+    const volScalarField dmdtf12(negPart(dmdtf));
+
+    dmdtUEqns[pair.phase1().index()] -=
+        dmdtf21*phase2.U() + fvm::SuSp(dmdtf12, phase1.URef());
+    dmdtPhiCoeffs[pair.phase1().index()] +=
+        fvc::interpolate(dmdtf21)*MRF.absolute(phase2.phi());
+
+    dmdtUEqns[pair.phase2().index()] +=
+        dmdtf12*phase1.U() - fvm::SuSp(-dmdtf21, phase2.URef());
+    dmdtPhiCoeffs[pair.phase2().index()] -=
+        fvc::interpolate(dmdtf12)*MRF.absolute(phase1.phi());
+}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragDiag.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragDiag.H
new file mode 100644
index 000000000..ce183e3bf
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragDiag.H
@@ -0,0 +1,71 @@
+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(phaseSystem::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 (phase.stationary()) continue;
+
+        if (partialElimination)
+        {
+            dragDiag[phase.index()] +=
+                K
+               *(
+                    1
+                  + rAs[otherPhase.index()]
+                   *(
+                        HUEqns[phase.index()].A()
+                      + ddtUEqns[phase.index()].A()
+                      + fvOptionsEqns[phase.index()].A()
+                    )
+                );
+        }
+        else
+        {
+            dragDiag[phase.index()] += K;
+        }
+    }
+}
+
+forAllConstIter(phaseSystem::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 (phase.stationary()) continue;
+
+        if (partialElimination)
+        {
+            dragDiagf[phase.index()] +=
+                Kf
+               *(
+                    1
+                  + rAfs[otherPhase.index()]
+                   *(
+                        fvc::interpolate(HUEqns[phase.index()].A())
+                      + fvc::interpolate(ddtUEqns[phase.index()].A())
+                      + fvc::interpolate(fvOptionsEqns[phase.index()].A())
+                    )
+                );
+
+            dragPCoefff[phase.index()] += Kf*rAfs[otherPhase.index()];
+        }
+        else
+        {
+            dragDiagf[phase.index()] += Kf;
+        }
+    }
+}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSource.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSource.H
new file mode 100644
index 000000000..7eace2f3d
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSource.H
@@ -0,0 +1,75 @@
+initZeroPtrList(dragSource, mesh, dimensionSet(1, -2, -2, 0, 0), "dragSource");
+
+{
+    // 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(movingPhases, movingPhasei)
+        {
+            const phaseModel& phase = movingPhases[movingPhasei];
+
+            HSum +=
+                HUEqns[phase.index()].H()
+              + ddtUEqns[phase.index()].H()
+              + fvOptionsEqns[phase.index()].H();
+        }
+    }
+
+    forAllConstIter(phaseSystem::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 (phase.stationary() || otherPhase.stationary()) continue;
+
+            if (partialElimination)
+            {
+                dragSource[phase.index()] += K*rAs[otherPhase.index()]*HSum;
+
+                forAll(movingPhases, movingPhasei)
+                {
+                    const phaseModel& remainingPhase =
+                        movingPhases[movingPhasei];
+
+                    if
+                    (
+                        (remainingPhase.index() == phase.index())
+                     || (remainingPhase.index() == otherPhase.index())
+                    ) continue;
+
+                    dragSource[phase.index()] -=
+                    (
+                        K*rAs[otherPhase.index()]
+                       *(
+                            HUEqns[remainingPhase.index()].A()
+                          + ddtUEqns[remainingPhase.index()].A()
+                          + fvOptionsEqns[remainingPhase.index()].A()
+                        )
+                       *phases[remainingPhase.index()].U()
+                    );
+                }
+            }
+            else
+            {
+                dragSource[phase.index()] += K*phases[otherPhase.index()].U();
+            }
+        }
+    }
+}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSourcef.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSourcef.H
new file mode 100644
index 000000000..5f24aa1e6
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSourcef.H
@@ -0,0 +1,84 @@
+initZeroPtrList(dragSourcef, mesh, dimensionSet(0, 3, -1, 0, 0), "dragSourcef");
+
+{
+    // Summed RHS of face-based momentum balance equation
+    surfaceScalarField HSumf
+    (
+        IOobject
+        (
+            "HSumf",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar(dimensionSet(1, 0, -2, 0, 0), Zero)
+    );
+
+    if (partialElimination)
+    {
+        forAll(movingPhases, movingPhasei)
+        {
+            const phaseModel& phase = movingPhases[movingPhasei];
+
+            HSumf +=
+                fvc::flux(HUEqns[phase.index()].H())
+              + ddtPhiCoeffs[phase.index()]
+              + fvc::flux(fvOptionsEqns[phase.index()].H());
+        }
+    }
+
+    forAllConstIter(phaseSystem::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 (phase.stationary() || otherPhase.stationary()) continue;
+
+            if (partialElimination)
+            {
+                dragSourcef[phase.index()]
+                    += Kf*rAfs[otherPhase.index()]*rHUfs[phase.index()]*HSumf;
+
+                forAll(movingPhases, movingPhasei)
+                {
+                    const phaseModel& remainingPhase =
+                        movingPhases[movingPhasei];
+
+                    if
+                    (
+                        (remainingPhase.index() == phase.index())
+                     || (remainingPhase.index() == otherPhase.index())
+                    ) continue;
+
+                    dragSourcef[phase.index()] -=
+                    (
+                        Kf*rAfs[otherPhase.index()]*rHUfs[phase.index()]
+                       *(
+                            fvc::interpolate(HUEqns[remainingPhase.index()].A())
+                          + fvc::interpolate
+                            (
+                                ddtUEqns[remainingPhase.index()].A()
+                            )
+                          + fvc::interpolate
+                            (
+                                fvOptionsEqns[remainingPhase.index()].A()
+                            )
+                        )
+                       *MRF.absolute(phases[remainingPhase.index()].phi())
+                    );
+                }
+            }
+            else
+            {
+                dragSourcef[phase.index()]
+                    += Kf*rHUfs[phase.index()]
+                      *MRF.absolute(phases[otherPhase.index()].phi());
+            }
+        }
+    }
+}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSourcef.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSourcef.H
new file mode 100644
index 000000000..2c674a319
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/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(phaseSystem::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(phaseSystem::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(phaseSystem::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/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/fvOptionsEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/fvOptionsEqn.H
new file mode 100644
index 000000000..bab37091e
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/fvOptionsEqn.H
@@ -0,0 +1,18 @@
+forAll(movingPhases, movingPhasei)
+{
+    phaseModel& phase = movingPhases[movingPhasei];
+    const volScalarField& alpha = phase;
+    volVectorField& UAlpha = phase.URef();
+    volScalarField& rhoAlpha = phase.thermoRef().rho();
+
+    fvOptionsEqns.set
+    (
+        phase.index(),
+        fvm::Sp(dimensionedScalar(dimDensity/dimTime, Zero), UAlpha)
+    );
+
+    if(fvOptions.appliesToField(UAlpha.name()))
+    {
+        fvOptionsEqns[phase.index()] -= fvOptions(alpha, rhoAlpha, UAlpha);
+    }
+}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/pEqn.H
new file mode 100644
index 000000000..7f2dc1e87
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/pEqn.H
@@ -0,0 +1,347 @@
+rho = fluid.rho();
+
+// Correct p_rgh for consistency with p and the updated densities due to
+// updated temperature
+p_rgh = p - rho*gh;
+
+// Correct fixed-flux BCs to be consistent with the velocity BCs
+forAll(fluid.movingPhases(), movingPhasei)
+{
+    phaseModel& phase = fluid.movingPhases()[movingPhasei];
+    MRF.correctBoundaryFlux(phase.U(), phase.phiRef());
+}
+
+{
+    // Mixture density gradient times gravitational acceleration
+    const surfaceScalarField ghSnGradRho
+    (
+        "ghSnGradRho",
+        ghf*fvc::snGrad(rho)*mesh.magSf()
+    );
+
+    // phigs WITH field relaxation of UAlpha
+    forAll(movingPhases, movingPhasei)
+    {
+        const phaseModel& phase = movingPhases[movingPhasei];
+
+        // Combined buoyancy and force fluxes
+        phigs.set
+        (
+            phase.index(),
+            alpharAUfs[phase.index()]
+           *(
+              - ghSnGradRho
+              + fluid.surfaceTension(phase)*mesh.magSf()
+            )
+        );
+    }
+}
+
+// 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/forcesSourcef.H"
+
+// Construct Total predicted flux
+surfaceScalarField phiHbyA
+(
+    IOobject
+    (
+        "phiHbyA",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(dimArea*dimVelocity, Zero)
+);
+
+forAll(movingPhases, movingPhasei)
+{
+    const phaseModel& phase = fluid.movingPhases()[movingPhasei];
+    const volScalarField& rhoAlpha = phase.rho();
+    const volVectorField& UAlpha = phase.U();
+
+    // Pseudo convected velocities
+    HbyAs.set
+    (
+        phase.index(),
+        constrainHbyA
+        (
+            rHUs[phase.index()]*HUEqns[phase.index()].H(),
+            UAlpha,
+            p_rgh
+        )
+    );
+
+    // Predicted fluxes for each phase
+    phiHbyAs.set
+    (
+        phase.index(),
+        constrainPhiHbyA
+        (
+            rAUfs[phase.index()]
+           *(
+                fvc::flux(HbyAs[phase.index()])
+              + rHUfs[phase.index()]
+               *(
+                    ddtPhiCoeffs[phase.index()]
+                  + dmdtPhiCoeffs[phase.index()]
+                  + alphafs[phase.index()]*(g & mesh.Sf())
+                   *(
+                        fvc::interpolate(rhoAlpha) - fvc::interpolate(rho)
+                    )
+                  + forcesSourcef[phase.index()]
+                  + fvc::flux(fvOptionsEqns[phase.index()].H())
+                )
+              + (1 - relaxU[phase.index()])/max(relaxU[phase.index()], small)
+               *MRF.absolute(phases[phase.index()].phi()().prevIter())
+              + dragSourcef[phase.index()]
+              + vmSourcef[phase.index()]
+            ),
+            UAlpha,
+            p_rgh
+        )
+      + phigs[phase.index()]
+    );
+
+    phiHbyA += alphafs[phase.index()]*phiHbyAs[phase.index()];
+}
+
+MRF.makeRelative(phiHbyA);
+
+// Update the fixedFluxPressure BCs to ensure flux consistency
+{
+    surfaceScalarField::Boundary phib(phi.boundaryField());
+    phib = 0;
+    forAll(phases, phasei)
+    {
+        const phaseModel& phase = phases[phasei];
+        phib += alphafs[phasei].boundaryField()*phase.phi()().boundaryField();
+    }
+
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryFieldRef(),
+        (
+            phiHbyA.boundaryField() - phib
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+}
+
+// Compressible pressure equations
+PtrList<fvScalarMatrix> pEqnComps(phases.size());
+PtrList<volScalarField> dmdts(fluid.dmdts());
+forAll(phases, phasei)
+{
+    phaseModel& phase = phases[phasei];
+    const volScalarField& alpha = phase;
+    volScalarField& rhoAlpha = phase.thermoRef().rho();
+
+    pEqnComps.set(phasei, new fvScalarMatrix(p_rgh, dimVolume/dimTime));
+    fvScalarMatrix& pEqnComp = pEqnComps[phasei];
+
+    // Density variation
+    if (!phase.isochoric() || !phase.pure())
+    {
+        pEqnComp +=
+            (
+                fvc::ddt(alpha, rhoAlpha) + fvc::div(phase.alphaRhoPhi())
+              - fvc::Sp(fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), rhoAlpha)
+            )/rhoAlpha;
+    }
+
+    // Compressibility
+    if (!phase.incompressible())
+    {
+        if (pimple.transonic())
+        {
+            const surfaceScalarField phid
+            (
+                IOobject::groupName("phid", phase.name()),
+                fvc::interpolate(phase.thermo().psi())*phase.phi()
+            );
+
+            pEqnComp +=
+                correction
+                (
+                    (alpha/rhoAlpha)*
+                    (
+                        phase.thermo().psi()*fvm::ddt(p_rgh)
+                      + fvm::div(phid, p_rgh)
+                      - fvm::Sp(fvc::div(phid), p_rgh)
+                    )
+                );
+
+            deleteDemandDrivenData
+            (
+                pEqnComps[phasei].faceFluxCorrectionPtr()
+            );
+
+            pEqnComps[phasei].relax();
+        }
+        else
+        {
+            pEqnComp +=
+                (alpha*phase.thermo().psi()/rhoAlpha)
+               *correction(fvm::ddt(p_rgh));
+        }
+    }
+
+    // Option sources
+    if (fvOptions.appliesToField(rhoAlpha.name()))
+    {
+        pEqnComp -= (fvOptions(alpha, rhoAlpha) & rhoAlpha)/rhoAlpha;
+    }
+
+    // Mass transfer
+    if (dmdts.set(phasei))
+    {
+        pEqnComp -= dmdts[phasei]/rhoAlpha;
+    }
+}
+
+// Cache p prior to solve for density update
+volScalarField p_rgh_0(p_rgh);
+
+// Iterate over the pressure equation to correct for non-orthogonality
+while (pimple.correctNonOrthogonal())
+{
+    // Construct the transport part of the pressure equation
+    fvScalarMatrix pEqnIncomp
+    (
+        fvc::div(phiHbyA)
+      - fvm::laplacian(rAUf, p_rgh)
+    );
+
+    // Solve
+    {
+        fvScalarMatrix pEqn(pEqnIncomp);
+
+        forAll(phases, phasei)
+        {
+            pEqn += pEqnComps[phasei];
+        }
+
+        if (fluid.incompressible())
+        {
+            pEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        }
+
+        pEqn.solve();
+    }
+
+    // Correct fluxes and velocities on last non-orthogonal iteration
+    if (pimple.finalNonOrthogonalIter())
+    {
+        phi = phiHbyA + pEqnIncomp.flux();
+
+        surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
+
+        forAll(movingPhases, movingPhasei)
+        {
+            phaseModel& phase = movingPhases[movingPhasei];
+
+            phase.phiRef() =
+                phiHbyAs[phase.index()]
+              + alpharAUfs[phase.index()]*mSfGradp;
+
+            // Set the phase dilatation rate
+            phase.divU(-pEqnComps[phase.index()] & p_rgh);
+        }
+
+        // Optionally relax pressure for velocity correction
+        p_rgh.relax();
+
+        mSfGradp = pEqnIncomp.flux()/rAUf;
+
+        // Momentum corrector
+        forAll(movingPhases, movingPhasei)
+        {
+            phaseModel& phase = movingPhases[movingPhasei];
+            const volScalarField& alpha = phase;
+            const volScalarField& rhoAlpha = phase.rho();
+            volVectorField& UAlpha = phase.URef();
+
+            if (faceMomentum)
+            {
+                UAlpha = fvc::reconstruct(phase.phi());
+            }
+            else
+            {
+                UAlpha =
+                    1.0
+                   /(
+                        ddtUEqns[phase.index()].A()
+                      + dmdtUEqns[phase.index()].A()
+                      + HUEqns[phase.index()].A()
+                       /max(relaxU[phase.index()], small)
+                      + dragDiag[phase.index()]
+                      + vmDiag[phase.index()]
+                      + fvOptionsEqns[phase.index()].A()
+                    )
+                   *(
+                        HUEqns[phase.index()].A()*HbyAs[phase.index()]
+                      + ddtUEqns[phase.index()].H()
+                      + dmdtUEqns[phase.index()].H()
+                      + dragSource[phase.index()]
+                      + alpha*(rhoAlpha - rho)*g
+                      + (1 - relaxU[phase.index()])
+                       /max(relaxU[phase.index()], small)
+                       *HUEqns[phase.index()].A()*UAlpha.prevIter()
+                      + vmSource[phase.index()]
+                      + fvOptionsEqns[phase.index()].H()
+                    )
+                  + fvc::reconstruct
+                    (
+                        phigs[phase.index()]
+                      + alpharAUfs[phase.index()]*mSfGradp
+                      + rHUfs[phase.index()]*rAUfs[phase.index()]
+                       *forcesSourcef[phase.index()]
+                    );
+            }
+
+            UAlpha.correctBoundaryConditions();
+            fvOptions.correct(UAlpha);
+
+            MRF.makeRelative(phase.phiRef());
+        }
+    }
+}
+
+// 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)
+    );
+}
+
+// Limit p_rgh
+p_rgh = p - rho*gh;
+
+// Update densities from change in p_rgh
+forAll(phases, phasei)
+{
+    phaseModel& phase = phases[phasei];
+    phase.thermoRef().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
+}
+
+// Correct p_rgh for consistency with p and the updated densities
+rho = fluid.rho();
+p_rgh = p - rho*gh;
+p_rgh.correctBoundaryConditions();
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/stabilizeHUEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/stabilizeHUEqn.H
new file mode 100644
index 000000000..22fb7747b
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/stabilizeHUEqn.H
@@ -0,0 +1,43 @@
+forAll(movingPhases, movingPhasei)
+{
+    phaseModel& phase = movingPhases[movingPhasei];
+    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[phase.index()] +=
+        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[phase.index()] +=
+        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/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmDiag.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmDiag.H
new file mode 100644
index 000000000..1a6f80659
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/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(phaseSystem::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/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSource.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSource.H
new file mode 100644
index 000000000..ddaabab31
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSource.H
@@ -0,0 +1,22 @@
+initZeroPtrList(vmSource, mesh, dimensionSet(1, -2, -2, 0, 0), "vmSource");
+
+forAllConstIter(phaseSystem::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/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSourcef.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSourcef.H
new file mode 100644
index 000000000..eed96da3b
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSourcef.H
@@ -0,0 +1,23 @@
+initZeroPtrList(vmSourcef, mesh, dimensionSet(0, 3, -1, 0, 0), "vmSourcef");
+
+forAllConstIter(phaseSystem::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
+            (
+                MRF.absolute(phase.phi()().oldTime())
+              + MRF.absolute(otherPhase.phi())
+              - MRF.absolute(otherPhase.phi()().oldTime())
+            );
+    }
+}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/Make/files b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/Make/files
index 24000605a..8850485af 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/Make/files
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/Make/files
@@ -22,4 +22,9 @@ saturationModels/polynomial/polynomial.C
 saturationModels/function1/function1.C
 saturationModels/constantSaturationConditions/constantSaturationConditions.C
 
+SGSdeltas/SGSdelta/SGSdelta.C
+SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.C
+SGSdeltas/maxDeltaxyz/maxDeltaxyz.C
+SGSdeltas/oneDimensionalDelta/oneDimensionalDelta.C
+
 LIB = $(FOAM_LIBBIN)/libeulerianInterfacialCompositionModels
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.C
new file mode 100644
index 000000000..fedaf01a2
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.C
@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | (C) 2019 Helmholtz-Zentrum Dresden-Rossendorf
+  ==  == ====== ====   ==  ==  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SGSdelta.H"
+#include "calculatedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(SGSdelta, 0);
+    defineRunTimeSelectionTable(SGSdelta, dictionary);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::SGSdelta::SGSdelta
+(
+    const word& name,
+    const fvMesh& mesh
+)
+:
+    mesh_(mesh),
+    delta_
+    (
+        IOobject
+        (
+            name,
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar(name, dimLength, small),
+        calculatedFvPatchScalarField::typeName
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::SGSdelta> Foam::SGSdelta::New
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+{
+    const word deltaType(dict.lookup("delta"));
+
+    Info<< "Selecting SGS delta type " << deltaType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(deltaType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown SGSdelta type "
+            << deltaType << nl << nl
+            << "Valid SGSdelta types are :" << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<SGSdelta>(cstrIter()(name, mesh, dict));
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.H
new file mode 100644
index 000000000..d8faefa03
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.H
@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | (C) 2019 Helmholtz-Zentrum Dresden-Rossendorf
+  ==  == ====== ====   ==  ==  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::SGSdelta
+
+Description
+    Abstract base class for SGS deltas
+
+See also
+    Foam::LESdelta
+
+SourceFiles
+    SGSdelta.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SGSdelta_H
+#define SGSdelta_H
+
+#include "volFields.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class SGSdelta Declaration
+\*---------------------------------------------------------------------------*/
+
+class SGSdelta
+{
+protected:
+
+    // Protected data
+
+        //- Reference to the mesh
+        const fvMesh& mesh_;
+
+        volScalarField delta_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("SGSdelta");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            SGSdelta,
+            dictionary,
+            (
+                const word& name,
+                const fvMesh& mesh,
+                const dictionary& dict
+            ),
+            (name, mesh, dict)
+        );
+
+
+    // Constructors
+
+        //- Construct from name and mesh
+        SGSdelta
+        (
+            const word& name,
+            const fvMesh& mesh
+        );
+
+        //- Disallow default bitwise copy construction
+        SGSdelta(const SGSdelta&) = delete;
+
+    // Selectors
+
+        //- Return a reference to the selected SGS delta
+        static autoPtr<SGSdelta> New
+        (
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~SGSdelta()
+    {}
+
+
+    // Member Functions
+
+        //- Read the SGSdelta dictionary
+        virtual void read(const dictionary&) = 0;
+
+        // Correct values
+        virtual void correct() = 0;
+
+
+    // Member Operators
+
+        void operator=(const SGSdelta&) = delete;
+
+        virtual operator const volScalarField&() const
+        {
+            return delta_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.C
new file mode 100644
index 000000000..c5ed71b21
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.C
@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | (C) 2019 Helmholtz-Zentrum Dresden-Rossendorf
+  ==  == ====== ====   ==  ==  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "cubeRootVolDelta.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace SGSModels
+{
+    defineTypeNameAndDebug(cubeRootVolDelta, 0);
+    addToRunTimeSelectionTable(SGSdelta, cubeRootVolDelta, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::SGSModels::cubeRootVolDelta::calcDelta()
+{
+    label nD = mesh_.nGeometricD();
+
+    if (nD == 3)
+    {
+        delta_.primitiveFieldRef() = deltaCoeff_*pow(mesh_.V(), 1.0/3.0);
+    }
+    else if (nD == 2)
+    {
+        WarningInFunction
+            << "Case is 2D, SGS is not strictly applicable\n"
+            << endl;
+
+        const Vector<label>& directions = mesh_.geometricD();
+
+        scalar thickness = 0.0;
+        for (direction dir=0; dir<directions.nComponents; dir++)
+        {
+            if (directions[dir] == -1)
+            {
+                thickness = mesh_.bounds().span()[dir];
+                break;
+            }
+        }
+
+        delta_.primitiveFieldRef() = deltaCoeff_*sqrt(mesh_.V()/thickness);
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Case is not 3D or 2D, this 'delta' method is not applicable"
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::SGSModels::cubeRootVolDelta::cubeRootVolDelta
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    SGSdelta(name, mesh),
+    deltaCoeff_
+    (
+        dict.optionalSubDict(type() + "Coeffs").lookupOrDefault<scalar>
+        (
+            "deltaCoeff",
+            1
+        )
+    )
+{
+    calcDelta();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::SGSModels::cubeRootVolDelta::read(const dictionary& dict)
+{
+    dict.optionalSubDict(type() + "Coeffs").readIfPresent<scalar>
+    (
+        "deltaCoeff",
+        deltaCoeff_
+    );
+
+    calcDelta();
+}
+
+
+void Foam::SGSModels::cubeRootVolDelta::correct()
+{
+    if (mesh_.changing())
+    {
+        calcDelta();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.H
new file mode 100644
index 000000000..4f8682c01
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.H
@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | (C) 2019 Helmholtz-Zentrum Dresden-Rossendorf
+  ==  == ====== ====   ==  ==  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::SGSModels::cubeRootVolDelta
+
+Description
+    Simple cube-root of cell volume delta used in SGS models.
+
+See also
+    Foam::LESModels::cubeRootVolDelta
+
+SourceFiles
+    cubeRootVolDelta.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cubeRootVolDelta_H
+#define cubeRootVolDelta_H
+
+#include "SGSdelta.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace SGSModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class cubeRootVolDelta Declaration
+\*---------------------------------------------------------------------------*/
+
+class cubeRootVolDelta
+:
+    public SGSdelta
+{
+    // Private Data
+
+        scalar deltaCoeff_;
+
+
+    // Private Member Functions
+
+        // Calculate the delta values
+        void calcDelta();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("cubeRootVol");
+
+
+    // Constructors
+
+        //- Construct from name, turbulenceModel and dictionary
+        cubeRootVolDelta
+        (
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary&
+        );
+
+        //- Disallow default bitwise copy construction
+        cubeRootVolDelta(const cubeRootVolDelta&) = delete;
+
+
+    //- Destructor
+    virtual ~cubeRootVolDelta()
+    {}
+
+
+    // Member Functions
+
+        //- Read the SGSdelta dictionary
+        virtual void read(const dictionary&);
+
+        // Correct values
+        virtual void correct();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const cubeRootVolDelta&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace SGSModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/maxDeltaxyz/maxDeltaxyz.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/maxDeltaxyz/maxDeltaxyz.C
new file mode 100644
index 000000000..831c73814
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/maxDeltaxyz/maxDeltaxyz.C
@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | (C) 2019 Helmholtz-Zentrum Dresden-Rossendorf
+  ==  == ====== ====   ==  ==  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "maxDeltaxyz.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace SGSModels
+{
+    defineTypeNameAndDebug(maxDeltaxyz, 0);
+    addToRunTimeSelectionTable(SGSdelta, maxDeltaxyz, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::SGSModels::maxDeltaxyz::calcDelta()
+{
+    label nD = mesh_.nGeometricD();
+
+    const cellList& cells = mesh_.cells();
+    scalarField hmax(cells.size());
+
+    forAll(cells,celli)
+    {
+        scalar deltaMaxTmp = 0.0;
+        const labelList& cFaces = mesh_.cells()[celli];
+        const point& centrevector = mesh_.cellCentres()[celli];
+
+        forAll(cFaces, cFacei)
+        {
+            label facei = cFaces[cFacei];
+            const point& facevector = mesh_.faceCentres()[facei];
+            scalar tmp = mag(facevector - centrevector);
+            if (tmp > deltaMaxTmp)
+            {
+                deltaMaxTmp = tmp;
+            }
+        }
+
+        hmax[celli] = deltaCoeff_*deltaMaxTmp;
+    }
+
+    if (nD == 3)
+    {
+        delta_.primitiveFieldRef() = hmax;
+    }
+    else if (nD == 2)
+    {
+        delta_.primitiveFieldRef() = hmax;
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Case is not 3D or 2D, this 'delta' method is not applicable"
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::SGSModels::maxDeltaxyz::maxDeltaxyz
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    SGSdelta(name, mesh),
+    deltaCoeff_
+    (
+        dict.optionalSubDict(type() + "Coeffs").lookupOrDefault<scalar>
+        (
+            "deltaCoeff",
+            1
+        )
+    )
+{
+    calcDelta();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::SGSModels::maxDeltaxyz::read(const dictionary& dict)
+{
+    const dictionary& coeffsDict(dict.optionalSubDict(type() + "Coeffs"));
+
+    coeffsDict.readIfPresent<scalar>("deltaCoeff", deltaCoeff_);
+
+    calcDelta();
+}
+
+
+void Foam::SGSModels::maxDeltaxyz::correct()
+{
+    if (mesh_.changing())
+    {
+        calcDelta();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/maxDeltaxyz/maxDeltaxyz.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/maxDeltaxyz/maxDeltaxyz.H
new file mode 100644
index 000000000..f73a7e49d
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/maxDeltaxyz/maxDeltaxyz.H
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | (C) 2019 Helmholtz-Zentrum Dresden-Rossendorf
+  ==  == ====== ====   ==  ==  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::SGSModels::maxDeltaxyz
+
+Description
+    Delta calculated by taking the maximum distance between the cell centre
+    and any face centre.  For a regular hex cell, the computed delta will
+    equate to half of the cell width; accordingly, the deltaCoeff model
+    coefficient should be set to 2 for this case.
+
+See also
+    Foam::LESModels::maxDeltaxyz
+
+SourceFiles
+    maxDeltaxyz.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef maxDeltaxyz_H
+#define maxDeltaxyz_H
+
+#include "SGSdelta.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace SGSModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class maxDeltaxyz Declaration
+\*---------------------------------------------------------------------------*/
+
+class maxDeltaxyz
+:
+    public SGSdelta
+{
+    // Private Data
+
+        //- Model coefficient
+        scalar deltaCoeff_;
+
+
+    // Private Member Functions
+
+        // Calculate the delta values
+        void calcDelta();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("maxDeltaxyz");
+
+
+    // Constructors
+
+        //- Construct from name, turbulenceModel and dictionary
+        maxDeltaxyz
+        (
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary&
+        );
+
+        //- Disallow default bitwise copy construction
+        maxDeltaxyz(const maxDeltaxyz&) = delete;
+
+
+    //- Destructor
+    virtual ~maxDeltaxyz()
+    {}
+
+
+    // Member Functions
+
+        //- Read the SGSdelta dictionary
+        virtual void read(const dictionary&);
+
+        // Correct values
+        virtual void correct();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const maxDeltaxyz&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace SGSModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/oneDimensionalDelta/oneDimensionalDelta.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/oneDimensionalDelta/oneDimensionalDelta.C
new file mode 100644
index 000000000..083a64360
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/oneDimensionalDelta/oneDimensionalDelta.C
@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | (C) 2020 Helmholtz-Zentrum Dresden-Rossendorf
+  ==  == ====== ====   ==  ==  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "oneDimensionalDelta.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace SGSModels
+{
+    defineTypeNameAndDebug(oneDimensionalDelta, 0);
+    addToRunTimeSelectionTable(SGSdelta, oneDimensionalDelta, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::SGSModels::oneDimensionalDelta::calcDelta()
+{
+    const label nD = mesh_.nGeometricD();
+
+    if (nD != 1)
+    {
+        FatalErrorInFunction
+            << "Case is not 1D, this 'delta' method is not applicable"
+            << exit(FatalError);
+    }
+
+    const Vector<label>& directions = mesh_.geometricD();
+
+    // Get non-empty direction
+    vector basevector(0, 0, 0);
+    for (direction dir=0; dir<directions.nComponents; dir++)
+    {
+        if (directions[dir] == 1)
+        {
+            basevector[dir] = 1;
+            break;
+        }
+    }
+
+    const cellList& cells = mesh_.cells();
+    scalarField hmax(cells.size());
+
+    forAll(cells,celli)
+    {
+        scalar minCoord = great;
+        scalar maxCoord = 0.0;
+        const labelList& cFaces = mesh_.cells()[celli];
+
+        forAll(cFaces, cFacei)
+        {
+            label facei = cFaces[cFacei];
+            const point& facevector = mesh_.faceCentres()[facei];
+            const scalar tmp = facevector & basevector;
+            if (tmp > maxCoord)
+            {
+                maxCoord = tmp;
+            }
+
+            if (tmp < minCoord)
+            {
+                minCoord = tmp;
+            }
+        }
+
+        hmax[celli] = maxCoord - minCoord;
+    }
+
+    delta_.primitiveFieldRef() = hmax;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::SGSModels::oneDimensionalDelta::oneDimensionalDelta
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    SGSdelta(name, mesh)
+{
+    calcDelta();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::SGSModels::oneDimensionalDelta::correct()
+{
+    if (mesh_.changing())
+    {
+        calcDelta();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/oneDimensionalDelta/oneDimensionalDelta.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/oneDimensionalDelta/oneDimensionalDelta.H
new file mode 100644
index 000000000..fc8704c77
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/oneDimensionalDelta/oneDimensionalDelta.H
@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | (C) 2020 Helmholtz-Zentrum Dresden-Rossendorf
+  ==  == ====== ====   ==  ==  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::SGSModels::oneDimensionalDelta
+
+Description
+    Return cell dimension in one-dimensional space.
+
+SourceFiles
+    oneDimensionalDelta.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef oneDimensionalDelta_H
+#define oneDimensionalDelta_H
+
+#include "SGSdelta.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace SGSModels
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class oneDimensionalDelta Declaration
+\*---------------------------------------------------------------------------*/
+
+class oneDimensionalDelta
+:
+    public SGSdelta
+{
+    // Private Member Functions
+
+        // Calculate the delta values
+        void calcDelta();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("oneDimensional");
+
+
+    // Constructors
+
+        //- Construct from name, turbulenceModel and dictionary
+        oneDimensionalDelta
+        (
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary&
+        );
+
+        //- Disallow default bitwise copy construction
+        oneDimensionalDelta(const oneDimensionalDelta&) = delete;
+
+
+    //- Destructor
+    virtual ~oneDimensionalDelta()
+    {}
+
+
+    // Member Functions
+
+        //- Read the SGSdelta dictionary
+        virtual void read(const dictionary&)
+        {}
+
+        // Correct values
+        virtual void correct();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const oneDimensionalDelta&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace SGSModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/Make/files b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/Make/files
index 0f4a8754d..2f63f5dfc 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/Make/files
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/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/swarmCorrectionNew.C
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.C
new file mode 100644
index 000000000..412f482d2
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/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/multiphaseEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.H
new file mode 100644
index 000000000..473c91ad9
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/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/multiphaseEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.C
new file mode 100644
index 000000000..aa1c8f26b
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/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/multiphaseEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.H
new file mode 100644
index 000000000..9250a4f39
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/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/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseSystems/multiphaseSystems.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseSystems/multiphaseSystems.C
index f7ccafa7b..e9d9d6490 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseSystems/multiphaseSystems.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseSystems/multiphaseSystems.C
@@ -31,6 +31,7 @@ License
 #include "InterfaceCompositionPhaseChangePhaseSystem.H"
 #include "PopulationBalancePhaseSystem.H"
 #include "ThermalPhaseChangePhaseSystem.H"
+#include "ResolvedPhaseSystem.H"
 
 #include "addToRunTimeSelectionTable.H"
 
@@ -43,7 +44,10 @@ namespace Foam
         <
             OneResistanceHeatTransferPhaseSystem
             <
-                MomentumTransferPhaseSystem<phaseSystem>
+                    MomentumTransferPhaseSystem
+                    <
+                        ResolvedPhaseSystem<phaseSystem>
+                    >
             >
         >
         basicMultiphaseSystem;
@@ -63,7 +67,10 @@ namespace Foam
             <
                 TwoResistanceHeatTransferPhaseSystem
                 <
-                    MomentumTransferPhaseSystem<phaseSystem>
+                        MomentumTransferPhaseSystem
+                        <
+                            ResolvedPhaseSystem<phaseSystem>
+                        >
                 >
             >
         >
@@ -84,7 +91,10 @@ namespace Foam
             <
                 TwoResistanceHeatTransferPhaseSystem
                 <
-                    MomentumTransferPhaseSystem<phaseSystem>
+                        MomentumTransferPhaseSystem
+                        <
+                            ResolvedPhaseSystem<phaseSystem>
+                        >
                 >
             >
         >
@@ -105,7 +115,10 @@ namespace Foam
             <
                 OneResistanceHeatTransferPhaseSystem
                 <
-                    MomentumTransferPhaseSystem<phaseSystem>
+                        MomentumTransferPhaseSystem
+                        <
+                            ResolvedPhaseSystem<phaseSystem>
+                        >
                 >
             >
         >
@@ -128,7 +141,10 @@ namespace Foam
                 <
                     TwoResistanceHeatTransferPhaseSystem
                     <
-                        MomentumTransferPhaseSystem<phaseSystem>
+                            MomentumTransferPhaseSystem
+                            <
+                                ResolvedPhaseSystem<phaseSystem>
+                            >
                     >
                 >
             >
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.C
new file mode 100644
index 000000000..d3b7e5dc3
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/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/multiphaseEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.H
new file mode 100644
index 000000000..116c9ceb2
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.H
@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | 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> '.
+
+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/multiphaseEulerFoam/phaseSystems/Make/files b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/Make/files
index e3c9be13a..214a1e564 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/Make/files
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/Make/files
@@ -86,7 +86,23 @@ BlendedInterfacialModel/blendingMethods/blendingMethod/blendingMethodNew.C
 BlendedInterfacialModel/blendingMethods/noBlending/noBlending.C
 BlendedInterfacialModel/blendingMethods/linear/linear.C
 BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.C
+BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.C
 
 alphaContactAngle/alphaContactAngleFvPatchScalarField.C
 
+PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModel.C
+PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModelNew.C
+PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/noInterfaceCompression/noInterfaceCompression.C
+PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/constantInterfaceCompression/constantInterfaceCompression.C
+
+PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModel.C
+PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModelNew.C
+PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/noSurfaceTensionForce/noSurfaceTensionForce.C
+PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/continuumSurfaceForce/continuumSurfaceForce.C
+
+PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethod.C
+PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethodNew.C
+PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/useBlendingMethod/useBlendingMethod.C
+PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/gradAlpha/gradAlpha.C
+
 LIB = $(FOAM_LIBBIN)/libphaseSystem
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index f8703293e..74c4026e9 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
@@ -193,6 +193,89 @@ MomentumTransferPhaseSystem
         "turbulentDispersion",
         turbulentDispersionModels_
     );
+    forAllConstIter
+    (
+        liftModelTable,
+        liftModels_,
+        liftModelIter
+    )
+    {
+        const phasePair& pair(this->phasePairs_[liftModelIter.key()]);
+
+        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()]
+        );
+
+        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()]
+        );
+
+        Ftfs_.insert
+        (
+            pair,
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("Ftf", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh()
+                ),
+                this->mesh(),
+                dimensionedScalar(turbulentDispersionModel::dimF*dimArea, 0)
+            )
+        );
+    }
 
     forAllConstIter
     (
@@ -385,6 +468,41 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer()
         }
     }
 
+    // Update the lift forces
+    forAllConstIter
+    (
+        liftModelTable,
+        liftModels_,
+        liftModelIter
+    )
+    {
+        *Flfs_[liftModelIter.key()] = liftModelIter()->Ff();
+    }
+
+    // Update the wall lubrication forces
+    forAllConstIter
+    (
+        wallLubricationModelTable,
+        wallLubricationModels_,
+        wallLubricationModelIter
+    )
+    {
+        *Fwfs_[wallLubricationModelIter.key()] =
+            wallLubricationModelIter()->Ff();
+    }
+
+    // Update the turbulent dispersion coefficients
+    forAllConstIter
+    (
+        turbulentDispersionModelTable,
+        turbulentDispersionModels_,
+        turbulentDispersionModelIter
+    )
+    {
+        *Ftfs_[turbulentDispersionModelIter.key()] =
+            turbulentDispersionModelIter()->Ff();
+    }
+
     return eqnsPtr;
 }
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
index ef82b54b3..811bd160f 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
@@ -134,8 +134,46 @@ private:
         > turbulentDispersionModelTable;
 
 
+protected:
+
+    // Protected typedefs
+
+        typedef HashPtrTable
+        <
+            surfaceScalarField,
+            phasePairKey,
+            phasePairKey::hash
+        > FlfTable;
+
+        typedef HashPtrTable
+        <
+            surfaceScalarField,
+            phasePairKey,
+            phasePairKey::hash
+        > FwfTable;
+
+        typedef HashPtrTable
+        <
+            surfaceScalarField,
+            phasePairKey,
+            phasePairKey::hash
+        > FtfTable;
+
+
+private:
+
+
     // Private Data
 
+        //- Face lift forces
+        FlfTable Flfs_;
+
+        //- Face wall lubrication forces
+        FwfTable Fwfs_;
+
+        //- Face Turbulent dispersion forces
+        FtfTable Ftfs_;
+
         //- Drag coefficients
         KdTable Kds_;
 
@@ -297,6 +335,48 @@ public:
 
         //- Read base phaseProperties dictionary
         virtual bool read();
+
+        //- 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 face lift force for all phase-pairs
+        virtual const FlfTable& Flfs() const
+        {
+            return Flfs_;
+        };
+
+        //- Return the face wall lubrication force for all phase-pairs
+        virtual const FwfTable& Fwfs() const
+        {
+            return Fwfs_;
+        };
+
+        //- Return the face turbulent dispersion force for all phase-pairs
+        virtual const FtfTable& Ftfs() const
+        {
+            return Ftfs_;
+        };
 };
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.H
index 0fcb1d034..2f11c11c5 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.H
@@ -145,6 +145,11 @@ public:
 
         //- Read base phaseProperties dictionary
         virtual bool read();
+
+        virtual const phaseSystem::dmdtfTable& dmdtfsPhaseTransfer() const
+        {
+            return dmdtfs_;
+        }
 };
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/ResolvedPhaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/ResolvedPhaseSystem.C
new file mode 100644
index 000000000..d390aef0d
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/ResolvedPhaseSystem.C
@@ -0,0 +1,220 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | 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 "ResolvedPhaseSystem.H"
+
+#include "interfaceCompressionModel.H"
+#include "surfaceTensionForceModel.H"
+#include "interfaceIndicatorMethod.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasePhaseSystem>
+Foam::ResolvedPhaseSystem<BasePhaseSystem>::ResolvedPhaseSystem
+(
+    const fvMesh& mesh
+)
+:
+    BasePhaseSystem(mesh)
+{
+    this->generatePairsAndSubModels
+    (
+        "interfaceCompression",
+        interfaceCompressionModels_
+    );
+
+    surfaceTensionForceModel_ = surfaceTensionForceModel::New(*this);
+
+    interfaceIndicatorMethod_ = interfaceIndicatorMethod::New(*this);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class BasePhaseSystem>
+Foam::ResolvedPhaseSystem<BasePhaseSystem>::~ResolvedPhaseSystem()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class BasePhaseSystem>
+Foam::tmp<Foam::surfaceScalarField>
+Foam::ResolvedPhaseSystem<BasePhaseSystem>::phir
+(
+    const phaseModel& phase1,
+    const phaseModel& phase2
+) const
+{
+    tmp<surfaceScalarField> phir(BasePhaseSystem::phir(phase1, phase2));
+
+    phasePairKey key(phase1.name(), phase2.name());
+
+    if (interfaceCompressionModels_.found(key))
+    {
+        if (interfaceCompressionModels_[key]().type() != "none")
+        {
+            phir.ref() += interfaceCompressionModels_[key]().phirc(phase1);
+        }
+    }
+
+    return phir;
+}
+
+
+template<class BasePhaseSystem>
+void Foam::ResolvedPhaseSystem<BasePhaseSystem>::addCompressionFlux
+(
+    const phaseModel& phase,
+    surfaceScalarField& alphaPhi
+) const
+{
+    BasePhaseSystem::addCompressionFlux(phase, alphaPhi);
+
+    if (!interfaceCompressionModels_.empty())
+    {
+        const volScalarField& alpha = phase;
+
+        forAll(this->phases(), phasei)
+        {
+            const phaseModel& phase2 = this->phases()[phasei];
+            const volScalarField& alpha2 = phase2;
+            const phasePairKey key(phase.name(), phase2.name());
+
+            if (&phase2 == &phase) continue;
+
+            interfaceCompressionModelTable::const_iterator cAlpha
+            (
+                interfaceCompressionModels_.find(key)
+            );
+
+            if (cAlpha != interfaceCompressionModels_.end())
+            {
+                const surfaceScalarField phir
+                (
+                    interfaceCompressionModels_[key]().phirc(phase)
+                );
+
+                const word phirScheme
+                (
+                    "div(phir,"
+                  + alpha2.name() + ',' + alpha.name()
+                  + ')'
+                );
+
+                alphaPhi += fvc::flux
+                (
+                    -fvc::flux(-phir, alpha2, phirScheme),
+                    alpha,
+                    phirScheme
+                );
+            }
+        }
+    }
+}
+
+
+template<class BasePhaseSystem>
+Foam::tmp<Foam::surfaceScalarField>
+Foam::ResolvedPhaseSystem<BasePhaseSystem>::surfaceTension
+(
+    const phaseModel& phase1
+) const
+{
+    tmp<surfaceScalarField> tSurfaceTension
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "surfaceTension",
+                this->mesh().time().timeName(),
+                this->mesh()
+            ),
+            this->mesh(),
+            dimensionedScalar
+            (
+                "surfaceTension",
+                dimensionSet(1, -2, -2, 0, 0),
+                0
+            )
+        )
+    );
+
+    forAll(this->phases(), phasej)
+    {
+        const phaseModel& phase2 = this->phases()[phasej];
+
+        if (&phase2 != &phase1)
+        {
+            phasePairKey key12(phase1.name(), phase2.name());
+
+            if (interfaceCompressionModels_.found(key12))
+            {
+                if (interfaceCompressionModels_[key12]().type() != "none")
+                {
+                    const phasePair& pair(this->phasePairs_[key12]);
+
+                    tSurfaceTension.ref() +=
+                        surfaceTensionForceModel_->surfaceTensionForcef(pair);
+                }
+            }
+        }
+    }
+
+    return tSurfaceTension;
+}
+
+
+template<class BasePhaseSystem>
+Foam::tmp<Foam::volScalarField>
+Foam::ResolvedPhaseSystem<BasePhaseSystem>::interfaceIndicator
+(
+    const phaseModel& phase
+) const
+{
+   return interfaceIndicatorMethod_->interfaceIndicator(phase);
+}
+
+
+template<class BasePhaseSystem>
+bool Foam::ResolvedPhaseSystem<BasePhaseSystem>::read()
+{
+    if (BasePhaseSystem::read())
+    {
+        bool readOK = true;
+
+        // Read models ...
+
+        return readOK;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/ResolvedPhaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/ResolvedPhaseSystem.H
new file mode 100644
index 000000000..718c98886
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/ResolvedPhaseSystem.H
@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | 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::ResolvedPhaseSystem
+
+Description
+    Class which adds functionality for modeling resolved interfaces such as
+    interface compression, surface tension and contact angle correction.
+
+SourceFiles
+    ResolvedPhaseSystem.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ResolvedPhaseSystem_H
+#define ResolvedPhaseSystem_H
+
+#include "phaseSystem.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class interfaceCompressionModel;
+class surfaceTensionForceModel;
+class interfaceIndicatorMethod;
+
+/*---------------------------------------------------------------------------*\
+                     Class ResolvedPhaseSystem Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasePhaseSystem>
+class ResolvedPhaseSystem
+:
+    public BasePhaseSystem
+{
+private:
+
+    // Private Typedefs
+
+        typedef HashTable
+        <
+            autoPtr<interfaceCompressionModel>,
+            phasePairKey,
+            phasePairKey::hash
+        > interfaceCompressionModelTable;
+
+
+
+    // Private Data
+
+        //- Interface compression models
+        interfaceCompressionModelTable interfaceCompressionModels_;
+
+        //- Surface tension force models
+        autoPtr<surfaceTensionForceModel> surfaceTensionForceModel_;
+
+        //- Resolved interface detection method
+        autoPtr<interfaceIndicatorMethod> interfaceIndicatorMethod_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh
+        ResolvedPhaseSystem(const fvMesh&);
+
+
+    //- Destructor
+    virtual ~ResolvedPhaseSystem();
+
+
+    // Member Functions
+
+        //- Return relative flux for interface compression
+        virtual tmp<surfaceScalarField> phir
+        (
+            const phaseModel& phase1,
+            const phaseModel& phase2
+        ) const;
+
+        //- Add interface compression flux
+        virtual void addCompressionFlux
+        (
+            const phaseModel& phase,
+            surfaceScalarField& alphaPhi
+        ) const;
+
+        //- Return the surface tension force
+        tmp<surfaceScalarField> surfaceTension(const phaseModel& phase) const;
+
+        //- Indicator field of the proximity of a resolved phase interface
+        //  Values are distributed from 0 (far) to 1 (near).
+        virtual tmp<volScalarField> interfaceIndicator
+        (
+            const phaseModel& phase
+        ) const;
+
+        //- Read base phaseProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "ResolvedPhaseSystem.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/constantInterfaceCompression/constantInterfaceCompression.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/constantInterfaceCompression/constantInterfaceCompression.C
new file mode 100644
index 000000000..d03146971
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/constantInterfaceCompression/constantInterfaceCompression.C
@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | 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 "constantInterfaceCompression.H"
+#include "phasePair.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace interfaceCompressionModels
+{
+    defineTypeNameAndDebug(constantInterfaceCompression, 0);
+    addToRunTimeSelectionTable
+    (
+        interfaceCompressionModel,
+        constantInterfaceCompression,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::interfaceCompressionModels::constantInterfaceCompression::
+constantInterfaceCompression
+(
+    const dictionary& dict,
+    const phasePair& pair
+)
+:
+    interfaceCompressionModel(dict, pair),
+    C_("C", dimless, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::interfaceCompressionModels::constantInterfaceCompression::
+~constantInterfaceCompression()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::interfaceCompressionModels::constantInterfaceCompression::C() const
+{
+    return
+        tmp<surfaceScalarField>
+        (
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    "zero",
+                    mesh_.time().timeName(),
+                    mesh_
+                ),
+                mesh_,
+                C_
+            )
+        );
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/constantInterfaceCompression/constantInterfaceCompression.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/constantInterfaceCompression/constantInterfaceCompression.H
new file mode 100644
index 000000000..ef2a7ad67
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/constantInterfaceCompression/constantInterfaceCompression.H
@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | 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::interfaceCompressionModels::constantInterfaceCompression
+
+Description
+    Interface compression with a constant, user-specified coefficient.
+
+SourceFiles
+    constantInterfaceCompression.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantInterfaceCompression_H
+#define constantInterfaceCompression_H
+
+#include "interfaceCompressionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class phasePair;
+
+namespace interfaceCompressionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                Class constantInterfaceCompression Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantInterfaceCompression
+:
+    public interfaceCompressionModel
+{
+    // Private Data
+
+        //- Constant interface compression coefficient
+        const dimensionedScalar C_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constantCoefficient");
+
+
+    // Constructors
+
+        //- Construct from a dictionary and a phase pair
+        constantInterfaceCompression
+        (
+            const dictionary& dict,
+            const phasePair& pair
+        );
+
+
+    //- Destructor
+    virtual ~constantInterfaceCompression();
+
+
+    // Member Functions
+
+        //- Interface compression coefficient
+        virtual tmp<surfaceScalarField> C() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace interfaceCompressionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModel.C
new file mode 100644
index 000000000..8b387da07
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModel.C
@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | 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 "interfaceCompressionModel.H"
+#include "phasePair.H"
+#include "phaseSystem.H"
+#include "fvcGrad.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(interfaceCompressionModel, 0);
+    defineRunTimeSelectionTable(interfaceCompressionModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::interfaceCompressionModel::interfaceCompressionModel
+(
+    const dictionary& dict,
+    const phasePair& pair
+)
+:
+    pair_(pair),
+    mesh_(pair_.phase1().mesh()),
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(mesh_.V()), 1.0/3.0)
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::interfaceCompressionModel::~interfaceCompressionModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::surfaceScalarField> Foam::interfaceCompressionModel::nHatf
+(
+    const phaseModel& phase
+) const
+{
+    const volScalarField& alpha1 = phase;
+    const volScalarField& alpha2 = pair_.otherPhase(phase);
+
+    surfaceVectorField gradAlphaf
+    (
+        fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1))
+      - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
+    );
+
+    return gradAlphaf/(mag(gradAlphaf) + deltaN_) & mesh_.Sf();
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::interfaceCompressionModel::phirc
+(
+    const phaseModel& phase
+) const
+{
+    surfaceScalarField phic
+    (
+        (
+            mag(phase.fluid().phi())
+          + mag(phase.phi() - pair_.otherPhase(phase).phi())
+        )
+       /mesh_.magSf()
+    );
+
+    return min(C()*phic, max(phic))*nHatf(phase);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModel.H
new file mode 100644
index 000000000..47191d7e1
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModel.H
@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | 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::interfaceCompressionModel
+
+Description
+    Base class for modeling the interface compression flux which is
+    added to the correction fluxes within multiphaseSystem::solveAlphas().
+
+SourceFiles
+    interfaceCompressionModel.C
+    interfaceCompressionModelNew.C
+    phaseSystem.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interfaceCompressionModel_H
+#define interfaceCompressionModel_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "volFields.H"
+#include "dictionary.H"
+#include "runTimeSelectionTables.H"
+
+namespace Foam
+{
+
+class phasePair;
+class phaseModel;
+
+/*---------------------------------------------------------------------------*\
+                  Class interfaceCompressionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class interfaceCompressionModel
+{
+protected:
+
+    // Protected data
+
+        //- Phase pair
+        const phasePair& pair_;
+
+        //- Reference to the mesh
+        const fvMesh& mesh_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("interfaceCompressionModel");
+
+
+    // Declare runtime construction
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            interfaceCompressionModel,
+            dictionary,
+            (
+                const dictionary& dict,
+                const phasePair& pair
+            ),
+            (dict, pair)
+        );
+
+
+    // Constructors
+
+        //- Construct from a dictionary and a phase pair
+        interfaceCompressionModel
+        (
+            const dictionary& dict,
+            const phasePair& pair
+        );
+
+
+    //- Destructor
+    virtual ~interfaceCompressionModel();
+
+
+    // Selectors
+
+        static autoPtr<interfaceCompressionModel> New
+        (
+            const dictionary& dict,
+            const phasePair& pair
+        );
+
+
+    // Member Functions
+
+        //- Interface compression coefficient
+        virtual tmp<surfaceScalarField> C() const = 0;
+
+        //- Face unit interface normal flux
+        virtual tmp<surfaceScalarField> nHatf(const phaseModel& phase) const;
+
+        //- Interface compression flux
+        virtual tmp<surfaceScalarField> phirc(const phaseModel& phase) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModelNew.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModelNew.C
new file mode 100644
index 000000000..be059471f
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModelNew.C
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | 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 "interfaceCompressionModel.H"
+#include "phasePair.H"
+
+// * * * * * * * * * * * * * * * * Selector  * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::interfaceCompressionModel>
+Foam::interfaceCompressionModel::New
+(
+    const dictionary& dict,
+    const phasePair& pair
+)
+{
+    word interfaceCompressionModelType(dict.lookup("type"));
+
+    Info<< "Selecting interfaceCompressionModel for "
+        << pair << ": " << interfaceCompressionModelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(interfaceCompressionModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown interfaceCompressionModelType type "
+            << interfaceCompressionModelType << endl << endl
+            << "Valid interfaceCompressionModel types are : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return cstrIter()(dict, pair);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/noInterfaceCompression/noInterfaceCompression.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/noInterfaceCompression/noInterfaceCompression.C
new file mode 100644
index 000000000..1d38b4560
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/noInterfaceCompression/noInterfaceCompression.C
@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | 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 "noInterfaceCompression.H"
+#include "phasePair.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace interfaceCompressionModels
+{
+    defineTypeNameAndDebug(noInterfaceCompression, 0);
+    addToRunTimeSelectionTable
+    (
+        interfaceCompressionModel,
+        noInterfaceCompression,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::interfaceCompressionModels::noInterfaceCompression::noInterfaceCompression
+(
+    const dictionary& dict,
+    const phasePair& pair
+)
+:
+    interfaceCompressionModel(dict, pair)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::interfaceCompressionModels::noInterfaceCompression::
+~noInterfaceCompression()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::interfaceCompressionModels::noInterfaceCompression::C() const
+{
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "C",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::interfaceCompressionModels::noInterfaceCompression::phir
+(
+    const phaseModel& phase
+) const
+{
+    return
+        tmp<surfaceScalarField>
+        (
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    "phir",
+                    mesh_.time().timeName(),
+                    mesh_
+                ),
+                mesh_,
+                dimensionedScalar(dimVolume/dimTime, Zero)
+            )
+        );
+}
+
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/noInterfaceCompression/noInterfaceCompression.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/noInterfaceCompression/noInterfaceCompression.H
new file mode 100644
index 000000000..abbc5e50e
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/noInterfaceCompression/noInterfaceCompression.H
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | 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::interfaceCompressionModels::noInterfaceCompression
+
+Description
+
+SourceFiles
+    noInterfaceCompression.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noInterfaceCompression_H
+#define noInterfaceCompression_H
+
+#include "interfaceCompressionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class phasePair;
+
+namespace interfaceCompressionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class noInterfaceCompression Declaration
+\*---------------------------------------------------------------------------*/
+
+class noInterfaceCompression
+:
+    public interfaceCompressionModel
+{
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+        //- Construct from a dictionary and a phase pair
+        noInterfaceCompression
+        (
+            const dictionary& dict,
+            const phasePair& pair
+        );
+
+
+    //- Destructor
+    virtual ~noInterfaceCompression();
+
+
+    // Member Functions
+
+        //- Interface compression coefficient
+        virtual tmp<surfaceScalarField> C() const;
+
+        //- Interface compression flux
+        virtual tmp<surfaceScalarField> phir(const phaseModel& phase) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace interfaceCompressionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/gradAlpha/gradAlpha.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/gradAlpha/gradAlpha.C
new file mode 100644
index 000000000..ec049a1a8
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/gradAlpha/gradAlpha.C
@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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 "addToRunTimeSelectionTable.H"
+#include "fvcGrad.H"
+#include "gradAlpha.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace interfaceIndicatorMethods
+{
+    defineTypeNameAndDebug(gradAlpha, 0);
+    addToRunTimeSelectionTable
+    (
+        interfaceIndicatorMethod,
+        gradAlpha,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::volVectorField>
+Foam::interfaceIndicatorMethods::gradAlpha::gradAlphaPair
+(
+    const volScalarField& phase1,
+    const volScalarField& phase2
+) const
+{
+    return phase2*fvc::grad(phase1) - phase1*fvc::grad(phase2);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::interfaceIndicatorMethods::gradAlpha::interfaceIndicator
+(
+    const phaseModel& phase1,
+    const phaseModel& phase2
+) const
+{
+    // Limit the indicator to 1
+    return min(mag(gradAlphaPair(phase1, phase2))*dx(), 1.0);
+}
+
+
+const Foam::volScalarField&
+Foam::interfaceIndicatorMethods::gradAlpha::dx() const
+{
+    return delta_();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::interfaceIndicatorMethods::gradAlpha::gradAlpha
+(
+    const phaseSystem& fluid
+)
+:
+    interfaceIndicatorMethod(fluid),
+    delta_
+    (
+        SGSdelta::New
+        (
+            "delta",
+            fluid.mesh(),
+            fluid.subDict(interfaceIndicatorMethod::typeName)
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::interfaceIndicatorMethods::gradAlpha::~gradAlpha()
+{}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/gradAlpha/gradAlpha.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/gradAlpha/gradAlpha.H
new file mode 100644
index 000000000..6a94ee5db
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/gradAlpha/gradAlpha.H
@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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::interfaceIndicatorMethods::gradAlpha
+
+Description
+    Determine the interface indicator field for a pair of continuous phases,
+    i and j:
+
+    \f[
+        \mathrm{RI}_{ij} =
+            \mathrm{max}(\mathrm{mag}((\nabla\alpha)_{ij}\Delta x), 1),
+    \f]
+
+    where
+
+    \f[
+        (\nabla\alpha)_{ij} = \alpha_j\nabla\alpha_i - \alpha_i\nabla\alpha_j,
+    \f]
+
+    and \f$\Delta x\f$ is the typical grid size calculated using
+    a SGSdelta model.
+
+Usage
+    \verbatim
+    interfaceIndicatorMethod
+    {
+        type    gradAlpha;
+        delta   cubeRootVol;
+    }
+    \endverbatim
+
+See also
+    Foam::SGSdelta
+
+SourceFiles
+    gradAlpha.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef gradAlpha_H
+#define gradAlpha_H
+
+#include "interfaceIndicatorMethod.H"
+#include "SGSdelta.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace interfaceIndicatorMethods
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class gradAlpha Declaration
+\*---------------------------------------------------------------------------*/
+
+class gradAlpha
+:
+    public interfaceIndicatorMethod
+{
+private:
+
+    // Private Data
+
+        // Local grid size approximation model
+        autoPtr<SGSdelta> delta_;
+
+    // Private Member Functions
+
+        //- Calculate interface gradient between two phases
+        tmp<volVectorField> gradAlphaPair
+        (
+            const volScalarField& phase1, const volScalarField& phase2
+        ) const;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Indicator field between a phase pair
+        virtual tmp<volScalarField> interfaceIndicator
+        (
+            const phaseModel& phase1, const phaseModel& phase2
+        ) const override;
+
+        //- Return local grid size approximation
+        const volScalarField& dx() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("gradAlpha");
+
+
+    // Constructors
+
+        //- Construct from a dictionary from phaseSystem
+        gradAlpha
+        (
+            const phaseSystem& fluid
+        );
+
+
+    //- Destructor
+    virtual ~gradAlpha();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace interfaceIndicatorMethods
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethod.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethod.C
new file mode 100644
index 000000000..ea5cfa374
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethod.C
@@ -0,0 +1,118 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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 "addToRunTimeSelectionTable.H"
+#include "fvcGrad.H"
+#include "interfaceIndicatorMethod.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(interfaceIndicatorMethod, 0);
+    defineRunTimeSelectionTable(interfaceIndicatorMethod, dictionary);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+//- Provides default implementation for pure virtual function
+Foam::tmp<Foam::volScalarField>
+Foam::interfaceIndicatorMethod::
+interfaceIndicator(const phaseModel& phase1, const phaseModel& phase2) const
+{
+    NotImplemented;
+    return tmp<volScalarField>(nullptr);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::interfaceIndicatorMethod::interfaceIndicatorMethod
+(
+    const phaseSystem& fluid
+)
+:
+    fluid_(fluid)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::interfaceIndicatorMethod::~interfaceIndicatorMethod()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::interfaceIndicatorMethod::
+interfaceIndicator(const phaseModel& phase) const
+{
+    // Consider only interfaces between continuous phases
+    if (!phase.continuous())
+    {
+        FatalErrorInFunction
+            << "Attempt to detect interface for non-continuous phase: "
+            << phase.name() << endl
+            << exit(FatalError);
+    }
+    const phaseSystem::phaseModelPartialList& phases(fluid_.continuousPhases());
+
+    tmp<volScalarField> ri
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("interfaceIndicator", phase.name()),
+                fluid_.mesh().time().timeName(),
+                fluid_.mesh()
+            ),
+            fluid_.mesh(),
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+
+    // Use maximum pair contribution
+    forAll(phases, phasei)
+    {
+        const phaseModel& phase2 = phases[phasei];
+
+        if (phase2 != phase)
+        {
+            ri.ref() = max(this->interfaceIndicator(phase, phase2), ri());
+        }
+    }
+
+    if (debug && phase.mesh().time().outputTime())
+    {
+        ri().write();
+    }
+    return ri;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethod.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethod.H
new file mode 100644
index 000000000..b77060552
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethod.H
@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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::interfaceIndicatorMethod
+
+Description
+    Base class for methods to track resolved interfaces in a multiphase system.
+
+    Determine the interface indicator field for a continuous phase i
+    by determining the maximum contribution from interface pairs between
+    phases i and j:
+
+    \f[
+        \mathrm{RI}_i = \mathrm{max}(\mathrm{RI}_{ij})\text{, where }i \neq j.
+    \f]
+
+SourceFiles
+    interfaceIndicatorMethod.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interfaceIndicatorMethod_H
+#define interfaceIndicatorMethod_H
+
+#include "dictionary.H"
+#include "runTimeSelectionTables.H"
+#include "phaseModel.H"
+#include "phaseSystem.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class interfaceIndicatorMethod Declaration
+\*---------------------------------------------------------------------------*/
+
+class interfaceIndicatorMethod
+{
+protected:
+
+    // Protected References
+
+        //- Private reference to the phaseSystem this method belongs to.
+        const phaseSystem& fluid_;
+
+    // Protected Member Functions
+
+        //- Indicator field between a phase pair
+        virtual tmp<volScalarField> interfaceIndicator
+        (
+            const phaseModel& phase1, const phaseModel& phase2
+        ) const = 0;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("interfaceIndicatorMethod");
+
+
+    // Declare runtime construction
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        interfaceIndicatorMethod,
+        dictionary,
+        (
+            const phaseSystem& fluid
+        ),
+        (fluid)
+    );
+
+
+    // Constructors
+
+        //- Construct from a dictionary from phaseSystem
+        interfaceIndicatorMethod
+        (
+            const phaseSystem& fluid
+        );
+
+
+    // Selector
+
+        static autoPtr<interfaceIndicatorMethod> New
+        (
+            const phaseSystem& fluid
+        );
+
+
+    //- Destructor
+    virtual ~interfaceIndicatorMethod();
+
+
+    // Member Functions
+
+        //- Indicator field of the proximity of a resolved phase interface
+        //  Values are distributed from 0 (far) to 1 (near).
+        virtual tmp<volScalarField> interfaceIndicator
+        (
+            const phaseModel& phase
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethodNew.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethodNew.C
new file mode 100644
index 000000000..30aa75c0a
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethodNew.C
@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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 "interfaceIndicatorMethod.H"
+#include "phaseSystem.H"
+
+// Helper class to handle error reporting
+namespace Foam
+{
+namespace interfaceIndicatorMethods
+{
+
+class undefined
+:  public interfaceIndicatorMethod
+{
+protected:
+
+    //- Handles error reporting for missing dict entry in phaseProperties.
+    virtual tmp<volScalarField> interfaceIndicator
+    (
+        const phaseModel& phase1, const phaseModel& phase2
+    ) const
+    {
+        // Use dictionary error reporting to handle the missing dict case
+        word selectedType
+        (
+            fluid_.subDict(interfaceIndicatorMethod::typeName).lookup("type")
+        );
+
+        // Handles other strange cases (which should not occur)
+        // Uses default "not implemented" from base class
+        return interfaceIndicatorMethod::interfaceIndicator(phase1, phase2);
+    }
+
+
+public:
+
+    //- Inherit all constructors from base
+    using interfaceIndicatorMethod::interfaceIndicatorMethod;
+};
+
+}
+}
+
+// * * * * * * * * * * * * * * * * Selector  * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::interfaceIndicatorMethod>
+Foam::interfaceIndicatorMethod::New
+(
+    const phaseSystem& fluid
+)
+{
+    // Select the base method if method information is missing.
+    const word methodDictName = Foam::interfaceIndicatorMethod::typeName;
+    const word typeKeyword = "type";
+    word selectedType;
+    if (fluid.found(methodDictName))
+    {
+        const dictionary& methodDict = fluid.subDict(methodDictName);
+        if (methodDict.found(typeKeyword))
+        {
+            selectedType = methodDict.lookup<word>(typeKeyword);
+            Info<< "Selecting interfaceIndicatorMethod: "
+                << selectedType << endl;
+
+            dictionaryConstructorTable::iterator cstrIter =
+                dictionaryConstructorTablePtr_->find(selectedType);
+
+            if (cstrIter == dictionaryConstructorTablePtr_->end())
+            {
+                FatalErrorInFunction
+                    << "Unknown " << methodDictName << " type "
+                    << selectedType << endl << endl
+                    << "Valid " << methodDictName << " types are : " << endl
+                    << dictionaryConstructorTablePtr_->sortedToc() << endl
+                    << exit(FatalError);
+            }
+
+            return cstrIter()(fluid);
+        }
+    }
+
+    // No user input - select the undefined method which lazily reports errors
+    // if client model code requests a member function of base method.
+    return
+        autoPtr<interfaceIndicatorMethod>
+        (
+            new interfaceIndicatorMethods::undefined(fluid)
+        );
+}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/useBlendingMethod/useBlendingMethod.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/useBlendingMethod/useBlendingMethod.C
new file mode 100644
index 000000000..1d8b52849
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/useBlendingMethod/useBlendingMethod.C
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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 "addToRunTimeSelectionTable.H"
+#include "fvcGrad.H"
+#include "useBlendingMethod.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace interfaceIndicatorMethods
+{
+    defineTypeNameAndDebug(useBlendingMethod, 0);
+    addToRunTimeSelectionTable
+    (
+        interfaceIndicatorMethod,
+        useBlendingMethod,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::interfaceIndicatorMethods::useBlendingMethod::
+interfaceIndicator(const phaseModel& phase1, const phaseModel& phase2) const
+{
+    return
+    (
+          scalar(1)
+        - blendingMethod_->f1(phase1, phase2)
+        - blendingMethod_->f2(phase1, phase2)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::interfaceIndicatorMethods::useBlendingMethod::useBlendingMethod
+(
+    const phaseSystem& fluid
+)
+:
+    interfaceIndicatorMethod(fluid)
+{
+    const dictionary& blendingMethodDict
+    (
+        fluid_.subDict(interfaceIndicatorMethod::typeName).subDict
+        (
+            "blendingMethod"
+
+        )
+    );
+
+    const word blendingName = useBlendingMethod::typeName;
+
+    const phaseSystem::phaseModelPartialList& continuousPhases =
+        fluid_.continuousPhases();
+
+    wordList phaseNames;
+    forAll(continuousPhases, cphasei)
+    {
+        phaseNames.append(continuousPhases[cphasei].name());
+    }
+
+    blendingMethod_ = blendingMethod::New
+    (
+        blendingName,
+        blendingMethodDict,
+        phaseNames
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::interfaceIndicatorMethods::useBlendingMethod::~useBlendingMethod()
+{}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/useBlendingMethod/useBlendingMethod.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/useBlendingMethod/useBlendingMethod.H
new file mode 100644
index 000000000..d6a229292
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/useBlendingMethod/useBlendingMethod.H
@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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::interfaceIndicatorMethods::useBlendingMethod
+
+Description
+    Determine the interface indicator field for a pair of continuous phases,
+    i and j:
+
+    \f[
+        \mathrm{RI}_{ij} = (1 - f_1 - f_2)_{ij},
+    \f]
+
+    where \f$f_1\f$ and \f$f_2\f$ are provided by the selected blendingMethod.
+
+Usage
+    Example setup for continuous air and water phases:
+    \verbatim
+    interfaceIndicatorMethod
+    {
+        type    useBlendingMethod;
+        blendingMethod
+        {
+            type                psiSurf;
+            critAlpha.air       0.3;
+            critAlpha.water     0.3;
+            deltaAlpha.air      0.05;
+            deltaAlpha.water    0.05;
+            f                   0.1;
+            n                   5;
+        }
+    }
+    \endverbatim
+
+See also
+    Foam::blendingMethod
+    Foam::blendingMethods::psiSurf
+
+SourceFiles
+    useBlendingMethod.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef useBlendingMethod_H
+#define useBlendingMethod_H
+
+#include "blendingMethod.H"
+#include "interfaceIndicatorMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace interfaceIndicatorMethods
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class useBlendingMethod Declaration
+\*---------------------------------------------------------------------------*/
+
+class useBlendingMethod
+:
+    public interfaceIndicatorMethod
+{
+private:
+
+    // Private Data
+
+        //- Blending method used for interface detection
+        autoPtr<blendingMethod> blendingMethod_;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Indicator field between a phase pair
+        virtual tmp<volScalarField> interfaceIndicator
+        (
+            const phaseModel& phase1, const phaseModel& phase2
+        ) const override;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("useBlendingMethod");
+
+
+    // Constructors
+
+        //- Construct from a dictionary from phaseSystem
+        useBlendingMethod
+        (
+            const phaseSystem& fluid
+        );
+
+
+    //- Destructor
+    virtual ~useBlendingMethod();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace interfaceIndicatorMethods
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/continuumSurfaceForce/continuumSurfaceForce.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/continuumSurfaceForce/continuumSurfaceForce.C
new file mode 100644
index 000000000..5fb7fd35c
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/continuumSurfaceForce/continuumSurfaceForce.C
@@ -0,0 +1,241 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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 "continuumSurfaceForce.H"
+#include "phasePair.H"
+#include "addToRunTimeSelectionTable.H"
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "unitConversion.H"
+#include "fvc.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceTensionForceModels
+{
+    defineTypeNameAndDebug(continuumSurfaceForce, 0);
+    addToRunTimeSelectionTable
+    (
+        surfaceTensionForceModel,
+        continuumSurfaceForce,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::surfaceVectorField>
+Foam::surfaceTensionForceModels::continuumSurfaceForce::nHatfv
+(
+    const phasePair& pair
+) const
+{
+    surfaceVectorField gradAlphaf
+    (
+        fvc::interpolate(pair.phase2())
+       *fvc::interpolate(fvc::grad(pair.phase1()))
+      - fvc::interpolate(pair.phase1())
+       *fvc::interpolate(fvc::grad(pair.phase2()))
+    );
+
+    return gradAlphaf/(mag(gradAlphaf) + deltaN_);
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::surfaceTensionForceModels::continuumSurfaceForce::nHatf
+(
+    const phasePair& pair
+) const
+{
+    return nHatfv(pair) & this->fluid_.mesh().Sf();
+}
+
+
+// Correction for the boundary condition on the unit normal nHat on
+// walls to produce the correct contact angle.
+
+// The dynamic contact angle is calculated from the component of the
+// velocity on the direction of the interface, parallel to the wall.
+
+void
+Foam::surfaceTensionForceModels::continuumSurfaceForce::correctContactAngle
+(
+    surfaceVectorField::Boundary& nHatb,
+    const phasePair& pair
+) const
+{
+    const volScalarField::Boundary& gbf = pair.phase1().boundaryField();
+
+    const fvBoundaryMesh& boundary = this->fluid_.mesh().boundary();
+
+    forAll(boundary, patchi)
+    {
+        if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
+        {
+            const alphaContactAngleFvPatchScalarField& acap =
+                refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
+
+            vectorField& nHatPatch = nHatb[patchi];
+
+            vectorField AfHatPatch
+            (
+                this->fluid_.mesh().Sf().boundaryField()[patchi]
+               /this->fluid_.mesh().magSf().boundaryField()[patchi]
+            );
+
+            alphaContactAngleFvPatchScalarField::thetaPropsTable::
+                const_iterator tp =
+                    acap.thetaProps()
+                   .find
+                    (
+                        phasePairKey
+                        (
+                            pair.phase1().name(),
+                            pair.phase2().name()
+                        )
+                    );
+
+            if (tp == acap.thetaProps().end())
+            {
+                FatalErrorInFunction
+                    << "Cannot find interface "
+                    << phasePairKey
+                       (
+                            pair.phase1().name(),
+                            pair.phase2().name()
+                       )
+                    << "\n    in table of theta properties for patch "
+                    << acap.patch().name()
+                    << exit(FatalError);
+            }
+
+            bool matched = (tp.key().first() == pair.phase1().name());
+
+            scalar theta0 = degToRad(tp().theta0(matched));
+            scalarField theta(boundary[patchi].size(), theta0);
+
+            scalar uTheta = tp().uTheta();
+
+            // Calculate the dynamic contact angle if required
+            if (uTheta > small)
+            {
+                scalar thetaA = degToRad(tp().thetaA(matched));
+                scalar thetaR = degToRad(tp().thetaR(matched));
+
+                // Calculated the component of the velocity parallel to the wall
+                vectorField Uwall
+                (
+                    pair.phase1().U()().boundaryField()[patchi].
+                    patchInternalField()
+                  - pair.phase1().U()().boundaryField()[patchi]
+                );
+                Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
+
+                // Find the direction of the interface parallel to the wall
+                vectorField nWall
+                (
+                    nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
+                );
+
+                // Normalise nWall
+                nWall /= (mag(nWall) + small);
+
+                // Calculate Uwall resolved normal to the interface parallel to
+                // the interface
+                scalarField uwall(nWall & Uwall);
+
+                theta += (thetaA - thetaR)*tanh(uwall/uTheta);
+            }
+
+
+            // Reset nHatPatch to correspond to the contact angle
+
+            scalarField a12(nHatPatch & AfHatPatch);
+
+            scalarField b1(cos(theta));
+
+            scalarField b2(nHatPatch.size());
+
+            forAll(b2, facei)
+            {
+                b2[facei] = cos(acos(a12[facei]) - theta[facei]);
+            }
+
+            scalarField det(1 - a12*a12);
+
+            scalarField a((b1 - a12*b2)/det);
+            scalarField b((b2 - a12*b1)/det);
+
+            nHatPatch = a*AfHatPatch + b*nHatPatch;
+
+            nHatPatch /= (mag(nHatPatch) + deltaN_.value());
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceTensionForceModels::continuumSurfaceForce::continuumSurfaceForce
+(
+    const phaseSystem& fluid
+)
+:
+    surfaceTensionForceModel(fluid),
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(this->fluid_.mesh().V()), 1.0/3.0)
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfaceTensionForceModels::continuumSurfaceForce::~continuumSurfaceForce()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::surfaceTensionForceModels::continuumSurfaceForce::K
+(
+    const phasePair& pair
+) const
+{
+    tmp<surfaceVectorField> tnHatfv = nHatfv(pair);
+
+    correctContactAngle(tnHatfv.ref().boundaryFieldRef(), pair);
+
+    return -fvc::div(tnHatfv & this->fluid_.mesh().Sf());
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/continuumSurfaceForce/continuumSurfaceForce.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/continuumSurfaceForce/continuumSurfaceForce.H
new file mode 100644
index 000000000..cb8a95b67
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/continuumSurfaceForce/continuumSurfaceForce.H
@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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::surfaceTensionForceModels::continuumSurfaceForce
+
+Description
+    Continuum surface force approach according to Brackbill et al. (1992)
+
+    Reference:
+    \verbatim
+        Brackbill, J. U., Kothe, D. B., Zemach, C. (1992).
+        A continuum method for modeling surface tension.
+        Journal of Computational Physics, 100, 335-354.
+    \endverbatim
+
+SourceFiles
+    continuumSurfaceForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef continuumSurfaceForce_H
+#define continuumSurfaceForce_H
+
+#include "surfaceTensionForceModel.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class phasePair;
+
+namespace surfaceTensionForceModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class continuumSurfaceForce Declaration
+\*---------------------------------------------------------------------------*/
+
+class continuumSurfaceForce
+:
+    public surfaceTensionForceModel
+{
+    // Private Data
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+
+    // Private Member Functions
+
+        //- Return the face-based interface normal vector
+        tmp<surfaceVectorField> nHatfv(const phasePair& pair) const;
+
+        //- Return the face-projection of interface normal vector
+        tmp<surfaceScalarField> nHatf(const phasePair& pair) const;
+
+        void correctContactAngle
+        (
+            surfaceVectorField::Boundary& nHatb,
+            const phasePair& pair
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("continuumSurfaceForce");
+
+
+    // Constructors
+
+        //- Construct from a dictionary and a phase pair
+        continuumSurfaceForce(const phaseSystem& fluid);
+
+
+    //- Destructor
+    virtual ~continuumSurfaceForce();
+
+
+    // Member Functions
+
+        //- Return the curvature of this pair
+        virtual tmp<volScalarField> K(const phasePair& pair) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceTensionForceModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/noSurfaceTensionForce/noSurfaceTensionForce.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/noSurfaceTensionForce/noSurfaceTensionForce.C
new file mode 100644
index 000000000..a457afe79
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/noSurfaceTensionForce/noSurfaceTensionForce.C
@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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 "noSurfaceTensionForce.H"
+#include "phasePair.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceTensionForceModels
+{
+    defineTypeNameAndDebug(noSurfaceTensionForce, 0);
+    addToRunTimeSelectionTable
+    (
+        surfaceTensionForceModel,
+        noSurfaceTensionForce,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceTensionForceModels::noSurfaceTensionForce::noSurfaceTensionForce
+(
+    const phaseSystem& fluid
+)
+:
+    surfaceTensionForceModel(fluid)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfaceTensionForceModels::noSurfaceTensionForce::~noSurfaceTensionForce()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::surfaceTensionForceModels::noSurfaceTensionForce::K
+(
+    const phasePair& pair
+) const
+{
+    NotImplemented;
+    return volScalarField::null();
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::surfaceTensionForceModels::noSurfaceTensionForce::surfaceTensionForcef
+(
+    const phasePair& pair
+) const
+{
+    return
+        tmp<surfaceScalarField>
+        (
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    "surfaceTensionForce",
+                    this->fluid_.mesh().time().timeName(),
+                    this->fluid_.mesh()
+                ),
+                this->fluid_.mesh(),
+                dimensionedScalar(dimensionSet(1, -2, -2, 0, 0), Zero)
+            )
+        );
+}
+
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/noSurfaceTensionForce/noSurfaceTensionForce.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/noSurfaceTensionForce/noSurfaceTensionForce.H
new file mode 100644
index 000000000..f5fa7f692
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/noSurfaceTensionForce/noSurfaceTensionForce.H
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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::surfaceTensionForceModels::noSurfaceTensionForce
+
+Description
+    Deactivated surface tension force.
+
+SourceFiles
+    noSurfaceTensionForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noSurfaceTensionForce_H
+#define noSurfaceTensionForce_H
+
+#include "surfaceTensionForceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class phasePair;
+
+namespace surfaceTensionForceModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class noSurfaceTensionForce Declarationc
+\*---------------------------------------------------------------------------*/
+
+class noSurfaceTensionForce
+:
+    public surfaceTensionForceModel
+{
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+        //- Construct from a dictionary and a phase pair
+        noSurfaceTensionForce(const phaseSystem& fluid);
+
+
+    //- Destructor
+    virtual ~noSurfaceTensionForce();
+
+
+    // Member Functions
+
+        //- Return the curvature of this pair
+        virtual tmp<volScalarField> K(const phasePair& pair) const;
+
+        //- Interface compression flux
+        virtual tmp<surfaceScalarField> surfaceTensionForcef
+        (
+            const phasePair& pair
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceTensionForceModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModel.C
new file mode 100644
index 000000000..64ae5eb80
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModel.C
@@ -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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceTensionForceModel.H"
+#include "phasePair.H"
+#include "phaseSystem.H"
+#include "fvcSnGrad.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(surfaceTensionForceModel, 0);
+    defineRunTimeSelectionTable(surfaceTensionForceModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::surfaceTensionForceModel::gradAlphaf(const phasePair& pair) const
+{
+    return
+        fvc::interpolate(pair.phase2())*fvc::snGrad(pair.phase1())
+      - fvc::interpolate(pair.phase1())*fvc::snGrad(pair.phase2());
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceTensionForceModel::surfaceTensionForceModel
+(
+    const phaseSystem& fluid
+)
+:
+    fluid_(fluid)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfaceTensionForceModel::~surfaceTensionForceModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::surfaceTensionForceModel::surfaceTensionForcef
+(
+    const phasePair& pair
+) const
+{
+    return
+        fvc::interpolate(fluid_.sigma(pair)*this->K(pair))
+       *this->gradAlphaf(pair);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModel.H
new file mode 100644
index 000000000..bc76b5c0b
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModel.H
@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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::surfaceTensionForceModel
+
+Description
+    Base class for modeling the surface tension force flux.
+
+SourceFiles
+    surfaceTensionForceModel.C
+    surfaceTensionForceModelNew.C
+    phaseSystem.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfaceTensionForceModel_H
+#define surfaceTensionForceModel_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "volFields.H"
+#include "dictionary.H"
+#include "runTimeSelectionTables.H"
+#include "phaseSystem.H"
+
+namespace Foam
+{
+
+class phasePair;
+class phaseModel;
+
+/*---------------------------------------------------------------------------*\
+                  Class surfaceTensionForceModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class surfaceTensionForceModel
+{
+protected:
+
+    // Protected data
+
+        //- Reference to the mesh
+        const phaseSystem& fluid_;
+
+
+    // Protected Member Functions
+
+        //- Return the face-based interface vector
+        tmp<surfaceScalarField> gradAlphaf(const phasePair& pair) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("surfaceTensionForceModel");
+
+
+    // Declare runtime construction
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            surfaceTensionForceModel,
+            dictionary,
+            (
+                const phaseSystem& fluid
+            ),
+            (fluid)
+        );
+
+
+    // Constructors
+
+        //- Construct from mesh
+        surfaceTensionForceModel(const phaseSystem& fluid);
+
+
+    //- Destructor
+    virtual ~surfaceTensionForceModel();
+
+
+    // Selectors
+
+        static autoPtr<surfaceTensionForceModel> New(const phaseSystem& fluid);
+
+
+    // Member Functions
+
+        //- Return the curvature of this pair
+        virtual tmp<volScalarField> K(const phasePair& pair) const = 0;
+
+        //- Surface tension force flux
+        virtual tmp<surfaceScalarField> surfaceTensionForcef
+        (
+            const phasePair& pair
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModelNew.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModelNew.C
new file mode 100644
index 000000000..872c2049b
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModelNew.C
@@ -0,0 +1,78 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 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 "continuumSurfaceForce.H"
+
+// * * * * * * * * * * * * * * * * Selector  * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::surfaceTensionForceModel>
+Foam::surfaceTensionForceModel::New
+(
+    const phaseSystem& fluid
+)
+{
+    // Select the base method if method information is missing.
+    const word methodDictName = Foam::surfaceTensionForceModel::typeName;
+    const word typeKeyword = "type";
+    word selectedType;
+
+    if (fluid.found(methodDictName))
+    {
+        const dictionary& methodDict = fluid.subDict(methodDictName);
+
+        if (methodDict.found(typeKeyword))
+        {
+            selectedType = methodDict.lookup<word>(typeKeyword);
+            Info<< "Selecting surfaceTensionForceModel: "
+                << selectedType<< endl;
+
+            dictionaryConstructorTable::iterator cstrIter =
+                dictionaryConstructorTablePtr_->find(selectedType);
+
+            if (cstrIter == dictionaryConstructorTablePtr_->end())
+            {
+                FatalErrorInFunction
+                    << "Unknown " << methodDictName << " type "
+                    << selectedType << endl << endl
+                    << "Valid " << methodDictName << " types are : " << endl
+                    << dictionaryConstructorTablePtr_->sortedToc() << endl
+                    << exit(FatalError);
+            }
+
+            return cstrIter()(fluid);
+        }
+    }
+
+    Info<< "Selecting surfaceTensionForceModel: continuumSurfaceForce" << endl;
+
+    return
+        autoPtr<surfaceTensionForceModel>
+        (
+            new surfaceTensionForceModels::continuumSurfaceForce(fluid)
+        );
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
index 0f4eab24f..2e326073d 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
@@ -468,6 +468,14 @@ Foam::MovingPhaseModel<BasePhaseModel>::pPrime() const
 }
 
 
+template<class BasePhaseModel>
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::MovingPhaseModel<BasePhaseModel>::divDevTau(volVectorField& U) const
+{
+    return turbulence_->divDevTau(U);
+}
+
+
 template<class BasePhaseModel>
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::MovingPhaseModel<BasePhaseModel>::divq(volScalarField& he) const
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
index 452e97334..a716e4034 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
@@ -219,6 +219,9 @@ public:
             //  (derivative of phase-pressure w.r.t. phase-fraction)
             virtual tmp<volScalarField> pPrime() const;
 
+            //- Return the turbulent term for the momentum equation
+            virtual tmp<fvVectorMatrix> divDevTau(volVectorField& U) const;
+
 
         // Thermophysical transport
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C
index 9a8ab37d3..837ea0d5e 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C
@@ -286,6 +286,21 @@ Foam::StationaryPhaseModel<BasePhaseModel>::pPrime() const
 }
 
 
+template<class BasePhaseModel>
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::StationaryPhaseModel<BasePhaseModel>::divDevTau
+(
+    volVectorField& U
+) const
+{
+    FatalErrorInFunction
+        << "Cannot access the turbulent term of a stationary phase"
+        << exit(FatalError);
+
+    return fvm::Su(U, U);
+}
+
+
 template<class BasePhaseModel>
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::StationaryPhaseModel<BasePhaseModel>::divq(volScalarField& he) const
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H
index f30bc27cd..d97080e6e 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H
@@ -146,6 +146,9 @@ public:
             //  (derivative of phase-pressure w.r.t. phase-fraction)
             virtual tmp<volScalarField> pPrime() const;
 
+            //- Return the turbulent term for the momentum equation
+            virtual tmp<fvVectorMatrix> divDevTau(volVectorField& U) const;
+
 
         // Thermophysical transport
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
index 40a3709dd..f3039f33b 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
@@ -85,7 +85,8 @@ 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)),
+    continuous_(fluid.subDict(phaseName).lookupOrDefault("continuous", false))
 {
     diameterModel_ = diameterModel::New(fluid.subDict(phaseName), *this);
 }
@@ -142,6 +143,12 @@ Foam::scalar Foam::phaseModel::alphaMax() const
 }
 
 
+bool Foam::phaseModel::continuous() const
+{
+    return continuous_;
+}
+
+
 Foam::tmp<Foam::volScalarField> Foam::phaseModel::d() const
 {
     return diameterModel_().d();
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
index 3bbbc3f1c..46963db11 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
@@ -75,6 +75,9 @@ class phaseModel
         //- Optional maximum phase-fraction (e.g. packing limit)
         scalar alphaMax_;
 
+        //- Morphology of phase (default is false, i.e. dispersed)
+        const Switch continuous_;
+
         //- Diameter model
         autoPtr<diameterModel> diameterModel_;
 
@@ -192,6 +195,9 @@ public:
         //- Return the maximum phase-fraction (e.g. packing limit)
         scalar alphaMax() const;
 
+        //- Return true if phase is continuous
+        bool continuous() const;
+
         //- Return the Sauter-mean diameter
         tmp<volScalarField> d() const;
 
@@ -349,6 +355,12 @@ public:
             //- Return the phase-pressure'
             //  (derivative of phase-pressure w.r.t. phase-fraction)
             virtual tmp<volScalarField> pPrime() const = 0;
+
+            //- Return the turbulent term for the momentum equation
+            virtual tmp<fvVectorMatrix> divDevTau
+            (
+                volVectorField& U
+            ) const = 0;
 };
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phasePair/phasePair/phasePair.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phasePair/phasePair/phasePair.C
index 3e4123a33..440d21d8c 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phasePair/phasePair/phasePair.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phasePair/phasePair/phasePair.C
@@ -104,7 +104,24 @@ Foam::word Foam::phasePair::otherName() const
 
 Foam::tmp<Foam::volScalarField> Foam::phasePair::rho() const
 {
-    return phase1()*phase1().rho() + phase2()*phase2().rho();
+    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
+        );
 }
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
index ef75b4a9c..09477b88d 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
@@ -194,164 +194,6 @@ void Foam::phaseSystem::setMixturePhi
 }
 
 
-Foam::tmp<Foam::surfaceVectorField> Foam::phaseSystem::nHatfv
-(
-    const volScalarField& alpha1,
-    const volScalarField& alpha2
-) const
-{
-    /*
-    // Cell gradient of alpha
-    volVectorField gradAlpha =
-        alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
-
-    // Interpolated face-gradient of alpha
-    surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
-    */
-
-    surfaceVectorField gradAlphaf
-    (
-        fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1))
-      - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
-    );
-
-    // Face unit interface normal
-    return gradAlphaf/(mag(gradAlphaf) + deltaN_);
-}
-
-
-Foam::tmp<Foam::surfaceScalarField> Foam::phaseSystem::nHatf
-(
-    const volScalarField& alpha1,
-    const volScalarField& alpha2
-) const
-{
-    // Face unit interface normal flux
-    return nHatfv(alpha1, alpha2) & mesh_.Sf();
-}
-
-
-void Foam::phaseSystem::correctContactAngle
-(
-    const phaseModel& phase1,
-    const phaseModel& phase2,
-    surfaceVectorField::Boundary& nHatb
-) const
-{
-    const volScalarField::Boundary& gbf
-        = phase1.boundaryField();
-
-    const fvBoundaryMesh& boundary = mesh_.boundary();
-
-    forAll(boundary, patchi)
-    {
-        if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
-        {
-            const alphaContactAngleFvPatchScalarField& acap =
-                refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
-
-            vectorField& nHatPatch = nHatb[patchi];
-
-            vectorField AfHatPatch
-            (
-                mesh_.Sf().boundaryField()[patchi]
-               /mesh_.magSf().boundaryField()[patchi]
-            );
-
-            alphaContactAngleFvPatchScalarField::thetaPropsTable::
-                const_iterator tp =
-                    acap.thetaProps()
-                   .find(phasePairKey(phase1.name(), phase2.name()));
-
-            if (tp == acap.thetaProps().end())
-            {
-                FatalErrorInFunction
-                    << "Cannot find interface "
-                    << phasePairKey(phase1.name(), phase2.name())
-                    << "\n    in table of theta properties for patch "
-                    << acap.patch().name()
-                    << exit(FatalError);
-            }
-
-            bool matched = (tp.key().first() == phase1.name());
-
-            scalar theta0 = degToRad(tp().theta0(matched));
-            scalarField theta(boundary[patchi].size(), theta0);
-
-            scalar uTheta = tp().uTheta();
-
-            // Calculate the dynamic contact angle if required
-            if (uTheta > small)
-            {
-                scalar thetaA = degToRad(tp().thetaA(matched));
-                scalar thetaR = degToRad(tp().thetaR(matched));
-
-                // Calculated the component of the velocity parallel to the wall
-                vectorField Uwall
-                (
-                    phase1.U()().boundaryField()[patchi].patchInternalField()
-                  - phase1.U()().boundaryField()[patchi]
-                );
-                Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
-
-                // Find the direction of the interface parallel to the wall
-                vectorField nWall
-                (
-                    nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
-                );
-
-                // Normalise nWall
-                nWall /= (mag(nWall) + small);
-
-                // Calculate Uwall resolved normal to the interface parallel to
-                // the interface
-                scalarField uwall(nWall & Uwall);
-
-                theta += (thetaA - thetaR)*tanh(uwall/uTheta);
-            }
-
-
-            // Reset nHatPatch to correspond to the contact angle
-
-            scalarField a12(nHatPatch & AfHatPatch);
-
-            scalarField b1(cos(theta));
-
-            scalarField b2(nHatPatch.size());
-
-            forAll(b2, facei)
-            {
-                b2[facei] = cos(acos(a12[facei]) - theta[facei]);
-            }
-
-            scalarField det(1 - a12*a12);
-
-            scalarField a((b1 - a12*b2)/det);
-            scalarField b((b2 - a12*b1)/det);
-
-            nHatPatch = a*AfHatPatch + b*nHatPatch;
-
-            nHatPatch /= (mag(nHatPatch) + deltaN_.value());
-        }
-    }
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::phaseSystem::K
-(
-    const phaseModel& phase1,
-    const phaseModel& phase2
-) const
-{
-    tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
-
-    correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
-
-    // Simple expression for curvature
-    return -fvc::div(tnHatfv & mesh_.Sf());
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::phaseSystem::phaseSystem
@@ -402,21 +244,14 @@ Foam::phaseSystem::phaseSystem
         dimensionedScalar(dimPressure/dimTime, 0)
     ),
 
-    MRF_(mesh_),
-
-    cAlphas_(lookupOrDefault("interfaceCompression", cAlphaTable())),
-
-    deltaN_
-    (
-        "deltaN",
-        1e-8/pow(average(mesh_.V()), 1.0/3.0)
-    )
+    MRF_(mesh_)
 {
     // Groupings
     label movingPhasei = 0;
     label stationaryPhasei = 0;
     label anisothermalPhasei = 0;
     label multiComponentPhasei = 0;
+    label continuousPhasei = 0;
     forAll(phaseModels_, phasei)
     {
         phaseModel& phase = phaseModels_[phasei];
@@ -424,16 +259,19 @@ Foam::phaseSystem::phaseSystem
         stationaryPhasei += phase.stationary();
         anisothermalPhasei += !phase.isothermal();
         multiComponentPhasei += !phase.pure();
+        continuousPhasei += phase.continuous();
     }
     movingPhaseModels_.resize(movingPhasei);
     stationaryPhaseModels_.resize(stationaryPhasei);
     anisothermalPhaseModels_.resize(anisothermalPhasei);
     multiComponentPhaseModels_.resize(multiComponentPhasei);
+    continuousPhaseModels_.resize(continuousPhasei);
 
     movingPhasei = 0;
     stationaryPhasei = 0;
     anisothermalPhasei = 0;
     multiComponentPhasei = 0;
+    continuousPhasei = 0;
     forAll(phaseModels_, phasei)
     {
         phaseModel& phase = phaseModels_[phasei];
@@ -453,6 +291,20 @@ Foam::phaseSystem::phaseSystem
         {
             multiComponentPhaseModels_.set(multiComponentPhasei++, &phase);
         }
+        if (phase.continuous())
+        {
+            continuousPhaseModels_.set(continuousPhasei++, &phase);
+        }
+    }
+
+    if (!continuousPhaseModels_.empty())
+    {
+        Info<< "Continuous phases:";
+        forAll(continuousPhaseModels_, phasei)
+        {
+            Info<< " " << continuousPhaseModels_[phasei].name();
+        }
+        Info<< endl;
     }
 
     // Write phi
@@ -700,47 +552,6 @@ bool Foam::phaseSystem::implicitPhasePressure() const
 }
 
 
-Foam::tmp<Foam::surfaceScalarField> Foam::phaseSystem::surfaceTension
-(
-    const phaseModel& phase1
-) const
-{
-    tmp<surfaceScalarField> tSurfaceTension
-    (
-        surfaceScalarField::New
-        (
-            "surfaceTension",
-            mesh_,
-            dimensionedScalar(dimensionSet(1, -2, -2, 0, 0), 0)
-        )
-    );
-
-    forAll(phases(), phasej)
-    {
-        const phaseModel& phase2 = phases()[phasej];
-
-        if (&phase2 != &phase1)
-        {
-            phasePairKey key12(phase1.name(), phase2.name());
-
-            cAlphaTable::const_iterator cAlpha(cAlphas_.find(key12));
-
-            if (cAlpha != cAlphas_.end())
-            {
-                tSurfaceTension.ref() +=
-                    fvc::interpolate(sigma(key12)*K(phase1, phase2))
-                   *(
-                        fvc::interpolate(phase2)*fvc::snGrad(phase1)
-                      - fvc::interpolate(phase1)*fvc::snGrad(phase2)
-                    );
-            }
-        }
-    }
-
-    return tSurfaceTension;
-}
-
-
 void Foam::phaseSystem::correct()
 {
     forAll(phaseModels_, phasei)
@@ -897,4 +708,80 @@ Foam::tmp<Foam::surfaceScalarField> Foam::byDt(const surfaceScalarField& sf)
 }
 
 
+Foam::tmp<Foam::surfaceScalarField>
+Foam::phaseSystem::phir
+(
+    const phaseModel& phase1,
+    const phaseModel& phase2
+) const
+{
+    return phase1.phi() - phase2.phi();
+}
+
+
+void Foam::phaseSystem::addCompressionFlux
+(
+    const phaseModel& phase,
+    surfaceScalarField& alphaPhi
+) const
+{}
+
+
+Foam::tmp<Foam::volScalarField> Foam::phaseSystem::rhoContinuous() const
+{
+    const label nPhases = movingPhaseModels_.size();
+
+    tmp<volScalarField> rhoCont
+    (
+        movingPhaseModels_[0].continuous()
+       *movingPhaseModels_[0]
+       *movingPhaseModels_[0].rho()
+    );
+
+    for (label phasei = 1; phasei < nPhases; ++phasei)
+    {
+        rhoCont.ref() +=
+            movingPhaseModels_[phasei].continuous()
+           *movingPhaseModels_[phasei]
+           *movingPhaseModels_[phasei].rho();
+    }
+
+    return rhoCont/max(sumAlphaContinuous(), small);
+}
+
+
+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;
+}
+
+
+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();
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
index 7d98e03f7..6067c1eeb 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
@@ -100,6 +100,68 @@ public:
             >
             dmidtfTable;
 
+        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
+            <
+                surfaceScalarField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            FlfTable;
+
+        typedef
+            HashPtrTable
+            <
+                surfaceScalarField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            FwfTable;
+
+        typedef
+            HashPtrTable
+            <
+                surfaceScalarField,
+                phasePairKey,
+                phasePairKey::hash
+            >
+            FtfTable;
 
 
 protected:
@@ -132,9 +194,6 @@ protected:
             >
             aspectRatioModelTable;
 
-        typedef HashTable<scalar, phasePairKey, phasePairKey::hash>
-            cAlphaTable;
-
 
     // Protected data
 
@@ -160,6 +219,9 @@ protected:
         //- Multi-component phase models
         phaseModelPartialList multiComponentPhaseModels_;
 
+        //- Continuous phase models
+        phaseModelPartialList continuousPhaseModels_;
+
         //- Phase pairs
         phasePairTable phasePairs_;
 
@@ -175,12 +237,6 @@ protected:
         //- Blending methods
         blendingMethodTable blendingMethods_;
 
-        //- Interface compression coefficients
-        cAlphaTable cAlphas_;
-
-        //- Stabilisation for normalisation of the interface normal
-        const dimensionedScalar deltaN_;
-
 
         // Sub Models
 
@@ -279,47 +335,6 @@ protected:
         );
 
 
-        // Functions required for interface compression
-
-        //- Normal to interface between two phases
-        //  Used for interface compression
-        tmp<surfaceVectorField> nHatfv
-        (
-            const volScalarField& alpha1,
-            const volScalarField& alpha2
-        ) const;
-
-        //- Normal to interface between two phases dotted with face areas
-        //  Used for interface compression
-        tmp<surfaceScalarField> nHatf
-        (
-            const volScalarField& alpha1,
-            const volScalarField& alpha2
-        ) const;
-
-        //- Correction for the boundary condition on the unit normal nHat on
-        //  walls to produce the correct contact angle.
-        //
-        //  The dynamic contact angle is calculated from the component of the
-        //  velocity on the direction of the interface, parallel to the wall.
-        //  Used for interface compression
-        void correctContactAngle
-        (
-            const phaseModel& alpha1,
-            const phaseModel& alpha2,
-            surfaceVectorField::Boundary& nHatb
-        ) const;
-
-
-        //- Curvature of interface between two phases
-        //  Used for interface compression
-        tmp<volScalarField> K
-        (
-            const phaseModel& alpha1,
-            const phaseModel& alpha2
-        ) const;
-
-
         // Functions required by twoPhaseSystem
 
             //- Return the drag coefficient for phase pair
@@ -413,6 +428,9 @@ public:
             //- Access the models for phases that have multiple species
             inline phaseModelPartialList& multiComponentPhases();
 
+            //- Return the models for phases that are continuous
+            inline const phaseModelPartialList& continuousPhases() const;
+
             //- Return the phase pairs
             inline const phasePairTable& phasePairs() const;
 
@@ -633,10 +651,11 @@ public:
             virtual autoPtr<specieTransferTable> specieTransfer() const = 0;
 
             //- Return the surface tension force
-            tmp<surfaceScalarField> surfaceTension
+            //  Remark: moved to ResolvedPhaseSystem
+            virtual tmp<surfaceScalarField> surfaceTension
             (
                 const phaseModel& phase
-            ) const;
+            ) const = 0;
 
 
         // Evolution
@@ -677,6 +696,62 @@ public:
 
             //- Read base phaseProperties dictionary
             virtual bool read();
+
+
+        // HZDR
+
+            //- Return relative flux
+            virtual tmp<surfaceScalarField> phir
+            (
+                const phaseModel& phase1,
+                const phaseModel& phase2
+            ) const;
+
+            //- Add interface compression flux
+            virtual void addCompressionFlux
+            (
+                const phaseModel& phase,
+                surfaceScalarField& alphaPhi
+            ) const;
+
+            //- Return the mixture density based on continuous phases only
+            tmp<volScalarField> rhoContinuous() const;
+
+            //- Return the total volume fraction of continuous phases
+            tmp<volScalarField> sumAlphaContinuous() const;
+
+            //- Indicator field of the proximity of a resolved phase interface
+            //  Values are distributed from 0 (far) to 1 (near).
+            //  Remark: implemented in ResolvedPhaseSystem
+            virtual tmp<volScalarField> interfaceIndicator
+            (
+                const phaseModel& phase
+            ) const = 0;
+
+            //- 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 face lift force for all phase-pairs
+            virtual const FlfTable& Flfs() const = 0;
+
+            //- Return the face wall lubrication force for all phase-pairs
+            virtual const FwfTable& Fwfs() const = 0;
+
+            //- Return the face turbulent dispersion force for all phase-pairs
+            virtual const FtfTable& Ftfs() const = 0;
+
+            //- Return the interfacial mass transfer rate for all phase-pairs
+            virtual const dmdtfTable& dmdtfsPhaseTransfer() const = 0;
 };
 
 
@@ -684,6 +759,7 @@ public:
 
 tmp<volScalarField> byDt(const volScalarField& vf);
 tmp<surfaceScalarField> byDt(const surfaceScalarField& sf);
+tmp<volVectorField> byDt(const volVectorField& vf);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
index 234f431c8..54e2806db 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
@@ -101,6 +101,13 @@ Foam::phaseSystem::multiComponentPhases()
 }
 
 
+inline const Foam::phaseSystem::phaseModelPartialList&
+Foam::phaseSystem::continuousPhases() const
+{
+    return continuousPhaseModels_;
+}
+
+
 inline const Foam::phaseSystem::phasePairTable&
 Foam::phaseSystem::phasePairs() const
 {
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemSolve.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemSolve.C
index 39bd635e1..941472ec4 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemSolve.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemSolve.C
@@ -300,27 +300,11 @@ void Foam::phaseSystem::solve
 
                         if (&phase2 == &phase) continue;
 
-                        surfaceScalarField phir(phase.phi() - phase2.phi());
-
-                        cAlphaTable::const_iterator cAlpha
-                        (
-                            cAlphas_.find
-                            (
-                                phasePairKey(phase.name(), phase2.name())
-                            )
-                        );
-
-                        if (cAlpha != cAlphas_.end())
-                        {
-                            surfaceScalarField phic
-                            (
-                                (mag(phi_) + mag(phir))/mesh_.magSf()
-                            );
-
-                            phir +=
-                                min(cAlpha()*phic, max(phic))
-                               *nHatf(alpha, alpha2);
-                        }
+                        // Calculation of relative flux outsourced into new
+                        // virtual function such that additional fluxes
+                        // (e.g. interface compression) can be added in separate
+                        // phaseSystems (e.g. ResolvedPhaseSytem).
+                        surfaceScalarField phir(this->phir(phase, phase2));
 
                         const word phirScheme
                         (
@@ -337,56 +321,9 @@ void Foam::phaseSystem::solve
                         );
                     }
                 }
-                else if (!cAlphas_.empty())
+                else
                 {
-                    forAll(phases(), phasei)
-                    {
-                        const phaseModel& phase2 = phases()[phasei];
-                        const volScalarField& alpha2 = phase2;
-
-                        if (&phase2 == &phase) continue;
-
-                        cAlphaTable::const_iterator cAlpha
-                        (
-                            cAlphas_.find
-                            (
-                                phasePairKey(phase.name(), phase2.name())
-                            )
-                        );
-
-                        if (cAlpha != cAlphas_.end())
-                        {
-                            const surfaceScalarField phir
-                            (
-                                phase.phi() - phase2.phi()
-                            );
-
-                            const surfaceScalarField phic
-                            (
-                                (mag(phi_) + mag(phir))/mesh_.magSf()
-                            );
-
-                            const surfaceScalarField phirc
-                            (
-                                min(cAlpha()*phic, max(phic))
-                               *nHatf(alpha, alpha2)
-                            );
-
-                            const word phirScheme
-                            (
-                                "div(phir,"
-                              + alpha2.name() + ',' + alpha.name()
-                              + ')'
-                            );
-
-                            alphaPhi += fvc::flux
-                            (
-                                -fvc::flux(-phirc, alpha2, phirScheme),
-                                alpha,
-                                phirScheme
-                            );
-                        }
-                    }
+                    this->addCompressionFlux(phase, alphaPhi);
                 }
 
                 if (alphaPhiDbyA0s.set(phase.index()))
@@ -602,6 +539,19 @@ void Foam::phaseSystem::solve
             }
         }
     }
+
+    if (!nAlphaCorr)
+    {
+        forAll(phases(), phasei)
+        {
+            phaseModel& phase = phases()[phasei];
+
+            phase.alphaPhiRef() = fvc::interpolate(phase)*phase.phi();
+
+            phase.alphaRhoPhiRef() =
+                fvc::interpolate(phase.rho())*phase.alphaPhi();
+        }
+    }
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/include/createPIDTimeControls.H b/src/finiteVolume/cfdTools/general/include/createPIDTimeControls.H
new file mode 100644
index 000000000..4e338cec1
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/include/createPIDTimeControls.H
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 2019 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/>.
+
+Global
+    readPIDTimeControls
+
+Description
+    Read the control parameters used by setPIDDeltaT
+
+\*---------------------------------------------------------------------------*/
+
+const dictionary& cDict(runTime.controlDict());
+bool adjustTimeStep(cDict.lookupOrDefault("adjustTimeStep", false));
+
+// Default values
+const scalar defaultTargetCo(great);
+const scalar defaultKP(-1e-2);
+const scalar defaultKI(-2e-2);
+const scalar defaultKD(-1e-6);
+
+const scalar defaultMaxCo(1.0);
+const scalar defaultMinDeltaT(small);
+const scalar defaultMaxDeltaT(great);
+
+// PID control variable and coefficients
+scalar targetCo(cDict.lookupOrDefault<scalar>("targetCo", defaultTargetCo));
+scalar kP(defaultKP);
+scalar kI(defaultKI);
+scalar kD(defaultKD);
+
+// Limiters
+scalar maxCo(cDict.lookupOrDefault<scalar>("maxCo", defaultMaxCo));
+scalar minDeltaT(cDict.lookupOrDefault<scalar>("minDeltaT", defaultMinDeltaT));
+scalar maxDeltaT(cDict.lookupOrDefault<scalar>("maxDeltaT", defaultMaxDeltaT));
+
+// Recent PID values
+scalar deltaCoSum(cDict.lookupOrDefault<scalar>("deltaCoSum", 0));
+scalar deltaCo0(cDict.lookupOrDefault<scalar>("deltaCo0", 0));
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/include/readPIDTimeControls.H b/src/finiteVolume/cfdTools/general/include/readPIDTimeControls.H
new file mode 100644
index 000000000..78c84617b
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/include/readPIDTimeControls.H
@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 2019 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/>.
+
+Global
+    readPIDTimeControls
+
+Description
+    Read the control parameters used by setPIDDeltaT
+
+\*---------------------------------------------------------------------------*/
+
+adjustTimeStep = cDict.lookupOrDefault("adjustTimeStep", false);
+
+if (adjustTimeStep)
+{
+    // PID control variable and coefficients
+    targetCo = cDict.lookupOrDefault<scalar>("targetCo", defaultTargetCo);
+    kP = cDict.lookupOrDefault<scalar>("kP", defaultKP);
+    kI = cDict.lookupOrDefault<scalar>("kI", defaultKI);
+    kD = cDict.lookupOrDefault<scalar>("kD", defaultKD);
+
+    // Limiters
+    maxCo = cDict.lookupOrDefault<scalar>("maxCo", defaultMaxCo);
+    minDeltaT = cDict.lookupOrDefault<scalar>("minDeltaT", defaultMinDeltaT);
+    maxDeltaT = cDict.lookupOrDefault<scalar>("maxDeltaT", defaultMaxDeltaT);
+
+    if (runTime.timeIndex() == 0)
+    {
+        if (targetCo < maxCo)
+        {
+            Info<< "PID controlled time step" << nl
+                << "  targetCo        = " << targetCo << nl
+                << "  kP              = " << kP << nl
+                << "  kI              = " << kI << nl
+                << "  kD              = " << kD << nl
+                << "  maxCo           = " << maxCo << nl
+                << "  minDeltaT       = " << minDeltaT << nl
+                << "  maxDeltaT       = " << maxDeltaT << nl
+                << endl;
+
+            // Warnings
+            if
+            (
+                cDict.lookupOrDefault<word>("writeControl", "runTime")
+             == "adjustableRunTime"
+            )
+            {
+                Warning
+                    << "adjustableRunTime may compromise PID controller "
+                    << "performance for time step." << nl << endl;
+            }
+
+            if (kP > 0 || kI > 0 || kD > 0)
+            {
+                Warning
+                    << "PID coefficients for time step control should be "
+                    << "negative." << endl;
+            }
+        }
+        else
+        {
+            Info<< "Controlled time step" << nl
+                << "  maxCo           = " << maxCo << nl
+                << "  minDeltaT       = " << minDeltaT << nl
+                << "  maxDeltaT       = " << maxDeltaT << nl
+                << endl;
+        }
+    }
+}
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/include/setInitialPIDDeltaT.H b/src/finiteVolume/cfdTools/general/include/setInitialPIDDeltaT.H
new file mode 100644
index 000000000..0ff26cdb1
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/include/setInitialPIDDeltaT.H
@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 2019 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/>.
+
+Global
+    setInitialPIDDeltaT
+
+Description
+    Set the initial timestep corresponding to the original timestep adjustment
+    algorithm in setDeltaT but only if it would reduce the timestep.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    if ((runTime.timeIndex() == 0) && (CoNum > small))
+    {
+        runTime.setDeltaT
+        (
+            max
+            (
+                min
+                (
+                    min(targetCo, maxCo)*runTime.deltaTValue()
+                   /max(CoNum, small),
+                    min(runTime.deltaTValue(), maxDeltaT)
+                ),
+                minDeltaT
+            )
+        );
+    }
+}
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/include/setPIDDeltaT.H b/src/finiteVolume/cfdTools/general/include/setPIDDeltaT.H
new file mode 100644
index 000000000..c5862d8c5
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/include/setPIDDeltaT.H
@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+  ==  == ====== ====   ====    |
+                    \\     ||  | Unsupported Contributions for OpenFOAM
+  ======   //   ||  || ===//   |
+  ||  ||  //    ||  // || \\   | Copyright (C) 2019 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/>.
+
+Global
+    setPIDDeltaT
+
+Description
+    Set time step size with PID controller to maintain target Courant number
+    (targetCo). Reduction of time-step is immediate, but increase is limited to
+    a fixed rate of increase (factor 1.2) to avoid unstable oscillations.
+    Additionally time step size is bounded by the value from the original time
+    step adjustment algorithm, which is controlled with the maximum Courant
+    number (maxCo).
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    // Standard time step control
+    const scalar maxDeltaTFact(maxCo/(CoNum + small));
+    const scalar deltaTFact
+    (
+        min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2)
+    );
+    const scalar stdDeltaT(deltaTFact*runTime.deltaTValue());
+
+    // PID time step control
+    scalar PIDDeltaT(great);
+    if (targetCo < maxCo)
+    {
+        const scalar deltaCo(CoNum - targetCo);
+
+        const scalar deltaTP(kP*deltaCo);
+
+        deltaCoSum += kI*deltaCo*runTime.deltaTValue();
+        const scalar deltaTI(deltaCoSum);
+        Info<< "deltaCoSum = " << deltaCoSum << endl;
+
+        const scalar deltaTD(kD*(deltaCo - deltaCo0)/runTime.deltaTValue());
+        deltaCo0 = deltaCo;
+
+        PIDDeltaT = runTime.deltaTValue() + deltaTP + deltaTI + deltaTD;
+    }
+
+    // Limit between minDeltaT and maxDeltaT
+    runTime.setDeltaT
+    (
+        max
+        (
+            min
+            (
+                min(stdDeltaT, PIDDeltaT),
+                maxDeltaT
+            ),
+            minDeltaT
+        )
+    );
+
+    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
+}
+
+// ************************************************************************* //
-- 
2.29.2

