driftFluxFoam: New version of settlingFoam including MULES and MRF/fvOptions support

2014-02-19 16:33:51 +00:00
parent cecf40bc81
commit 9fbfc7da57
67 changed files with 4515 additions and 0 deletions
--- a/applications/solvers/multiphase/driftFluxFoam/Make/files
+++ b/applications/solvers/multiphase/driftFluxFoam/Make/files
@ -0,0 +1,3 @@
+driftFluxFoam.C
+
+EXE = $(FOAM_APPBIN)/driftFluxFoam
--- a/applications/solvers/multiphase/driftFluxFoam/Make/options
+++ b/applications/solvers/multiphase/driftFluxFoam/Make/options
@ -0,0 +1,11 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions
--- a/applications/solvers/multiphase/driftFluxFoam/UEqn.H
+++ b/applications/solvers/multiphase/driftFluxFoam/UEqn.H
@ -0,0 +1,38 @@
+    // Solve the Momentum equation
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(rhoPhi, U)
+      + fvc::div
+        (
+            (alpha/(scalar(1.001) - alpha))*((rhoc*rhod)/rho)*Vdj*Vdj,
+            "div(phiVdj,Vdj)"
+        )
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev2(T(fvc::grad(U))))
+     ==
+        fvOptions(rho, U)
+    );
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+          ==
+            fvc::reconstruct
+            (
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            )
+        );
+
+        fvOptions.correct(U);
+    }
--- a/applications/solvers/multiphase/driftFluxFoam/alphaControls.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaControls.H
@ -0,0 +1,4 @@
+const dictionary& alphaControls = mesh.solverDict(alpha.name());
+
+label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
+label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
--- a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
@ -0,0 +1,31 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+    {
+        phiAlpha =
+        (
+            fvc::flux
+            (
+                phi,
+                alpha,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                phir,
+                alpha,
+                alpharScheme
+            )
+        );
+
+        MULES::explicitSolve(alpha, phi, phiAlpha, 1, 0);
+    }
+
+    Info<< "Phase-1 volume fraction = "
+        << alpha.weightedAverage(mesh.Vsc()).value()
+        << "  Min(alpha) = " << min(alpha).value()
+        << "  Max(alpha) = " << max(alpha).value()
+        << endl;
+}
--- a/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
@ -0,0 +1,73 @@
+{
+    surfaceScalarField phiAlpha
+    (
+        IOobject
+        (
+            "phiAlpha",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("0", phi.dimensions(), 0)
+    );
+
+    surfaceScalarField phir
+    (
+        rhoc*(mesh.Sf() & fvc::interpolate(Vdj/rho))
+    );
+
+    if (nAlphaSubCycles > 1)
+    {
+        dimensionedScalar totalDeltaT = runTime.deltaT();
+        surfaceScalarField phiAlphaSum
+        (
+            IOobject
+            (
+                "phiAlphaSum",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar("0", phi.dimensions(), 0)
+        );
+
+        for
+        (
+            subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
+        )
+        {
+            #include "alphaEqn.H"
+            phiAlphaSum += (runTime.deltaT()/totalDeltaT)*phiAlpha;
+        }
+
+        phiAlpha = phiAlphaSum;
+    }
+    else
+    {
+        #include "alphaEqn.H"
+    }
+
+    // Apply the diffusion term separately to allow implicit solution
+    // and boundedness of the explicit advection
+    {
+        fvScalarMatrix alphaEqn
+        (
+            fvm::ddt(alpha) - fvc::ddt(alpha)
+          - fvm::laplacian(mut/rho, alpha)
+        );
+
+        alphaEqn.solve();
+
+        phiAlpha += alphaEqn.flux();
+    }
+
+    Info<< "Phase-1 volume fraction = "
+        << alpha.weightedAverage(mesh.Vsc()).value()
+        << "  Min(alpha) = " << min(alpha).value()
+        << "  Max(alpha) = " << max(alpha).value()
+        << endl;
+
+    rhoPhi = phiAlpha*(rhod - rhoc) + phi*rhoc;
+    rho == alpha*rhod + (scalar(1) - alpha)*rhoc;
+}
--- a/applications/solvers/multiphase/driftFluxFoam/calcVdj.H
+++ b/applications/solvers/multiphase/driftFluxFoam/calcVdj.H
@ -0,0 +1,20 @@
+if (VdjModel == "general")
+{
+    Vdj = V0*
+    (
+        exp(-a*max(alpha - alphaMin, scalar(0)))
+      - exp(-a1*max(alpha - alphaMin, scalar(0)))
+    );
+}
+else if (VdjModel == "simple")
+{
+    Vdj = V0*pow(10.0, -a*max(alpha, scalar(0)));
+}
+else
+{
+    FatalErrorIn(args.executable())
+        << "Unknown VdjModel : " << VdjModel
+        << abort(FatalError);
+}
+
+Vdj.correctBoundaryConditions();
--- a/applications/solvers/multiphase/driftFluxFoam/correctViscosity.H
+++ b/applications/solvers/multiphase/driftFluxFoam/correctViscosity.H
@ -0,0 +1,39 @@
+{
+    mul = muc +
+        plasticViscosity
+        (
+            plasticViscosityCoeff,
+            plasticViscosityExponent,
+            alpha
+        );
+
+    if (BinghamPlastic)
+    {
+        volScalarField tauy = yieldStress
+        (
+            yieldStressCoeff,
+            yieldStressExponent,
+            yieldStressOffset,
+            alpha
+        );
+
+        mul =
+            tauy/
+            (
+                mag(fvc::grad(U))
+              + 1.0e-4*
+                (
+                    tauy
+                  + dimensionedScalar
+                    (
+                        "deltaTauy",
+                        tauy.dimensions(),
+                        1.0e-15
+                    )
+                )/mul
+            )
+          + mul;
+    }
+
+    mul = min(mul, muMax);
+}
--- a/applications/solvers/multiphase/driftFluxFoam/createFields.H
+++ b/applications/solvers/multiphase/driftFluxFoam/createFields.H
@ -0,0 +1,384 @@
+    Info<< "Reading field p_rgh\n" << endl;
+    volScalarField p_rgh
+    (
+        IOobject
+        (
+            "p_rgh",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field alpha\n" << endl;
+    volScalarField alpha
+    (
+        IOobject
+        (
+            "alpha",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "createPhi.H"
+
+
+    Info<< "Reading transportProperties\n" << endl;
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+
+    dimensionedScalar rhoc(transportProperties.lookup("rhoc"));
+
+    dimensionedScalar rhod(transportProperties.lookup("rhod"));
+
+    dimensionedScalar muc(transportProperties.lookup("muc"));
+    dimensionedScalar muMax(transportProperties.lookup("muMax"));
+
+    dimensionedScalar plasticViscosityCoeff
+    (
+        transportProperties.lookup("plasticViscosityCoeff")
+    );
+
+    dimensionedScalar plasticViscosityExponent
+    (
+        transportProperties.lookup("plasticViscosityExponent")
+    );
+
+    dimensionedScalar yieldStressCoeff
+    (
+        transportProperties.lookup("yieldStressCoeff")
+    );
+
+    dimensionedScalar yieldStressExponent
+    (
+        transportProperties.lookup("yieldStressExponent")
+    );
+
+    dimensionedScalar yieldStressOffset
+    (
+        transportProperties.lookup("yieldStressOffset")
+    );
+
+    Switch BinghamPlastic(transportProperties.lookup("BinghamPlastic"));
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        alpha*rhod + (scalar(1) - alpha)*rhoc
+    );
+    rho.oldTime();
+
+    // Mass flux
+    surfaceScalarField rhoPhi
+    (
+        IOobject
+        (
+            "rhoPhi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        fvc::interpolate(rho)*phi
+    );
+
+    Info<< "Calculating field mul\n" << endl;
+    volScalarField mul
+    (
+        IOobject
+        (
+            "mul",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        muc
+      + plasticViscosity
+        (
+            plasticViscosityCoeff,
+            plasticViscosityExponent,
+            alpha
+        )
+    );
+
+
+    Info<< "Initialising field Vdj\n" << endl;
+    volVectorField Vdj
+    (
+        IOobject
+        (
+            "Vdj",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("0.0", U.dimensions(), vector::zero),
+        U.boundaryField().types()
+    );
+
+
+    Info<< "Selecting Drift-Flux model " << endl;
+
+    const word VdjModel(transportProperties.lookup("VdjModel"));
+
+    Info<< tab << VdjModel << " selected\n" << endl;
+
+    const dictionary& VdjModelCoeffs
+    (
+        transportProperties.subDict(VdjModel + "Coeffs")
+    );
+
+    dimensionedVector V0(VdjModelCoeffs.lookup("V0"));
+
+    dimensionedScalar a(VdjModelCoeffs.lookup("a"));
+
+    dimensionedScalar a1(VdjModelCoeffs.lookup("a1"));
+
+    dimensionedScalar alphaMin(VdjModelCoeffs.lookup("alphaMin"));
+
+
+    IOdictionary RASProperties
+    (
+        IOobject
+        (
+            "RASProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+
+    Switch turbulence(RASProperties.lookup("turbulence"));
+
+    dictionary kEpsilonDict(RASProperties.subDictPtr("kEpsilonCoeffs"));
+
+    dimensionedScalar Cmu
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "Cmu",
+            kEpsilonDict,
+            0.09
+        )
+    );
+
+    dimensionedScalar C1
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "C1",
+            kEpsilonDict,
+            1.44
+        )
+    );
+
+    dimensionedScalar C2
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "C2",
+            kEpsilonDict,
+            1.92
+        )
+    );
+
+    dimensionedScalar C3
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "C3",
+            kEpsilonDict,
+            0.85
+        )
+    );
+
+    dimensionedScalar sigmak
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "sigmak",
+            kEpsilonDict,
+            1.0
+        )
+    );
+
+    dimensionedScalar sigmaEps
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "sigmaEps",
+            kEpsilonDict,
+            1.3
+        )
+    );
+
+    dictionary wallFunctionDict(RASProperties.subDictPtr("wallFunctionCoeffs"));
+
+    dimensionedScalar kappa
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "kappa",
+            wallFunctionDict,
+            0.41
+        )
+    );
+
+    dimensionedScalar E
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "E",
+            wallFunctionDict,
+            9.8
+        )
+    );
+
+    if (RASProperties.lookupOrDefault("printCoeffs", false))
+    {
+        Info<< "kEpsilonCoeffs" << kEpsilonDict << nl
+            << "wallFunctionCoeffs" << wallFunctionDict << endl;
+    }
+
+
+    nearWallDist y(mesh);
+
+    Info<< "Reading field k\n" << endl;
+    volScalarField k
+    (
+        IOobject
+        (
+            "k",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field epsilon\n" << endl;
+    volScalarField epsilon
+    (
+        IOobject
+        (
+            "epsilon",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Calculating field mut\n" << endl;
+    volScalarField mut
+    (
+        IOobject
+        (
+            "mut",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        Cmu*rho*sqr(k)/epsilon
+    );
+
+
+    Info<< "Calculating field muEff\n" << endl;
+    volScalarField muEff
+    (
+        IOobject
+        (
+            "muEff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mut + mul
+    );
+
+
+    Info<< "Calculating field (g.h)f\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("gh", g & mesh.Cf());
+
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        p_rgh + rho*gh
+    );
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell
+    (
+        p,
+        p_rgh,
+        mesh.solutionDict().subDict("PIMPLE"),
+        pRefCell,
+        pRefValue
+    );
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
--- a/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
+++ b/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 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
+    driftFluxFoam
+
+Description
+    Solver for 2 incompressible fluids using the mixture approach with the
+    drift-flux approximation for relative motion of the phases.
+
+    Used for simulating the settling of the dispersed phase and other similar
+    separation problems.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "CMULES.H"
+#include "subCycle.H"
+#include "nearWallDist.H"
+#include "wallFvPatch.H"
+#include "bound.H"
+#include "Switch.H"
+#include "plasticViscosity.H"
+#include "yieldStress.H"
+#include "pimpleControl.H"
+#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    pimpleControl pimple(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "alphaControls.H"
+
+            #include "calcVdj.H"
+            #include "alphaEqnSubCycle.H"
+            #include "correctViscosity.H"
+
+            #include "UEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                #include "kEpsilon.H"
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H
+++ b/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H
@ -0,0 +1,78 @@
+if (turbulence)
+{
+    if (mesh.changing())
+    {
+        y.correct();
+    }
+
+    dimensionedScalar k0("k0", k.dimensions(), 0);
+    dimensionedScalar kMin("kMin", k.dimensions(), SMALL);
+    dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), 0);
+    dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL);
+
+    volScalarField divU(fvc::div(rhoPhi/fvc::interpolate(rho)));
+
+    tmp<volTensorField> tgradU = fvc::grad(U);
+    volScalarField G(mut*(tgradU() && dev(twoSymm(tgradU()))));
+    tgradU.clear();
+
+    volScalarField Gcoef
+    (
+        Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin)
+    );
+
+    #include "wallFunctions.H"
+
+    // Dissipation equation
+    fvScalarMatrix epsEqn
+    (
+        fvm::ddt(rho, epsilon)
+      + fvm::div(rhoPhi, epsilon)
+      - fvm::laplacian
+        (
+            mut/sigmaEps + muc, epsilon,
+            "laplacian(DepsilonEff,epsilon)"
+        )
+     ==
+        C1*G*epsilon/(k + kMin)
+      - fvm::SuSp(C1*(1.0 - C3)*Gcoef + (2.0/3.0*C1)*rho*divU, epsilon)
+      - fvm::Sp(C2*rho*epsilon/(k + kMin), epsilon)
+    );
+
+    #include "wallDissipation.H"
+
+    epsEqn.relax();
+    epsEqn.solve();
+
+    bound(epsilon, epsilon0);
+
+
+    // Turbulent kinetic energy equation
+    fvScalarMatrix kEqn
+    (
+        fvm::ddt(rho, k)
+      + fvm::div(rhoPhi, k)
+      - fvm::laplacian
+        (
+            mut/sigmak + muc, k,
+            "laplacian(DkEff,k)"
+        )
+     ==
+        G
+      - fvm::SuSp(Gcoef + 2.0/3.0*rho*divU, k)
+      - fvm::Sp(rho*epsilon/(k + kMin), k)
+    );
+
+    kEqn.relax();
+    kEqn.solve();
+
+    bound(k, k0);
+
+
+    //- Re-calculate viscosity
+    mut = rho*Cmu*sqr(k)/(epsilon + epsilonMin);
+
+    #include "wallViscosity.H"
+}
+
+muEff = mut + mul;
--- a/applications/solvers/multiphase/driftFluxFoam/pEqn.H
+++ b/applications/solvers/multiphase/driftFluxFoam/pEqn.H
@ -0,0 +1,73 @@
+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
+    );
+    adjustPhi(phiHbyA, U, p_rgh);
+    fvOptions.makeRelative(phiHbyA);
+
+    surfaceScalarField phig
+    (
+        (
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
+
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_rghEqn.flux();
+
+            p_rgh.relax();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+        }
+    }
+
+    #include "continuityErrs.H"
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+}
--- a/applications/solvers/multiphase/driftFluxFoam/plasticViscosity.H
+++ b/applications/solvers/multiphase/driftFluxFoam/plasticViscosity.H
@ -0,0 +1,21 @@
+volScalarField plasticViscosity
+(
+    const dimensionedScalar& plasticViscosityCoeff,
+    const dimensionedScalar& plasticViscosityExponent,
+    const volScalarField& Alpha
+)
+{
+    tmp<volScalarField> tfld
+    (
+         plasticViscosityCoeff*
+         (
+             pow
+             (
+                 10.0,
+                 plasticViscosityExponent*Alpha + SMALL
+             ) - scalar(1)
+         )
+    );
+
+    return tfld();
+}
--- a/applications/solvers/multiphase/driftFluxFoam/wallDissipation.H
+++ b/applications/solvers/multiphase/driftFluxFoam/wallDissipation.H
@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    wallDissipation
+
+Description
+    Set wall dissipation in the epsilon matrix
+
+\*---------------------------------------------------------------------------*/
+
+{
+    const fvPatchList& patches = mesh.boundary();
+
+    forAll(patches, patchi)
+    {
+        const fvPatch& p = patches[patchi];
+
+        if (isA<wallFvPatch>(p))
+        {
+            epsEqn.setValues
+            (
+                p.faceCells(),
+                epsilon.boundaryField()[patchi].patchInternalField()
+            );
+        }
+    }
+}
+
+// ************************************************************************* //
--- a/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H
+++ b/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H
@ -0,0 +1,85 @@
+{
+    labelList cellBoundaryFaceCount(epsilon.size(), 0);
+
+    const scalar Cmu25 = ::pow(Cmu.value(), 0.25);
+    const scalar Cmu75 = ::pow(Cmu.value(), 0.75);
+    const scalar kappa_ = kappa.value();
+    const scalar muc_ = muc.value();
+
+    const fvPatchList& patches = mesh.boundary();
+
+    //- Initialise the near-wall P field to zero
+    forAll(patches, patchi)
+    {
+        const fvPatch& curPatch = patches[patchi];
+
+        if (isA<wallFvPatch>(curPatch))
+        {
+            forAll(curPatch, facei)
+            {
+                label faceCelli = curPatch.faceCells()[facei];
+
+                epsilon[faceCelli] = 0.0;
+                G[faceCelli] = 0.0;
+            }
+        }
+    }
+
+    //- Accumulate the wall face contributions to epsilon and G
+    //  Increment cellBoundaryFaceCount for each face for averaging
+    forAll(patches, patchi)
+    {
+        const fvPatch& curPatch = patches[patchi];
+
+        if (isA<wallFvPatch>(curPatch))
+        {
+            const scalarField& mutw = mut.boundaryField()[patchi];
+
+            scalarField magFaceGradU
+            (
+                mag(U.boundaryField()[patchi].snGrad())
+            );
+
+            forAll(curPatch, facei)
+            {
+                label faceCelli = curPatch.faceCells()[facei];
+
+                // For corner cells (with two boundary or more faces),
+                // epsilon and G in the near-wall cell are calculated
+                // as an average
+
+                cellBoundaryFaceCount[faceCelli]++;
+
+                epsilon[faceCelli] +=
+                    Cmu75*::pow(k[faceCelli], 1.5)
+                   /(kappa_*y[patchi][facei]);
+
+                G[faceCelli] +=
+                    (mutw[facei] + muc_)
+                   *magFaceGradU[facei]
+                   *Cmu25*::sqrt(k[faceCelli])
+                   /(kappa_*y[patchi][facei]);
+            }
+        }
+    }
+
+
+    // perform the averaging
+
+    forAll(patches, patchi)
+    {
+        const fvPatch& curPatch = patches[patchi];
+
+        if (isA<wallFvPatch>(curPatch))
+        {
+            forAll(curPatch, facei)
+            {
+                label faceCelli = curPatch.faceCells()[facei];
+
+                epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
+                G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
+                cellBoundaryFaceCount[faceCelli] = 1;
+            }
+        }
+    }
+}
--- a/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H
+++ b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H
@ -0,0 +1,38 @@
+{
+    const scalar Cmu25 = ::pow(Cmu.value(), 0.25);
+    const scalar kappa_ = kappa.value();
+    const scalar E_ = E.value();
+    const scalar muc_ = muc.value();
+    const scalar nuc_ = muc_/rhoc.value();
+
+    const fvPatchList& patches = mesh.boundary();
+
+    forAll(patches, patchi)
+    {
+        const fvPatch& curPatch = patches[patchi];
+
+        if (isA<wallFvPatch>(curPatch))
+        {
+            scalarField& mutw = mut.boundaryField()[patchi];
+
+            forAll(curPatch, facei)
+            {
+                label faceCelli = curPatch.faceCells()[facei];
+
+                scalar yPlus =
+                    Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
+                   /nuc_;
+
+                if (yPlus > 11.6)
+                {
+                    mutw[facei] =
+                        muc_*(yPlus*kappa_/::log(E_*yPlus) - 1);
+                }
+                else
+                {
+                    mutw[facei] = 0.0;
+                }
+            }
+        }
+    }
+}
--- a/applications/solvers/multiphase/driftFluxFoam/yieldStress.H
+++ b/applications/solvers/multiphase/driftFluxFoam/yieldStress.H
@ -0,0 +1,27 @@
+volScalarField yieldStress
+(
+    const dimensionedScalar& yieldStressCoeff,
+    const dimensionedScalar& yieldStressExponent,
+    const dimensionedScalar& yieldStressOffset,
+    const volScalarField& alpha
+)
+{
+    tmp<volScalarField> tfld
+    (
+        yieldStressCoeff*
+        (
+            pow
+            (
+                10.0,
+                yieldStressExponent*(max(alpha, scalar(0)) + yieldStressOffset)
+            )
+          - pow
+            (
+                10.0,
+                yieldStressExponent*yieldStressOffset
+            )
+        )
+    );
+
+    return tfld();
+}
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/0/U
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/0/U
@ -0,0 +1,58 @@
+/*--------------------------------*- 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      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0.0191 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
+    }
+
+    bottomWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    endWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    top
+    {
+        type            slip;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/0/alpha
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/0/alpha
@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.001;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+    }
+
+    bottomWall
+    {
+        type            zeroGradient;
+    }
+
+    endWall
+    {
+        type            zeroGradient;
+    }
+
+    top
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/0/epsilon
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/0/epsilon
@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 1.50919e-06;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+    }
+
+    bottomWall
+    {
+        type            zeroGradient;
+    }
+
+    endWall
+    {
+        type            zeroGradient;
+    }
+
+    top
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/0/k
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/0/k
@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.00015;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+    }
+
+    bottomWall
+    {
+        type            zeroGradient;
+    }
+
+    endWall
+    {
+        type            zeroGradient;
+    }
+
+    top
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/0/p_rgh
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/0/p_rgh
@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    bottomWall
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    endWall
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    top
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/RASProperties
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/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      on;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/g
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/g
@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 -9.81 0);
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/polyMesh/blockMeshDict
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/polyMesh/blockMeshDict
@ -0,0 +1,105 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    (0 0 -0.1)
+    (8.65 0 -0.1)
+    (8.65 0.1 -0.1)
+    (0 0.1 -0.1)
+    (8.65 1 -0.1)
+    (0 1 -0.1)
+    (0 0 0.1)
+    (8.65 0 0.1)
+    (8.65 0.1 0.1)
+    (0 0.1 0.1)
+    (8.65 1 0.1)
+    (0 1 0.1)
+);
+
+blocks
+(
+    hex (0 1 2 3 6 7 8 9) (200 4 1) simpleGrading (1 1 1)
+    hex (3 2 4 5 9 8 10 11) (200 36 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 6 9 3)
+            (3 9 11 5)
+        );
+    }
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (1 2 8 7)
+        );
+    }
+    bottomWall
+    {
+        type wall;
+        faces
+        (
+            (0 1 7 6)
+        );
+    }
+    endWall
+    {
+        type wall;
+        faces
+        (
+            (2 4 10 8)
+        );
+    }
+    top
+    {
+        type patch;
+        faces
+        (
+            (5 11 10 4)
+        );
+    }
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (0 3 2 1)
+            (6 7 8 9)
+            (3 5 4 2)
+            (9 8 10 11)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/polyMesh/boundary
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/polyMesh/boundary
@ -0,0 +1,61 @@
+/*--------------------------------*- 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      binary;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+6
+(
+    inlet
+    {
+        type            patch;
+        nFaces          40;
+        startFace       15760;
+    }
+    outlet
+    {
+        type            patch;
+        nFaces          4;
+        startFace       15800;
+    }
+    bottomWall
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          200;
+        startFace       15804;
+    }
+    endWall
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          36;
+        startFace       16004;
+    }
+    top
+    {
+        type            patch;
+        nFaces          200;
+        startFace       16040;
+    }
+    frontAndBack
+    {
+        type            empty;
+        inGroups        1(empty);
+        nFaces          16000;
+        startFace       16240;
+    }
+)
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/transportProperties
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/transportProperties
@ -0,0 +1,57 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+muc             muc [ 1 -1 -1 0 0 0 0 ] 0.00178;
+
+plasticViscosityCoeff plasticViscosityCoeff [ 1 -1 -1 0 0 0 0 ] 0.00023143;
+
+plasticViscosityExponent plasticViscosityExponent [ 0 0 0 0 0 0 0 ] 179.26;
+
+BinghamPlastic  on;
+
+yieldStressCoeff yieldStressCoeff [ 1 -1 -2 0 0 0 0 ] 0.00042189;
+
+yieldStressExponent yieldStressExponent [ 0 0 0 0 0 0 0 ] 1050.8;
+
+yieldStressOffset yieldStressOffset [ 0 0 0 0 0 0 0 ] 0;
+
+muMax           muMax [ 1 -1 -1 0 0 0 0 ] 10.0;
+
+rhoc            rhoc [ 1 -3 0 0 0 0 0 ] 996;
+
+rhod            rhod [ 1 -3 0 0 0 0 0 ] 1996;
+
+VdjModel        simple;
+
+simpleCoeffs
+{
+    V0              V0 [ 0 1 -1 0 0 0 0 ] ( 0 -0.002198 0 );
+    a               a [ 0 0 0 0 0 0 0 ] 285.84;
+    a1              a1 [ 0 0 0 0 0 0 0 ] 0;
+    alphaMin        alphaMin [ 0 0 0 0 0 0 0 ] 0;
+}
+
+generalCoeffs
+{
+    V0              V0 [ 0 1 -1 0 0 0 0 ] ( 0 -0.0018 0 );
+    a               a [ 0 0 0 0 0 0 0 ] 1e-05;
+    a1              a1 [ 0 0 0 0 0 0 0 ] 0.1;
+    alphaMin        alphaMin [ 0 0 0 0 0 0 0 ] 2e-05;
+}
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/system/controlDict
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/system/controlDict
@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     driftFluxFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         6400;
+
+deltaT          0.1;
+
+writeControl    adjustableRunTime;
+
+writeInterval   50;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  on;
+
+maxCo           1;
+
+maxDeltaT       1;
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSchemes
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSchemes
@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default             none;
+
+    div(rhoPhi,U)       Gauss linearUpwind grad(U);
+    div(phiVdj,Vdj)     Gauss linear;
+    div(phi,alpha)      Gauss vanLeer;
+    div(phirb,alpha)    Gauss vanLeer;
+    div(rhoPhi,k)       Gauss limitedLinear 1;
+    div(rhoPhi,epsilon) Gauss limitedLinear 1;
+
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p_rgh;
+    alpha;
+}
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSolution
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSolution
@ -0,0 +1,86 @@
+/*--------------------------------*- 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
+{
+    "alpha.*"
+    {
+        nAlphaCorr      1;
+        nAlphaSubCycles 4;
+
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-6;
+        relTol          0;
+        minIter         1;
+    }
+
+    p_rgh
+    {
+        solver          GAMG;
+        tolerance       1e-7;
+        relTol          0.01;
+        smoother        GaussSeidel;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        relTol          0;
+    }
+
+    "(U|k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-7;
+        relTol          0.1;
+        minIter         1;
+    }
+
+    "(U|k|epsilon)Final"
+    {
+        $k;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        "U.*"           1;
+        "k.*"           1;
+        "epsilon.*"     1;
+    }
+}
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/U
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/U
@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+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;
+    }
+}
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/alpha
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/alpha
@ -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      alpha;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.001;
+
+boundaryField
+{
+    rotor
+    {
+        type            zeroGradient;
+    }
+
+    stator
+    {
+        type            zeroGradient;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/epsilon
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/epsilon
@ -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      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 7e-3;
+
+boundaryField
+{
+    rotor
+    {
+        type            zeroGradient;
+    }
+
+    stator
+    {
+        type            zeroGradient;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/k
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/k
@ -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      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.01;
+
+boundaryField
+{
+    rotor
+    {
+        type            zeroGradient;
+    }
+
+    stator
+    {
+        type            zeroGradient;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/p_rgh
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/p_rgh
@ -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      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    rotor
+    {
+        type            zeroGradient;
+    }
+
+    stator
+    {
+        type            zeroGradient;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/Allrun
+++ b/tutorials/multiphase/driftFluxFoam/ras/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
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/RASProperties
+++ b/tutorials/multiphase/driftFluxFoam/ras/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      on;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/g
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/g
@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 0 0);
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/polyMesh/blockMeshDict.m4
+++ b/tutorials/multiphase/driftFluxFoam/ras/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)
+    )
+);
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/polyMesh/boundary
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/polyMesh/boundary
@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    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;
+        inGroups        1(wall);
+        nFaces          192;
+        startFace       5952;
+    }
+    stator
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          192;
+        startFace       6144;
+    }
+    front
+    {
+        type            empty;
+        inGroups        1(empty);
+        nFaces          3072;
+        startFace       6336;
+    }
+    back
+    {
+        type            empty;
+        inGroups        1(empty);
+        nFaces          3072;
+        startFace       9408;
+    }
+)
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/transportProperties
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/transportProperties
@ -0,0 +1,57 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+muc             muc [ 1 -1 -1 0 0 0 0 ] 0.00178;
+
+plasticViscosityCoeff plasticViscosityCoeff [ 1 -1 -1 0 0 0 0 ] 0.00023143;
+
+plasticViscosityExponent plasticViscosityExponent [ 0 0 0 0 0 0 0 ] 179.26;
+
+BinghamPlastic  on;
+
+yieldStressCoeff yieldStressCoeff [ 1 -1 -2 0 0 0 0 ] 0.00042189;
+
+yieldStressExponent yieldStressExponent [ 0 0 0 0 0 0 0 ] 1050.8;
+
+yieldStressOffset yieldStressOffset [ 0 0 0 0 0 0 0 ] 0;
+
+muMax           muMax [ 1 -1 -1 0 0 0 0 ] 10.0;
+
+rhoc            rhoc [ 1 -3 0 0 0 0 0 ] 996;
+
+rhod            rhod [ 1 -3 0 0 0 0 0 ] 1996;
+
+VdjModel        simple;
+
+simpleCoeffs
+{
+    V0              V0 [ 0 1 -1 0 0 0 0 ] ( 0 -0.002198 0 );
+    a               a [ 0 0 0 0 0 0 0 ] 285.84;
+    a1              a1 [ 0 0 0 0 0 0 0 ] 0;
+    alphaMin        alphaMin [ 0 0 0 0 0 0 0 ] 0;
+}
+
+generalCoeffs
+{
+    V0              V0 [ 0 1 -1 0 0 0 0 ] ( 0 -0.0018 0 );
+    a               a [ 0 0 0 0 0 0 0 ] 1e-05;
+    a1              a1 [ 0 0 0 0 0 0 0 ] 0.1;
+    alphaMin        alphaMin [ 0 0 0 0 0 0 0 ] 2e-05;
+}
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/makeMesh
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/makeMesh
@ -0,0 +1,6 @@
+#!/bin/sh
+
+m4 < constant/polyMesh/blockMeshDict.m4 > constant/polyMesh/blockMeshDict
+blockMesh
+topoSet
+setsToZones -noFlipMap
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/controlDict
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/controlDict
@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     driftFluxFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         10;
+
+deltaT          1e-3;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.25;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           0.5;
+maxAlphaCo      0.5;
+
+maxDeltaT       0.05;
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvOptions
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvOptions
@ -0,0 +1,37 @@
+/*--------------------------------*- 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      fvOptions;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+MRF1
+{
+    type            MRFSource;
+    active          yes;
+    selectionMode   cellZone;
+    cellZone        rotor;
+
+    MRFSourceCoeffs
+    {
+        // Fixed patches (by default they 'move' with the MRF zone)
+        nonRotatingPatches ();
+
+        origin    (0 0 0);
+        axis      (0 0 1);
+        omega     constant 6.2831853;
+    }
+}
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSchemes
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSchemes
@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default             none;
+
+    div(rhoPhi,U)       Gauss linearUpwind grad(U);
+    div(phiVdj,Vdj)     Gauss linear;
+    div(phi,alpha)      Gauss vanLeer;
+    div(phirb,alpha)    Gauss vanLeer;
+    div(rhoPhi,k)       Gauss limitedLinear 1;
+    div(rhoPhi,epsilon) Gauss limitedLinear 1;
+
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p_rgh           ;
+    alpha           ;
+}
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSolution
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSolution
@ -0,0 +1,83 @@
+/*--------------------------------*- 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
+{
+    "alpha.*"
+    {
+        nAlphaCorr      1;
+        nAlphaSubCycles 3;
+
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-6;
+        relTol          0;
+    }
+
+    "rho.*"
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-8;
+        relTol          0;
+    }
+
+    p_rgh
+    {
+        solver          GAMG;
+        tolerance       1e-6;
+        relTol          0.05;
+        smoother        GaussSeidel;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        relTol          0;
+    }
+
+    "(U|k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-6;
+        relTol          0.1;
+    }
+
+    "(U|k|epsilon)Final"
+    {
+        $k;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    momentumPredictor yes;
+    nCorrectors     3;
+    nNonOrthogonalCorrectors 0;
+
+    pRefCell        0;
+    pRefValue       0;
+}
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/setFieldsDict
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/setFieldsDict
@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha.water 0
+);
+
+regions
+(
+    boxToCell
+    {
+        box (0 0 -1) (1 1 1);
+        fieldValues
+        (
+            volScalarFieldValue alpha.water 1
+        );
+    }
+);
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/topoSetDict
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/topoSetDict
@ -0,0 +1,32 @@
+/*--------------------------------*- 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      topoSetDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    {
+        name    rotor;
+        type    cellSet;
+        action  new;
+        source  zoneToCell;
+        sourceInfo
+        {
+            name rotor;
+        }
+    }
+);
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/U
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/U
@ -0,0 +1,144 @@
+/*--------------------------------*- 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      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    SYMP3
+    {
+        type            slip;
+    }
+
+    INLE1
+    {
+        type            fixedValue;
+        value           uniform (0.1315 0 0);
+    }
+
+    OUTL9
+    {
+        type            fixedValue;
+        value           uniform (0 0.0177 0);
+    }
+
+    OUTL10
+    {
+        type            fixedValue;
+        value           uniform (0 0.0177 0);
+    }
+
+    OUTL11
+    {
+        type            fixedValue;
+        value           uniform (0 0.0177 0);
+    }
+
+    OUTL12
+    {
+        type            fixedValue;
+        value           uniform (0 0.0177 0);
+    }
+
+    WALL6
+    {
+        type            translatingWallVelocity;
+        U               (-0.003 0 0);
+        value           uniform (-0.003 0 0);
+    }
+
+    WALL8
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    WALL61
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    WALL62
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    WALL63
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    WALL64
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    WALL65
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    WALL66
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    WALL67
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    WALL68
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    WALL69
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    WALL7
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    WALL70
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    OUTL15
+    {
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/alpha
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/alpha
@ -0,0 +1,126 @@
+/*--------------------------------*- 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      alpha;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    SYMP3
+    {
+        type            zeroGradient;
+    }
+
+    INLE1
+    {
+        type            fixedValue;
+        value           uniform 0.002;
+    }
+
+    OUTL9
+    {
+        type            zeroGradient;
+    }
+
+    OUTL10
+    {
+        type            zeroGradient;
+    }
+
+    OUTL11
+    {
+        type            zeroGradient;
+    }
+
+    OUTL12
+    {
+        type            zeroGradient;
+    }
+
+    WALL6
+    {
+        type            zeroGradient;
+    }
+
+    WALL8
+    {
+        type            zeroGradient;
+    }
+
+    WALL61
+    {
+        type            zeroGradient;
+    }
+
+    WALL62
+    {
+        type            zeroGradient;
+    }
+
+    WALL63
+    {
+        type            zeroGradient;
+    }
+
+    WALL64
+    {
+        type            zeroGradient;
+    }
+
+    WALL65
+    {
+        type            zeroGradient;
+    }
+
+    WALL66
+    {
+        type            zeroGradient;
+    }
+
+    WALL67
+    {
+        type            zeroGradient;
+    }
+
+    WALL68
+    {
+        type            zeroGradient;
+    }
+
+    WALL69
+    {
+        type            zeroGradient;
+    }
+
+    WALL7
+    {
+        type            zeroGradient;
+    }
+
+    WALL70
+    {
+        type            zeroGradient;
+    }
+
+    OUTL15
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/epsilon
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/epsilon
@ -0,0 +1,126 @@
+/*--------------------------------*- 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      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 1.973e-07;
+
+boundaryField
+{
+    SYMP3
+    {
+        type            zeroGradient;
+    }
+
+    INLE1
+    {
+        type            fixedValue;
+        value           uniform 1.973e-07;
+    }
+
+    OUTL9
+    {
+        type            zeroGradient;
+    }
+
+    OUTL10
+    {
+        type            zeroGradient;
+    }
+
+    OUTL11
+    {
+        type            zeroGradient;
+    }
+
+    OUTL12
+    {
+        type            zeroGradient;
+    }
+
+    WALL6
+    {
+        type            zeroGradient;
+    }
+
+    WALL8
+    {
+        type            zeroGradient;
+    }
+
+    WALL61
+    {
+        type            zeroGradient;
+    }
+
+    WALL62
+    {
+        type            zeroGradient;
+    }
+
+    WALL63
+    {
+        type            zeroGradient;
+    }
+
+    WALL64
+    {
+        type            zeroGradient;
+    }
+
+    WALL65
+    {
+        type            zeroGradient;
+    }
+
+    WALL66
+    {
+        type            zeroGradient;
+    }
+
+    WALL67
+    {
+        type            zeroGradient;
+    }
+
+    WALL68
+    {
+        type            zeroGradient;
+    }
+
+    WALL69
+    {
+        type            zeroGradient;
+    }
+
+    WALL7
+    {
+        type            zeroGradient;
+    }
+
+    WALL70
+    {
+        type            zeroGradient;
+    }
+
+    OUTL15
+    {
+        type            inletOutlet;
+        inletValue      uniform 1.973e-07;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/k
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/k
@ -0,0 +1,126 @@
+/*--------------------------------*- 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      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.000259;
+
+boundaryField
+{
+    SYMP3
+    {
+        type            zeroGradient;
+    }
+
+    INLE1
+    {
+        type            fixedValue;
+        value           uniform 0.000259;
+    }
+
+    OUTL9
+    {
+        type            zeroGradient;
+    }
+
+    OUTL10
+    {
+        type            zeroGradient;
+    }
+
+    OUTL11
+    {
+        type            zeroGradient;
+    }
+
+    OUTL12
+    {
+        type            zeroGradient;
+    }
+
+    WALL6
+    {
+        type            zeroGradient;
+    }
+
+    WALL8
+    {
+        type            zeroGradient;
+    }
+
+    WALL61
+    {
+        type            zeroGradient;
+    }
+
+    WALL62
+    {
+        type            zeroGradient;
+    }
+
+    WALL63
+    {
+        type            zeroGradient;
+    }
+
+    WALL64
+    {
+        type            zeroGradient;
+    }
+
+    WALL65
+    {
+        type            zeroGradient;
+    }
+
+    WALL66
+    {
+        type            zeroGradient;
+    }
+
+    WALL67
+    {
+        type            zeroGradient;
+    }
+
+    WALL68
+    {
+        type            zeroGradient;
+    }
+
+    WALL69
+    {
+        type            zeroGradient;
+    }
+
+    WALL7
+    {
+        type            zeroGradient;
+    }
+
+    WALL70
+    {
+        type            zeroGradient;
+    }
+
+    OUTL15
+    {
+        type            inletOutlet;
+        inletValue      uniform 0.000259;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/p_rgh
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/p_rgh
@ -0,0 +1,144 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    SYMP3
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    INLE1
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    OUTL9
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    OUTL10
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    OUTL11
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    OUTL12
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL6
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL8
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL61
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL62
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL63
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL64
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL65
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL66
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL67
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL68
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL69
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL7
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    WALL70
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+
+    OUTL15
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+}
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/Allclean
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/Allclean
@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Clean time directories only
+
+rm -rf *[1-9]*
+rm -f log.* 2>/dev/null
+
+# ----------------------------------------------------------------- end-of-file
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/Allrun
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/Allrun
@ -0,0 +1,11 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+application=`getApplication`
+
+runApplication $application
+
+# ----------------------------------------------------------------- end-of-file
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/RASProperties
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/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      on;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/g
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/g
@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 -9.81 0 );
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/polyMesh/boundary
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/polyMesh/boundary
@ -0,0 +1,160 @@
+/*--------------------------------*- 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;
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+(
+    SYMP3
+    {
+        type            patch;
+        startFace       53708;
+        nFaces          3191;
+    }
+
+    INLE1
+    {
+        type            patch;
+        startFace       56899;
+        nFaces          45;
+    }
+
+    OUTL9
+    {
+        type            patch;
+        startFace       56944;
+        nFaces          11;
+    }
+
+    OUTL10
+    {
+        type            patch;
+        startFace       56955;
+        nFaces          11;
+    }
+
+    OUTL11
+    {
+        type            patch;
+        startFace       56966;
+        nFaces          11;
+    }
+
+    OUTL12
+    {
+        type            patch;
+        startFace       56977;
+        nFaces          11;
+    }
+
+    WALL6
+    {
+        type            wall;
+        startFace       56988;
+        nFaces          1474;
+    }
+
+    WALL8
+    {
+        type            wall;
+        startFace       58462;
+        nFaces          1640;
+    }
+
+    WALL61
+    {
+        type            wall;
+        startFace       60102;
+        nFaces          44;
+    }
+
+    WALL62
+    {
+        type            wall;
+        startFace       60146;
+        nFaces          78;
+    }
+
+    WALL63
+    {
+        type            wall;
+        startFace       60224;
+        nFaces          255;
+    }
+
+    WALL64
+    {
+        type            wall;
+        startFace       60479;
+        nFaces          157;
+    }
+
+    WALL65
+    {
+        type            wall;
+        startFace       60636;
+        nFaces          35;
+    }
+
+    WALL66
+    {
+        type            wall;
+        startFace       60671;
+        nFaces          60;
+    }
+
+    WALL67
+    {
+        type            wall;
+        startFace       60731;
+        nFaces          55;
+    }
+
+    WALL68
+    {
+        type            wall;
+        startFace       60786;
+        nFaces          116;
+    }
+
+    WALL69
+    {
+        type            wall;
+        startFace       60902;
+        nFaces          24;
+    }
+
+    WALL7
+    {
+        type            wall;
+        startFace       60926;
+        nFaces          121;
+    }
+
+    WALL70
+    {
+        type            wall;
+        startFace       61047;
+        nFaces          182;
+    }
+
+    OUTL15
+    {
+        type            patch;
+        startFace       61229;
+        nFaces          14;
+    }
+
+)
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/polyMesh/cells.gz
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/polyMesh/cells.gz
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/polyMesh/faces.gz
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/polyMesh/faces.gz
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/polyMesh/neighbour.gz
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/polyMesh/neighbour.gz
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/polyMesh/owner.gz
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/polyMesh/owner.gz
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/polyMesh/points.gz
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/polyMesh/points.gz
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/transportProperties
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/transportProperties
@ -0,0 +1,57 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+muc             muc [ 1 -1 -1 0 0 0 0 ] 0.00178;
+
+plasticViscosityCoeff plasticViscosityCoeff [ 1 -1 -1 0 0 0 0 ] 0.00023143;
+
+plasticViscosityExponent plasticViscosityExponent [ 0 0 0 0 0 0 0 ] 0.17926;
+
+BinghamPlastic  on;
+
+yieldStressCoeff yieldStressCoeff [ 1 -1 -2 0 0 0 0 ] 5.5469e-07;
+
+yieldStressExponent yieldStressExponent [ 0 0 0 0 0 0 0 ] 95.25;
+
+yieldStressOffset yieldStressOffset [ 0 0 0 0 0 0 0 ] 0;
+
+muMax           muMax [ 1 -1 -1 0 0 0 0 ] 10.0;
+
+rhoc            rhoc [ 1 -3 0 0 0 0 0 ] 1000;
+
+rhod            rhod [ 1 -3 0 0 0 0 0 ] 1042;
+
+VdjModel        simple;
+
+simpleCoeffs
+{
+    V0              V0 [ 0 1 -1 0 0 0 0 ] ( 0 -0.002198 0 );
+    a               a [ 0 0 0 0 0 0 0 ] 8.84;
+    a1              a1 [ 0 0 0 0 0 0 0 ] 0;
+    alphaMin        alphaMin [ 0 0 0 0 0 0 0 ] 0;
+}
+
+generalCoeffs
+{
+    V0              V0 [ 0 1 -1 0 0 0 0 ] ( 0 -0.0018 0 );
+    a               a [ 0 0 0 0 0 0 0 ] 1e-05;
+    a1              a1 [ 0 0 0 0 0 0 0 ] 0.1;
+    alphaMin        alphaMin [ 0 0 0 0 0 0 0 ] 2e-05;
+}
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/controlDict
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/controlDict
@ -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       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     driftFluxFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         8000;
+
+deltaT          0.1;
+
+writeControl    runTime;
+
+writeInterval   50;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression compressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSchemes
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSchemes
@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+
+    div(rhoPhi,U)       Gauss linearUpwind grad(U);
+    div(phiVdj,Vdj)     Gauss linear;
+    div(phi,alpha)      Gauss vanLeer;
+    div(phirb,alpha)    Gauss vanLeer;
+    div(rhoPhi,k)       Gauss limitedLinear 1;
+    div(rhoPhi,epsilon) Gauss limitedLinear 1;
+
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p_rgh;
+    alpha;
+}
+
+
+// ************************************************************************* //
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSolution
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSolution
@ -0,0 +1,84 @@
+/*--------------------------------*- 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
+{
+    "alpha.*"
+    {
+        nAlphaCorr      1;
+        nAlphaSubCycles 4;
+
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-6;
+        relTol          0;
+        minIter         1;
+    }
+
+    p_rgh
+    {
+        solver          GAMG;
+        tolerance       1e-7;
+        relTol          0.01;
+        smoother        GaussSeidel;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        relTol          0;
+    }
+
+    "(U|k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-7;
+        relTol          0.1;
+        minIter         1;
+    }
+
+    "(U|k|epsilon)Final"
+    {
+        $k;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        ".*"       1;
+    }
+}
+
+
+// ************************************************************************* //