From 0b240922e9ab95a44c85e9e3f76f6bb9e5f750eb Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 19 Mar 2012 16:14:58 +0000
Subject: [PATCH] twoPhaseEulerFoam: updated alphaEqn to use MULES

---
 .../multiphase/twoPhaseEulerFoam/alphaEqn.H   | 66 ++++++++++++++-----
 .../readTwoPhaseEulerFoamControls.H           |  2 -
 .../twoPhaseEulerFoam/twoPhaseEulerFoam.C     |  2 +
 3 files changed, 51 insertions(+), 19 deletions(-)

diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
index 52f09e4793..34fdc24afc 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
@@ -1,6 +1,9 @@
 {
-    word scheme("div(phi,alpha1)");
-    word schemer("div(phir,alpha1)");
+    label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
+    label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
+
+    word alphaScheme("div(phi,alpha1)");
+    word alpharScheme("div(phir,alpha1)");
 
     surfaceScalarField phic("phic", phi);
     surfaceScalarField phir("phir", phi1 - phi2);
@@ -15,14 +18,39 @@
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
     {
-        fvScalarMatrix alpha1Eqn
+        for
         (
-             fvm::ddt(alpha1)
-           + fvm::div(phic, alpha1, scheme)
-           + fvm::div(-fvc::flux(-phir, alpha2, schemer), alpha1, schemer)
-        );
+            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
+        )
+        {
+            surfaceScalarField phiAlpha
+            (
+                fvc::flux
+                (
+                    phic,
+                    alpha1,
+                    alphaScheme
+                )
+              + fvc::flux
+                (
+                    -fvc::flux(-phir, alpha2, alpharScheme),
+                    alpha1,
+                    alpharScheme
+                )
+            );
 
-        if (g0.value() > 0.0)
+            MULES::explicitSolve
+            (
+                alpha1,
+                phi,
+                phiAlpha,
+                (g0.value() > 0 ? alphaMax : 1),
+                0
+            );
+        }
+
+        if (g0.value() > 0)
         {
             ppMagf = rAU1f*fvc::interpolate
             (
@@ -30,19 +58,23 @@
                *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
             );
 
-            alpha1Eqn -= fvm::laplacian
+            fvScalarMatrix alpha1Eqn
             (
-                (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
-                alpha1,
-                "laplacian(alpha1PpMag,alpha1)"
+                fvm::ddt(alpha1) - fvc::ddt(alpha1)
+              - fvm::laplacian
+                (
+                    (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
+                    alpha1,
+                    "laplacian(alpha1PpMag,alpha1)"
+                )
             );
+
+            alpha1Eqn.relax();
+            alpha1Eqn.solve();
+
+            #include "packingLimiter.H"
         }
 
-        alpha1Eqn.relax();
-        alpha1Eqn.solve();
-
-        #include "packingLimiter.H"
-
         alpha2 = scalar(1) - alpha1;
 
         Info<< "Dispersed phase volume fraction = "
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
index a345c4e53b..4e6847debe 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
@@ -1,5 +1,3 @@
     #include "readTimeControls.H"
 
-    int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
-
     Switch correctAlpha(pimple.dict().lookup("correctAlpha"));
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index 357c7e413a..d9da532aef 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -31,6 +31,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "MULES.H"
+#include "subCycle.H"
 #include "nearWallDist.H"
 #include "wallFvPatch.H"
 #include "Switch.H"