From 177b54ae9437ced1cd8d50abb97167ba7a9a708f Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 24 Feb 2014 17:10:20 +0000
Subject: [PATCH] driftFluxFoam: Rationalise phase-fraction handling and
 implement semi-implicit MULES

---
 .../solvers/multiphase/driftFluxFoam/UEqn.H   |   2 +-
 .../multiphase/driftFluxFoam/alphaControls.H  |   4 -
 .../multiphase/driftFluxFoam/alphaEqn.H       | 102 ++++++++++++++++--
 .../driftFluxFoam/alphaEqnSubCycle.H          |  32 +++---
 .../multiphase/driftFluxFoam/calcVdj.H        |   6 +-
 .../driftFluxFoam/correctViscosity.H          |   4 +-
 .../multiphase/driftFluxFoam/createFields.H   |  20 ++--
 .../multiphase/driftFluxFoam/driftFluxFoam.C  |   2 +
 .../multiphase/driftFluxFoam/wallViscosity.H  |   2 +-
 .../multiphase/driftFluxFoam/yieldStress.H    |   4 +-
 .../ras/dahl/0/{alpha => alpha1}              |   2 +-
 .../ras/dahl/constant/transportProperties     |   5 +-
 .../driftFluxFoam/ras/dahl/system/controlDict |   2 +-
 .../driftFluxFoam/ras/dahl/system/fvSchemes   |   4 +-
 .../driftFluxFoam/ras/dahl/system/fvSolution  |  32 ++++--
 .../ras/mixerVessel2D/0/{alpha => alpha1}     |   2 +-
 .../constant/transportProperties              |   5 +-
 .../ras/mixerVessel2D/system/controlDict      |   3 +-
 .../ras/mixerVessel2D/system/fvSchemes        |   6 +-
 .../ras/mixerVessel2D/system/fvSolution       |  22 +++-
 .../multiphase/driftFluxFoam/ras/tank3D/0/U   |   4 +-
 .../ras/tank3D/0/{alpha => alpha1}            |   2 +-
 .../ras/tank3D/constant/transportProperties   |   5 +-
 .../ras/tank3D/system/controlDict             |   6 +-
 .../driftFluxFoam/ras/tank3D/system/fvSchemes |   4 +-
 .../ras/tank3D/system/fvSolution              |  25 +++--
 26 files changed, 215 insertions(+), 92 deletions(-)
 delete mode 100644 applications/solvers/multiphase/driftFluxFoam/alphaControls.H
 rename tutorials/multiphase/driftFluxFoam/ras/dahl/0/{alpha => alpha1} (98%)
 rename tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/{alpha => alpha1} (97%)
 rename tutorials/multiphase/driftFluxFoam/ras/tank3D/0/{alpha => alpha1} (98%)

diff --git a/applications/solvers/multiphase/driftFluxFoam/UEqn.H b/applications/solvers/multiphase/driftFluxFoam/UEqn.H
index b6d87564861..010417a0f90 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 5a3d10bd71d..00000000000
--- 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 89c139dfd08..817a0c91242 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 d0c08ba9f23..f134702c681 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 = "
-        << alpha.weightedAverage(mesh.Vsc()).value()
-        << "  Min(alpha) = " << min(alpha).value()
-        << "  Max(alpha) = " << max(alpha).value()
-        << endl;
+        Info<< "Phase-1 volume fraction = "
+            << alpha1.weightedAverage(mesh.Vsc()).value()
+            << "  Min(alpha1) = " << min(alpha1).value()
+            << "  Max(alpha1) = " << max(alpha1).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 7a921ef7084..cf02fafdd5c 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 0bf086a4f42..38c6f804d02 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 950235cfe55..01917ba26d1 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 d0aafc6eab3..00e80087b8c 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 d63f7e79b1f..ac87a026ece 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 eda55d2ecb9..8d3fd8ee76c 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 807802b1d15..650cc6d14c8 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 ee877fce583..e2079db0db3 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 04cdc807cdb..67e02247c9c 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 9f169f3c341..98e6fa16c96 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 c9120401361..b78a4821e45 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 eb8492e6ca9..22f0efce92c 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 ee877fce583..e2079db0db3 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 25ba33a488f..5f9bee13499 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 4fb64f81ae8..98e6fa16c96 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 3f5cff1ea0b..3fa09bb4f50 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 1c824278e69..da99a43dd39 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 67aa1fa666a..1a3bd64e205 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 ac45d4a30b2..8c92706eaa1 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 93e66698b0e..dd41af85ba9 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 0d47c3816e8..8809a091d24 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 403f82dc409..f6f1b18e531 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;
-- 
GitLab