diff --git a/applications/solvers/multiphase/driftFluxFoam/UEqn.H b/applications/solvers/multiphase/driftFluxFoam/UEqn.H
index b6d8756486..010417a0f9 100644
--- a/applications/solvers/multiphase/driftFluxFoam/UEqn.H
+++ b/applications/solvers/multiphase/driftFluxFoam/UEqn.H
@@ -6,7 +6,7 @@
       + fvm::div(rhoPhi, U)
       + fvc::div
         (
-            (alpha/(scalar(1.001) - alpha))*((rhoc*rhod)/rho)*Vdj*Vdj,
+            (alpha1/(scalar(1.001) - alpha1))*((rho2*rho1)/rho)*Vdj*Vdj,
             "div(phiVdj,Vdj)"
         )
       - fvm::laplacian(muEff, U)
diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaControls.H b/applications/solvers/multiphase/driftFluxFoam/alphaControls.H
deleted file mode 100644
index 5a3d10bd71..0000000000
--- a/applications/solvers/multiphase/driftFluxFoam/alphaControls.H
+++ /dev/null
@@ -1,4 +0,0 @@
-const dictionary& alphaControls = mesh.solverDict(alpha.name());
-
-label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
-label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
index 89c139dfd0..817a0c9124 100644
--- a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
@@ -2,30 +2,118 @@
     word alphaScheme("div(phi,alpha)");
     word alpharScheme("div(phirb,alpha)");
 
+    if (MULESCorr)
+    {
+        fvScalarMatrix alpha1Eqn
+        (
+            #ifdef LTSSOLVE
+            fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
+            #else
+            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+            #endif
+          + fv::gaussConvectionScheme<scalar>
+            (
+                mesh,
+                phi,
+                upwind<scalar>(mesh, phi)
+            ).fvmDiv(phi, alpha1)
+          - fv::gaussLaplacianScheme<scalar, scalar>
+            (
+                mesh,
+                linear<scalar>(mesh),
+                fv::uncorrectedSnGrad<scalar>(mesh)
+            ).fvmLaplacian(fvc::interpolate(mut/rho), alpha1)
+        );
+
+        alpha1Eqn.solve();
+
+        Info<< "Phase-1 volume fraction = "
+            << alpha1.weightedAverage(mesh.Vsc()).value()
+            << "  Min(alpha1) = " << min(alpha1).value()
+            << "  Max(alpha1) = " << max(alpha1).value()
+            << endl;
+
+        tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
+        phiAlpha = tphiAlphaUD();
+
+        if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
+        {
+            Info<< "Applying the previous iteration correction flux" << endl;
+            #ifdef LTSSOLVE
+            MULES::LTScorrect(alpha1, phiAlpha, tphiAlphaCorr0(), 1, 0);
+            #else
+            MULES::correct(alpha1, phiAlpha, tphiAlphaCorr0(), 1, 0);
+            #endif
+
+            phiAlpha += tphiAlphaCorr0();
+        }
+
+        // Cache the upwind-flux
+        tphiAlphaCorr0 = tphiAlphaUD;
+    }
+
     for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
     {
-        phiAlpha =
+        tmp<surfaceScalarField> tphiAlphaUn
         (
             fvc::flux
             (
                 phi,
-                alpha,
+                alpha1,
                 alphaScheme
             )
           + fvc::flux
             (
                 phir,
-                alpha,
+                alpha1,
                 alpharScheme
             )
         );
 
-        MULES::explicitSolve(alpha, phi, phiAlpha, 1, 0);
+        if (MULESCorr)
+        {
+            tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - phiAlpha);
+            volScalarField alpha10(alpha1);
+
+            #ifdef LTSSOLVE
+            MULES::LTScorrect(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
+            #else
+            MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
+            #endif
+
+            // Under-relax the correction for all but the 1st corrector
+            if (aCorr == 0)
+            {
+                phiAlpha += tphiAlphaCorr();
+            }
+            else
+            {
+                alpha1 = 0.5*alpha1 + 0.5*alpha10;
+                phiAlpha += 0.5*tphiAlphaCorr();
+            }
+        }
+        else
+        {
+            phiAlpha = tphiAlphaUn;
+
+            #ifdef LTSSOLVE
+            MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
+            #else
+            MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
+            #endif
+        }
     }
 
