diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
index 8dbb489c8d2d0fad19b6e8035a64b71be2767152..fa56de00bb40474dc995a83b92f6a8e2a7997894 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
@@ -1,5 +1,4 @@
 label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-
 label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
 if (nAlphaSubCycles > 1)
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index 1c42c7b91a64104651b4b2948ed7ee8b4bc823be..9c8c2d4f18c2bbdf4267a19af53d7d2ecb690f37 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -132,6 +132,17 @@
         phasei++;
     }
 
+    phasei = 0;
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+
+        mrfZones.relativeFlux(phase.phi().oldTime());
+        mrfZones.relativeFlux(phase.phi());
+
+        phasei++;
+    }
+
     surfaceScalarField Dp
     (
         IOobject
@@ -185,7 +196,6 @@
                 phase.phi() =
                     phiHbyAs[phasei]
                   + rAlphaAUfs[phasei]*mSfGradp/phase.rho();
-                mrfZones.relativeFlux(phase.phi().oldTime());
                 phi += alphafs[phasei]*phase.phi();
 
                 phasei++;
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean b/applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean
index 63db39ff0566cf06b8259c8c5040c97d50c172e5..cc138bc068e6882e24eb995668e00a0fc18d3e58 100755
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean
@@ -6,6 +6,5 @@ wclean libso phaseModel
 wclean libso interfacialModels
 wclean libso kineticTheoryModels
 wclean
-wclean MRFtwoPhaseEulerFoam
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake b/applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake
index faf438d0bddd640f77b0ad233a16266f3d8bf80e..29294d166a947be8f5f391168d9e4f67faeb8718 100755
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake
@@ -6,6 +6,5 @@ wmake libso phaseModel
 wmake libso interfacialModels
 wmake libso kineticTheoryModels
 wmake
-wmake MRFtwoPhaseEulerFoam
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C
deleted file mode 100644
index 357c7e413ab889f1a6845ceeaec7fdf46f655c93..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C
+++ /dev/null
@@ -1,119 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    twoPhaseEulerFoam
-
-Description
-    Solver for a system of 2 incompressible fluid phases with one phase
-    dispersed, e.g. gas bubbles in a liquid or solid particles in a gas.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "nearWallDist.H"
-#include "wallFvPatch.H"
-#include "Switch.H"
-
-#include "IFstream.H"
-#include "OFstream.H"
-
-#include "dragModel.H"
-#include "phaseModel.H"
-#include "kineticTheoryModel.H"
-
-#include "pimpleControl.H"
-#include "MRFZones.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "setRootCase.H"
-
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "readGravitationalAcceleration.H"
-    #include "createFields.H"
-    #include "readPPProperties.H"
-    #include "initContinuityErrs.H"
-    #include "createMRFZones.H"
-    #include "readTimeControls.H"
-    #include "CourantNo.H"
-    #include "setInitialDeltaT.H"
-
-    pimpleControl pimple(mesh);
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.run())
-    {
-        #include "readTwoPhaseEulerFoamControls.H"
-        #include "CourantNos.H"
-        #include "setDeltaT.H"
-
-        runTime++;
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            #include "alphaEqn.H"
-            #include "liftDragCoeffs.H"
-            #include "UEqns.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                #include "pEqn.H"
-
-                if (correctAlpha && !pimple.finalIter())
-                {
-                    #include "alphaEqn.H"
-                }
-            }
-
-            #include "DDtU.H"
-
-            if (pimple.turbCorr())
-            {
-                #include "kEpsilon.H"
-            }
-        }
-
-        #include "write.H"
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files
deleted file mode 100644
index 45960722aefce40efa346d1723b98b35c5651da0..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-MRFTwoPhaseEulerFoam.C
-
-EXE = $(FOAM_APPBIN)/MRFTwoPhaseEulerFoam
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options
deleted file mode 100644
index b9b19059da672e0822af3f61d0314ab3988ebd39..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options
+++ /dev/null
@@ -1,17 +0,0 @@
-EXE_INC = \
-    -I.. \
-    -I../../bubbleFoam \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-    -IturbulenceModel \
-    -I../kineticTheoryModels/lnInclude \
-    -I../interfacialModels/lnInclude \
-    -I../phaseModel/lnInclude \
-
-EXE_LIBS = \
-    -lEulerianInterfacialModels \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lincompressibleTransportModels \
-    -lphaseModel \
-    -lkineticTheoryModel
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H
deleted file mode 100644
index 0c0cc1543a7a5aad203013d0ed97b15505121f36..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H
+++ /dev/null
@@ -1,99 +0,0 @@
-fvVectorMatrix U1Eqn(U1, U1.dimensions()*dimVol/dimTime);
-fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
-
-{
-    {
-        volTensorField gradU1T(T(fvc::grad(U1)));
-
-        if (kineticTheory.on())
-        {
-            kineticTheory.solve(gradU1T);
-            nuEff1 = kineticTheory.mu1()/rho1;
-        }
-        else // If not using kinetic theory is using Ct model
-        {
-            nuEff1 = sqr(Ct)*nut2 + nu1;
-        }
-
-        volTensorField Rc1
-        (
-            "Rc1",
-            (((2.0/3.0)*I)*nuEff1)*tr(gradU1T) - nuEff1*gradU1T
-        );
-
-        if (kineticTheory.on())
-        {
-            Rc1 -= ((kineticTheory.lambda()/rho1)*tr(gradU1T))*tensor(I);
-        }
-
-        surfaceScalarField phiR1
-        (
-            -fvc::interpolate(nuEff1)*mesh.magSf()*fvc::snGrad(alpha1)
-            /fvc::interpolate(alpha1 + scalar(0.001))
-        );
-
-        U1Eqn =
-        (
-            (scalar(1) + Cvm*rho2*alpha2/rho1)*
-            (
-                fvm::ddt(U1)
-              + fvm::div(phi1, U1, "div(phi1,U1)")
-              - fvm::Sp(fvc::div(phi1), U1)
-            )
-
-          - fvm::laplacian(nuEff1, U1)
-          + fvc::div(Rc1)
-
-          + fvm::div(phiR1, U1, "div(phi1,U1)")
-          - fvm::Sp(fvc::div(phiR1), U1)
-          + (fvc::grad(alpha1)/(fvc::average(alpha1) + scalar(0.001)) & Rc1)
-         ==
-        //  g                          // Buoyancy term transfered to p-equation
-          - fvm::Sp(alpha2/rho1*K, U1)
-        //+ alpha2/rho1*K*U2           // Explicit drag transfered to p-equation
-          - alpha2/rho1*(liftCoeff - Cvm*rho2*DDtU2)
-        );
-        mrfZones.addCoriolis(U1Eqn);
-        U1Eqn.relax();
-    }
-
-    {
-        volTensorField gradU2T(T(fvc::grad(U2)));
-        volTensorField Rc2
-        (
-            "Rc2",
-            (((2.0/3.0)*I)*nuEff2)*tr(gradU2T) - nuEff2*gradU2T
-        );
-
-        surfaceScalarField phiR2
-        (
-            -fvc::interpolate(nuEff2)*mesh.magSf()*fvc::snGrad(alpha2)
-            /fvc::interpolate(alpha2 + scalar(0.001))
-        );
-
-        U2Eqn =
-        (
-            (scalar(1) + Cvm*rho2*alpha1/rho2)*
-            (
-                fvm::ddt(U2)
-              + fvm::div(phi2, U2, "div(phi2,U2)")
-              - fvm::Sp(fvc::div(phi2), U2)
-            )
-
-          - fvm::laplacian(nuEff2, U2)
-          + fvc::div(Rc2)
-
-          + fvm::div(phiR2, U2, "div(phi2,U2)")
-          - fvm::Sp(fvc::div(phiR2), U2)
-
-          + (fvc::grad(alpha2)/(fvc::average(alpha2) + scalar(0.001)) & Rc2)
-         ==
-        //  g                          // Buoyancy term transfered to p-equation
-          - fvm::Sp(alpha1/rho2*K, U2)
-        //+ alpha1/rho2*K*U1           // Explicit drag transfered to p-equation
-          + alpha1/rho2*(liftCoeff + Cvm*rho2*DDtU1)
-        );
-        mrfZones.addCoriolis(U2Eqn);
-        U2Eqn.relax();
-    }
-}
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H
deleted file mode 100644
index 348cc847d71411e7e4c2c4ada254244c01ef59c6..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H
+++ /dev/null
@@ -1,118 +0,0 @@
-{
-    surfaceScalarField alpha1f(fvc::interpolate(alpha1));
-    surfaceScalarField alpha2f(scalar(1) - alpha1f);
-
-    volScalarField rAU1(1.0/U1Eqn.A());
-    volScalarField rAU2(1.0/U2Eqn.A());
-
-    rAU1f = fvc::interpolate(rAU1);
-    surfaceScalarField rAU2f(fvc::interpolate(rAU2));
-
-    volVectorField HbyA1("HbyA1", U1);
-    HbyA1 = rAU1*U1Eqn.H();
-
-    volVectorField HbyA2("HbyA2", U2);
-    HbyA2 = rAU2*U2Eqn.H();
-
-    mrfZones.absoluteFlux(phi1.oldTime());
-    mrfZones.absoluteFlux(phi1);
-
-    mrfZones.absoluteFlux(phi2.oldTime());
-    mrfZones.absoluteFlux(phi2);
-
-    surfaceScalarField phiDrag1
-    (
-        fvc::interpolate(alpha2/rho1*K*rAU1)*phi2 + rAU1f*(g & mesh.Sf())
-    );
-
-    if (g0.value() > 0.0)
-    {
-        phiDrag1 -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
-    }
-
-    if (kineticTheory.on())
-    {
-        phiDrag1 -= rAU1f*fvc::snGrad(kineticTheory.pa()/rho1)*mesh.magSf();
-    }
-
-
-    surfaceScalarField phiDrag2
-    (
-        fvc::interpolate(alpha1/rho2*K*rAU2)*phi1 + rAU2f*(g & mesh.Sf())
-    );
-
-    // Fix for gravity on outlet boundary.
-    forAll(p.boundaryField(), patchi)
-    {
-        if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
-        {
-            phiDrag1.boundaryField()[patchi] = 0.0;
-            phiDrag2.boundaryField()[patchi] = 0.0;
-        }
-    }
-
-    surfaceScalarField phiHbyA1
-    (
-        "phiHbyA1",
-        (fvc::interpolate(HbyA1) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU1, U1, phi1)
-      + phiDrag1
-    );
-    mrfZones.relativeFlux(phiHbyA1);
-
-    surfaceScalarField phiHbyA2
-    (
-        "phiHbyA2",
-        (fvc::interpolate(HbyA2) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU2, U2, phi2)
-      + phiDrag2
-    );
-    mrfZones.relativeFlux(phiHbyA2);
-
-    mrfZones.relativeFlux(phi1.oldTime());
-    mrfZones.relativeFlux(phi1);
-    mrfZones.relativeFlux(phi2.oldTime());
-    mrfZones.relativeFlux(phi2);
-
-    surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
-
-    surfaceScalarField Dp
-    (
-        "Dp",
-        alpha1f*rAU1f/rho1 + alpha2f*rAU2f/rho2
-    );
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
-        );
-
-        pEqn.setReference(pRefCell, pRefValue);
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            surfaceScalarField SfGradp(pEqn.flux()/Dp);
-
-            phi1 = phiHbyA1 - rAU1f*SfGradp/rho1;
-            phi2 = phiHbyA2 - rAU2f*SfGradp/rho2;
-            phi = alpha1f*phi1 + alpha2f*phi2;
-
-            p.relax();
-            SfGradp = pEqn.flux()/Dp;
-
-            U1 = HbyA1 + fvc::reconstruct(phiDrag1 - rAU1f*SfGradp/rho1);
-            U1.correctBoundaryConditions();
-
-            U2 = HbyA2 + fvc::reconstruct(phiDrag2 - rAU2f*SfGradp/rho2);
-            U2.correctBoundaryConditions();
-
-            U = alpha1*U1 + alpha2*U2;
-        }
-    }
-}
-
-#include "continuityErrs.H"
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
index 8d09ccd3ef0c74db4cc9ef096eca61d0747ddb77..0c0cc1543a7a5aad203013d0ed97b15505121f36 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
@@ -53,7 +53,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         //+ alpha2/rho1*K*U2           // Explicit drag transfered to p-equation
           - alpha2/rho1*(liftCoeff - Cvm*rho2*DDtU2)
         );
