From 314c9475ff7cc5f92e923a14a26437e64d438e6a Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Wed, 24 Oct 2012 17:25:23 +0100
Subject: [PATCH] twoLiquidMixingFoam: Converted to using MULES and sub-cycling
 for alpha1

Operator-splitting is used for the diffusion
---
 .../twoLiquidMixingFoam/alphaDiffusionEqn.H   | 18 +++++++++++++
 .../twoLiquidMixingFoam/alphaEqnSubCycle.H    | 26 +++++++++++++++++++
 2 files changed, 44 insertions(+)
 create mode 100644 applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H
 create mode 100644 applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H

diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H
new file mode 100644
index 00000000000..139646e43c6
--- /dev/null
+++ b/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;
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
new file mode 100644
index 00000000000..6c82f940499
--- /dev/null
+++ b/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;
-- 
GitLab