+    if (alphaApplyPrevCorr && MULESCorr)
+    {
+        tphiAlphaCorr0 = phiAlpha - tphiAlphaCorr0;
+    }
+
+    alpha2 = 1.0 - alpha1;
+
     Info<< "Phase-1 volume fraction = "
-        << alpha.weightedAverage(mesh.Vsc()).value()
-        << "  Min(alpha) = " << min(alpha).value()
-        << "  Max(alpha) = " << max(alpha).value()
+        << alpha1.weightedAverage(mesh.Vsc()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
         << endl;
 }
diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
index d0c08ba9f2..f134702c68 100644
--- a/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
@@ -13,7 +13,7 @@
 
     surfaceScalarField phir
     (
-        rhoc*(mesh.Sf() & fvc::interpolate(Vdj/rho))
+        rho2*(mesh.Sf() & fvc::interpolate(Vdj/rho))
     );
 
     if (nAlphaSubCycles > 1)
@@ -33,7 +33,7 @@
 
         for
         (
-            subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
+            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
             !(++alphaSubCycle).end();
         )
         {
@@ -50,24 +50,26 @@
 
     // Apply the diffusion term separately to allow implicit solution
     // and boundedness of the explicit advection
+    if (!MULESCorr)
     {
-        fvScalarMatrix alphaEqn
+        fvScalarMatrix alpha1Eqn
         (
-            fvm::ddt(alpha) - fvc::ddt(alpha)
-          - fvm::laplacian(mut/rho, alpha)
+            fvm::ddt(alpha1) - fvc::ddt(alpha1)
+          - fvm::laplacian(mut/rho, alpha1)
         );
 
-        alphaEqn.solve();
+        alpha1Eqn.solve(mesh.solver("alpha1Diffusion"));
 
-        phiAlpha += alphaEqn.flux();
+        phiAlpha += alpha1Eqn.flux();
+        alpha2 = 1.0 - alpha1;
+
+        Info<< "Phase-1 volume fraction = "
+            << alpha1.weightedAverage(mesh.Vsc()).value()
+            << "  Min(alpha1) = " << min(alpha1).value()
+            << "  Max(alpha1) = " << max(alpha1).value()
+            << endl;
     }
 
-    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;
+    rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
+    rho == alpha1*rho1 + alpha2*rho2;
 }
diff --git a/applications/solvers/multiphase/driftFluxFoam/calcVdj.H b/applications/solvers/multiphase/driftFluxFoam/calcVdj.H
index 7a921ef708..cf02fafdd5 100644
--- a/applications/solvers/multiphase/driftFluxFoam/calcVdj.H
+++ b/applications/solvers/multiphase/driftFluxFoam/calcVdj.H
@@ -2,13 +2,13 @@ if (VdjModel == "general")
 {
     Vdj = V0*
     (
-        exp(-a*max(alpha - alphaMin, scalar(0)))
-      - exp(-a1*max(alpha - alphaMin, scalar(0)))
+        exp(-a*max(alpha1 - alphaMin, scalar(0)))
+      - exp(-a1*max(alpha1 - alphaMin, scalar(0)))
     );
 }
 else if (VdjModel == "simple")
 {
-    Vdj = V0*pow(10.0, -a*max(alpha, scalar(0)));
+    Vdj = V0*pow(10.0, -a*max(alpha1, scalar(0)));
 }
 else
 {
diff --git a/applications/solvers/multiphase/driftFluxFoam/correctViscosity.H b/applications/solvers/multiphase/driftFluxFoam/correctViscosity.H
index 0bf086a4f4..38c6f804d0 100644
--- a/applications/solvers/multiphase/driftFluxFoam/correctViscosity.H
+++ b/applications/solvers/multiphase/driftFluxFoam/correctViscosity.H
@@ -4,7 +4,7 @@
         (
             plasticViscosityCoeff,
             plasticViscosityExponent,
-            alpha
+            alpha1
         );
 
     if (BinghamPlastic)
@@ -14,7 +14,7 @@
             yieldStressCoeff,
             yieldStressExponent,
             yieldStressOffset,
-            alpha
+            alpha1
         );
 
         mul =
diff --git a/applications/solvers/multiphase/driftFluxFoam/createFields.H b/applications/solvers/multiphase/driftFluxFoam/createFields.H
index 950235cfe5..01917ba26d 100644
--- a/applications/solvers/multiphase/driftFluxFoam/createFields.H
+++ b/applications/solvers/multiphase/driftFluxFoam/createFields.H
@@ -12,12 +12,12 @@
         mesh
     );
 
