twoLiquidMixingFoam: Converted to using MULES and sub-cycling for alpha1

Operator-splitting is used for the diffusion

applications/solvers/multiphase/twoLiquidMixingFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I../interFoam \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \

applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H

@@ -1,22 +1,22 @@
 {
-    fvScalarMatrix alpha1Eqn
+    word alphaScheme("div(phi,alpha1)");
+
+    surfaceScalarField phiAlpha
     (
-        fvm::ddt(alpha1)
-      + fvm::div(phi, alpha1)
-      - fvm::laplacian
+        fvc::flux
         (
-            Dab + alphatab*turbulence->nut(), alpha1,
-            "laplacian(Dab,alpha1)"
+            phi,
+            alpha1,
+            alphaScheme
         )
     );
 
-    alpha1Eqn.solve();
+    MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
 
-    rhoPhi = alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2;
-    rho = alpha1*rho1 + (scalar(1) - alpha1)*rho2;
+    rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
 
-    Info<< "Phase 1 volume fraction = "
-        << alpha1.weightedAverage(mesh.V()).value()
+    Info<< "Phase-1 volume fraction = "
+        << alpha1.weightedAverage(mesh.Vsc()).value()
         << "  Min(alpha1) = " << min(alpha1).value()
         << "  Max(alpha1) = " << max(alpha1).value()
         << endl;

applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C

@@ -32,6 +32,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "MULES.H"
+#include "subCycle.H"
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
@@ -60,19 +62,21 @@ int main(int argc, char *argv[])
     {
         #include "readTimeControls.H"
         #include "CourantNo.H"
+        #include "alphaCourantNo.H"
         #include "setDeltaT.H"
 
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        twoPhaseProperties.correct();
+
+        #include "alphaEqnSubCycle.H"
+        #include "alphaDiffusionEqn.H"
+
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
-            twoPhaseProperties.correct();
-
-            #include "alphaEqn.H"
-
             #include "UEqn.H"
 
             // --- Pressure corrector loop

tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/controlDict

@@ -48,6 +48,7 @@ runTimeModifiable yes;
 adjustTimeStep  on;
 
 maxCo           0.5;
+maxAlphaCo      0.5;
 maxDeltaT       1;
 
 // ************************************************************************* //

tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSolution

@@ -21,7 +21,7 @@ solvers
     {
         solver          smoothSolver;
         smoother        GaussSeidel;
-        tolerance       1e-7;
+        tolerance       1e-9;
         relTol          0;
         nSweeps         1;
     }
@@ -66,6 +66,7 @@ PIMPLE
     nOuterCorrectors    1;
     nCorrectors         2;
     nNonOrthogonalCorrectors 0;
+    nAlphaSubCycles     2;
     pRefValue           0;
     pRefPoint           (0.1 0.1 1);
 }