driftFluxFoam: Rationalise phase-fraction handling and implement semi-implicit MULES

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)

applications/solvers/multiphase/driftFluxFoam/alphaControls.H

@@ -1,4 +0,0 @@
-const dictionary& alphaControls = mesh.solverDict(alpha.name());
-
-label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
-label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));

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()
+        << alpha1.weightedAverage(mesh.Vsc()).value()
-        << "  Min(alpha) = " << min(alpha).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
-        << "  Max(alpha) = " << max(alpha).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
         << endl;
 }

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::ddt(alpha1) - fvc::ddt(alpha1)
-          - fvm::laplacian(mut/rho, alpha)
+          - 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 = "
+    rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
-        << alpha.weightedAverage(mesh.Vsc()).value()
+    rho == alpha1*rho1 + alpha2*rho2;
-        << "  Min(alpha) = " << min(alpha).value()
-        << "  Max(alpha) = " << max(alpha).value()
-        << endl;
-
-    rhoPhi = phiAlpha*(rhod - rhoc) + phi*rhoc;
-    rho == alpha*rhod + (scalar(1) - alpha)*rhoc;
 }

applications/solvers/multiphase/driftFluxFoam/calcVdj.H

@@ -2,13 +2,13 @@ if (VdjModel == "general")
 {
     Vdj = V0*
     (
-        exp(-a*max(alpha - alphaMin, scalar(0)))
+        exp(-a*max(alpha1 - alphaMin, scalar(0)))
-      - exp(-a1*max(alpha - 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
 {

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 =

applications/solvers/multiphase/driftFluxFoam/createFields.H

@@ -12,12 +12,12 @@
         mesh
     );
 
-    Info<< "Reading field alpha\n" << endl;
+    Info<< "Reading field alpha1\n" << endl;
-    volScalarField alpha
+    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 rho1(transportProperties.lookup("rho1"));
-
+    dimensionedScalar rho2(transportProperties.lookup("rho2"));
-    dimensionedScalar rhod(transportProperties.lookup("rhod"));
 
     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;

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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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();
 

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
             (

tutorials/multiphase/driftFluxFoam/ras/dahl/0/alpha1

@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha;
+    object      alpha1;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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;
+rho1            rho1 [ 1 -3 0 0 0 0 0 ] 1996;
-
+rho2            rho2 [ 1 -3 0 0 0 0 0 ] 996;
-rhod            rhod [ 1 -3 0 0 0 0 0 ] 1996;
 
 VdjModel        simple;
 

tutorials/multiphase/driftFluxFoam/ras/dahl/system/controlDict

@@ -47,7 +47,7 @@ runTimeModifiable yes;
 
 adjustTimeStep  on;
 
-maxCo           1;
+maxCo           5;
 
 maxDeltaT       1;
 

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;
 }
 
 

tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSolution

@@ -17,10 +17,14 @@ FoamFile
 
 solvers
 {
-    "alpha.*"
+    "alpha1.*"
     {
-        nAlphaCorr      1;
+        nAlphaCorr      2;
-        nAlphaSubCycles 4;
+        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;
+        ".*"           1;
-        "k.*"           1;
-        "epsilon.*"     1;
     }
 }
 

tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/0/alpha1

@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha;
+    object      alpha1;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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;
+rho1            rho1 [ 1 -3 0 0 0 0 0 ] 1996;
-
+rho2            rho2 [ 1 -3 0 0 0 0 0 ] 996;
-rhod            rhod [ 1 -3 0 0 0 0 0 ] 1996;
 
 VdjModel        simple;
 

tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/controlDict

@@ -47,8 +47,7 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.5;
+maxCo           1;
-maxAlphaCo      0.5;
 
 maxDeltaT       0.05;
 

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           ;
+    p_rgh;
-    alpha           ;
+    alpha1;
 }
 
 

tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSolution

@@ -17,15 +17,29 @@ FoamFile
 
 solvers
 {
-    "alpha.*"
+    "alpha1.*"
     {
-        nAlphaCorr      1;
+        nAlphaCorr      2;
-        nAlphaSubCycles 3;
+        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;
 

tutorials/multiphase/driftFluxFoam/ras/tank3D/0/U

@@ -136,8 +136,8 @@ boundaryField
 
     OUTL15
     {
-        type            inletOutlet;
+        type            pressureInletOutletVelocity;
-        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
     }
 }
 

tutorials/multiphase/driftFluxFoam/ras/tank3D/0/alpha1

@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha;
+    object      alpha1;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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;
+rho1            rho1 [ 1 -3 0 0 0 0 0 ] 1042;
-
+rho2            rho2 [ 1 -3 0 0 0 0 0 ] 1000;
-rhod            rhod [ 1 -3 0 0 0 0 0 ] 1042;
 
 VdjModel        simple;
 

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;
 

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;
 }
 
 

tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSolution

@@ -17,10 +17,14 @@ FoamFile
 
 solvers
 {
-    "alpha.*"
+    "alpha1.*"
     {
-        nAlphaCorr      1;
+        nAlphaCorr      2;
-        nAlphaSubCycles 4;
+        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;