-    Info<< "Reading field alpha\n" << endl;
-    volScalarField alpha
+    Info<< "Reading field alpha1\n" << endl;
+    volScalarField alpha1
     (
         IOobject
         (
-            "alpha",
+            "alpha1",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -26,6 +26,8 @@
         mesh
     );
 
+    volScalarField alpha2("alpha2", scalar(1) - alpha1);
+
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -58,9 +60,8 @@
     );
 
 
-    dimensionedScalar rhoc(transportProperties.lookup("rhoc"));
-
-    dimensionedScalar rhod(transportProperties.lookup("rhod"));
+    dimensionedScalar rho1(transportProperties.lookup("rho1"));
+    dimensionedScalar rho2(transportProperties.lookup("rho2"));
 
     dimensionedScalar muc(transportProperties.lookup("muc"));
     dimensionedScalar muMax(transportProperties.lookup("muMax"));
@@ -102,7 +103,7 @@
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        alpha*rhod + (scalar(1) - alpha)*rhoc
+        alpha1*rho1 + alpha2*rho2
     );
     rho.oldTime();
 
@@ -136,7 +137,7 @@
         (
             plasticViscosityCoeff,
             plasticViscosityExponent,
-            alpha
+            alpha1
         )
     );
 
@@ -382,3 +383,6 @@
         );
         p_rgh = p - rho*gh;
     }
+
+
+    tmp<surfaceScalarField> tphiAlphaCorr0;
diff --git a/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C b/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
index d0aafc6eab..00e80087b8 100644
--- a/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
+++ b/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
@@ -45,6 +45,8 @@ Description
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
 #include "fixedFluxPressureFvPatchScalarField.H"
+#include "gaussLaplacianScheme.H"
+#include "uncorrectedSnGrad.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H
index d63f7e79b1..ac87a026ec 100644
--- a/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H
+++ b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H
@@ -3,7 +3,7 @@
     const scalar kappa_ = kappa.value();
     const scalar E_ = E.value();
     const scalar muc_ = muc.value();
-    const scalar nuc_ = muc_/rhoc.value();
+    const scalar nuc_ = muc_/rho2.value();
 
     const fvPatchList& patches = mesh.boundary();
 
diff --git a/applications/solvers/multiphase/driftFluxFoam/yieldStress.H b/applications/solvers/multiphase/driftFluxFoam/yieldStress.H
index eda55d2ecb..8d3fd8ee76 100644
--- a/applications/solvers/multiphase/driftFluxFoam/yieldStress.H
+++ b/applications/solvers/multiphase/driftFluxFoam/yieldStress.H
@@ -3,7 +3,7 @@ volScalarField yieldStress
     const dimensionedScalar& yieldStressCoeff,
     const dimensionedScalar& yieldStressExponent,
     const dimensionedScalar& yieldStressOffset,
