diff --git a/applications/solvers/incompressible/porousSimpleFoam/Make/files b/applications/solvers/incompressible/porousSimpleFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..e72f7ed78c75d9aec66277a58f9807aee0073d75
--- /dev/null
+++ b/applications/solvers/incompressible/porousSimpleFoam/Make/files
@@ -0,0 +1,3 @@
+porousSimpleFoam.C
+
+EXE = $(FOAM_APPBIN)/porousSimpleFoam
diff --git a/applications/solvers/incompressible/porousSimpleFoam/Make/options b/applications/solvers/incompressible/porousSimpleFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..f490133ba67a4764a473d9a150b2b06f839a3790
--- /dev/null
+++ b/applications/solvers/incompressible/porousSimpleFoam/Make/options
@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I../simpleFoam \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lincompressibleRASModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools
diff --git a/applications/solvers/incompressible/porousSimpleFoam/UEqn.H b/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..c1925c7708c8fb821515092efcc420b22f816ad1
--- /dev/null
+++ b/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
@@ -0,0 +1,47 @@
+    // Construct the Momentum equation
+
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::div(phi, U)
+      - fvm::Sp(fvc::div(phi), U)
+      + turbulence->divDevReff(U)
+    );
+
+    UEqn().relax();
+
+    // Include the porous media resistance and solve the momentum equation
+    // either implicit in the tensorial resistance or transport using by
+    // including the spherical part of the resistance in the momentum diagonal
+
+    tmp<volScalarField> trAU;
+    tmp<volTensorField> trTU;
+
+    if (pressureImplicitPorosity)
+    {
+        tmp<volTensorField> tTU = tensor(I)*UEqn().A();
+        pZones.addResistance(UEqn(), tTU());
+        trTU = inv(tTU());
+        trTU().rename("rAU");
+
+        volVectorField gradp = fvc::grad(p);
+
+        for (int UCorr=0; UCorr<nUCorr; UCorr++)
+        {
+            U = trTU() & (UEqn().H() - gradp);
+        }
+        U.correctBoundaryConditions();
+    }
+    else
+    {
+        pZones.addResistance(UEqn());
+
+        eqnResidual = solve
+        (
+            UEqn() == -fvc::grad(p)
+        ). initialResidual();
+
+        maxResidual = max(eqnResidual, maxResidual);
+
+        trAU = 1.0/UEqn().A();
+        trAU().rename("rAU");
+    }
diff --git a/applications/solvers/incompressible/porousSimpleFoam/createFields.H b/applications/solvers/incompressible/porousSimpleFoam/createFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..04d57d0571472ea1ee00c9a8ec1e24171eefe640
--- /dev/null
+++ b/applications/solvers/incompressible/porousSimpleFoam/createFields.H
@@ -0,0 +1,64 @@
+    Info << "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            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"
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::RASModel> turbulence
+    (
+        incompressible::RASModel::New(U, phi, laminarTransport)
+    );
+
+
+    porousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);
+
+    int nUCorr = 0;
+    if (pZones.size())
+    {
+        // nUCorrectors for pressureImplicitPorosity
+        if (mesh.solutionDict().subDict("SIMPLE").found("nUCorrectors"))
+        {
+            nUCorr = readInt
+            (
+                mesh.solutionDict().subDict("SIMPLE").lookup("nUCorrectors")
+            );
+        }
+
+        if (nUCorr > 0)
+        {
+            pressureImplicitPorosity = true;
+        }
+    }
diff --git a/applications/solvers/incompressible/porousSimpleFoam/pEqn.H b/applications/solvers/incompressible/porousSimpleFoam/pEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..e6797cd8603b8032e1c901583b6ddb15992b7ddc
--- /dev/null
+++ b/applications/solvers/incompressible/porousSimpleFoam/pEqn.H
@@ -0,0 +1,59 @@
+if (pressureImplicitPorosity)
+{
+    U = trTU()&UEqn().H();
+}
+else
+{
+    U = trAU()*UEqn().H();
+}
+
+UEqn.clear();
+phi = fvc::interpolate(U) & mesh.Sf();
+adjustPhi(phi, U, p);
+
+for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+{
+    tmp<fvScalarMatrix> tpEqn;
+
+    if (pressureImplicitPorosity)
+    {
+        tpEqn = (fvm::laplacian(trTU(), p) == fvc::div(phi));
+    }
+    else
+    {
+        tpEqn = (fvm::laplacian(trAU(), p) == fvc::div(phi));
+    }
+
+    tpEqn().setReference(pRefCell, pRefValue);
+    // retain the residual from the first iteration
+    if (nonOrth == 0)
+    {
+        eqnResidual = tpEqn().solve().initialResidual();
+        maxResidual = max(eqnResidual, maxResidual);
+    }
+    else
+    {
+        tpEqn().solve();
+    }
+
+    if (nonOrth == nNonOrthCorr)
+    {
+        phi -= tpEqn().flux();
+    }
+}
+
+#include "continuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+if (pressureImplicitPorosity)
+{
+    U -= trTU()&fvc::grad(p);
+}
+else
+{
+    U -= trAU()*fvc::grad(p);
+}
+
+U.correctBoundaryConditions();
diff --git a/applications/solvers/incompressible/porousSimpleFoam/porousSimpleFoam.C b/applications/solvers/incompressible/porousSimpleFoam/porousSimpleFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..b78f1688fe9bacd942529d1edc031e7fce981eff
--- /dev/null
+++ b/applications/solvers/incompressible/porousSimpleFoam/porousSimpleFoam.C
@@ -0,0 +1,85 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    porousSimpleFoam
+
+Description
+    Steady-state solver for incompressible, turbulent flow with
+    implicit or explicit porosity treatment
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "RASModel.H"
+#include "porousZones.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "readSIMPLEControls.H"
+        #include "initConvergenceCheck.H"
+
+        p.storePrevIter();
+
+        // Pressure-velocity SIMPLE corrector
+        {
+            #include "UEqn.H"
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        #include "convergenceCheck.H"
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interMixingFoam/Make/files b/applications/solvers/multiphase/interMixingFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..488cd77e56dca265db6dab1b98933df8fe9ae72a
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/Make/files
@@ -0,0 +1,6 @@
+incompressibleThreePhaseMixture/threePhaseMixture.C
+threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
+interMixingFoam.C
+
+EXE = $(FOAM_APPBIN)/interMixingFoam
+
diff --git a/applications/solvers/multiphase/interMixingFoam/Make/options b/applications/solvers/multiphase/interMixingFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..d8e4da2313e4343d0c6200c7fe9c7e8d2c090d48
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/Make/options
@@ -0,0 +1,17 @@
+INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam
+
+EXE_INC = \
+    -I$(INTERFOAM) \
+    -IincompressibleThreePhaseMixture \
+    -IthreePhaseInterfaceProperties \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/transportModels
+
+EXE_LIBS = \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume
diff --git a/applications/solvers/multiphase/interMixingFoam/alphaEqns.H b/applications/solvers/multiphase/interMixingFoam/alphaEqns.H
new file mode 100644
index 0000000000000000000000000000000000000000..f9bad0a7058577487b0cb6e4fd7bb14f1e9ee24d
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/alphaEqns.H
@@ -0,0 +1,164 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    surfaceScalarField phir
+    (
+        IOobject
+        (
+            "phir",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        interface.cAlpha()*mag(phi/mesh.magSf())*interface.nHatf()
+    );
+
+    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
+    {
+        // Create the limiter to be used for all phase-fractions
+        scalarField allLambda(mesh.nFaces(), 1.0);
+
+        // Split the limiter into a surfaceScalarField
+        slicedSurfaceScalarField lambda
+        (
+            IOobject
+            (
+                "lambda",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimless,
+            allLambda
+        );
+
+
+        // Create the complete convection flux for alpha1
+        surfaceScalarField phiAlpha1 =
+            fvc::flux
+            (
+                phi,
+                alpha1,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha1,
+                alpharScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(-phir, alpha3, alpharScheme),
+                alpha1,
+                alpharScheme
+            );
+
+        // Create the bounded (upwind) flux for alpha1
+        surfaceScalarField phiAlpha1BD =
+            upwind<scalar>(mesh, phi).flux(alpha1);
+
+        // Calculate the flux correction for alpha1
+        phiAlpha1 -= phiAlpha1BD;
+
+        // Calculate the limiter for alpha1
+        MULES::limiter
+        (
+            allLambda,
+            oneField(),
+            alpha1,
+            phiAlpha1BD,
+            phiAlpha1,
+            zeroField(),
+            zeroField(),
+            1,
+            0,
+            3
+        );
+
+        // Create the complete flux for alpha2
+        surfaceScalarField phiAlpha2 =
+            fvc::flux
+            (
+                phi,
+                alpha2,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(phir, alpha1, alpharScheme),
+                alpha2,
+                alpharScheme
+            );
+
+        // Create the bounded (upwind) flux for alpha2
+        surfaceScalarField phiAlpha2BD =
+            upwind<scalar>(mesh, phi).flux(alpha2);
+
+        // Calculate the flux correction for alpha2
+        phiAlpha2 -= phiAlpha2BD;
+
+        // Further limit the limiter for alpha2
+        MULES::limiter
+        (
+            allLambda,
+            oneField(),
+            alpha2,
+            phiAlpha2BD,
+            phiAlpha2,
+            zeroField(),
+            zeroField(),
+            1,
+            0,
+            3
+        );
+
+        // Construct the limited fluxes
+        phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1;
+        phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2;
+
+        // Solve for alpha1
+        solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1));
+
+        // Create the diffusion coefficients for alpha2<->alpha3
+        volScalarField Dc23 = D23*max(alpha3, scalar(0))*pos(alpha2);
+        volScalarField Dc32 = D23*max(alpha2, scalar(0))*pos(alpha3);
+
+        // Add the diffusive flux for alpha3->alpha2
+        phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
+
+        // Solve for alpha2
+        fvScalarMatrix alpha2Eqn
+        (
+            fvm::ddt(alpha2)
+          + fvc::div(phiAlpha2)
+          - fvm::laplacian(Dc23 + Dc32, alpha2)
+        );
+        alpha2Eqn.solve();
+
+        // Construct the complete mass flux
+        rhoPhi =
+              phiAlpha1*(rho1 - rho3)
+            + (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3)
+            + phi*rho3;
+
+        alpha3 = 1.0 - alpha1 - alpha2;
+    }
+
+    Info<< "Air phase volume fraction = "
+        << alpha1.weightedAverage(mesh.V()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
+        << endl;
+
+    Info<< "Liquid phase volume fraction = "
+        << alpha2.weightedAverage(mesh.V()).value()
+        << "  Min(alpha2) = " << min(alpha2).value()
+        << "  Max(alpha2) = " << max(alpha2).value()
+        << endl;
+}
diff --git a/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H
new file mode 100644
index 0000000000000000000000000000000000000000..765087a183baf3db62ec9b4dcf382d667b6209f3
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H
@@ -0,0 +1,43 @@
+label nAlphaCorr
+(
+    readLabel(piso.lookup("nAlphaCorr"))
+);
+
+label nAlphaSubCycles
+(
+    readLabel(piso.lookup("nAlphaSubCycles"))
+);
+
+if (nAlphaSubCycles > 1)
+{
+    surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+#       include "alphaEqns.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+#       include "alphaEqns.H"
+}
+
+interface.correct();
+
+{
+    volScalarField rhoNew = alpha1*rho1 + alpha2*rho2 + alpha3*rho3;
+
+    //solve(fvm::ddt(rho) + fvc::div(rhoPhi));
+    //Info<< "density error = "
+    //    << max((mag(rho - rhoNew)/mag(rhoNew))().internalField()) << endl;
+
+    rho == rhoNew;
+}
diff --git a/applications/solvers/multiphase/interMixingFoam/createFields.H b/applications/solvers/multiphase/interMixingFoam/createFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..9615e6c18645286de8cbdf87ba5c1a127b251236
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/createFields.H
@@ -0,0 +1,132 @@
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field alpha1\n" << endl;
+    volScalarField alpha1
+    (
+        IOobject
+        (
+            "alpha1",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    Info<< "Reading field alpha2\n" << endl;
+    volScalarField alpha2
+    (
+        IOobject
+        (
+            "alpha2",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    Info<< "Reading field alpha3\n" << endl;
+    volScalarField alpha3
+    (
+        IOobject
+        (
+            "alpha3",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    alpha3 == 1.0 - alpha1 - alpha2;
+
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+    threePhaseMixture threePhaseProperties(U, phi);
+
+    const dimensionedScalar& rho1 = threePhaseProperties.rho1();
+    const dimensionedScalar& rho2 = threePhaseProperties.rho2();
+    const dimensionedScalar& rho3 = threePhaseProperties.rho3();
+
+    dimensionedScalar D23(threePhaseProperties.lookup("D23"));
+
+    // Need to store rho for ddt(rho, U)
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT
+        ),
+        alpha1*rho1 + alpha2*rho2 + alpha3*rho3,
+        alpha1.boundaryField().types()
+    );
+    rho.oldTime();
+
+
+    // Mass flux
+    // Initialisation does not matter because rhoPhi is reset after the
+    // alpha solution before it is used in the U equation.
+    surfaceScalarField rhoPhi
+    (
+        IOobject
+        (
+            "rho*phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        rho1*phi
+    );
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+
+
+    // Construct interface from alpha distribution
+    threePhaseInterfaceProperties interface(threePhaseProperties);
+
+
+    // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, threePhaseProperties)
+    );
diff --git a/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C
new file mode 100644
index 0000000000000000000000000000000000000000..f27e17923f952e5fafec44cf6536241071ad9743
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C
@@ -0,0 +1,204 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Class
+    threePhaseMixture
+
+\*---------------------------------------------------------------------------*/
+
+#include "threePhaseMixture.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceFields.H"
+#include "fvc.H"
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+//- Calculate and return the laminar viscosity
+void Foam::threePhaseMixture::calcNu()
+{
+    // Average kinematic viscosity calculated from dynamic viscosity
+    nu_ = mu()/(alpha1_*rho1_ + alpha2_*rho2_ + alpha3_*rho3_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::threePhaseMixture::threePhaseMixture
+(
+    const volVectorField& U,
+    const surfaceScalarField& phi
+)
+:
+    transportModel(U, phi),
+
+    phase1Name_("phase1"),
+    phase2Name_("phase2"),
+    phase3Name_("phase3"),
+
+    nuModel1_
+    (
+        viscosityModel::New
+        (
+            "nu1",
+            subDict(phase1Name_),
+            U,
+            phi
+        )
+    ),
+    nuModel2_
+    (
+        viscosityModel::New
+        (
+            "nu2",
+            subDict(phase2Name_),
+            U,
+            phi
+        )
+    ),
+    nuModel3_
+    (
+        viscosityModel::New
+        (
+            "nu3",
+            subDict(phase2Name_),
+            U,
+            phi
+        )
+    ),
+
+    rho1_(nuModel1_->viscosityProperties().lookup("rho")),
+    rho2_(nuModel2_->viscosityProperties().lookup("rho")),
+    rho3_(nuModel3_->viscosityProperties().lookup("rho")),
+
+    U_(U),
+    phi_(phi),
+
+    alpha1_(U_.db().lookupObject<const volScalarField> ("alpha1")),
+    alpha2_(U_.db().lookupObject<const volScalarField> ("alpha2")),
+    alpha3_(U_.db().lookupObject<const volScalarField> ("alpha3")),
+
+    nu_
+    (
+        IOobject
+        (
+            "nu",
+            U_.time().timeName(),
+            U_.db()
+        ),
+        U_.mesh(),
+        dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0),
+        calculatedFvPatchScalarField::typeName
+    )
+{
+    calcNu();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::threePhaseMixture::mu() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            "mu",
+            alpha1_*rho1_*nuModel1_->nu()
+          + alpha2_*rho2_*nuModel2_->nu()
+          + alpha3_*rho3_*nuModel3_->nu()
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::muf() const
+{
+    surfaceScalarField alpha1f = fvc::interpolate(alpha1_);
+    surfaceScalarField alpha2f = fvc::interpolate(alpha2_);
+    surfaceScalarField alpha3f = fvc::interpolate(alpha3_);
+
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            "mu",
+            alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
+          + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu())
+          + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu())
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::nuf() const
+{
+    surfaceScalarField alpha1f = fvc::interpolate(alpha1_);
+    surfaceScalarField alpha2f = fvc::interpolate(alpha2_);
+    surfaceScalarField alpha3f = fvc::interpolate(alpha3_);
+
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            "nu",
+            (
+                alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
+              + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu())
+              + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu())
+            )/(alpha1f*rho1_ + alpha2f*rho2_ + alpha3f*rho3_)
+        )
+    );
+}
+
+
+bool Foam::threePhaseMixture::read()
+{
+    if (transportModel::read())
+    {
+        if
+        (
+            nuModel1_().read(*this)
+         && nuModel2_().read(*this)
+         && nuModel3_().read(*this)
+        )
+        {
+            nuModel1_->viscosityProperties().lookup("rho") >> rho1_;
+            nuModel2_->viscosityProperties().lookup("rho") >> rho2_;
+            nuModel3_->viscosityProperties().lookup("rho") >> rho3_;
+
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H
new file mode 100644
index 0000000000000000000000000000000000000000..2712b5d0cbfc18500f47c9b8db5661a5c0721fa6
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H
@@ -0,0 +1,197 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Class
+    threePhaseMixture
+
+Description
+
+SourceFiles
+    threePhaseMixture.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef threePhaseMixture_H
+#define threePhaseMixture_H
+
+#include "incompressible/transportModel/transportModel.H"
+#include "incompressible/viscosityModels/viscosityModel/viscosityModel.H"
+#include "dimensionedScalar.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class threePhaseMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+class threePhaseMixture
+:
+    public transportModel
+{
+    // Private data
+
+        word phase1Name_;
+        word phase2Name_;
+        word phase3Name_;
+
+        autoPtr<viscosityModel> nuModel1_;
+        autoPtr<viscosityModel> nuModel2_;
+        autoPtr<viscosityModel> nuModel3_;
+
+        dimensionedScalar rho1_;
+        dimensionedScalar rho2_;
+        dimensionedScalar rho3_;
+
+        const volVectorField& U_;
+        const surfaceScalarField& phi_;
+
+        const volScalarField& alpha1_;
+        const volScalarField& alpha2_;
+        const volScalarField& alpha3_;
+
+        volScalarField nu_;
+
+
+    // Private Member Functions
+
+        //- Calculate and return the laminar viscosity
+        void calcNu();
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        threePhaseMixture
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+
+    // Destructor
+
+        ~threePhaseMixture()
+        {}
+
+
+    // Member Functions
+
+        //- Return const-access to phase1 viscosityModel
+        const viscosityModel& nuModel1() const
+        {
+            return nuModel1_();
+        }
+
+        //- Return const-access to phase2 viscosityModel
+        const viscosityModel& nuModel2() const
+        {
+            return nuModel2_();
+        }
+
+        //- Return const-access to phase3 viscosityModel
+        const viscosityModel& nuModel3() const
+        {
+            return nuModel3_();
+        }
+
+        //- Return const-access to phase1 density
+        const dimensionedScalar& rho1() const
+        {
+            return rho1_;
+        }
+
+        //- Return const-access to phase2 density
+        const dimensionedScalar& rho2() const
+        {
+            return rho2_;
+        };
+
+        //- Return const-access to phase3 density
+        const dimensionedScalar& rho3() const
+        {
+            return rho3_;
+        };
+
+        const volScalarField& alpha1() const
+        {
+            return alpha1_;
+        }
+
+        const volScalarField& alpha2() const
+        {
+            return alpha2_;
+        }
+
+        const volScalarField& alpha3() const
+        {
+            return alpha3_;
+        }
+
+        //- Return the velocity
+        const volVectorField& U() const
+        {
+            return U_;
+        }
+
+        //- Return the dynamic laminar viscosity
+        tmp<volScalarField> mu() const;
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> muf() const;
+
+        //- Return the kinematic laminar viscosity
+        tmp<volScalarField> nu() const
+        {
+            return nu_;
+        }
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> nuf() const;
+
+        //- Correct the laminar viscosity
+        void correct()
+        {
+            calcNu();
+        }
+
+        //- Read base transportProperties dictionary
+        bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..a18df8d12a59ee7d3914109374467a5a10e14a80
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Application
+    interMixingFoam
+
+Description
+    Solver for 3 incompressible fluids, two of which are miscible,
+    using a VOF method to capture the interface.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "MULES.H"
+#include "subCycle.H"
+#include "threePhaseMixture.H"
+#include "threePhaseInterfaceProperties.H"
+#include "turbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+    #include "correctPhi.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readPISOControls.H"
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        threePhaseProperties.correct();
+
+        #include "alphaEqnsSubCycle.H"
+
+        #define twoPhaseProperties threePhaseProperties
+        #include "UEqn.H"
+
+        // --- PISO loop
+        for (int corr=0; corr<nCorr; corr++)
+        {
+            #include "pEqn.H"
+        }
+
+        #include "continuityErrs.H"
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "\n end \n";
+
+    return(0);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
new file mode 100644
index 0000000000000000000000000000000000000000..1ab10bd74348be3faafce4e5db2b66c57c5f7a05
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
@@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Application
+    threePhaseInterfaceProperties
+
+Description
+    Properties to aid interFoam :
+    1. Correct the alpha boundary condition for dynamic contact angle.
+    2. Calculate interface curvature.
+
+\*---------------------------------------------------------------------------*/
+
+#include "threePhaseInterfaceProperties.H"
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "mathematicalConstants.H"
+#include "surfaceInterpolate.H"
+#include "fvcDiv.H"
+#include "fvcGrad.H"
+
+// * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
+
+const Foam::scalar Foam::threePhaseInterfaceProperties::convertToRad =
+    Foam::constant::mathematical::pi/180.0;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// 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::threePhaseInterfaceProperties::correctContactAngle
+(
+    surfaceVectorField::GeometricBoundaryField& nHatb
+) const
+{
+    const volScalarField::GeometricBoundaryField& alpha1 =
+        mixture_.alpha1().boundaryField();
+    const volScalarField::GeometricBoundaryField& alpha2 =
+        mixture_.alpha2().boundaryField();
+    const volScalarField::GeometricBoundaryField& alpha3 =
+        mixture_.alpha3().boundaryField();
+    const volVectorField::GeometricBoundaryField& U =
+        mixture_.U().boundaryField();
+
+    const fvMesh& mesh = mixture_.U().mesh();
+    const fvBoundaryMesh& boundary = mesh.boundary();
+
+    forAll(boundary, patchi)
+    {
+        if (isA<alphaContactAngleFvPatchScalarField>(alpha1[patchi]))
+        {
+            const alphaContactAngleFvPatchScalarField& a2cap =
+                refCast<const alphaContactAngleFvPatchScalarField>
+                (alpha2[patchi]);
+
+            const alphaContactAngleFvPatchScalarField& a3cap =
+                refCast<const alphaContactAngleFvPatchScalarField>
+                (alpha3[patchi]);
+
+            scalarField twoPhaseAlpha2 = max(a2cap, scalar(0));
+            scalarField twoPhaseAlpha3 = max(a3cap, scalar(0));
+
+            scalarField sumTwoPhaseAlpha =
+                twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL;
+
+            twoPhaseAlpha2 /= sumTwoPhaseAlpha;
+            twoPhaseAlpha3 /= sumTwoPhaseAlpha;
+
+            fvsPatchVectorField& nHatp = nHatb[patchi];
+
+            scalarField theta =
+                convertToRad
+               *(
+                   twoPhaseAlpha2*(180 - a2cap.theta(U[patchi], nHatp))
+                 + twoPhaseAlpha3*(180 - a3cap.theta(U[patchi], nHatp))
+               );
+
+            vectorField nf = boundary[patchi].nf();
+
+            // Reset nHatPatch to correspond to the contact angle
+
+            scalarField a12 = nHatp & nf;
+
+            scalarField b1 = cos(theta);
+
+            scalarField b2(nHatp.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;
+
+            nHatp = a*nf + b*nHatp;
+
+            nHatp /= (mag(nHatp) + deltaN_.value());
+        }
+    }
+}
+
+
+void Foam::threePhaseInterfaceProperties::calculateK()
+{
+    const volScalarField& alpha1 = mixture_.alpha1();
+
+    const fvMesh& mesh = alpha1.mesh();
+    const surfaceVectorField& Sf = mesh.Sf();
+
+    // Cell gradient of alpha
+    volVectorField gradAlpha = fvc::grad(alpha1);
+
+    // Interpolated face-gradient of alpha
+    surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
+
+    // Face unit interface normal
+    surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_);
+    correctContactAngle(nHatfv.boundaryField());
+
+    // Face unit interface normal flux
+    nHatf_ = nHatfv & Sf;
+
+    // Simple expression for curvature
+    K_ = -fvc::div(nHatf_);
+
+    // Complex expression for curvature.
+    // Correction is formally zero but numerically non-zero.
+    //volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
+    //nHat.boundaryField() = nHatfv.boundaryField();
+    //K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
+(
+    const threePhaseMixture& mixture
+)
+:
+    mixture_(mixture),
+    cAlpha_
+    (
+        readScalar
+        (
+            mixture.U().mesh().solutionDict().subDict("PISO").lookup("cAlpha")
+        )
+    ),
+    sigma12_(mixture.lookup("sigma12")),
+    sigma13_(mixture.lookup("sigma13")),
+
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(mixture.U().mesh().V()), 1.0/3.0)
+    ),
+
+    nHatf_
+    (
+        (
+            fvc::interpolate(fvc::grad(mixture.alpha1()))
+           /(mag(fvc::interpolate(fvc::grad(mixture.alpha1()))) + deltaN_)
+        ) & mixture.alpha1().mesh().Sf()
+    ),
+
+    K_
+    (
+        IOobject
+        (
+            "K",
+            mixture.alpha1().time().timeName(),
+            mixture.alpha1().mesh()
+        ),
+        -fvc::div(nHatf_)
+    )
+{
+    calculateK();
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H
new file mode 100644
index 0000000000000000000000000000000000000000..3a402feeedf37ebda540220c8773b76b4cb4f6aa
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H
@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Class
+    threePhaseInterfaceProperties
+
+Description
+    Properties to aid interFoam :
+    1. Correct the alpha boundary condition for dynamic contact angle.
+    2. Calculate interface curvature.
+
+SourceFiles
+    threePhaseInterfaceProperties.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef threePhaseInterfaceProperties_H
+#define threePhaseInterfaceProperties_H
+
+#include "threePhaseMixture.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Class threePhaseInterfaceProperties Declaration
+\*---------------------------------------------------------------------------*/
+
+class threePhaseInterfaceProperties
+{
+    // Private data
+
+        const threePhaseMixture& mixture_;
+
+        //- Compression coefficient
+        scalar cAlpha_;
+
+        //- Surface tension 1-2
+        dimensionedScalar sigma12_;
+
+        //- Surface tension 1-3
+        dimensionedScalar sigma13_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+        surfaceScalarField nHatf_;
+        volScalarField K_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct and assignment
+        threePhaseInterfaceProperties(const threePhaseInterfaceProperties&);
+        void operator=(const threePhaseInterfaceProperties&);
+
+        //- Correction for the boundary condition on the unit normal nHat on
+        //  walls to produce the correct contact dynamic angle
+        //  calculated from the component of U parallel to the wall
+        void correctContactAngle
+        (
+            surfaceVectorField::GeometricBoundaryField& nHat
+        ) const;
+
+        //- Re-calculate the interface curvature
+        void calculateK();
+
+
+public:
+
+    //- Conversion factor for degrees into radians
+    static const scalar convertToRad;
+
+
+    // Constructors
+
+        //- Construct from volume fraction field alpha and IOdictionary
+        threePhaseInterfaceProperties(const threePhaseMixture& mixture);
+
+
+    // Member Functions
+
+        scalar cAlpha() const
+        {
+            return cAlpha_;
+        }
+
+        const dimensionedScalar& deltaN() const
+        {
+            return deltaN_;
+        }
+
+        const surfaceScalarField& nHatf() const
+        {
+            return nHatf_;
+        }
+
+        const volScalarField& K() const
+        {
+            return K_;
+        }
+
+        tmp<volScalarField> sigma() const
+        {
+            volScalarField limitedAlpha2 = max(mixture_.alpha2(), scalar(0));
+            volScalarField limitedAlpha3 = max(mixture_.alpha3(), scalar(0));
+
+            return
+                (limitedAlpha2*sigma12_ + limitedAlpha3*sigma13_)
+               /(limitedAlpha2 + limitedAlpha3 + SMALL);
+        }
+
+        tmp<volScalarField> sigmaK() const
+        {
+            return sigma()*K_;
+        }
+
+        void correct()
+        {
+            calculateK();
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/0 b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/0
new file mode 120000
index 0000000000000000000000000000000000000000..f1c4a884b51ae4171513426c3b7f5f585d0d889e
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/0
@@ -0,0 +1 @@
+../angledDuctImplicit/0
\ No newline at end of file
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/constant b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/constant
new file mode 120000
index 0000000000000000000000000000000000000000..28205c782b36471b118c5be2948a56e345961ad4
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/constant
@@ -0,0 +1 @@
+../angledDuctImplicit/constant
\ No newline at end of file
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/controlDict b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..9ffe4063cf457efe9e2a167ad38617dc42f6dd8b
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/controlDict
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     rhoPorousSimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         200;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   10;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+graphFormat     raw;
+
+runTimeModifiable yes;
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSchemes b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..a70c35e1bf81ab94f439453dfac9b9349f09bcae
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSchemes
@@ -0,0 +1,64 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(U)         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+divSchemes
+{
+    div(phi,U)      Gauss upwind;
+    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div(phi,epsilon) Gauss upwind;
+    div(phi,k)      Gauss upwind;
+}
+
+laplacianSchemes
+{
+    laplacian(nuEff,U) Gauss linear corrected;
+    laplacian(rAU,p) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(1,p)  Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSolution b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..b7d4dbd573caabf8c8e753563360f72f2c29af0b
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSolution
@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.05;
+        smoother        GaussSeidel;
+        cacheAgglomeration off;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        nSweeps         2;
+        tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    "(k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        nSweeps         2;
+        tolerance       1e-07;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    p               0.3;
+    U               0.7;
+    k               0.9;
+    epsilon         0.9;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/T b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/T
new file mode 100644
index 0000000000000000000000000000000000000000..e1df94398773b013813d4f9f37769c8587889a1c
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/T
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 293;
+
+boundaryField
+{
+
+    front
+    {
+        type            zeroGradient;
+    }
+    back
+    {
+        type            zeroGradient;
+    }
+    wall
+    {
+        type            zeroGradient;
+    }
+    porosityWall
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        value           $internalField;
+        inletValue      $internalField;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/U b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/U
new file mode 100644
index 0000000000000000000000000000000000000000..88f509425865e94093e76cca468b5839e437ab77
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/U
@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    front
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    back
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    wall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    porosityWall
+    {
+        type            slip;
+        value           uniform (0 0 0);
+    }
+    inlet
+    {
+        type            flowRateInletVelocity;
+        flowRate        0.1;
+        value           uniform (0 0 0);
+    }
+    outlet
+    {
+        type            inletOutlet;
+        value           uniform (0 0 0);
+        inletValue      uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/epsilon b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/epsilon
new file mode 100644
index 0000000000000000000000000000000000000000..d43279e3d7b721b6783cb6f10ccc0c884fde38c0
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/epsilon
@@ -0,0 +1,64 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -3 0 0 0 0 ];
+
+internalField   uniform 200;
+
+boundaryField
+{
+    front
+    {
+        type            epsilonWallFunction;
+        value           uniform 200;
+    }
+
+    back
+    {
+        type            epsilonWallFunction;
+        value           uniform 200;
+    }
+
+    wall
+    {
+        type            epsilonWallFunction;
+        value           uniform 200;
+    }
+
+    porosityWall
+    {
+        type            epsilonWallFunction;
+        value           uniform 200;
+    }
+
+    inlet
+    {
+        type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.005;
+        value           uniform 200;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 200;
+        value           uniform 200;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/k b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/k
new file mode 100644
index 0000000000000000000000000000000000000000..8a7e20d154faf52c2bf200cd7c0bf3e66f42d817
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/k
@@ -0,0 +1,64 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -2 0 0 0 0 ];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    front
+    {
+        type            kqRWallFunction;
+        value           uniform 1;
+    }
+
+    back
+    {
+        type            kqRWallFunction;
+        value           uniform 1;
+    }
+
+    wall
+    {
+        type            kqRWallFunction;
+        value           uniform 1;
+    }
+
+    porosityWall
+    {
+        type            kqRWallFunction;
+        value           uniform 1;
+    }
+
+    inlet
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.05;
+        value           uniform 1;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/nut b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/nut
new file mode 100644
index 0000000000000000000000000000000000000000..de1c9477df5f7ccb984f523c27282804fe92337f
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/nut
@@ -0,0 +1,62 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    front
+    {
+        type            nutWallFunction;
+        value           uniform 0;
+    }
+
+    back
+    {
+        type            nutWallFunction;
+        value           uniform 0;
+    }
+
+    wall
+    {
+        type            nutWallFunction;
+        value           uniform 0;
+    }
+
+    porosityWall
+    {
+        type            nutWallFunction;
+        value           uniform 0;
+    }
+
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/p b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..4794b84b189f7020d3f1bc09e5d43bbf914d7b45
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/p
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    front
+    {
+        type            zeroGradient;
+    }
+    back
+    {
+        type            zeroGradient;
+    }
+    wall
+    {
+        type            zeroGradient;
+    }
+    porosityWall
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/RASProperties b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/RASProperties
new file mode 100644
index 0000000000000000000000000000000000000000..48cf724e269cffbb00c1dc8bcf837a8a6af5543d
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/RASProperties
@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "constant";
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel        kEpsilon;
+
+turbulence      on;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..0438819b2509cf399923c3a2ce92ea0999735aea
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict
@@ -0,0 +1,123 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// block definition for a porosity with an angled inlet/outlet
+// the porosity is not aligned with the main axes
+//
+                             
+convertToMeters 0.001;
+
+vertices
+(
+    // inlet region
+    ( -150  0  -25 )  // pt 0 (in1b) 
+    ( -150 35.35533906  -25 ) // pt 1 (in2b) 
+    ( -150  0  25 )  // pt 2 (in1f) 
+    ( -150 35.35533906  25 ) // pt 3 (in2f) 
+
+    // join inlet->outlet
+    (  0 0  -25 )    // pt 4 (join1b) 
+    ( -35.35533906   35.35533906  -25 ) // pt 5 (join2b) 
+    (  0 0  25 )    // pt 6 (join1f) 
+    ( -35.35533906   35.35533906  25 ) // pt 7 (join2f) 
+
+    // porosity ends ->outlet
+    ( 70.71067812 70.71067812  -25 )  // pt 8 (poro1b) 
+    ( 35.35533906 106.06601718  -25 )  // pt 9 (poro2b) 
+    ( 70.71067812 70.71067812  25 )  // pt 10 (poro1f) 
+    ( 35.35533906 106.06601718  25 )  // pt 11 (poro2f) 
+
+    // outlet
+    ( 141.42135624 141.42135624 -25 ) // pt 12 (out1b) 
+    ( 106.06601718 176.7766953 -25 ) // pt 13 (out2b) 
+    ( 141.42135624 141.42135624 25 ) // pt 14 (out1f) 
+    ( 106.06601718 176.7766953 25 ) // pt 15 (out2f) 
+);
+
+blocks
+(
+    // inlet block
+    hex (0 4 5 1 2 6 7 3)
+    inlet ( 15 20 20 ) simpleGrading (1 1 1)
+
+    // porosity block
+    hex (4 8 9 5 6 10 11 7)
+    porosity ( 20 20 20 ) simpleGrading (1 1 1)
+
+    // outlet block
+    hex (8 12 13 9 10 14 15 11)
+    outlet ( 20 20 20 )  simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+patches
+(
+    // is there no way of defining all my 'defaultFaces' to be 'wall'?
+    wall front
+    (
+    // inlet block
+    (2 6 7 3)
+    // outlet block
+    (10 14 15 11)
+    )
+
+    wall back
+    (
+    // inlet block
+    (1 5 4 0)
+    // outlet block
+    (9 13 12 8)
+    )
+
+    wall wall
+    (
+    // inlet block
+    (2 0 4 6)
+    (7 5 1 3)
+    // outlet block
+    (10 8 12 14)
+    (15 13 9 11)
+    )
+
+    wall porosityWall
+    (
+    // porosity block
+    (6 10 11 7)
+    // porosity block
+    (5 9 8 4)
+    // porosity block
+    (6 4 8 10)
+    (11 9 5 7)
+    )
+
+    patch inlet
+    (
+    (3 1 0 2)
+    )
+
+    patch outlet
+    (
+    (15 13 12 14)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict.m4 b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict.m4
new file mode 100644
index 0000000000000000000000000000000000000000..6d6d0669392b2396160799c4747d2e9f6c2c54d2
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict.m4
@@ -0,0 +1,165 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    `format'      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// block definition for a porosity with an angled inlet/outlet
+// the porosity is not aligned with the main axes
+//
+dnl> -----------------------------------------------------------------
+dnl> <STANDARD DEFINTIONS>
+dnl>
+changecom(//)changequote([,]) dnl>
+define(calc, [esyscmd(perl -e 'print ($1)')]) dnl>
+define(VCOUNT, 0)  dnl>
+define(vlabel, [[// ]pt VCOUNT ($1) define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))])  dnl>
+dnl>
+define(hex2D, hex ($1b $2b $3b $4b $1f $2f $3f $4f)) dnl>
+define(quad2D, ($1f $1b $2b $2f))  dnl>
+define(frontQuad, ($1f $2f $3f $4f)) dnl>
+define(backQuad, ($4b $3b $2b $1b)) dnl>
+dnl>
+dnl> </STANDARD DEFINTIONS>
+dnl> -----------------------------------------------------------------
+dnl>
+define(ncells, 20) dnl>
+define(ninlet, 15) dnl>
+define(nporo, 20) dnl>
+define(noutlet, 20) dnl>
+dnl>
+define(x0,0) dnl>
+define(y0,0) dnl>
+define(y0,0) dnl>
+define(Cos,0.7071067812) dnl>   == cos(45)
+define(Sin,0.7071067812) dnl>   == sin(45)
+dnl>
+define(width,50) dnl>
+define(zBack,calc(-width/2)) dnl>
+define(zFront,calc(width/2)) dnl>
+define(leninlet,150)dnl>
+define(lenporo,100)dnl>
+define(lenoutlet,100)dnl>
+dnl>
+define(xhyp,calc(Sin*width)) dnl>
+define(yhyp,calc(Cos*width)) dnl>
+define(xinlet,leninlet)dnl>
+define(xporo,calc(Cos*lenporo)) dnl>
+define(yporo,calc(Sin*lenporo)) dnl>
+define(xoutlet,calc(xporo + Cos*lenoutlet)) dnl>
+define(youtlet,calc(yporo + Sin*lenoutlet)) dnl>
+dnl>
+
+convertToMeters 0.001;
+
+vertices
+(
+    // inlet region
+    ( -xinlet  y0  zBack )  vlabel(in1b)
+    ( -xinlet yhyp  zBack ) vlabel(in2b)
+    ( -xinlet  y0  zFront )  vlabel(in1f)
+    ( -xinlet yhyp  zFront ) vlabel(in2f)
+
+    // join inlet->outlet
+    (  x0 y0  zBack )    vlabel(join1b)
+    ( -xhyp   yhyp  zBack ) vlabel(join2b)
+    (  x0 y0  zFront )    vlabel(join1f)
+    ( -xhyp   yhyp  zFront ) vlabel(join2f)
+
+    // porosity ends ->outlet
+    ( xporo yporo  zBack )  vlabel(poro1b)
+    ( calc(xporo - xhyp) calc(yporo + yhyp)  zBack )  vlabel(poro2b)
+    ( xporo yporo  zFront )  vlabel(poro1f)
+    ( calc(xporo - xhyp) calc(yporo + yhyp)  zFront )  vlabel(poro2f)
+
+    // outlet
+    ( xoutlet youtlet zBack ) vlabel(out1b)
+    ( calc(xoutlet - xhyp) calc(youtlet + yhyp) zBack ) vlabel(out2b)
+    ( xoutlet youtlet zFront ) vlabel(out1f)
+    ( calc(xoutlet - xhyp) calc(youtlet + yhyp) zFront ) vlabel(out2f)
+);
+
+blocks
+(
+    // inlet block
+    hex2D(in1, join1, join2, in2)
+    inlet ( ninlet ncells ncells ) simpleGrading (1 1 1)
+
+    // porosity block
+    hex2D(join1, poro1, poro2, join2)
+    porosity ( nporo ncells ncells ) simpleGrading (1 1 1)
+
+    // outlet block
+    hex2D(poro1, out1, out2, poro2)
+    outlet ( noutlet ncells ncells )  simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+patches
+(
+    // is there no way of defining all my 'defaultFaces' to be 'wall'?
+    wall front
+    (
+    // inlet block
+    frontQuad(in1, join1, join2, in2)
+    // outlet block
+    frontQuad(poro1, out1, out2, poro2)
+    )
+
+    wall back
+    (
+    // inlet block
+    backQuad(in1, join1, join2, in2)
+    // outlet block
+    backQuad(poro1, out1, out2, poro2)
+    )
+
+    wall wall
+    (
+    // inlet block
+    quad2D(in1, join1)
+    quad2D(join2, in2)
+    // outlet block
+    quad2D(poro1, out1)
+    quad2D(out2, poro2)
+    )
+
+    wall porosityWall
+    (
+    // porosity block
+    frontQuad(join1, poro1, poro2, join2)
+    // porosity block
+    backQuad(join1, poro1, poro2, join2)
+    // porosity block
+    quad2D(join1, poro1)
+    quad2D(poro2, join2)
+    )
+
+    patch inlet
+    (
+    quad2D(in2, in1)
+    )
+
+    patch outlet
+    (
+    quad2D(out2, out1)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
new file mode 100644
index 0000000000000000000000000000000000000000..948cb99e4312759b00458320d02423017267fd81
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+6
+(
+    front
+    {
+        type            wall;
+        nFaces          700;
+        startFace       63400;
+    }
+    back
+    {
+        type            wall;
+        nFaces          700;
+        startFace       64100;
+    }
+    wall
+    {
+        type            wall;
+        nFaces          1400;
+        startFace       64800;
+    }
+    porosityWall
+    {
+        type            wall;
+        nFaces          1600;
+        startFace       66200;
+    }
+    inlet
+    {
+        type            patch;
+        nFaces          400;
+        startFace       67800;
+    }
+    outlet
+    {
+        type            patch;
+        nFaces          400;
+        startFace       68200;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/porousZones b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/porousZones
new file mode 100644
index 0000000000000000000000000000000000000000..634799837eaf7009528df07d027deafc7c6cc1f4
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/porousZones
@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      porousZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+1
+(
+    porosity
+    {
+        coordinateSystem
+        {
+            e1  (0.70710678 0.70710678 0);
+            e2  (0 0 1);
+        }
+
+        Darcy
+        {
+            d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
+            f   f [0 -1 0 0 0 0 0] (0 0 0);
+        }
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/transportProperties b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..7f3bedaae59eebf8d32dd99ac2f9c8ff8adf77a6
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/transportProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              nu [0 2 -1 0 0 0 0] 1.5e-05;
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/controlDict b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..88d62b98cc2648330cce6426b6043397caed194d
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/controlDict
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     rhoPorousSimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         100;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   10;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+graphFormat     raw;
+
+runTimeModifiable yes;
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSchemes b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..a70c35e1bf81ab94f439453dfac9b9349f09bcae
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSchemes
@@ -0,0 +1,64 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(U)         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+divSchemes
+{
+    div(phi,U)      Gauss upwind;
+    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div(phi,epsilon) Gauss upwind;
+    div(phi,k)      Gauss upwind;
+}
+
+laplacianSchemes
+{
+    laplacian(nuEff,U) Gauss linear corrected;
+    laplacian(rAU,p) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(1,p)  Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSolution b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..f5f9feafa8b313dea68bfc6599515ae1f6b7b6ec
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSolution
@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.05;
+        smoother        GaussSeidel;
+        cacheAgglomeration off;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    "(k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        nSweeps         2;
+        tolerance       1e-07;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nUCorrectors    2;
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    p               0.3;
+    U               0.7;
+    k               0.9;
+    epsilon         0.9;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/U b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/U
new file mode 100644
index 0000000000000000000000000000000000000000..26852141c2105c1e29071bcfc24a91c920937a02
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/U
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://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/interMixingFoam/laminar/damBreak/0/alpha1.org b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha1.org
new file mode 100644
index 0000000000000000000000000000000000000000..f0bdefa8cddb1ae29fe7f3c675625c288f9a509e
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha1.org
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha1;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+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 1;
+        value           uniform 1;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha2.org b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha2.org
new file mode 100644
index 0000000000000000000000000000000000000000..122e71dacedafedd83ea01b45380bc03d789950a
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha2.org
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha2;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+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/interMixingFoam/laminar/damBreak/0/alpha3.org b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha3.org
new file mode 100644
index 0000000000000000000000000000000000000000..f3d69ec08266a13c1c1f9ce2a80b2f1bdb6655cd
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha3.org
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha3;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+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/interMixingFoam/laminar/damBreak/0/p b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..b77bfe911838f761d5818a97593e1a4652949d22
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/p
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://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 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            buoyantPressure;
+        value           uniform 0;
+    }
+
+    rightWall
+    {
+        type            buoyantPressure;
+        value           uniform 0;
+    }
+
+    lowerWall
+    {
+        type            buoyantPressure;
+        value           uniform 0;
+    }
+
+    atmosphere
+    {
+        type            totalPressure;
+        p0              uniform 0;
+        U               U;
+        phi             phi;
+        rho             rho;
+        psi             none;
+        gamma           1;
+        value           uniform 0;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/Allclean b/tutorials/multiphase/interMixingFoam/laminar/damBreak/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..cff4dbd23191f6a4ca343bf6898b0762ed0cce6d
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/Allclean
@@ -0,0 +1,4 @@
+#!/bin/sh
+
+foamCleanTutorials cases
+rm -rf 0/alpha[1-3].gz
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/Allrun b/tutorials/multiphase/interMixingFoam/laminar/damBreak/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..15f6fafd90d19159e7c83d7192668c297d0e6798
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/Allrun
@@ -0,0 +1,10 @@
+#!/bin/sh
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+cp 0/alpha1.org 0/alpha1
+cp 0/alpha2.org 0/alpha2
+cp 0/alpha3.org 0/alpha3
+runApplication setFields
+runApplication interMixingFoam
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/dynamicMeshDict b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/dynamicMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..8fa634496563aaf871d8be56a7d7a1159758692c
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/dynamicMeshDict
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   staticFvMesh;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/g b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..27d4d324881b52b2a37d45b4cfbed46ed94d8051
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    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/interMixingFoam/laminar/damBreak/constant/polyMesh/blockMeshDict b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..a2553926d290a80076f1bf827324fd256d933b53
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/polyMesh/blockMeshDict
@@ -0,0 +1,92 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://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           
+(
+);
+
+patches         
+(
+    wall leftWall 
+    (
+        (0 12 16 4)
+        (4 16 20 8)
+    )
+    wall rightWall 
+    (
+        (7 19 15 3)
+        (11 23 19 7)
+    )
+    wall lowerWall 
+    (
+        (0 1 13 12)
+        (1 5 17 13)
+        (5 6 18 17)
+        (2 14 18 6)
+        (2 3 15 14)
+    )
+    patch atmosphere 
+    (
+        (8 20 21 9)
+        (9 21 22 10)
+        (10 22 23 11)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/polyMesh/boundary b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/polyMesh/boundary
new file mode 100644
index 0000000000000000000000000000000000000000..cfd3f8c4d431b8e5a2ab2b5035b4d4f6747f6b36
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/polyMesh/boundary
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6.x                                 |
+|   \\  /    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;
+        nFaces          4536;
+        startFace       4640;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/transportProperties b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..3cada9d4de2381a5996e8ef87393adffd7f10119
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/transportProperties
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Air
+phase1
+{
+    transportModel  Newtonian;
+    nu              nu [0 2 -1 0 0 0 0] 1.48e-05;
+    rho             rho [1 -3 0 0 0 0 0] 1;
+}
+
+// Other Liquid
+phase2
+{
+    transportModel  Newtonian;
+    nu              nu [0 2 -1 0 0 0 0]  1e-6;
+    rho             rho [1 -3 0 0 0 0 0] 1010;
+}
+
+// Water
+phase3
+{
+    transportModel  Newtonian;
+    nu              nu [0 2 -1 0 0 0 0]  1e-6;
+    rho             rho [1 -3 0 0 0 0 0] 1000;
+}
+
+// Surface tension coefficients
+sigma12           sigma12 [1 0 -2 0 0 0 0] 0.05;
+sigma13           sigma13 [1 0 -2 0 0 0 0] 0.04;
+
+// Diffusivity between miscible phases
+D23               D23   [0 2 -1 0 0 0 0]  3e-09;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/turbulenceProperties b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..9cfc50a3d927937cf2ce96c64d3c9d6bd3144e2c
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    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/interMixingFoam/laminar/damBreak/system/controlDict b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..a6f556508b1043f8de1653b692df0c37a7050d82
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/controlDict
@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     interFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         2;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.05;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression compressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           0.5;
+
+maxDeltaT       1;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/decomposeParDict b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..f1482468d421307958abf604cda46138a6feee91
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/decomposeParDict
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    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;
+}
+
+metisCoeffs
+{
+    processorWeights ( 1 1 1 1 );
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..bc9a0699584078169b5ede5cde5d4a94eac8c86d
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    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 limitedLinearV 1;
+    div(phi,alpha)  Gauss vanLeer;
+    div(phirb,alpha) Gauss interfaceCompression;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p;
+    pcorr;
+    "alpha.";
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSolution b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..f2543e14e1d63795f0d15e7deb675f27b1affb02
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSolution
@@ -0,0 +1,73 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "alpha."
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       1e-06;
+        relTol          0;
+        nSweeps         1;
+    }
+
+    pcorr
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-10;
+        relTol          0;
+    }
+
+    p
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-07;
+        relTol          0.05;
+    }
+
+    pFinal
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-07;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-06;
+        relTol          0;
+    }
+}
+
+PISO
+{
+    momentumPredictor no;
+    nCorrectors     3;
+    nNonOrthogonalCorrectors 0;
+    nAlphaCorr      1;
+    nAlphaSubCycles 2;
+    cAlpha          1;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/setFieldsDict b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/setFieldsDict
new file mode 100644
index 0000000000000000000000000000000000000000..51ad3170f569e195dd9a8613a81f53657a277b86
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/setFieldsDict
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha1 0
+    volScalarFieldValue alpha2 1
+    volScalarFieldValue alpha3 0
+);
+
+regions
+(
+    boxToCell
+    {
+        box (0 0 -1) (0.1461 0.292 1);
+        fieldValues
+        (
+            volScalarFieldValue alpha1 0
+            volScalarFieldValue alpha2 0
+            volScalarFieldValue alpha3 1
+        );
+    }
+    boxToCell
+    {
+        box (0.1461 0.05 -1) (1 1 1);
+        fieldValues
+        (
+            volScalarFieldValue alpha1 1
+            volScalarFieldValue alpha2 0
+            volScalarFieldValue alpha3 0
+        );
+    }
+);
+
+
+// ************************************************************************* //