From 25c4f31bfdb96e08de7608a8f3b352364e105b0c Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Fri, 4 Oct 2013 08:58:41 +0100
Subject: [PATCH] compressibleMultiphaseInterFoam: hack implementation of
 compressible multiphaseInterFoam Needs to be consolidated with
 multiphaseInterFoam with thermal and compressibility effects made run-time
 selectable

---
 .../compressibleMultiphaseInterFoam/Allwclean |    8 +
 .../compressibleMultiphaseInterFoam/Allwmake  |    8 +
 .../Make/files                                |    3 +
 .../Make/options                              |   17 +
 .../compressibleMultiphaseInterFoam/TEqn.H    |   17 +
 .../compressibleMultiphaseInterFoam/UEqn.H    |   27 +
 .../alphaCourantNo.H                          |   57 +
 .../compressibleMultiphaseInterFoam.C         |  111 ++
 .../createFields.H                            |   71 ++
 .../multiphaseMixtureThermo/Make/files        |    5 +
 .../multiphaseMixtureThermo/Make/options      |    8 +
 .../alphaContactAngleFvPatchScalarField.C     |  146 +++
 .../alphaContactAngleFvPatchScalarField.H     |  215 ++++
 .../multiphaseMixtureThermo.C                 | 1129 +++++++++++++++++
 .../multiphaseMixtureThermo.H                 |  434 +++++++
 .../phaseModel/phaseModel.C                   |   95 ++
 .../phaseModel/phaseModel.H                   |  156 +++
 .../compressibleMultiphaseInterFoam/pEqn.H    |  142 +++
 .../laminar/damBreak4phase/0.org/T            |   52 +
 .../laminar/damBreak4phase/0.org/U            |   51 +
 .../laminar/damBreak4phase/0.org/alpha.air    |   79 ++
 .../damBreak4phase/0.org/alpha.mercury        |   49 +
 .../laminar/damBreak4phase/0.org/alpha.oil    |   49 +
 .../laminar/damBreak4phase/0.org/alpha.water  |   49 +
 .../laminar/damBreak4phase/0.org/alphas       |   47 +
 .../laminar/damBreak4phase/0.org/p            |   53 +
 .../laminar/damBreak4phase/0.org/p_rgh        |   59 +
 .../laminar/damBreak4phase/Allclean           |   11 +
 .../laminar/damBreak4phase/Allrun             |   17 +
 .../laminar/damBreak4phase/constant/g         |   22 +
 .../constant/polyMesh/blockMeshDict           |  108 ++
 .../damBreak4phase/constant/polyMesh/boundary |   53 +
 .../constant/thermophysicalProperties         |   33 +
 .../constant/thermophysicalProperties.air     |   49 +
 .../constant/thermophysicalProperties.mercury |   54 +
 .../constant/thermophysicalProperties.oil     |   54 +
 .../constant/thermophysicalProperties.water   |   54 +
 .../constant/turbulenceProperties             |   21 +
 .../laminar/damBreak4phase/system/controlDict |   56 +
 .../damBreak4phase/system/decomposeParDict    |   45 +
 .../laminar/damBreak4phase/system/fvSchemes   |   64 +
 .../laminar/damBreak4phase/system/fvSolution  |  120 ++
 .../damBreak4phase/system/setFieldsDict       |   65 +
 43 files changed, 3963 insertions(+)
 create mode 100755 applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean
 create mode 100755 applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/alphaCourantNo.H
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/files
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.H
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H
 create mode 100644 applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/T
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/U
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.air
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.mercury
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.oil
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.water
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alphas
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p_rgh
 create mode 100755 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allclean
 create mode 100755 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allrun
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/g
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/blockMeshDict
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/boundary
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.air
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.mercury
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.oil
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.water
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/turbulenceProperties
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/controlDict
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/decomposeParDict
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution
 create mode 100644 tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/setFieldsDict

diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean
new file mode 100755
index 00000000000..8098c051604
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wclean libso multiphaseMixtureThermo
+wclean
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake
new file mode 100755
index 00000000000..dd002ee0698
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake libso multiphaseMixtureThermo
+wmake
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files
new file mode 100644
index 00000000000..8e6a4cb0de4
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files
@@ -0,0 +1,3 @@
+compressibleMultiphaseInterFoam.C
+
+EXE = $(FOAM_APPBIN)/compressibleMultiphaseInterFoam
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
new file mode 100644
index 00000000000..0f47e6979ed
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
@@ -0,0 +1,17 @@
+EXE_INC = \
+    -ImultiphaseMixtureThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lmultiphaseMixtureThermo \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -linterfaceProperties \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lfiniteVolume
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H
new file mode 100644
index 00000000000..a7ee3a7bb4b
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H
@@ -0,0 +1,17 @@
+{
+    fvScalarMatrix TEqn
+    (
+        fvm::ddt(rho, T)
+      + fvm::div(multiphaseProperties.rhoPhi(), T)
+      - fvm::laplacian(multiphaseProperties.alphaEff(turbulence->mut()), T)
+      + (
+            fvc::div(fvc::absolute(phi, U), p)
+          + fvc::ddt(rho, K) + fvc::div(multiphaseProperties.rhoPhi(), K)
+        )*multiphaseProperties.rCv()
+    );
+
+    TEqn.relax();
+    TEqn.solve();
+
+    multiphaseProperties.correct();
+}
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H
new file mode 100644
index 00000000000..38cfebde6f3
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H
@@ -0,0 +1,27 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(multiphaseProperties.rhoPhi(), U)
+      + turbulence->divDevRhoReff(U)
+    );
+
+    UEqn.relax();
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    multiphaseProperties.surfaceTensionForce()
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                ) * mesh.magSf()
+            )
+        );
+
+        K = 0.5*magSqr(U);
+    }
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/alphaCourantNo.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/alphaCourantNo.H
new file mode 100644
index 00000000000..1fd68ffb790
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/alphaCourantNo.H
@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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 maxAlphaCo
+(
+    readScalar(runTime.controlDict().lookup("maxAlphaCo"))
+);
+
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    scalarField sumPhi
+    (
+        multiphaseProperties.nearInterface()().internalField()
+      * fvc::surfaceSum(mag(phi))().internalField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C
new file mode 100644
index 00000000000..09018b603f0
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C
@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+Application
+    compressibleMultiphaseInterFoam
+
+Description
+    Solver for n compressible, non-isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach.
+
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    Turbulence modelling is generic, i.e.  laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "multiphaseMixtureThermo.H"
+#include "turbulenceModel.H"
+#include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "readTimeControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            multiphaseProperties.solve();
+
+            solve(fvm::ddt(rho) + fvc::div(multiphaseProperties.rhoPhi()));
+
+            #include "UEqn.H"
+            #include "TEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        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/compressibleMultiphaseInterFoam/createFields.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
new file mode 100644
index 00000000000..5c7f723a875
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
@@ -0,0 +1,71 @@
+    Info<< "Reading field p_rgh\n" << endl;
+    volScalarField p_rgh
+    (
+        IOobject
+        (
+            "p_rgh",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "createPhi.H"
+
+    Info<< "Constructing multiphaseMixtureThermo\n" << endl;
+    multiphaseMixtureThermo multiphaseProperties(U, phi);
+
+    volScalarField& p = multiphaseProperties.p();
+    volScalarField& T = multiphaseProperties.T();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT
+        ),
+        multiphaseProperties.rho()
+    );
+
+    Info<< max(rho) << min(rho);
+
+    dimensionedScalar pMin(multiphaseProperties.lookup("pMin"));
+
+
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    // Construct compressible turbulence model
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            multiphaseProperties.rhoPhi(),
+            multiphaseProperties
+        )
+    );
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/files b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/files
new file mode 100644
index 00000000000..93e6eb9e3dd
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/files
@@ -0,0 +1,5 @@
+phaseModel/phaseModel.C
+alphaContactAngle/alphaContactAngleFvPatchScalarField.C
+multiphaseMixtureThermo.C
+
+LIB = $(FOAM_LIBBIN)/libmultiphaseMixtureThermo
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options
new file mode 100644
index 00000000000..eab8cce15dd
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options
@@ -0,0 +1,8 @@
+EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+LIB_LIBS = \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lfiniteVolume
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C
new file mode 100644
index 00000000000..6872ae0321e
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C
@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+alphaContactAngleFvPatchScalarField::interfaceThetaProps::interfaceThetaProps
+(
+    Istream& is
+)
+:
+    theta0_(readScalar(is)),
+    uTheta_(readScalar(is)),
+    thetaA_(readScalar(is)),
+    thetaR_(readScalar(is))
+{}
+
+
+Istream& operator>>
+(
+    Istream& is,
+    alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_;
+    return is;
+}
+
+
+Ostream& operator<<
+(
+    Ostream& os,
+    const alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    os  << tp.theta0_ << token::SPACE
+        << tp.uTheta_ << token::SPACE
+        << tp.thetaA_ << token::SPACE
+        << tp.thetaR_;
+
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(p, iF)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, p, iF, mapper),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    zeroGradientFvPatchScalarField(p, iF),
+    thetaProps_(dict.lookup("thetaProperties"))
+{
+    evaluate();
+}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, iF),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void alphaContactAngleFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("thetaProperties")
+        << thetaProps_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    alphaContactAngleFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.H
new file mode 100644
index 00000000000..53ce61d24f3
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.H
@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::alphaContactAngleFvPatchScalarField
+
+Description
+    Contact-angle boundary condition for multi-phase interface-capturing
+    simulations.  Used in conjuction with multiphaseMixture.
+
+SourceFiles
+    alphaContactAngleFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef alphaContactAngleFvPatchScalarField_H
+#define alphaContactAngleFvPatchScalarField_H
+
+#include "zeroGradientFvPatchFields.H"
+#include "multiphaseMixtureThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class alphaContactAngleFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class alphaContactAngleFvPatchScalarField
+:
+    public zeroGradientFvPatchScalarField
+{
+public:
+
+    class interfaceThetaProps
+    {
+        //- Equilibrium contact angle
+        scalar theta0_;
+
+        //- Dynamic contact angle velocity scale
+        scalar uTheta_;
+
+        //- Limiting advancing contact angle
+        scalar thetaA_;
+
+        //- Limiting receeding contact angle
+        scalar thetaR_;
+
+
+    public:
+
+        // Constructors
+            interfaceThetaProps()
+            {}
+
+            interfaceThetaProps(Istream&);
+
+
+        // Member functions
+
+            //- Return the equilibrium contact angle theta0
+            scalar theta0(bool matched=true) const
+            {
+                if (matched) return theta0_;
+                else return 180.0 - theta0_;
+            }
+
+            //- Return the dynamic contact angle velocity scale
+            scalar uTheta() const
+            {
+                return uTheta_;
+            }
+
+            //- Return the limiting advancing contact angle
+            scalar thetaA(bool matched=true) const
+            {
+                if (matched) return thetaA_;
+                else return 180.0 - thetaA_;
+            }
+
+            //- Return the limiting receeding contact angle
+            scalar thetaR(bool matched=true) const
+            {
+                if (matched) return thetaR_;
+                else return 180.0 - thetaR_;
+            }
+
+
+        // IO functions
+
+            friend Istream& operator>>(Istream&, interfaceThetaProps&);
+            friend Ostream& operator<<(Ostream&, const interfaceThetaProps&);
+    };
+
+    typedef HashTable
+    <
+        interfaceThetaProps,
+        multiphaseMixtureThermo::interfacePair,
+        multiphaseMixtureThermo::interfacePair::hash
+    > thetaPropsTable;
+
+
+private:
+
+    // Private data
+
+        thetaPropsTable thetaProps_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("alphaContactAngle");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given alphaContactAngleFvPatchScalarField
+        //  onto a new patch
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Return the contact angle properties
+        const thetaPropsTable& thetaProps() const
+        {
+            return thetaProps_;
+        }
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
new file mode 100644
index 00000000000..16395505643
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
@@ -0,0 +1,1129 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "multiphaseMixtureThermo.H"
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "Time.H"
+#include "subCycle.H"
+#include "MULES.H"
+#include "fvcDiv.H"
+#include "fvcGrad.H"
+#include "fvcSnGrad.H"
+#include "fvcFlux.H"
+#include "fvcMeshPhi.H"
+#include "surfaceInterpolate.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(multiphaseMixtureThermo, 0);
+}
+
+
+const Foam::scalar Foam::multiphaseMixtureThermo::convertToRad =
+    Foam::constant::mathematical::pi/180.0;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::multiphaseMixtureThermo::calcAlphas()
+{
+    scalar level = 0.0;
+    alphas_ == 0.0;
+
+    forAllIter(PtrDictionary<phaseModel>, phases_, phase)
+    {
+        alphas_ += level*phase();
+        level += 1.0;
+    }
+
+    alphas_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::multiphaseMixtureThermo::multiphaseMixtureThermo
+(
+    const volVectorField& U,
+    const surfaceScalarField& phi
+)
+:
+    psiThermo(U.mesh(), word::null),
+    phases_(lookup("phases"), phaseModel::iNew(p_, T_)),
+
+    mesh_(U.mesh()),
+    U_(U),
+    phi_(phi),
+
+    rhoPhi_
+    (
+        IOobject
+        (
+            "rho*phi",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("rho*phi", dimMass/dimTime, 0.0)
+    ),
+
+    alphas_
+    (
+        IOobject
+        (
+            "alphas",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("alphas", dimless, 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+
+    sigmas_(lookup("sigmas")),
+    dimSigma_(1, 0, -2, 0, 0),
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(mesh_.V()), 1.0/3.0)
+    )
+{
+    calcAlphas();
+    alphas_.write();
+    correct();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::multiphaseMixtureThermo::correct()
+{
+    forAllIter(PtrDictionary<phaseModel>, phases_, phasei)
+    {
+        phasei().correct();
+    }
+
+    PtrDictionary<phaseModel>::iterator phasei = phases_.begin();
+
+    psi_ = phasei()*phasei().thermo().psi();
+    mu_ = phasei()*phasei().thermo().mu();
+    alpha_ = phasei()*phasei().thermo().alpha();
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        psi_ += phasei()*phasei().thermo().psi();
+        mu_ += phasei()*phasei().thermo().mu();
+        alpha_ += phasei()*phasei().thermo().alpha();
+    }
+}
+
+
+void Foam::multiphaseMixtureThermo::correctRho(const volScalarField& dp)
+{
+    forAllIter(PtrDictionary<phaseModel>, phases_, phasei)
+    {
+        phasei().thermo().rho() +=  phasei().thermo().psi()*dp;
+    }
+}
+
+
+bool Foam::multiphaseMixtureThermo::incompressible() const
+{
+    bool ico = true;
+
+    forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
+    {
+        ico &= phase().thermo().incompressible();
+    }
+
+    return ico;
+}
+
+
+bool Foam::multiphaseMixtureThermo::isochoric() const
+{
+    bool iso = true;
+
+    forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
+    {
+        iso &= phase().thermo().incompressible();
+    }
+
+    return iso;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::he
+(
+    const volScalarField& p,
+    const volScalarField& T
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<volScalarField> the(phasei()*phasei().thermo().he(p, T));
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        the() += phasei()*phasei().thermo().he(p, T);
+    }
+
+    return the;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::he
+(
+    const scalarField& p,
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<scalarField> the
+    (
+        scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells)
+    );
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        the() += scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells);
+    }
+
+    return the;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::he
+(
+    const scalarField& p,
+    const scalarField& T,
+    const label patchi
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<scalarField> the
+    (
+        phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi)
+    );
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        the() +=
+            phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi);
+    }
+
+    return the;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::hc() const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<volScalarField> thc(phasei()*phasei().thermo().hc());
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        thc() += phasei()*phasei().thermo().hc();
+    }
+
+    return thc;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::THE
+(
+    const scalarField& h,
+    const scalarField& p,
+    const scalarField& T0,
+    const labelList& cells
+) const
+{
+    notImplemented("multiphaseMixtureThermo::THE(...)");
+    return T0;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::THE
+(
+    const scalarField& h,
+    const scalarField& p,
+    const scalarField& T0,
+    const label patchi
+) const
+{
+    notImplemented("multiphaseMixtureThermo::THE(...)");
+    return T0;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::rho() const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<volScalarField> trho(phasei()*phasei().thermo().rho());
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        trho() += phasei()*phasei().thermo().rho();
+    }
+
+    return trho;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cp() const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<volScalarField> tCp(phasei()*phasei().thermo().Cp());
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tCp() += phasei()*phasei().thermo().Cp();
+    }
+
+    return tCp;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cp
+(
+    const scalarField& p,
+    const scalarField& T,
+    const label patchi
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<scalarField> tCp
+    (
+        phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi)
+    );
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tCp() +=
+            phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi);
+    }
+
+    return tCp;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cv() const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<volScalarField> tCv(phasei()*phasei().thermo().Cv());
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tCv() += phasei()*phasei().thermo().Cv();
+    }
+
+    return tCv;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cv
+(
+    const scalarField& p,
+    const scalarField& T,
+    const label patchi
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<scalarField> tCv
+    (
+        phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi)
+    );
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tCv() +=
+            phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi);
+    }
+
+    return tCv;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::gamma() const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<volScalarField> tgamma(phasei()*phasei().thermo().gamma());
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tgamma() += phasei()*phasei().thermo().gamma();
+    }
+
+    return tgamma;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::gamma
+(
+    const scalarField& p,
+    const scalarField& T,
+    const label patchi
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<scalarField> tgamma
+    (
+        phasei().boundaryField()[patchi]*phasei().thermo().gamma(p, T, patchi)
+    );
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tgamma() +=
+            phasei().boundaryField()[patchi]
+           *phasei().thermo().gamma(p, T, patchi);
+    }
+
+    return tgamma;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cpv() const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<volScalarField> tCpv(phasei()*phasei().thermo().Cpv());
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tCpv() += phasei()*phasei().thermo().Cpv();
+    }
+
+    return tCpv;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cpv
+(
+    const scalarField& p,
+    const scalarField& T,
+    const label patchi
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<scalarField> tCpv
+    (
+        phasei().boundaryField()[patchi]*phasei().thermo().Cpv(p, T, patchi)
+    );
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tCpv() +=
+            phasei().boundaryField()[patchi]
+           *phasei().thermo().Cpv(p, T, patchi);
+    }
+
+    return tCpv;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::CpByCpv() const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<volScalarField> tCpByCpv(phasei()*phasei().thermo().CpByCpv());
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tCpByCpv() += phasei()*phasei().thermo().CpByCpv();
+    }
+
+    return tCpByCpv;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::CpByCpv
+(
+    const scalarField& p,
+    const scalarField& T,
+    const label patchi
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<scalarField> tCpByCpv
+    (
+        phasei().boundaryField()[patchi]*phasei().thermo().CpByCpv(p, T, patchi)
+    );
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tCpByCpv() +=
+            phasei().boundaryField()[patchi]
+           *phasei().thermo().CpByCpv(p, T, patchi);
+    }
+
+    return tCpByCpv;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::kappa() const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<volScalarField> tkappa(phasei()*phasei().thermo().kappa());
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tkappa() += phasei()*phasei().thermo().kappa();
+    }
+
+    return tkappa;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::kappa
+(
+    const label patchi
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<scalarField> tkappa
+    (
+        phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi)
+    );
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tkappa() +=
+            phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi);
+    }
+
+    return tkappa;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::kappaEff
+(
+    const volScalarField& alphat
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<volScalarField> tkappaEff(phasei()*phasei().thermo().kappaEff(alphat));
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tkappaEff() += phasei()*phasei().thermo().kappaEff(alphat);
+    }
+
+    return tkappaEff;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::kappaEff
+(
+    const scalarField& alphat,
+    const label patchi
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<scalarField> tkappaEff
+    (
+        phasei().boundaryField()[patchi]
+       *phasei().thermo().kappaEff(alphat, patchi)
+    );
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        tkappaEff() +=
+            phasei().boundaryField()[patchi]
+           *phasei().thermo().kappaEff(alphat, patchi);
+    }
+
+    return tkappaEff;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::alphaEff
+(
+    const volScalarField& alphat
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<volScalarField> talphaEff(phasei()*phasei().thermo().alphaEff(alphat));
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        talphaEff() += phasei()*phasei().thermo().alphaEff(alphat);
+    }
+
+    return talphaEff;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::alphaEff
+(
+    const scalarField& alphat,
+    const label patchi
+) const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<scalarField> talphaEff
+    (
+        phasei().boundaryField()[patchi]
+       *phasei().thermo().alphaEff(alphat, patchi)
+    );
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        talphaEff() +=
+            phasei().boundaryField()[patchi]
+           *phasei().thermo().alphaEff(alphat, patchi);
+    }
+
+    return talphaEff;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::rCv() const
+{
+    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+
+    tmp<volScalarField> trCv(phasei()/phasei().thermo().Cv());
+
+    for (++phasei; phasei != phases_.end(); ++phasei)
+    {
+        trCv() += phasei()/phasei().thermo().Cv();
+    }
+
+    return trCv;
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::multiphaseMixtureThermo::surfaceTensionForce() const
+{
+    tmp<surfaceScalarField> tstf
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "surfaceTensionForce",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar
+            (
+                "surfaceTensionForce",
+                dimensionSet(1, -2, -2, 0, 0),
+                0.0
+            )
+        )
+    );
+
+    surfaceScalarField& stf = tstf();
+
+    forAllConstIter(PtrDictionary<phaseModel>, phases_, phase1)
+    {
+        const phaseModel& alpha1 = phase1();
+
+        PtrDictionary<phaseModel>::const_iterator phase2 = phase1;
+        ++phase2;
+
+        for (; phase2 != phases_.end(); ++phase2)
+        {
+            const phaseModel& alpha2 = phase2();
+
+            sigmaTable::const_iterator sigma =
+                sigmas_.find(interfacePair(alpha1, alpha2));
+
+            if (sigma == sigmas_.end())
+            {
+                FatalErrorIn
+                (
+                    "multiphaseMixtureThermo::surfaceTensionForce() const"
+                )   << "Cannot find interface " << interfacePair(alpha1, alpha2)
+                    << " in list of sigma values"
+                    << exit(FatalError);
+            }
+
+            stf += dimensionedScalar("sigma", dimSigma_, sigma())
+               *fvc::interpolate(K(alpha1, alpha2))*
+                (
+                    fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
+                  - fvc::interpolate(alpha1)*fvc::snGrad(alpha2)
+                );
+        }
+    }
+
+    return tstf;
+}
+
+
+void Foam::multiphaseMixtureThermo::solve()
+{
+    const Time& runTime = mesh_.time();
+
+    const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE");
+
+    label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
+
+    scalar cAlpha(readScalar(pimpleDict.lookup("cAlpha")));
+
+
+    volScalarField& alpha = phases_.first();
+
+    if (nAlphaSubCycles > 1)
+    {
+        surfaceScalarField rhoPhiSum(0.0*rhoPhi_);
+        dimensionedScalar totalDeltaT = runTime.deltaT();
+
+        for
+        (
+            subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
+        )
+        {
+            solveAlphas(cAlpha);
+            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
+        }
+
+        rhoPhi_ = rhoPhiSum;
+    }
+    else
+    {
+        solveAlphas(cAlpha);
+    }
+}
+
+
+Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixtureThermo::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::multiphaseMixtureThermo::nHatf
+(
+    const volScalarField& alpha1,
+    const volScalarField& alpha2
+) const
+{
+    // Face unit interface normal flux
+    return nHatfv(alpha1, alpha2) & 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::multiphaseMixtureThermo::correctContactAngle
+(
+    const phaseModel& alpha1,
+    const phaseModel& alpha2,
+    surfaceVectorField::GeometricBoundaryField& nHatb
+) const
+{
+    const volScalarField::GeometricBoundaryField& gbf
+        = alpha1.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(interfacePair(alpha1, alpha2));
+
+            if (tp == acap.thetaProps().end())
+            {
+                FatalErrorIn
+                (
+                    "multiphaseMixtureThermo::correctContactAngle"
+                    "(const phaseModel& alpha1, const phaseModel& alpha2, "
+                    "fvPatchVectorFieldField& nHatb) const"
+                )   << "Cannot find interface " << interfacePair(alpha1, alpha2)
+                    << "\n    in table of theta properties for patch "
+                    << acap.patch().name()
+                    << exit(FatalError);
+            }
+
+            bool matched = (tp.key().first() == alpha1.name());
+
+            scalar theta0 = convertToRad*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 = convertToRad*tp().thetaA(matched);
+                scalar thetaR = convertToRad*tp().thetaR(matched);
+
+                // Calculated the component of the velocity parallel to the wall
+                vectorField Uwall
+                (
+                    U_.boundaryField()[patchi].patchInternalField()
+                  - 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.0 - 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::multiphaseMixtureThermo::K
+(
+    const phaseModel& alpha1,
+    const phaseModel& alpha2
+) const
+{
+    tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
+
+    correctContactAngle(alpha1, alpha2, tnHatfv().boundaryField());
+
+    // Simple expression for curvature
+    return -fvc::div(tnHatfv & mesh_.Sf());
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixtureThermo::nearInterface() const
+{
+    tmp<volScalarField> tnearInt
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "nearInterface",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar("nearInterface", dimless, 0.0)
+        )
+    );
+
+    forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
+    {
+        tnearInt() = max(tnearInt(), pos(phase() - 0.01)*pos(0.99 - phase()));
+    }
+
+    return tnearInt;
+}
+
+
+void Foam::multiphaseMixtureThermo::solveAlphas
+(
+    const scalar cAlpha
+)
+{
+    static label nSolves=-1;
+    nSolves++;
+
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    surfaceScalarField phic(mag(phi_/mesh_.magSf()));
+    phic = min(cAlpha*phic, max(phic));
+
+    PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
+    int phasei = 0;
+
+    forAllIter(PtrDictionary<phaseModel>, phases_, phase)
+    {
+        phaseModel& alpha = phase();
+
+        phiAlphaCorrs.set
+        (
+            phasei,
+            new surfaceScalarField
+            (
+                fvc::flux
+                (
+                    phi_,
+                    alpha,
+                    alphaScheme
+                )
+            )
+        );
+
+        surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
+
+        forAllIter(PtrDictionary<phaseModel>, phases_, phase2)
+        {
+            phaseModel& alpha2 = phase2();
+
+            if (&alpha2 == &alpha) continue;
+
+            surfaceScalarField phir(phic*nHatf(alpha, alpha2));
+
+            phiAlphaCorr += fvc::flux
+            (
+                -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha,
+                alpharScheme
+            );
+        }
+
+        MULES::limit
+        (
+            geometricOneField(),
+            alpha,
+            phi_,
+            phiAlphaCorr,
+            zeroField(),
+            zeroField(),
+            1,
+            0,
+            3,
+            true
+        );
+
+        phasei++;
+    }
+
+    MULES::limitSum(phiAlphaCorrs);
+
+    rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0);
+
+    volScalarField sumAlpha
+    (
+        IOobject
+        (
+            "sumAlpha",
+            mesh_.time().timeName(),
+            mesh_
+        ),
+        mesh_,
+        dimensionedScalar("sumAlpha", dimless, 0)
+    );
+
+
+    volScalarField divU(fvc::div(fvc::absolute(phi_, U_)));
+
+
+    phasei = 0;
+
+    forAllIter(PtrDictionary<phaseModel>, phases_, phase)
+    {
+        phaseModel& alpha = phase();
+
+        surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
+        phiAlpha += upwind<scalar>(mesh_, phi_).flux(alpha);
+
+        volScalarField::DimensionedInternalField Sp
+        (
+            IOobject
+            (
+                "Sp",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar("Sp", alpha.dgdt().dimensions(), 0.0)
+        );
+
+        volScalarField::DimensionedInternalField Su
+        (
+            IOobject
+            (
+                "Su",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            // Divergence term is handled explicitly to be
+            // consistent with the explicit transport solution
+            divU*min(alpha, scalar(1))
+        );
+
+        {
+            const scalarField& dgdt = alpha.dgdt();
+
+            forAll(dgdt, celli)
+            {
+                if (dgdt[celli] < 0.0 && alpha[celli] > 0.0)
+                {
+                    Sp[celli] += dgdt[celli]*alpha[celli];
+                    Su[celli] -= dgdt[celli]*alpha[celli];
+                }
+                else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0)
+                {
+                    Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]);
+                }
+            }
+        }
+
+        forAllConstIter(PtrDictionary<phaseModel>, phases_, phase2)
+        {
+            const phaseModel& alpha2 = phase2();
+
+            if (&alpha2 == &alpha) continue;
+
+            const scalarField& dgdt2 = alpha2.dgdt();
+
+            forAll(dgdt2, celli)
+            {
+                if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0)
+                {
+                    Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]);
+                    Su[celli] += dgdt2[celli]*alpha[celli];
+                }
+                else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0)
+                {
+                    Sp[celli] += dgdt2[celli]*alpha2[celli];
+                }
+            }
+        }
+
+        MULES::explicitSolve
+        (
+            geometricOneField(),
+            alpha,
+            phiAlpha,
+            Sp,
+            Su
+        );
+
+        rhoPhi_ += fvc::interpolate(alpha.thermo().rho())*phiAlpha;
+
+        Info<< alpha.name() << " volume fraction, min, max = "
+            << alpha.weightedAverage(mesh_.V()).value()
+            << ' ' << min(alpha).value()
+            << ' ' << max(alpha).value()
+            << endl;
+
+        sumAlpha += alpha;
+
+        phasei++;
+    }
+
+    Info<< "Phase-sum volume fraction, min, max = "
+        << sumAlpha.weightedAverage(mesh_.V()).value()
+        << ' ' << min(sumAlpha).value()
+        << ' ' << max(sumAlpha).value()
+        << endl;
+
+    calcAlphas();
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H
new file mode 100644
index 00000000000..546048894f4
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H
@@ -0,0 +1,434 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::multiphaseMixtureThermo
+
+Description
+
+SourceFiles
+    multiphaseMixtureThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseMixtureThermo_H
+#define multiphaseMixtureThermo_H
+
+#include "phaseModel.H"
+#include "PtrDictionary.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "rhoThermo.H"
+#include "psiThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class multiphaseMixtureThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseMixtureThermo
+:
+    public psiThermo
+{
+public:
+
+    class interfacePair
+    :
+        public Pair<word>
+    {
+    public:
+
+        class hash
+        :
+            public Hash<interfacePair>
+        {
+        public:
+
+            hash()
+            {}
+
+            label operator()(const interfacePair& key) const
+            {
+                return word::hash()(key.first()) + word::hash()(key.second());
+            }
+        };
+
+
+        // Constructors
+
+            interfacePair()
+            {}
+
+            interfacePair(const word& alpha1Name, const word& alpha2Name)
+            :
+                Pair<word>(alpha1Name, alpha2Name)
+            {}
+
+            interfacePair(const phaseModel& alpha1, const phaseModel& alpha2)
+            :
+                Pair<word>(alpha1.name(), alpha2.name())
+            {}
+
+
+        // Friend Operators
+
+            friend bool operator==
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return
+                (
+                    ((a.first() == b.first()) && (a.second() == b.second()))
+                 || ((a.first() == b.second()) && (a.second() == b.first()))
+                );
+            }
+
+            friend bool operator!=
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return (!(a == b));
+            }
+    };
+
+
+private:
+
+    // Private data
+
+        //- Dictionary of phases
+        PtrDictionary<phaseModel> phases_;
+
+        const fvMesh& mesh_;
+        const volVectorField& U_;
+        const surfaceScalarField& phi_;
+
+        surfaceScalarField rhoPhi_;
+
+        volScalarField alphas_;
+
+        typedef HashTable<scalar, interfacePair, interfacePair::hash>
+            sigmaTable;
+
+        sigmaTable sigmas_;
+        dimensionSet dimSigma_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+        //- Conversion factor for degrees into radians
+        static const scalar convertToRad;
+
+
+    // Private member functions
+
+        void calcAlphas();
+
+        void solveAlphas(const scalar cAlpha);
+
+        tmp<surfaceVectorField> nHatfv
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        tmp<surfaceScalarField> nHatf
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        void correctContactAngle
+        (
+            const phaseModel& alpha1,
+            const phaseModel& alpha2,
+            surfaceVectorField::GeometricBoundaryField& nHatb
+        ) const;
+
+        tmp<volScalarField> K
+        (
+            const phaseModel& alpha1,
+            const phaseModel& alpha2
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiphaseMixtureThermo");
+
+
+    // Constructors
+
+        //- Construct from components
+        multiphaseMixtureThermo
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+
+    //- Destructor
+    virtual ~multiphaseMixtureThermo()
+    {}
+
+
+    // Member Functions
+
+        //- Return the phases
+        const PtrDictionary<phaseModel>& phases() const
+        {
+            return phases_;
+        }
+
+        //- Return non-const access to the phases
+        PtrDictionary<phaseModel>& phases()
+        {
+            return phases_;
+        }
+
+        //- Return the velocity
+        const volVectorField& U() const
+        {
+            return U_;
+        }
+
+        //- Return the volumetric flux
+        const surfaceScalarField& phi() const
+        {
+            return phi_;
+        }
+
+        const surfaceScalarField& rhoPhi() const
+        {
+            return rhoPhi_;
+        }
+
+        //- Update properties
+        virtual void correct();
+
+        //- Update densities for given pressure change
+        void correctRho(const volScalarField& dp);
+
+        //- Return true if the equation of state is incompressible
+        //  i.e. rho != f(p)
+        virtual bool incompressible() const;
+
+        //- Return true if the equation of state is isochoric
+        //  i.e. rho = const
+        virtual bool isochoric() const;
+
+
+        // Access to thermodynamic state variables
+
+            //- Enthalpy/Internal energy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& he()
+            {
+                notImplemented("multiphaseMixtureThermo::he()");
+                return phases_[0]->thermo().he();
+            }
+
+            //- Enthalpy/Internal energy [J/kg]
+            virtual const volScalarField& he() const
+            {
+                notImplemented("multiphaseMixtureThermo::he() const");
+                return phases_[0]->thermo().he();
+            }
+
+            //- Enthalpy/Internal energy
+            //  for given pressure and temperature [J/kg]
+            virtual tmp<volScalarField> he
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            ) const;
+
+            //- Enthalpy/Internal energy for cell-set [J/kg]
+            virtual tmp<scalarField> he
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
+            //- Enthalpy/Internal energy for patch [J/kg]
+            virtual tmp<scalarField> he
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Chemical enthalpy [J/kg]
+            virtual tmp<volScalarField> hc() const;
+
+            //- Temperature from enthalpy/internal energy for cell-set
+            virtual tmp<scalarField> THE
+            (
+                const scalarField& h,
+                const scalarField& p,
+                const scalarField& T0,      // starting temperature
+                const labelList& cells
+            ) const;
+
+            //- Temperature from enthalpy/internal energy for patch
+            virtual tmp<scalarField> THE
+            (
+                const scalarField& h,
+                const scalarField& p,
+                const scalarField& T0,      // starting temperature
+                const label patchi
+            ) const;
+
+
+        // Fields derived from thermodynamic state variables
+
+            //- Density [kg/m^3]
+            virtual tmp<volScalarField> rho() const;
+
+            //- Heat capacity at constant pressure [J/kg/K]
+            virtual tmp<volScalarField> Cp() const;
+
+            //- Heat capacity at constant pressure for patch [J/kg/K]
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant volume [J/kg/K]
+            virtual tmp<volScalarField> Cv() const;
+
+            //- Heat capacity at constant volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cv
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- gamma = Cp/Cv []
+            virtual tmp<volScalarField> gamma() const;
+
+            //- gamma = Cp/Cv for patch []
+            virtual tmp<scalarField> gamma
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure/volume [J/kg/K]
+            virtual tmp<volScalarField> Cpv() const;
+
+            //- Heat capacity at constant pressure/volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cpv
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity ratio []
+            virtual tmp<volScalarField> CpByCpv() const;
+
+            //- Heat capacity ratio for patch []
+            virtual tmp<scalarField> CpByCpv
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+
+        // Fields derived from transport state variables
+
+            //- Thermal diffusivity for temperature of mixture [J/m/s/K]
+            virtual tmp<volScalarField> kappa() const;
+
+            //- Thermal diffusivity of mixture for patch [J/m/s/K]
+            virtual tmp<scalarField> kappa
+            (
+                const label patchi
+            ) const;
+
+            //- Effective thermal diffusivity of mixture [J/m/s/K]
+            virtual tmp<volScalarField> kappaEff
+            (
+                const volScalarField& alphat
+            ) const;
+
+            //- Effective thermal diffusivity of mixture for patch [J/m/s/K]
+            virtual tmp<scalarField> kappaEff
+            (
+                const scalarField& alphat,
+                const label patchi
+            ) const;
+
+            //- Effective thermal diffusivity of mixture [J/m/s/K]
+            virtual tmp<volScalarField> alphaEff
+            (
+                const volScalarField& alphat
+            ) const;
+
+            //- Effective thermal diffusivity of mixture for patch [J/m/s/K]
+            virtual tmp<scalarField> alphaEff
+            (
+                const scalarField& alphat,
+                const label patchi
+            ) const;
+
+
+        //- Return the phase-averaged reciprocal Cv
+        tmp<volScalarField> rCv() const;
+
+        tmp<surfaceScalarField> surfaceTensionForce() const;
+
+        //- Indicator of the proximity of the interface
+        //  Field values are 1 near and 0 away for the interface.
+        tmp<volScalarField> nearInterface() const;
+
+        //- Solve for the mixture phase-fractions
+        void solve();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C
new file mode 100644
index 00000000000..1559f25a483
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C
@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "phaseModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phaseModel::phaseModel
+(
+    const word& phaseName,
+    const volScalarField& p,
+    const volScalarField& T
+)
+:
+    volScalarField
+    (
+        IOobject
+        (
+            IOobject::groupName("alpha", phaseName),
+            p.mesh().time().timeName(),
+            p.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        p.mesh()
+    ),
+    name_(phaseName),
+    p_(p),
+    T_(T),
+    thermo_(NULL),
+    dgdt_
+    (
+        IOobject
+        (
+            IOobject::groupName("dgdt", phaseName),
+            p.mesh().time().timeName(),
+            p.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        p.mesh(),
+        dimensionedScalar("0", dimless/dimTime, 0)
+    )
+{
+    {
+        volScalarField Tp(IOobject::groupName("T", phaseName), T);
+        Tp.write();
+    }
+
+    thermo_ = rhoThermo::New(p.mesh(), phaseName);
+    thermo_->validate(phaseName, "e");
+
+    correct();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::phaseModel> Foam::phaseModel::clone() const
+{
+    notImplemented("phaseModel::clone() const");
+    return autoPtr<phaseModel>(NULL);
+}
+
+
+void Foam::phaseModel::correct()
+{
+    thermo_->he() = thermo_->he(p_, T_);
+    thermo_->correct();
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H
new file mode 100644
index 00000000000..66d0ac8d63a
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H
@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::phaseModel
+
+Description
+    Single incompressible phase derived from the phase-fraction.
+    Used as part of the multiPhaseMixture for interface-capturing multi-phase
+    simulations.
+
+SourceFiles
+    phaseModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phaseModel_H
+#define phaseModel_H
+
+#include "rhoThermo.H"
+#include "volFields.H"
+#include "dictionaryEntry.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class phaseModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class phaseModel
+:
+    public volScalarField
+{
+    // Private data
+
+        word name_;
+        const volScalarField& p_;
+        const volScalarField& T_;
+        autoPtr<rhoThermo> thermo_;
+        volScalarField dgdt_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        phaseModel
+        (
+            const word& phaseName,
+            const volScalarField& p,
+            const volScalarField& T
+        );
+
+        //- Return clone
+        autoPtr<phaseModel> clone() const;
+
+        //- Return a pointer to a new phaseModel created on freestore
+        //  from Istream
+        class iNew
+        {
+            const volScalarField& p_;
+            const volScalarField& T_;
+
+        public:
+
+            iNew
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            )
+            :
+                p_(p),
+                T_(T)
+            {}
+
+            autoPtr<phaseModel> operator()(Istream& is) const
+            {
+                return autoPtr<phaseModel>(new phaseModel(is, p_, T_));
+            }
+        };
+
+
+    // Member Functions
+
+        const word& name() const
+        {
+            return name_;
+        }
+
+        const word& keyword() const
+        {
+            return name();
+        }
+
+        //- Return const-access to phase rhoThermo
+        const rhoThermo& thermo() const
+        {
+            return thermo_();
+        }
+
+        //- Return access to phase rhoThermo
+        rhoThermo& thermo()
+        {
+            return thermo_();
+        }
+
+        //- Return const-access to phase divergence
+        const volScalarField& dgdt() const
+        {
+            return dgdt_;
+        }
+
+        //- Return access to phase divergence
+        volScalarField& dgdt()
+        {
+            return dgdt_;
+        }
+
+        void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
new file mode 100644
index 00000000000..e42a54aabda
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
@@ -0,0 +1,142 @@
+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
+    );
+
+    surfaceScalarField phig
+    (
+        (
+            multiphaseProperties.surfaceTensionForce()
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
+
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
+    PtrList<fvScalarMatrix> p_rghEqnComps(multiphaseProperties.phases().size());
+
+    label phasei = 0;
+    forAllConstIter
+    (
+        PtrDictionary<phaseModel>,
+        multiphaseProperties.phases(),
+        phase
+    )
+    {
+        const rhoThermo& thermo = phase().thermo();
+        const volScalarField& rho = thermo.rho()();
+
+        p_rghEqnComps.set
+        (
+            phasei,
+            (
+                fvc::ddt(rho) + thermo.psi()*correction(fvm::ddt(p_rgh))
+              + fvc::div(phi, rho) - fvc::Sp(fvc::div(phi), rho)
+            ).ptr()
+        );
+
+        phasei++;
+    }
+
+    // Cache p_rgh prior to solve for density update
+    volScalarField p_rgh_0(p_rgh);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqnIncomp
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rAUf, p_rgh)
+        );
+
+        tmp<fvScalarMatrix> p_rghEqnComp;
+
+        phasei = 0;
+        forAllConstIter
+        (
+            PtrDictionary<phaseModel>,
+            multiphaseProperties.phases(),
+            phase
+        )
+        {
+            tmp<fvScalarMatrix> hmm
+            (
+                (max(phase(), scalar(0))/phase().thermo().rho())
+               *p_rghEqnComps[phasei]
+            );
+
+            if (phasei == 0)
+            {
+                p_rghEqnComp = hmm;
+            }
+            else
+            {
+                p_rghEqnComp() += hmm;
+            }
+
+            phasei++;
+        }
+
+        solve
+        (
+            p_rghEqnComp
+          + p_rghEqnIncomp,
+            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+        );
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            // p = max(p_rgh + multiphaseProperties.rho()*gh, pMin);
+            // p_rgh = p - multiphaseProperties.rho()*gh;
+
+            phasei = 0;
+            forAllIter
+            (
+                PtrDictionary<phaseModel>,
+                multiphaseProperties.phases(),
+                phase
+            )
+            {
+                phase().dgdt() =
+                    pos(phase())
+                  *(p_rghEqnComps[phasei] & p_rgh)/phase().thermo().rho();
+            }
+
+            phi = phiHbyA + p_rghEqnIncomp.flux();
+
+            U = HbyA
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+            U.correctBoundaryConditions();
+        }
+    }
+
+    p = max(p_rgh + multiphaseProperties.rho()*gh, pMin);
+
+    // Update densities from change in p_rgh
+    multiphaseProperties.correctRho(p_rgh - p_rgh_0);
+    rho = multiphaseProperties.rho();
+
+    K = 0.5*magSqr(U);
+
+    Info<< "max(U) " << max(mag(U)).value() << endl;
+    Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
+}
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/T b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/T
new file mode 100644
index 00000000000..1dc67cd149c
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/T
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 293;
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform 293;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        phi             rho*phi;
+        inletValue      $internalField;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/U b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/U
new file mode 100644
index 00000000000..7ea3a0c3232
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/U
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.air b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.air
new file mode 100644
index 00000000000..be8307fdf27
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.air
@@ -0,0 +1,79 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            alphaContactAngle;
+        thetaProperties
+        (
+            ( water air ) 90 0 0 0
+            ( oil air ) 90 0 0 0
+            ( mercury air ) 90 0 0 0
+            ( water oil ) 90 0 0 0
+            ( water mercury ) 90 0 0 0
+            ( oil mercury ) 90 0 0 0
+        );
+        value           uniform 0;
+    }
+    rightWall
+    {
+        type            alphaContactAngle;
+        thetaProperties
+        (
+            ( water air ) 90 0 0 0
+            ( oil air ) 90 0 0 0
+            ( mercury air ) 90 0 0 0
+            ( water oil ) 90 0 0 0
+            ( water mercury ) 90 0 0 0
+            ( oil mercury ) 90 0 0 0
+        );
+        value           uniform 1;
+    }
+    lowerWall
+    {
+        type            alphaContactAngle;
+        thetaProperties
+        (
+            ( water air ) 90 0 0 0
+            ( oil air ) 90 0 0 0
+            ( mercury air ) 90 0 0 0
+            ( water oil ) 90 0 0 0
+            ( water mercury ) 90 0 0 0
+            ( oil mercury ) 90 0 0 0
+        );
+        value           uniform 0;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.mercury b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.mercury
new file mode 100644
index 00000000000..d224de9509a
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.mercury
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.mercury;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+    rightWall
+    {
+        type            zeroGradient;
+    }
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.oil b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.oil
new file mode 100644
index 00000000000..bfcff63aedd
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.oil
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.oil;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+    rightWall
+    {
+        type            zeroGradient;
+    }
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.water b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.water
new file mode 100644
index 00000000000..cf96bb9d9f3
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.water
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+    rightWall
+    {
+        type            zeroGradient;
+    }
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alphas b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alphas
new file mode 100644
index 00000000000..10e3ea7a072
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alphas
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+    rightWall
+    {
+        type            zeroGradient;
+    }
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            zeroGradient;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p
new file mode 100644
index 00000000000..26edbd338c0
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    leftWall
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    rightWall
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    lowerWall
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    atmosphere
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p_rgh b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p_rgh
new file mode 100644
index 00000000000..16835050d0f
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p_rgh
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+
+    rightWall
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+
+    lowerWall
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+
+    atmosphere
+    {
+        type            totalPressure;
+        p0              uniform 1e5;
+        U               U;
+        phi             phi;
+        rho             rho;
+        psi             none;
+        gamma           1;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allclean b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allclean
new file mode 100755
index 00000000000..cd78b26599c
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allclean
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+
+\rm -rf 0
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allrun b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allrun
new file mode 100755
index 00000000000..887344985c0
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allrun
@@ -0,0 +1,17 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# Set application name
+application=`getApplication`
+
+\rm -rf 0
+cp -r 0.org 0
+
+runApplication blockMesh
+runApplication setFields
+runApplication $application
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/g b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/g
new file mode 100644
index 00000000000..e0ac2653b5b
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 -9.81 0 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/blockMeshDict b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/blockMeshDict
new file mode 100644
index 00000000000..11344be6ac8
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/blockMeshDict
@@ -0,0 +1,108 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.146;
+
+vertices
+(
+    (0 0 0)
+    (2 0 0)
+    (2.16438 0 0)
+    (4 0 0)
+    (0 0.32876 0)
+    (2 0.32876 0)
+    (2.16438 0.32876 0)
+    (4 0.32876 0)
+    (0 4 0)
+    (2 4 0)
+    (2.16438 4 0)
+    (4 4 0)
+    (0 0 0.1)
+    (2 0 0.1)
+    (2.16438 0 0.1)
+    (4 0 0.1)
+    (0 0.32876 0.1)
+    (2 0.32876 0.1)
+    (2.16438 0.32876 0.1)
+    (4 0.32876 0.1)
+    (0 4 0.1)
+    (2 4 0.1)
+    (2.16438 4 0.1)
+    (4 4 0.1)
+);
+
+blocks
+(
+    hex (0 1 5 4 12 13 17 16) (23 8 1) simpleGrading (1 1 1)
+    hex (2 3 7 6 14 15 19 18) (19 8 1) simpleGrading (1 1 1)
+    hex (4 5 9 8 16 17 21 20) (23 42 1) simpleGrading (1 1 1)
+    hex (5 6 10 9 17 18 22 21) (4 42 1) simpleGrading (1 1 1)
+    hex (6 7 11 10 18 19 23 22) (19 42 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    leftWall
+    {
+        type wall;
+        faces
+        (
+            (0 12 16 4)
+            (4 16 20 8)
+        );
+    }
+    rightWall
+    {
+        type wall;
+        faces
+        (
+            (7 19 15 3)
+            (11 23 19 7)
+        );
+    }
+    lowerWall
+    {
+        type wall;
+        faces
+        (
+            (0 1 13 12)
+            (1 5 17 13)
+            (5 6 18 17)
+            (2 14 18 6)
+            (2 3 15 14)
+        );
+    }
+    atmosphere
+    {
+        type patch;
+        faces
+        (
+            (8 20 21 9)
+            (9 21 22 10)
+            (10 22 23 11)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/boundary b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/boundary
new file mode 100644
index 00000000000..1b4dbb60aae
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/boundary
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+5
+(
+    leftWall
+    {
+        type            wall;
+        nFaces          50;
+        startFace       4432;
+    }
+    rightWall
+    {
+        type            wall;
+        nFaces          50;
+        startFace       4482;
+    }
+    lowerWall
+    {
+        type            wall;
+        nFaces          62;
+        startFace       4532;
+    }
+    atmosphere
+    {
+        type            patch;
+        nFaces          46;
+        startFace       4594;
+    }
+    defaultFaces
+    {
+        type            empty;
+        inGroups        1(empty);
+        nFaces          4536;
+        startFace       4640;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties
new file mode 100644
index 00000000000..129c3a278ec
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties
@@ -0,0 +1,33 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water oil mercury air);
+
+pMin            pMin [ 1 -1 -2 0 0 0 0 ] 10000;
+
+sigmas
+(
+    (air water) 0.07
+    (air oil) 0.07
+    (air mercury) 0.07
+    (water oil) 0.07
+    (water mercury) 0.07
+    (oil mercury) 0.07
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.air b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.air
new file mode 100644
index 00000000000..befc0aeae44
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.air
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   28.9;
+    }
+    thermodynamics
+    {
+        Cp          1007;
+        Hf          0;
+    }
+    transport
+    {
+        mu          1.84e-05;
+        Pr          0.7;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.mercury b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.mercury
new file mode 100644
index 00000000000..e90070ef131
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.mercury
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.mercury;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectFluid;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   200.59;
+    }
+    equationOfState
+    {
+        R           6818;
+        rho0        13529;
+    }
+    thermodynamics
+    {
+        Cp          139;
+        Hf          0;
+    }
+    transport
+    {
+        mu          1.522e-3;
+        Pr          0.022;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.oil b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.oil
new file mode 100644
index 00000000000..0bcdc33f4cb
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.oil
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.oil;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectFluid;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   100.21;
+    }
+    equationOfState
+    {
+        R           3564;
+        rho0        684;
+    }
+    thermodynamics
+    {
+        Cp          2240;
+        Hf          0;
+    }
+    transport
+    {
+        mu          3.76e-4;
+        Pr          6;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.water b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.water
new file mode 100644
index 00000000000..91e7adc381b
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.water
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectFluid;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   18.0;
+    }
+    equationOfState
+    {
+        R           7255;
+        rho0        1027;
+    }
+    thermodynamics
+    {
+        Cp          4195;
+        Hf          0;
+    }
+    transport
+    {
+        mu          3.645e-4;
+        Pr          2.289;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/turbulenceProperties b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/turbulenceProperties
new file mode 100644
index 00000000000..c2c3b28a1b4
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/controlDict b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/controlDict
new file mode 100644
index 00000000000..c2993ed4911
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/controlDict
@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     compressibleMultiphaseInterFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         10;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.05;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           0.5;
+maxAlphaCo      0.5;
+
+maxDeltaT       1;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/decomposeParDict b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/decomposeParDict
new file mode 100644
index 00000000000..e7e490bf74b
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 4;
+
+method          simple;
+
+simpleCoeffs
+{
+    n               ( 2 2 1 );
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               ( 1 1 1 );
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes
new file mode 100644
index 00000000000..9f3fac1dce5
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes
@@ -0,0 +1,64 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    div(rho*phi,U)  Gauss upwind;
+    div(phi,alpha)  Gauss vanLeer;
+    div(phirb,alpha) Gauss interfaceCompression;
+    "div\(phi,.*rho.*\)" Gauss upwind;
+    div(rho*phi,T)  Gauss upwind;
+    div(rho*phi,K)  Gauss upwind;
+    div(phi,p)      Gauss upwind;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    pcorr;
+    p_rgh;
+    "alpha.*";
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution
new file mode 100644
index 00000000000..6f88113f81b
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution
@@ -0,0 +1,120 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    pcorr
+    {
+        solver          PCG;
+        preconditioner
+        {
+            preconditioner  GAMG;
+            tolerance       1e-05;
+            relTol          0;
+            smoother        GaussSeidel;
+            nPreSweeps      0;
+            nPostSweeps     2;
+            nFinestSweeps   2;
+            cacheAgglomeration off;
+            nCellsInCoarsestLevel 10;
+            agglomerator    faceAreaPair;
+            mergeLevels     2;
+        }
+        tolerance       1e-05;
+        relTol          0;
+        maxIter         100;
+    }
+
+    ".*(rho|rhoFinal)"
+    {
+        solver          diagonal;
+    }
+
+    p_rgh
+    {
+        solver          GAMG;
+        tolerance       1e-07;
+        relTol          0.05;
+        smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        nFinestSweeps   2;
+        cacheAgglomeration on;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    p_rghFinal
+    {
+        solver          PCG;
+        preconditioner
+        {
+            preconditioner  GAMG;
+            tolerance       1e-07;
+            relTol          0;
+            nVcycles        2;
+            smoother        GaussSeidel;
+            nPreSweeps      0;
+            nPostSweeps     2;
+            nFinestSweeps   2;
+            cacheAgglomeration on;
+            nCellsInCoarsestLevel 10;
+            agglomerator    faceAreaPair;
+            mergeLevels     1;
+        }
+        tolerance       1e-07;
+        relTol          0;
+        maxIter         20;
+    }
+
+    "(U|T|k|B|nuTilda)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-08;
+        relTol          0.1;
+        nSweeps         1;
+    }
+
+    "(U|T|k|B|nuTilda)Final"
+    {
+        $U;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    nAlphaSubCycles 4;
+    cAlpha          2;
+}
+
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        "U.*"           1;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/setFieldsDict b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/setFieldsDict
new file mode 100644
index 00000000000..0a139a98ecd
--- /dev/null
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/setFieldsDict
@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha.air 1
+    volScalarFieldValue alpha.water 0
+    volScalarFieldValue alpha.oil 0
+    volScalarFieldValue alpha.mercury 0
+    volVectorFieldValue U ( 0 0 0 )
+);
+
+regions
+(
+    boxToCell
+    {
+        box ( 0 0 -1 ) ( 0.1461 0.292 1 );
+        fieldValues
+        (
+            volScalarFieldValue alpha.water 1
+            volScalarFieldValue alpha.oil 0
+            volScalarFieldValue alpha.mercury 0
+            volScalarFieldValue alpha.air 0
+        );
+    }
+    boxToCell
+    {
+        box ( 0.1461 0 -1 ) ( 0.2922 0.292 1 );
+        fieldValues
+        (
+            volScalarFieldValue alpha.water 0
+            volScalarFieldValue alpha.oil 1
+            volScalarFieldValue alpha.mercury 0
+            volScalarFieldValue alpha.air 0
+        );
+    }
+    boxToCell
+    {
+        box ( 0 0 -1 ) ( 0.1461 0.1 1 );
+        fieldValues
+        (
+            volScalarFieldValue alpha.water 0
+            volScalarFieldValue alpha.oil 0
+            volScalarFieldValue alpha.mercury 1
+            volScalarFieldValue alpha.air 0
+        );
+    }
+);
+
+
+// ************************************************************************* //
-- 
GitLab