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

Operator-splitting is used for the diffusion

applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H

@@ -0,0 +1,18 @@
+{
+    fvScalarMatrix alpha1Eqn
+    (
+        fvm::ddt(alpha1)
+      - fvc::ddt(alpha1)
+      - fvm::laplacian
+        (
+            volScalarField("Dab", Dab + alphatab*turbulence->nut()),
+            alpha1
+        )
+    );
+
+    alpha1Eqn.solve();
+
+    rhoPhi += alpha1Eqn.flux()*(rho1 - rho2);
+}
+
+rho = alpha1*rho1 + (scalar(1) - alpha1)*rho2;

applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H

@@ -0,0 +1,26 @@
+label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
+label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
+
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalarField rhoPhiSum(0.0*rhoPhi);
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+        #include "alphaEqn.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+    #include "alphaEqn.H"
+}
+
+rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;