-
+        mrfZones.addCoriolis(U1Eqn);
         U1Eqn.relax();
     }
 
@@ -93,7 +93,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         //+ alpha1/rho2*K*U1           // Explicit drag transfered to p-equation
           + alpha1/rho2*(liftCoeff + Cvm*rho2*DDtU1)
         );
-
+        mrfZones.addCoriolis(U2Eqn);
         U2Eqn.relax();
     }
 }
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
index 52f09e4793943737b24ed0f70805bdcce5ce5adc..34fdc24afc532d56c0fd0ddf0c08479237d2cdb9 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
@@ -1,6 +1,9 @@
 {
-    word scheme("div(phi,alpha1)");
-    word schemer("div(phir,alpha1)");
+    label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
+    label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
+
+    word alphaScheme("div(phi,alpha1)");
+    word alpharScheme("div(phir,alpha1)");
 
     surfaceScalarField phic("phic", phi);
     surfaceScalarField phir("phir", phi1 - phi2);
@@ -15,14 +18,39 @@
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
     {
-        fvScalarMatrix alpha1Eqn
+        for
         (
-             fvm::ddt(alpha1)
-           + fvm::div(phic, alpha1, scheme)
-           + fvm::div(-fvc::flux(-phir, alpha2, schemer), alpha1, schemer)
-        );
+            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
+        )
+        {
+            surfaceScalarField phiAlpha
+            (
+                fvc::flux
+                (
+                    phic,
+                    alpha1,
+                    alphaScheme
+                )
+              + fvc::flux
+                (
+                    -fvc::flux(-phir, alpha2, alpharScheme),
+                    alpha1,
+                    alpharScheme
+                )
+            );
+
+            MULES::explicitSolve
+            (
+                alpha1,
+                phi,
+                phiAlpha,
+                (g0.value() > 0 ? alphaMax : 1),
+                0
+            );
+        }
 
-        if (g0.value() > 0.0)
+        if (g0.value() > 0)
         {
             ppMagf = rAU1f*fvc::interpolate
             (
@@ -30,18 +58,22 @@
                *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
             );
 
-            alpha1Eqn -= fvm::laplacian
+            fvScalarMatrix alpha1Eqn
             (
-                (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
-                alpha1,
-                "laplacian(alpha1PpMag,alpha1)"
+                fvm::ddt(alpha1) - fvc::ddt(alpha1)
+              - fvm::laplacian
+                (
+                    (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
+                    alpha1,
+                    "laplacian(alpha1PpMag,alpha1)"
+                )
             );
-        }
 
-        alpha1Eqn.relax();
-        alpha1Eqn.solve();
+            alpha1Eqn.relax();
+            alpha1Eqn.solve();
 
-        #include "packingLimiter.H"
+            #include "packingLimiter.H"
+        }
 
         alpha2 = scalar(1) - alpha1;
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
index b71c7849e5b456591f19ef931ae4f102b58fc54d..314d9b55ea6d8a8edc0b90eb1ea11a2fad806df2 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
@@ -93,16 +93,22 @@
 
     dimensionedScalar Cvm
     (
+        "Cvm",
+        dimless,
         transportProperties.lookup("Cvm")
     );
 
     dimensionedScalar Cl
     (
+        "Cl",
+        dimless,
         transportProperties.lookup("Cl")
     );
 
     dimensionedScalar Ct
     (
+        "Ct",
+        dimless,
         transportProperties.lookup("Ct")
     );
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createMRFZones.H
similarity index 100%
rename from applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H
rename to applications/solvers/multiphase/twoPhaseEulerFoam/createMRFZones.H
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 91f8302d7ff1a3c27bf971068c89815a8a94db7e..348cc847d71411e7e4c2c4ada254244c01ef59c6 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -14,6 +14,12 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
+    mrfZones.absoluteFlux(phi1.oldTime());
+    mrfZones.absoluteFlux(phi1);
+
+    mrfZones.absoluteFlux(phi2.oldTime());
+    mrfZones.absoluteFlux(phi2);
+
     surfaceScalarField phiDrag1
     (
         fvc::interpolate(alpha2/rho1*K*rAU1)*phi2 + rAU1f*(g & mesh.Sf())
@@ -29,6 +35,7 @@
         phiDrag1 -= rAU1f*fvc::snGrad(kineticTheory.pa()/rho1)*mesh.magSf();
     }
 
+
     surfaceScalarField phiDrag2
     (
         fvc::interpolate(alpha1/rho2*K*rAU2)*phi1 + rAU2f*(g & mesh.Sf())
@@ -51,6 +58,7 @@
       + fvc::ddtPhiCorr(rAU1, U1, phi1)
       + phiDrag1
     );
+    mrfZones.relativeFlux(phiHbyA1);
 
     surfaceScalarField phiHbyA2
     (
@@ -59,6 +67,12 @@
       + fvc::ddtPhiCorr(rAU2, U2, phi2)
       + phiDrag2
     );
+    mrfZones.relativeFlux(phiHbyA2);
+
+    mrfZones.relativeFlux(phi1.oldTime());
+    mrfZones.relativeFlux(phi1);
+    mrfZones.relativeFlux(phi2.oldTime());
+    mrfZones.relativeFlux(phi2);
 
     surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
 
@@ -83,14 +97,8 @@
         {
             surfaceScalarField SfGradp(pEqn.flux()/Dp);
 
-            phi1.boundaryField() ==
-                (fvc::interpolate(U1) & mesh.Sf())().boundaryField();
             phi1 = phiHbyA1 - rAU1f*SfGradp/rho1;
-
-            phi2.boundaryField() ==
-                (fvc::interpolate(U2) & mesh.Sf())().boundaryField();
             phi2 = phiHbyA2 - rAU2f*SfGradp/rho2;
-
             phi = alpha1f*phi1 + alpha2f*phi2;
 
             p.relax();
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C
index d81ad8e703143ef3916e25ea55f5bfb7a69a6b2b..aebc7ba6cb01645ed9f638c4cb80501bfe8b0f3b 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C
@@ -43,14 +43,20 @@ Foam::phaseModel::phaseModel
     name_(phaseName),
     d_
     (
+        "d",
+        dimLength,
         dict_.lookup("d")
     ),
     nu_
     (
+        "nu",
+        dimensionSet(0, 2, -1, 0, 0),
         dict_.lookup("nu")
     ),
     rho_
     (
+        "rho",
+        dimDensity,
         dict_.lookup("rho")
     ),
     U_
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
index a345c4e53bc97e5b3ad14ddaf653805a78a75571..4e6847debe35e252aa3b696c93bb43d5f79d1213 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
@@ -1,5 +1,3 @@
     #include "readTimeControls.H"
 
-    int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
-
     Switch correctAlpha(pimple.dict().lookup("correctAlpha"));
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index 19810a15f19c7db8839497c21dd12b054655af28..d9da532aef570a4e206b7680531c72e7fff571f8 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "MULES.H"
+#include "subCycle.H"
 #include "nearWallDist.H"
 #include "wallFvPatch.H"
 #include "Switch.H"
@@ -43,6 +45,7 @@ Description
 #include "kineticTheoryModel.H"
 
 #include "pimpleControl.H"
+#include "MRFZones.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,6 +59,7 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "readPPProperties.H"
     #include "initContinuityErrs.H"
+    #include "createMRFZones.H"
     #include "readTimeControls.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
diff --git a/src/OpenFOAM/algorithms/subCycle/subCycle.H b/src/OpenFOAM/algorithms/subCycle/subCycle.H
index 7ead8235051bb033ffe23888798dfb20bcf0603a..a4a1ea1844c506b86ebeeba260e9895005405d70 100644
--- a/src/OpenFOAM/algorithms/subCycle/subCycle.H
+++ b/src/OpenFOAM/algorithms/subCycle/subCycle.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,13 +40,11 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class subCycle Declaration
+                          Class subCycleField Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class GeometricField>
-class subCycle
-:
-    public subCycleTime
+class subCycleField
 {
     // Private data
 
@@ -57,6 +55,41 @@ class subCycle
         GeometricField  gf0_;
 
 
+public:
+
+    // Constructors
+
+        //- Construct field and number of sub-cycles
+        subCycleField(GeometricField& gf)
+        :
+            gf_(gf),
+            gf0_(gf.oldTime())
+        {}
+
+
+    //- Destructor
+    ~subCycleField()
+    {
+        // Correct the time index of the field to correspond to the global time
+        gf_.timeIndex() = gf_.time().timeIndex();
+
+        // Reset the old-time field value
+        gf_.oldTime() = gf0_;
+        gf_.oldTime().timeIndex() = gf0_.timeIndex();
+    }
+};
+
+
+/*---------------------------------------------------------------------------*\
+                          Class subCycle Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class GeometricField>
+class subCycle
+:
+    public subCycleField<GeometricField>,
+    public subCycleTime
+{
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -73,9 +106,9 @@ public:
         //- Construct field and number of sub-cycles
         subCycle(GeometricField& gf, const label nSubCycles)
         :
-            subCycleTime(const_cast<Time&>(gf.time()), nSubCycles),
-            gf_(gf),
-            gf0_(gf.oldTime())
+
+            subCycleField<GeometricField>(gf),
+            subCycleTime(const_cast<Time&>(gf.time()), nSubCycles)
         {}
 
 
@@ -84,13 +117,6 @@ public:
     {
         // End the subCycleTime, which restores the time state
         endSubCycle();
-
-        // Correct the time index of the field to correspond to the global time
-        gf_.timeIndex() = gf_.time().timeIndex();
-
-        // Reset the old-time field value
-        gf_.oldTime() = gf0_;
-        gf_.oldTime().timeIndex() = gf_.time().timeIndex();
     }
 };
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed/constant/MRFZones b/tutorials/multiphase/twoPhaseEulerFoam/bed/constant/MRFZones
new file mode 100644
index 0000000000000000000000000000000000000000..0034a58c2293992131cec848f2dcb0f8ad50a0d4
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed/constant/MRFZones
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      MRFZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+0()
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed/constant/transportProperties b/tutorials/multiphase/twoPhaseEulerFoam/bed/constant/transportProperties
index ef9a0fccf7ebf3a1cae087fd9ceab3414a7f50b2..6d05c03c615194a5d93c828ac477e6e029df5161 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bed/constant/transportProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed/constant/transportProperties
@@ -17,25 +17,22 @@ FoamFile
 
 phase1
 {
-    rho             rho [ 1 -3 0 0 0 ] 2640;
-    nu              nu [ 0 2 -1 0 0 ] 1e-06;
-    d               d [ 0 1 0 0 0 0 0 ] 0.00048;
+    rho             2640;
+    nu              1e-06;
+    d               0.00048;
 }
 
 phase2
 {
-    rho             rho [ 1 -3 0 0 0 ] 1.28;
-    nu              nu [ 0 2 -1 0 0 ] 1.328e-05;
-    d               d [ 0 1 0 0 0 0 0 ] 1;
+    rho             1.28;
+    nu              1.328e-05;
+    d               1;
 }
 
-Cvm             Cvm [ 0 0 0 0 0 ] 0;
+Cvm             0;
 
-Cl              Cl [ 0 0 0 0 0 ] 0;
-
-Ct              Ct [ 0 0 0 0 0 ] 0;
-
-alpha1Alpha      alpha1Alpha [ 0 0 0 0 0 ] 0;
+Cl              0;
 
+Ct              0;
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/bed/system/fvSolution
index bdb02666e34b4bd55546ef3a66204ff11f3f10e5..5800e2a4b4dc1f19677781a8a5487b289bb2296e 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bed/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed/system/fvSolution
@@ -49,25 +49,21 @@ solvers
 
     "(k|epsilon)Final"
     {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-05;
+        $k;
         relTol          0;
     }
 
     alpha
     {
-        solver          PBiCG;
-        preconditioner  DILU;
+        solver          PCG;
+        preconditioner  DIC;
         tolerance       1e-10;
         relTol          0.1;
     }
 
     alpha1Final
     {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-10;
+        $alpha;
         relTol          0;
     }
 }
@@ -77,6 +73,7 @@ PIMPLE
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
     nAlphaCorr      2;
+    nAlphaSubCycles 3;
     correctAlpha    yes;
     pRefCell        0;
     pRefValue       0;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/MRFZones b/tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/MRFZones
new file mode 100644
index 0000000000000000000000000000000000000000..0034a58c2293992131cec848f2dcb0f8ad50a0d4
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/MRFZones
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      MRFZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+0()
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/transportProperties b/tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/transportProperties
index 859bd94263b9a0feb42e07c489ff4d1a48b79beb..1afd6e92c52bc8aa1afe3bf34f2b086f99c886a0 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/transportProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/transportProperties
@@ -17,25 +17,20 @@ FoamFile
 
 phase1
 {
-    rho             rho [ 1 -3 0 0 0 ] 2500;
-    nu              nu [ 0 2 -1 0 0 ] 1e-06;
-    d               d [ 0 1 0 0 0 0 0 ] 0.0003;
+    rho             2500;
+    nu              1e-06;
+    d               0.0003;
 }
 
 phase2
 {
-    rho             rho [ 1 -3 0 0 0 ] 1.2;
-    nu              nu [ 0 2 -1 0 0 ] 1.5e-05;
-    d               d [ 0 1 0 0 0 0 0 ] 1;
+    rho             1.2;
+    nu              1.5e-05;
+    d               1;
 }
 
-Cvm             Cvm [ 0 0 0 0 0 ] 0;
-
-Cl              Cl [ 0 0 0 0 0 ] 0;
-
-Ct              Ct [ 0 0 0 0 0 ] 0;
-
-alpha1Alpha      alpha1Alpha [ 0 0 0 0 0 ] 0;
-
+Cvm             0;
+Cl              0;
+Ct              0;
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed2/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/bed2/system/fvSolution
index bdb02666e34b4bd55546ef3a66204ff11f3f10e5..2e3cdfa4aaee926062296b0c2d40cb381f41ebcb 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bed2/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed2/system/fvSolution
@@ -76,7 +76,8 @@ PIMPLE
 {
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
-    nAlphaCorr      2;
+    nAlphaCorr      1;
+    nAlphaSubCycles 3;
     correctAlpha    yes;
     pRefCell        0;
     pRefValue       0;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/MRFZones b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/MRFZones
new file mode 100644
index 0000000000000000000000000000000000000000..0034a58c2293992131cec848f2dcb0f8ad50a0d4
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/MRFZones
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      MRFZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+0()
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/transportProperties b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/transportProperties
index 276c521ad3c505013c1cddb3b71b1d7074dbca50..e013730cb5cce44d7638b78954289a6b2a79d261 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/transportProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/transportProperties
@@ -17,23 +17,21 @@ FoamFile
 
 phase1
 {
-    nu              nu [ 0 2 -1 0 0 0 0 ] 1.6e-05;
-    rho             rho [ 1 -3 0 0 0 0 0 ] 1;
-    d               d [ 0 1 0 0 0 0 0 ] 0.003;
+    nu              1.6e-05;
+    rho             1;
+    d               0.003;
 }
 
 phase2
 {
-    nu              nu [ 0 2 -1 0 0 0 0 ] 1e-06;
-    rho             rho [ 1 -3 0 0 0 0 0 ] 1000;
-    d               d [ 0 1 0 0 0 0 0 ] 0.0001;
+    nu              1e-06;
+    rho             1000;
+    d               0.0001;
 }
 
-Cvm             Cvm [ 0 0 0 0 0 0 0 ] 0.5;
-
-Cl              Cl [ 0 0 0 0 0 0 0 ] 0;
-
-Ct              Ct [ 0 0 0 0 0 0 0 ] 1;
+Cvm             0.5;
+Cl              0;
+Ct              1;
 
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/system/fvSolution
index d0491f6d65a4f3c8b19e4823866f87a171501921..76a95c8327cadc21206927893a86f255c8cd1c8c 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/system/fvSolution
@@ -76,7 +76,8 @@ PIMPLE
 {
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
-    nAlphaCorr      2;
+    nAlphaCorr      1;
+    nAlphaSubCycles 2;
     correctAlpha    yes;
     pRefCell        0;
     pRefValue       0;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/Theta b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/Theta
new file mode 100644
index 0000000000000000000000000000000000000000..b18d202df894c0276dd7f546fd3820cd154623e4
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/Theta
@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Theta;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    rotor
+    {
+        type            zeroGradient;
+    }
+
+    stator
+    {
+        type            zeroGradient;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/U1 b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/U1
new file mode 100644
index 0000000000000000000000000000000000000000..7a88b384782d81444c1dbc8643f955a4720536da
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/U1
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U1;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    rotor
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    stator
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/U2 b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/U2
new file mode 100644
index 0000000000000000000000000000000000000000..b93b8d870de47ad7c0f257b0301895c8c8ed05d9
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/U2
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U2;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    rotor
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    stator
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/alpha1 b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/alpha1
new file mode 100644
index 0000000000000000000000000000000000000000..637ebed0e893d74b6cff62103ebf3c627ec682b9
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/alpha1
@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha1;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.5;
+
+boundaryField
+{
+    rotor
+    {
+        type            zeroGradient;
+    }
+
+    stator
+    {
+        type            zeroGradient;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/epsilon b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/epsilon
new file mode 100644
index 0000000000000000000000000000000000000000..3f507afed8ea6ba738e2fc7772a1fe96b57609de
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/epsilon
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -3 0 0 0 0 ];
+
+internalField   uniform 20;
+
+boundaryField
+{
+    rotor
+    {
+        type            zeroGradient;
+        value           uniform 20;
+    }
+
+    stator
+    {
+        type            zeroGradient;
+        value           uniform 20;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/k b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/k
new file mode 100644
index 0000000000000000000000000000000000000000..4704b49c838eb9f1f6c32dad5d2440ad7e7989f6
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/k
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -2 0 0 0 0 ];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    rotor
+    {
+        type            zeroGradient;
+        value           uniform 0;
+    }
+
+    stator
+    {
+        type            zeroGradient;
+        value           uniform 0;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/p b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..9c4de0da512d78d3fcd7499ef5621c328e3d56e0
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/0/p
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    rotor
+    {
+        type            multiphaseFixedFluxPressure;
+        value           $internalField;
+    }
+
+    stator
+    {
+        type            multiphaseFixedFluxPressure;
+        value           $internalField;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/Allrun b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..a718eaf7e662bad616179bd183caa9fd72de8b79
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/Allrun
@@ -0,0 +1,12 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+application=`getApplication`
+
+runApplication ./makeMesh
+runApplication $application
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/MRFZones b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/MRFZones
new file mode 100644
index 0000000000000000000000000000000000000000..25c3311d0b7ff552c41b5844854959e33685d807
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/MRFZones
@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      MRFZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+1
+(
+    rotor
+    {
+        // Fixed patches (by default they 'move' with the MRF zone)
+        nonRotatingPatches ();
+
+        origin    origin [0 1 0 0 0 0 0]  (0 0 0);
+        axis      axis   [0 0 0 0 0 0 0]  (0 0 1);
+        omega     omega  [0 0 -1 0 0 0 0] 10.472;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/RASProperties b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/RASProperties
new file mode 100644
index 0000000000000000000000000000000000000000..41b54318fe41ab6593fd868ca5c080769a1988c1
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/RASProperties
@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel        kEpsilon;
+
+turbulence      off;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/g b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..508d65849430f8e5abf4b12d7baa53d70521a1c3
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 0 0);
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/interfacialProperties b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/interfacialProperties
new file mode 100644
index 0000000000000000000000000000000000000000..41159fe9def1d68030d4d4c65aacca7947e0610b
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/interfacialProperties
@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      interfacialProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dragModel1          SchillerNaumann;
+dragModel2          SchillerNaumann;
+
+heatTransferModel1  RanzMarshall;
+heatTransferModel2  RanzMarshall;
+
+dispersedPhase      both;
+dragPhase           blended;
+
+residualSlip        1e-2;
+minInterfaceAlpha   1e-3;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/kineticTheoryProperties b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/kineticTheoryProperties
new file mode 100644
index 0000000000000000000000000000000000000000..53521fbae309deb8893eed645ceaf0776fb9786a
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/kineticTheoryProperties
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      kineticTheoryProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+kineticTheory   off;
+
+equilibrium     on;
+
+e               e [ 0 0 0 0 0 0 0 ] 0.9;
+
+alphaMax        alphaMax [ 0 0 0 0 0 0 0 ] 0.6;
+
+alphaMinFriction alphaMinFriction [ 0 0 0 0 0 0 0 ] 0.5;
+
+Fr              Fr [ 1 -1 -2 0 0 0 0 ] 0.05;
+
+eta             eta [ 0 0 0 0 0 0 0 ] 2;
+
+p               p [ 0 0 0 0 0 0 0 ] 5;
+
+phi             phi [ 0 0 0 0 0 0 0 ] 25;
+
+viscosityModel  Syamlal;
+
+conductivityModel HrenyaSinclair;
+
+granularPressureModel Lun;
+
+frictionalStressModel JohnsonJackson;
+
+radialModel     Gidaspow;
+
+HrenyaSinclairCoeffs
+{
+    L               L [ 0 1 0 0 0 0 0 ] 0.0005;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/polyMesh/blockMeshDict.m4 b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/polyMesh/blockMeshDict.m4
new file mode 100644
index 0000000000000000000000000000000000000000..a93868498ba68d11b25b0875ff69205f62f822ed
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/polyMesh/blockMeshDict.m4
@@ -0,0 +1,818 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    `format'      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// General macros to create 2D/extruded-2D meshes
+
+changecom(//)changequote([,])
+define(calc, [esyscmd(perl -e 'print ($1)')])
+define(VCOUNT, 0)
+define(vlabel, [[// ]Vertex $1 = VCOUNT define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))])
+define(pi, 3.14159265)
+
+define(hex2D, hex ($1b $2b $3b $4b $1t $2t $3t $4t))
+define(quad2D, ($1b $2b $2t $1t))
+define(frontQuad, ($1t $2t $3t $4t))
+define(backQuad, ($1b $4b $3b $2b))
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.1;
+
+// Hub radius
+define(r, 0.2)
+
+// Impeller-tip radius
+define(rb, 0.5)
+
+// Baffle-tip radius
+define(Rb, 0.7)
+
+// Tank radius
+define(R, 1)
+
+// MRF region radius
+define(ri, calc(0.5*(rb + Rb)))
+
+// Thickness of 2D slab
+define(z, 0.1)
+
+// Base z
+define(Zb, 0)
+
+// Top z
+define(Zt, calc(Zb + z))
+
+// Number of cells radially between hub and impeller tip
+define(Nr, 12)
+
+// Number of cells radially in each of the two regions between
+// impeller and baffle tips
+define(Ni, 4)
+
+// Number of cells radially between baffle tip and tank
+define(NR, 12)
+
+// Number of cells azimuthally in each of the 8 blocks
+define(Na, 12)
+
+// Number of cells in the thickness of the slab
+define(Nz, 1)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+define(vert, (x$1$2 y$1$2 $3))
+define(evert, (ex$1$2 ey$1$2 $3))
+
+define(a0, 0)
+define(a1, -45)
+define(a2, -90)
+define(a3, -135)
+define(a4, 180)
+define(a5, 135)
+define(a6, 90)
+define(a7, 45)
+
+define(ea0, -22.5)
+define(ea1, -67.5)
+define(ea2, -112.5)
+define(ea3, -157.5)
+define(ea4, 157.5)
+define(ea5, 112.5)
+define(ea6, 67.5)
+define(ea7, 22.5)
+
+define(ca0, calc(cos((pi/180)*a0)))
+define(ca1, calc(cos((pi/180)*a1)))
+define(ca2, calc(cos((pi/180)*a2)))
+define(ca3, calc(cos((pi/180)*a3)))
+define(ca4, calc(cos((pi/180)*a4)))
+define(ca5, calc(cos((pi/180)*a5)))
+define(ca6, calc(cos((pi/180)*a6)))
+define(ca7, calc(cos((pi/180)*a7)))
+
+define(sa0, calc(sin((pi/180)*a0)))
+define(sa1, calc(sin((pi/180)*a1)))
+define(sa2, calc(sin((pi/180)*a2)))
+define(sa3, calc(sin((pi/180)*a3)))
+define(sa4, calc(sin((pi/180)*a4)))
+define(sa5, calc(sin((pi/180)*a5)))
+define(sa6, calc(sin((pi/180)*a6)))
+define(sa7, calc(sin((pi/180)*a7)))
+
+define(cea0, calc(cos((pi/180)*ea0)))
+define(cea1, calc(cos((pi/180)*ea1)))
+define(cea2, calc(cos((pi/180)*ea2)))
+define(cea3, calc(cos((pi/180)*ea3)))
+define(cea4, calc(cos((pi/180)*ea4)))
+define(cea5, calc(cos((pi/180)*ea5)))
+define(cea6, calc(cos((pi/180)*ea6)))
+define(cea7, calc(cos((pi/180)*ea7)))
+
+define(sea0, calc(sin((pi/180)*ea0)))
+define(sea1, calc(sin((pi/180)*ea1)))
+define(sea2, calc(sin((pi/180)*ea2)))
+define(sea3, calc(sin((pi/180)*ea3)))
+define(sea4, calc(sin((pi/180)*ea4)))
+define(sea5, calc(sin((pi/180)*ea5)))
+define(sea6, calc(sin((pi/180)*ea6)))
+define(sea7, calc(sin((pi/180)*ea7)))
+
+define(x00, calc(r*ca0))
+define(x01, calc(r*ca1))
+define(x02, calc(r*ca2))
+define(x03, calc(r*ca3))
+define(x04, calc(r*ca4))
+define(x05, calc(r*ca5))
+define(x06, calc(r*ca6))
+define(x07, calc(r*ca7))
+
+define(x10, calc(rb*ca0))
+define(x11, calc(rb*ca1))
+define(x12, calc(rb*ca2))
+define(x13, calc(rb*ca3))
+define(x14, calc(rb*ca4))
+define(x15, calc(rb*ca5))
+define(x16, calc(rb*ca6))
+define(x17, calc(rb*ca7))
+
+define(x20, calc(ri*ca0))
+define(x21, calc(ri*ca1))
+define(x22, calc(ri*ca2))
+define(x23, calc(ri*ca3))
+define(x24, calc(ri*ca4))
+define(x25, calc(ri*ca5))
+define(x26, calc(ri*ca6))
+define(x27, calc(ri*ca7))
+
+define(x30, calc(Rb*ca0))
+define(x31, calc(Rb*ca1))
+define(x32, calc(Rb*ca2))
+define(x33, calc(Rb*ca3))
+define(x34, calc(Rb*ca4))
+define(x35, calc(Rb*ca5))
+define(x36, calc(Rb*ca6))
+define(x37, calc(Rb*ca7))
+
+define(x40, calc(R*ca0))
+define(x41, calc(R*ca1))
+define(x42, calc(R*ca2))
+define(x43, calc(R*ca3))
+define(x44, calc(R*ca4))
+define(x45, calc(R*ca5))
+define(x46, calc(R*ca6))
+define(x47, calc(R*ca7))
+
+define(y00, calc(r*sa0))
+define(y01, calc(r*sa1))
+define(y02, calc(r*sa2))
+define(y03, calc(r*sa3))
+define(y04, calc(r*sa4))
+define(y05, calc(r*sa5))
+define(y06, calc(r*sa6))
+define(y07, calc(r*sa7))
+
+define(y10, calc(rb*sa0))
+define(y11, calc(rb*sa1))
+define(y12, calc(rb*sa2))
+define(y13, calc(rb*sa3))
+define(y14, calc(rb*sa4))
+define(y15, calc(rb*sa5))
+define(y16, calc(rb*sa6))
+define(y17, calc(rb*sa7))
+
+define(y20, calc(ri*sa0))
+define(y21, calc(ri*sa1))
+define(y22, calc(ri*sa2))
+define(y23, calc(ri*sa3))
+define(y24, calc(ri*sa4))
+define(y25, calc(ri*sa5))
+define(y26, calc(ri*sa6))
+define(y27, calc(ri*sa7))
+
+define(y30, calc(Rb*sa0))
+define(y31, calc(Rb*sa1))
+define(y32, calc(Rb*sa2))
+define(y33, calc(Rb*sa3))
+define(y34, calc(Rb*sa4))
+define(y35, calc(Rb*sa5))
+define(y36, calc(Rb*sa6))
+define(y37, calc(Rb*sa7))
+
+define(y40, calc(R*sa0))
+define(y41, calc(R*sa1))
+define(y42, calc(R*sa2))
+define(y43, calc(R*sa3))
+define(y44, calc(R*sa4))
+define(y45, calc(R*sa5))
+define(y46, calc(R*sa6))
+define(y47, calc(R*sa7))
+
+define(ex00, calc(r*cea0))
+define(ex01, calc(r*cea1))
+define(ex02, calc(r*cea2))
+define(ex03, calc(r*cea3))
+define(ex04, calc(r*cea4))
+define(ex05, calc(r*cea5))
+define(ex06, calc(r*cea6))
+define(ex07, calc(r*cea7))
+
+define(ex10, calc(rb*cea0))
+define(ex11, calc(rb*cea1))
+define(ex12, calc(rb*cea2))
+define(ex13, calc(rb*cea3))
+define(ex14, calc(rb*cea4))
+define(ex15, calc(rb*cea5))
+define(ex16, calc(rb*cea6))
+define(ex17, calc(rb*cea7))
+
+define(ex20, calc(ri*cea0))
+define(ex21, calc(ri*cea1))
+define(ex22, calc(ri*cea2))
+define(ex23, calc(ri*cea3))
+define(ex24, calc(ri*cea4))
+define(ex25, calc(ri*cea5))
+define(ex26, calc(ri*cea6))
+define(ex27, calc(ri*cea7))
+
+define(ex30, calc(Rb*cea0))
+define(ex31, calc(Rb*cea1))
+define(ex32, calc(Rb*cea2))
+define(ex33, calc(Rb*cea3))
+define(ex34, calc(Rb*cea4))
+define(ex35, calc(Rb*cea5))
+define(ex36, calc(Rb*cea6))
+define(ex37, calc(Rb*cea7))
+
+define(ex40, calc(R*cea0))
+define(ex41, calc(R*cea1))
+define(ex42, calc(R*cea2))
+define(ex43, calc(R*cea3))
+define(ex44, calc(R*cea4))
+define(ex45, calc(R*cea5))
+define(ex46, calc(R*cea6))
+define(ex47, calc(R*cea7))
+
+define(ey00, calc(r*sea0))
+define(ey01, calc(r*sea1))
+define(ey02, calc(r*sea2))
+define(ey03, calc(r*sea3))
+define(ey04, calc(r*sea4))
+define(ey05, calc(r*sea5))
+define(ey06, calc(r*sea6))
+define(ey07, calc(r*sea7))
+
+define(ey10, calc(rb*sea0))
+define(ey11, calc(rb*sea1))
+define(ey12, calc(rb*sea2))
+define(ey13, calc(rb*sea3))
+define(ey14, calc(rb*sea4))
+define(ey15, calc(rb*sea5))
+define(ey16, calc(rb*sea6))
+define(ey17, calc(rb*sea7))
+
+define(ey20, calc(ri*sea0))
+define(ey21, calc(ri*sea1))
+define(ey22, calc(ri*sea2))
+define(ey23, calc(ri*sea3))
+define(ey24, calc(ri*sea4))
+define(ey25, calc(ri*sea5))
+define(ey26, calc(ri*sea6))
+define(ey27, calc(ri*sea7))
+
+define(ey30, calc(Rb*sea0))
+define(ey31, calc(Rb*sea1))
+define(ey32, calc(Rb*sea2))
+define(ey33, calc(Rb*sea3))
+define(ey34, calc(Rb*sea4))
+define(ey35, calc(Rb*sea5))
+define(ey36, calc(Rb*sea6))
+define(ey37, calc(Rb*sea7))
+
+define(ey40, calc(R*sea0))
+define(ey41, calc(R*sea1))
+define(ey42, calc(R*sea2))
+define(ey43, calc(R*sea3))
+define(ey44, calc(R*sea4))
+define(ey45, calc(R*sea5))
+define(ey46, calc(R*sea6))
+define(ey47, calc(R*sea7))
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+vertices
+(
+    vert(0, 0, Zb) vlabel(r0b)
+    vert(0, 0, Zb) vlabel(r0sb)
+    vert(0, 1, Zb) vlabel(r1b)
+    vert(0, 2, Zb) vlabel(r2b)
+    vert(0, 2, Zb) vlabel(r2sb)
+    vert(0, 3, Zb) vlabel(r3b)
+    vert(0, 4, Zb) vlabel(r4b)
+    vert(0, 4, Zb) vlabel(r4sb)
+    vert(0, 5, Zb) vlabel(r5b)
+    vert(0, 6, Zb) vlabel(r6b)
+    vert(0, 6, Zb) vlabel(r6sb)
+    vert(0, 7, Zb) vlabel(r7b)
+
+    vert(1, 0, Zb) vlabel(rb0b)
+    vert(1, 1, Zb) vlabel(rb1b)
+    vert(1, 2, Zb) vlabel(rb2b)
+    vert(1, 3, Zb) vlabel(rb3b)
+    vert(1, 4, Zb) vlabel(rb4b)
+    vert(1, 5, Zb) vlabel(rb5b)
+    vert(1, 6, Zb) vlabel(rb6b)
+    vert(1, 7, Zb) vlabel(rb7b)
+
+    vert(2, 0, Zb) vlabel(ri0b)
+    vert(2, 1, Zb) vlabel(ri1b)
+    vert(2, 2, Zb) vlabel(ri2b)
+    vert(2, 3, Zb) vlabel(ri3b)
+    vert(2, 4, Zb) vlabel(ri4b)
+    vert(2, 5, Zb) vlabel(ri5b)
+    vert(2, 6, Zb) vlabel(ri6b)
+    vert(2, 7, Zb) vlabel(ri7b)
+
+    vert(3, 0, Zb) vlabel(Rb0b)
+    vert(3, 1, Zb) vlabel(Rb1b)
+    vert(3, 2, Zb) vlabel(Rb2b)
+    vert(3, 3, Zb) vlabel(Rb3b)
+    vert(3, 4, Zb) vlabel(Rb4b)
+    vert(3, 5, Zb) vlabel(Rb5b)
+    vert(3, 6, Zb) vlabel(Rb6b)
+    vert(3, 7, Zb) vlabel(Rb7b)
+
+    vert(4, 0, Zb) vlabel(R0b)
+    vert(4, 1, Zb) vlabel(R1b)
+    vert(4, 1, Zb) vlabel(R1sb)
+    vert(4, 2, Zb) vlabel(R2b)
+    vert(4, 3, Zb) vlabel(R3b)
+    vert(4, 3, Zb) vlabel(R3sb)
+    vert(4, 4, Zb) vlabel(R4b)
+    vert(4, 5, Zb) vlabel(R5b)
+    vert(4, 5, Zb) vlabel(R5sb)
+    vert(4, 6, Zb) vlabel(R6b)
+    vert(4, 7, Zb) vlabel(R7b)
+    vert(4, 7, Zb) vlabel(R7sb)
+
+    vert(0, 0, Zt) vlabel(r0t)
+    vert(0, 0, Zt) vlabel(r0st)
+    vert(0, 1, Zt) vlabel(r1t)
+    vert(0, 2, Zt) vlabel(r2t)
+    vert(0, 2, Zt) vlabel(r2st)
+    vert(0, 3, Zt) vlabel(r3t)
+    vert(0, 4, Zt) vlabel(r4t)
+    vert(0, 4, Zt) vlabel(r4st)
+    vert(0, 5, Zt) vlabel(r5t)
+    vert(0, 6, Zt) vlabel(r6t)
+    vert(0, 6, Zt) vlabel(r6st)
+    vert(0, 7, Zt) vlabel(r7t)
+
+    vert(1, 0, Zt) vlabel(rb0t)
+    vert(1, 1, Zt) vlabel(rb1t)
+    vert(1, 2, Zt) vlabel(rb2t)
+    vert(1, 3, Zt) vlabel(rb3t)
+    vert(1, 4, Zt) vlabel(rb4t)
+    vert(1, 5, Zt) vlabel(rb5t)
+    vert(1, 6, Zt) vlabel(rb6t)
+    vert(1, 7, Zt) vlabel(rb7t)
+
+    vert(2, 0, Zt) vlabel(ri0t)
+    vert(2, 1, Zt) vlabel(ri1t)
+    vert(2, 2, Zt) vlabel(ri2t)
+    vert(2, 3, Zt) vlabel(ri3t)
+    vert(2, 4, Zt) vlabel(ri4t)
+    vert(2, 5, Zt) vlabel(ri5t)
+    vert(2, 6, Zt) vlabel(ri6t)
+    vert(2, 7, Zt) vlabel(ri7t)
+
+    vert(3, 0, Zt) vlabel(Rb0t)
+    vert(3, 1, Zt) vlabel(Rb1t)
+    vert(3, 2, Zt) vlabel(Rb2t)
+    vert(3, 3, Zt) vlabel(Rb3t)
+    vert(3, 4, Zt) vlabel(Rb4t)
+    vert(3, 5, Zt) vlabel(Rb5t)
+    vert(3, 6, Zt) vlabel(Rb6t)
+    vert(3, 7, Zt) vlabel(Rb7t)
+
+    vert(4, 0, Zt) vlabel(R0t)
+    vert(4, 1, Zt) vlabel(R1t)
+    vert(4, 1, Zt) vlabel(R1st)
+    vert(4, 2, Zt) vlabel(R2t)
+    vert(4, 3, Zt) vlabel(R3t)
+    vert(4, 3, Zt) vlabel(R3st)
+    vert(4, 4, Zt) vlabel(R4t)
+    vert(4, 5, Zt) vlabel(R5t)
+    vert(4, 5, Zt) vlabel(R5st)
+    vert(4, 6, Zt) vlabel(R6t)
+    vert(4, 7, Zt) vlabel(R7t)
+    vert(4, 7, Zt) vlabel(R7st)
+);
+
+blocks
+(
+    // block0
+    hex2D(r0, r1, rb1, rb0)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex2D(r1, r2s, rb2, rb1)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex2D(r2, r3, rb3, rb2)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block3
+    hex2D(r3, r4s, rb4, rb3)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block4
+    hex2D(r4, r5, rb5, rb4)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block5
+    hex2D(r5, r6s, rb6, rb5)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block6
+    hex2D(r6, r7, rb7, rb6)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block7
+    hex2D(r7, r0s, rb0, rb7)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block0
+    hex2D(rb0, rb1, ri1, ri0)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex2D(rb1, rb2, ri2, ri1)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex2D(rb2, rb3, ri3, ri2)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block3
+    hex2D(rb3, rb4, ri4, ri3)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block4
+    hex2D(rb4, rb5, ri5, ri4)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block5
+    hex2D(rb5, rb6, ri6, ri5)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block6
+    hex2D(rb6, rb7, ri7, ri6)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block7
+    hex2D(rb7, rb0, ri0, ri7)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block0
+    hex2D(ri0, ri1, Rb1, Rb0)
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex2D(ri1, ri2, Rb2, Rb1)
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex2D(ri2, ri3, Rb3, Rb2)
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block3
+    hex2D(ri3, ri4, Rb4, Rb3)
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block4
+    hex2D(ri4, ri5, Rb5, Rb4)
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block5
+    hex2D(ri5, ri6, Rb6, Rb5)
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block6
+    hex2D(ri6, ri7, Rb7, Rb6)
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block7
+    hex2D(ri7, ri0, Rb0, Rb7)
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block0
+    hex2D(Rb0, Rb1, R1s, R0)
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex2D(Rb1, Rb2, R2, R1)
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex2D(Rb2, Rb3, R3s, R2)
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block3
+    hex2D(Rb3, Rb4, R4, R3)
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block4
+    hex2D(Rb4, Rb5, R5s, R4)
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block5
+    hex2D(Rb5, Rb6, R6, R5)
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block6
+    hex2D(Rb6, Rb7, R7s, R6)
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block7
+    hex2D(Rb7, Rb0, R0, R7)
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+);
+
+edges
+(
+    arc r0b r1b evert(0, 0, Zb)
+    arc r1b r2sb evert(0, 1, Zb)
+    arc r2b r3b evert(0, 2, Zb)
+    arc r3b r4sb evert(0, 3, Zb)
+    arc r4b r5b evert(0, 4, Zb)
+    arc r5b r6sb evert(0, 5, Zb)
+    arc r6b r7b evert(0, 6, Zb)
+    arc r7b r0sb evert(0, 7, Zb)
+
+    arc rb0b rb1b evert(1, 0, Zb)
+    arc rb1b rb2b evert(1, 1, Zb)
+    arc rb2b rb3b evert(1, 2, Zb)
+    arc rb3b rb4b evert(1, 3, Zb)
+    arc rb4b rb5b evert(1, 4, Zb)
+    arc rb5b rb6b evert(1, 5, Zb)
+    arc rb6b rb7b evert(1, 6, Zb)
+    arc rb7b rb0b evert(1, 7, Zb)
+
+    arc ri0b ri1b evert(2, 0, Zb)
+    arc ri1b ri2b evert(2, 1, Zb)
+    arc ri2b ri3b evert(2, 2, Zb)
+    arc ri3b ri4b evert(2, 3, Zb)
+    arc ri4b ri5b evert(2, 4, Zb)
+    arc ri5b ri6b evert(2, 5, Zb)
+    arc ri6b ri7b evert(2, 6, Zb)
+    arc ri7b ri0b evert(2, 7, Zb)
+
+    arc Rb0b Rb1b evert(3, 0, Zb)
+    arc Rb1b Rb2b evert(3, 1, Zb)
+    arc Rb2b Rb3b evert(3, 2, Zb)
+    arc Rb3b Rb4b evert(3, 3, Zb)
+    arc Rb4b Rb5b evert(3, 4, Zb)
+    arc Rb5b Rb6b evert(3, 5, Zb)
+    arc Rb6b Rb7b evert(3, 6, Zb)
+    arc Rb7b Rb0b evert(3, 7, Zb)
+
+    arc R0b R1sb evert(4, 0, Zb)
+    arc R1b R2b evert(4, 1, Zb)
+    arc R2b R3sb evert(4, 2, Zb)
+    arc R3b R4b evert(4, 3, Zb)
+    arc R4b R5sb evert(4, 4, Zb)
+    arc R5b R6b evert(4, 5, Zb)
+    arc R6b R7sb evert(4, 6, Zb)
+    arc R7b R0b evert(4, 7, Zb)
+
+    arc r0t r1t evert(0, 0, Zt)
+    arc r1t r2st evert(0, 1, Zt)
+    arc r2t r3t evert(0, 2, Zt)
+    arc r3t r4st evert(0, 3, Zt)
+    arc r4t r5t evert(0, 4, Zt)
+    arc r5t r6st evert(0, 5, Zt)
+    arc r6t r7t evert(0, 6, Zt)
+    arc r7t r0st evert(0, 7, Zt)
+
+    arc rb0t rb1t evert(1, 0, Zt)
+    arc rb1t rb2t evert(1, 1, Zt)
+    arc rb2t rb3t evert(1, 2, Zt)
+    arc rb3t rb4t evert(1, 3, Zt)
+    arc rb4t rb5t evert(1, 4, Zt)
+    arc rb5t rb6t evert(1, 5, Zt)
+    arc rb6t rb7t evert(1, 6, Zt)
+    arc rb7t rb0t evert(1, 7, Zt)
+
+    arc ri0t ri1t evert(2, 0, Zt)
+    arc ri1t ri2t evert(2, 1, Zt)
+    arc ri2t ri3t evert(2, 2, Zt)
+    arc ri3t ri4t evert(2, 3, Zt)
+    arc ri4t ri5t evert(2, 4, Zt)
+    arc ri5t ri6t evert(2, 5, Zt)
+    arc ri6t ri7t evert(2, 6, Zt)
+    arc ri7t ri0t evert(2, 7, Zt)
+
+    arc Rb0t Rb1t evert(3, 0, Zt)
+    arc Rb1t Rb2t evert(3, 1, Zt)
+    arc Rb2t Rb3t evert(3, 2, Zt)
+    arc Rb3t Rb4t evert(3, 3, Zt)
+    arc Rb4t Rb5t evert(3, 4, Zt)
+    arc Rb5t Rb6t evert(3, 5, Zt)
+    arc Rb6t Rb7t evert(3, 6, Zt)
+    arc Rb7t Rb0t evert(3, 7, Zt)
+
+    arc R0t R1st evert(4, 0, Zt)
+    arc R1t R2t evert(4, 1, Zt)
+    arc R2t R3st evert(4, 2, Zt)
+    arc R3t R4t evert(4, 3, Zt)
+    arc R4t R5st evert(4, 4, Zt)
+    arc R5t R6t evert(4, 5, Zt)
+    arc R6t R7st evert(4, 6, Zt)
+    arc R7t R0t evert(4, 7, Zt)
+);
+
+patches
+(
+    wall rotor
+    (
+        quad2D(r0, r1)
+        quad2D(r1, r2s)
+        quad2D(r2, r3)
+        quad2D(r3, r4s)
+        quad2D(r4, r5)
+        quad2D(r5, r6s)
+        quad2D(r6, r7)
+        quad2D(r7, r0s)
+
+        quad2D(r0, rb0)
+        quad2D(r0s, rb0)
+
+        quad2D(r2, rb2)
+        quad2D(r2s, rb2)
+
+        quad2D(r4, rb4)
+        quad2D(r4s, rb4)
+
+        quad2D(r6, rb6)
+        quad2D(r6s, rb6)
+    )
+
+    wall stator
+    (
+        quad2D(R0, R1s)
+        quad2D(R1, R2)
+        quad2D(R2, R3s)
+        quad2D(R3, R4)
+        quad2D(R4, R5s)
+        quad2D(R5, R6)
+        quad2D(R6, R7s)
+        quad2D(R7, R0)
+
+        quad2D(R1, Rb1)
+        quad2D(R1s, Rb1)
+
+        quad2D(R3, Rb3)
+        quad2D(R3s, Rb3)
+
+        quad2D(R5, Rb5)
+        quad2D(R5s, Rb5)
+
+        quad2D(R7, Rb7)
+        quad2D(R7s, Rb7)
+    )
+
+    empty front
+    (
+        frontQuad(r0, r1, rb1, rb0)
+        frontQuad(r1, r2s, rb2, rb1)
+        frontQuad(r2, r3, rb3, rb2)
+        frontQuad(r3, r4s, rb4, rb3)
+        frontQuad(r4, r5, rb5, rb4)
+        frontQuad(r5, r6s, rb6, rb5)
+        frontQuad(r6, r7, rb7, rb6)
+        frontQuad(r7, r0s, rb0, rb7)
+        frontQuad(rb0, rb1, ri1, ri0)
+        frontQuad(rb1, rb2, ri2, ri1)
+        frontQuad(rb2, rb3, ri3, ri2)
+        frontQuad(rb3, rb4, ri4, ri3)
+        frontQuad(rb4, rb5, ri5, ri4)
+        frontQuad(rb5, rb6, ri6, ri5)
+        frontQuad(rb6, rb7, ri7, ri6)
+        frontQuad(rb7, rb0, ri0, ri7)
+        frontQuad(ri0, ri1, Rb1, Rb0)
+        frontQuad(ri1, ri2, Rb2, Rb1)
+        frontQuad(ri2, ri3, Rb3, Rb2)
+        frontQuad(ri3, ri4, Rb4, Rb3)
+        frontQuad(ri4, ri5, Rb5, Rb4)
+        frontQuad(ri5, ri6, Rb6, Rb5)
+        frontQuad(ri6, ri7, Rb7, Rb6)
+        frontQuad(ri7, ri0, Rb0, Rb7)
+        frontQuad(Rb0, Rb1, R1s, R0)
+        frontQuad(Rb1, Rb2, R2, R1)
+        frontQuad(Rb2, Rb3, R3s, R2)
+        frontQuad(Rb3, Rb4, R4, R3)
+        frontQuad(Rb4, Rb5, R5s, R4)
+        frontQuad(Rb5, Rb6, R6, R5)
+        frontQuad(Rb6, Rb7, R7s, R6)
+        frontQuad(Rb7, Rb0, R0, R7)
+    )
+
+    empty back
+    (
+        backQuad(r0, r1, rb1, rb0)
+        backQuad(r1, r2s, rb2, rb1)
+        backQuad(r2, r3, rb3, rb2)
+        backQuad(r3, r4s, rb4, rb3)
+        backQuad(r4, r5, rb5, rb4)
+        backQuad(r5, r6s, rb6, rb5)
+        backQuad(r6, r7, rb7, rb6)
+        backQuad(r7, r0s, rb0, rb7)
+        backQuad(rb0, rb1, ri1, ri0)
+        backQuad(rb1, rb2, ri2, ri1)
+        backQuad(rb2, rb3, ri3, ri2)
+        backQuad(rb3, rb4, ri4, ri3)
+        backQuad(rb4, rb5, ri5, ri4)
+        backQuad(rb5, rb6, ri6, ri5)
+        backQuad(rb6, rb7, ri7, ri6)
+        backQuad(rb7, rb0, ri0, ri7)
+        backQuad(ri0, ri1, Rb1, Rb0)
+        backQuad(ri1, ri2, Rb2, Rb1)
+        backQuad(ri2, ri3, Rb3, Rb2)
+        backQuad(ri3, ri4, Rb4, Rb3)
+        backQuad(ri4, ri5, Rb5, Rb4)
+        backQuad(ri5, ri6, Rb6, Rb5)
+        backQuad(ri6, ri7, Rb7, Rb6)
+        backQuad(ri7, ri0, Rb0, Rb7)
+        backQuad(Rb0, Rb1, R1s, R0)
+        backQuad(Rb1, Rb2, R2, R1)
+        backQuad(Rb2, Rb3, R3s, R2)
+        backQuad(Rb3, Rb4, R4, R3)
+        backQuad(Rb4, Rb5, R5s, R4)
+        backQuad(Rb5, Rb6, R6, R5)
+        backQuad(Rb6, Rb7, R7s, R6)
+        backQuad(Rb7, Rb0, R0, R7)
+    )
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/polyMesh/boundary b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/polyMesh/boundary
new file mode 100644
index 0000000000000000000000000000000000000000..292f25b806357d9df75c7731f74dee0ec0aa3a40
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/polyMesh/boundary
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+4
+(
+    rotor
+    {
+        type            wall;
+        nFaces          192;
+        startFace       5952;
+    }
+    stator
+    {
+        type            wall;
+        nFaces          192;
+        startFace       6144;
+    }
+    front
+    {
+        type            empty;
+        nFaces          3072;
+        startFace       6336;
+    }
+    back
+    {
+        type            empty;
+        nFaces          3072;
+        startFace       9408;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/ppProperties b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/ppProperties
new file mode 100644
index 0000000000000000000000000000000000000000..eddbfd7ad7910a8a85af9c310b1424eb18d3eb1e
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/ppProperties
@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      ppProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+preAlphaExp     300;
+
+expMax          10;
+
+alphaMax        0.6;
+
+g0              g0 [ 1 -1 -2 0 0 0 0 ] 0;
+
+packingLimiter  off;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/transportProperties b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..9ea9ed21895c06facd069c9507e3de6663ff40bf
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/constant/transportProperties
@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phase1
+{
+    rho0            0;
+    rho             0.88;
+    R               287;
+    Cp              1007;
+    nu              2.46e-05;
+    d               4e-3;
+
+    kappa           2.63e-2;
+    diameterModel   isothermal;
+    isothermalCoeffs
+    {
+        d0              3e-3;
+        p0              1e5;
+    }
+}
+
+phase2
+{
+    rho             733;
+    rho0            733;
+    R               1e10;
+    Cp              4195;
+    nu              2.73e-6;
+    d               1e-4;
+
+    kappa           0.668;
+    diameterModel constant;
+    constantCoeffs
+    {
+        d               1e-4;
+    }
+}
+
+// Virtual-mass ceofficient
+Cvm             0.5;
+
+// Lift coefficient
+Cl               0;
+
+// Dispersed-phase turbulence coefficient
+Ct               0.2;
+
+// Minimum allowable pressure
+pMin            10000;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/makeMesh b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/makeMesh
new file mode 100755
index 0000000000000000000000000000000000000000..8ef4993fdebc0faac8bb6c69426059aba4b3faac
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/makeMesh
@@ -0,0 +1,6 @@
+#!/bin/sh
+
+m4 < constant/polyMesh/blockMeshDict.m4 > constant/polyMesh/blockMeshDict
+blockMesh
+topoSet
+setsToZones -noFlipMap
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/system/controlDict b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..7b2bf5a816c8966dabb1ec9ec14fefd77ddb630e
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/system/controlDict
@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     twoPhaseEulerFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         20;
+
+deltaT          1e-4;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           0.5;
+
+maxDeltaT       1;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..22f964670702bc54b4221f5d022ce0f0bf942bd8
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/system/fvSchemes
@@ -0,0 +1,62 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default             none;
+
+    div(phi,alpha1)         Gauss limitedLinear01 1;
+    div(phir,alpha1)        Gauss limitedLinear01 1;
+    "div\(phi.,U.\)"        Gauss limitedLinearV 1;
+    "div\(Rc.\)"            Gauss linear;
+    "div\(phi.,(k|epsilon)\)"  Gauss limitedLinear 1;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+    alpha1;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..69208cbe950050e6f399b103952793cc7037bdf1
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/system/fvSolution
@@ -0,0 +1,103 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.05;
+        smoother        GaussSeidel;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    pFinal
+    {
+        $p;
+        tolerance       1e-08;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        nSweeps         2;
+        tolerance       1e-07;
+        relTol          0.1;
+    }
+
+    UFinal
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-07;
+        relTol          0;
+    }
+
+    "alpha.*"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-10;
+        relTol          0;
+    }
+
+    "(k|epsilon|Theta|T).*"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-06;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nOuterCorrectors 1;
+    nCorrectors     3;
+    nNonOrthogonalCorrectors 0;
+    nAlphaCorr      1;
+    nAlphaSubCycles 3;
+    correctAlpha    yes;
+    pRefCell        0;
+    pRefValue       0;
+}
+
+relaxationFactors
+{
+    fields
+    {
+        p                1;
+    }
+    equations
+    {
+        "U.*"            1;
+        "alpha.*"        1;
+        "Theta.*"        1;
+        "k.*"            1;
+        "epsilon.*"      1;
+        "T.*"            1;
+    }
+}
+
+
+// ************************************************************************* //