-    const volScalarField& alpha
+    const volScalarField& alpha1
 )
 {
     tmp<volScalarField> tfld
@@ -13,7 +13,7 @@ volScalarField yieldStress
             pow
             (
                 10.0,
-                yieldStressExponent*(max(alpha, scalar(0)) + yieldStressOffset)
+                yieldStressExponent*(max(alpha1, scalar(0)) + yieldStressOffset)
             )
           - pow
             (
diff --git a/tutorials/multiphase/driftFluxFoam/ras/dahl/0/alpha b/tutorials/multiphase/driftFluxFoam/ras/dahl/0/alpha1
similarity index 98%
rename from tutorials/multiphase/driftFluxFoam/ras/dahl/0/alpha
rename to tutorials/multiphase/driftFluxFoam/ras/dahl/0/alpha1
index 807802b1d1..650cc6d14c 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/0/alpha
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/0/alpha1
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha;
+    object      alpha1;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/transportProperties b/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/transportProperties
index ee877fce58..e2079db0db 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/transportProperties
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/constant/transportProperties
@@ -31,9 +31,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;
+rho1            rho1 [ 1 -3 0 0 0 0 0 ] 1996;
+rho2            rho2 [ 1 -3 0 0 0 0 0 ] 996;
 
 VdjModel        simple;
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/dahl/system/controlDict b/tutorials/multiphase/driftFluxFoam/ras/dahl/system/controlDict
index 04cdc807cd..67e02247c9 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/system/controlDict
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/system/controlDict
@@ -47,7 +47,7 @@ runTimeModifiable yes;
 
 adjustTimeStep  on;
 
-maxCo           1;
+maxCo           5;
 
 maxDeltaT       1;
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSchemes b/tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSchemes
index 9f169f3c34..98e6fa16c9 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSchemes
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSchemes
@@ -32,7 +32,7 @@ divSchemes
     div(rhoPhi,U)       Gauss linearUpwind grad(U);
     div(phiVdj,Vdj)     Gauss linear;
     div(phi,alpha)      Gauss vanLeer;
-    div(phirb,alpha)    Gauss vanLeer;
+    div(phirb,alpha)    Gauss linear;
     div(rhoPhi,k)       Gauss limitedLinear 1;
     div(rhoPhi,epsilon) Gauss limitedLinear 1;
 
@@ -58,7 +58,7 @@ fluxRequired
 {
     default         no;
     p_rgh;
-    alpha;
+    alpha1;
 }
 
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSolution b/tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSolution
index c912040136..b78a4821e4 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSolution
+++ b/tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSolution
@@ -17,10 +17,14 @@ FoamFile
 
 solvers
 {
-    "alpha.*"
+    "alpha1.*"
     {
-        nAlphaCorr      1;
-        nAlphaSubCycles 4;
+        nAlphaCorr      2;
+        nAlphaSubCycles 1;
+
+        MULESCorr       yes;
+        nLimiterIter    3;
+        alphaApplyPrevCorr  yes;
 
         solver          smoothSolver;
         smoother        symGaussSeidel;
@@ -29,6 +33,15 @@ solvers
         minIter         1;
     }
 
+    alpha1Diffusion
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-6;
+        relTol          0;
+        minIter         1;
+    }
+
     p_rgh
     {
         solver          GAMG;
@@ -49,7 +62,8 @@ solvers
 
     "(U|k|epsilon)"
     {
-        solver          smoothSolver;
+        solver          PBiCG;
+        preconditioner  DILU;
         smoother        symGaussSeidel;
         tolerance       1e-7;
         relTol          0.1;
@@ -65,20 +79,16 @@ solvers
 
 PIMPLE
 {
-    nCorrectors     2;
+    momentumPredictor no;
+    nCorrectors     3;
     nNonOrthogonalCorrectors 0;
 }
 
 relaxationFactors
 {
-    fields
-    {
-    }
     equations
     {
-        "U.*"           1;
-        "k.*"           1;
-        "epsilon.*"     1;
+        ".*"           1;
     }
 }
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/alpha b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/alpha1
similarity index 97%
rename from tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/alpha
rename to tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/alpha1
index eb8492e6ca..22f0efce92 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/alpha
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/alpha1
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha;
+    object      alpha1;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/transportProperties b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/transportProperties
index ee877fce58..e2079db0db 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/transportProperties
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/transportProperties
@@ -31,9 +31,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;
+rho1            rho1 [ 1 -3 0 0 0 0 0 ] 1996;
+rho2            rho2 [ 1 -3 0 0 0 0 0 ] 996;
 
 VdjModel        simple;
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/controlDict b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/controlDict
index 25ba33a488..5f9bee1349 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/controlDict
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/controlDict
@@ -47,8 +47,7 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.5;
-maxAlphaCo      0.5;
+maxCo           1;
 
 maxDeltaT       0.05;
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSchemes b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSchemes
index 4fb64f81ae..98e6fa16c9 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSchemes
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSchemes
@@ -32,7 +32,7 @@ divSchemes
     div(rhoPhi,U)       Gauss linearUpwind grad(U);
     div(phiVdj,Vdj)     Gauss linear;
     div(phi,alpha)      Gauss vanLeer;
-    div(phirb,alpha)    Gauss vanLeer;
+    div(phirb,alpha)    Gauss linear;
     div(rhoPhi,k)       Gauss limitedLinear 1;
     div(rhoPhi,epsilon) Gauss limitedLinear 1;
 
@@ -57,8 +57,8 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p_rgh           ;
-    alpha           ;
+    p_rgh;
+    alpha1;
 }
 
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSolution b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSolution
index 3f5cff1ea0..3fa09bb4f5 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSolution
+++ b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSolution
@@ -17,15 +17,29 @@ FoamFile
 
 solvers
 {
-    "alpha.*"
+    "alpha1.*"
     {
-        nAlphaCorr      1;
-        nAlphaSubCycles 3;
+        nAlphaCorr      2;
+        nAlphaSubCycles 1;
+
+        MULESCorr       yes;
+        nLimiterIter    3;
+        alphaApplyPrevCorr  yes;
 
         solver          smoothSolver;
         smoother        symGaussSeidel;
         tolerance       1e-6;
         relTol          0;
+        minIter         1;
+    }
+
+    alpha1Diffusion
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-6;
+        relTol          0;
+        minIter         1;
     }
 
     "rho.*"
@@ -71,7 +85,7 @@ solvers
 
 PIMPLE
 {
-    momentumPredictor yes;
+    momentumPredictor no;
     nCorrectors     3;
     nNonOrthogonalCorrectors 0;
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/U b/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/U
index 1c824278e6..da99a43dd3 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/U
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/U
@@ -136,8 +136,8 @@ boundaryField
 
     OUTL15
     {
-        type            inletOutlet;
-        inletValue      uniform (0 0 0);
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
     }
 }
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/alpha b/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/alpha1
similarity index 98%
rename from tutorials/multiphase/driftFluxFoam/ras/tank3D/0/alpha
rename to tutorials/multiphase/driftFluxFoam/ras/tank3D/0/alpha1
index 67aa1fa666..1a3bd64e20 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/alpha
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/0/alpha1
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha;
+    object      alpha1;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/transportProperties b/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/transportProperties
index ac45d4a30b..8c92706eaa 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/transportProperties
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/transportProperties
@@ -31,9 +31,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 ] 1000;
-
-rhod            rhod [ 1 -3 0 0 0 0 0 ] 1042;
+rho1            rho1 [ 1 -3 0 0 0 0 0 ] 1042;
+rho2            rho2 [ 1 -3 0 0 0 0 0 ] 1000;
 
 VdjModel        simple;
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/controlDict b/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/controlDict
index 93e66698b0..dd41af85ba 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/controlDict
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/controlDict
@@ -25,7 +25,7 @@ stopAt          endTime;
 
 endTime         8000;
 
-deltaT          0.1;
+deltaT          0.5;
 
 writeControl    runTime;
 
@@ -33,11 +33,11 @@ writeInterval   50;
 
 purgeWrite      0;
 
-writeFormat     ascii;
+writeFormat     binary;
 
 writePrecision  6;
 
-writeCompression compressed;
+writeCompression uncompressed;
 
 timeFormat      general;
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSchemes b/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSchemes
index 0d47c3816e..8809a091d2 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSchemes
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSchemes
@@ -32,7 +32,7 @@ divSchemes
     div(rhoPhi,U)       Gauss linearUpwind grad(U);
     div(phiVdj,Vdj)     Gauss linear;
     div(phi,alpha)      Gauss vanLeer;
-    div(phirb,alpha)    Gauss vanLeer;
+    div(phirb,alpha)    Gauss linear;
     div(rhoPhi,k)       Gauss limitedLinear 1;
     div(rhoPhi,epsilon) Gauss limitedLinear 1;
 
@@ -58,7 +58,7 @@ fluxRequired
 {
     default         no;
     p_rgh;
-    alpha;
+    alpha1;
 }
 
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSolution b/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSolution
index 403f82dc40..f6f1b18e53 100644
--- a/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSolution
+++ b/tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSolution
@@ -17,10 +17,14 @@ FoamFile
 
 solvers
 {
-    "alpha.*"
+    "alpha1.*"
     {
-        nAlphaCorr      1;
-        nAlphaSubCycles 4;
+        nAlphaCorr      2;
+        nAlphaSubCycles 1;
+
+        MULESCorr       yes;
+        nLimiterIter    3;
+        alphaApplyPrevCorr  yes;
 
         solver          smoothSolver;
         smoother        symGaussSeidel;
@@ -29,6 +33,15 @@ solvers
         minIter         1;
     }
 
+    alpha1Diffusion
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-6;
+        relTol          0;
+        minIter         1;
+    }
+
     p_rgh
     {
         solver          GAMG;
@@ -65,15 +78,13 @@ solvers
 
 PIMPLE
 {
-    nCorrectors     2;
+    momentumPredictor no;
+    nCorrectors     3;
     nNonOrthogonalCorrectors 0;
 }
 
 relaxationFactors
 {
-    fields
-    {
-    }
     equations
     {
         ".*"       1;