From 59f9f9ebd15e269ded39a1b80af7512023c4241c Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Wed, 6 Feb 2013 10:52:49 +0000
Subject: [PATCH] twoPhaseMixture: Add alpha2

---
 .../cavitatingFoam/cavitatingDyMFoam/pEqn.H      |  8 ++++----
 .../multiphase/cavitatingFoam/continuityErrs.H   |  2 +-
 .../multiphase/cavitatingFoam/createFields.H     |  6 ++++--
 .../solvers/multiphase/cavitatingFoam/gammaPsi.H |  1 +
 .../solvers/multiphase/cavitatingFoam/pEqn.H     | 12 ++++++------
 .../multiphase/interFoam/LTSInterFoam/alphaEqn.H |  3 ++-
 .../interFoam/LTSInterFoam/alphaEqnSubCycle.H    |  2 +-
 .../solvers/multiphase/interFoam/alphaEqn.H      |  3 ++-
 .../multiphase/interFoam/alphaEqnSubCycle.H      |  2 +-
 .../solvers/multiphase/interFoam/createFields.H  |  3 ++-
 .../twoLiquidMixingFoam/alphaDiffusionEqn.H      |  3 ++-
 .../twoLiquidMixingFoam/alphaEqnSubCycle.H       |  2 +-
 .../twoLiquidMixingFoam/createFields.H           |  3 ++-
 .../twoPhaseMixture.C                            | 14 +++++++++++++-
 .../twoPhaseMixture.H                            | 16 +++++++++++++++-
 15 files changed, 57 insertions(+), 23 deletions(-)

diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index a4e9e0b40f6..3c4cbb0d32e 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@@ -4,8 +4,8 @@
         p =
         (
             rho
-          - (1.0 - gamma)*rhol0
-          - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+          - gamma2*rhol0
+          - ((gamma*psiv + gamma2*psil) - psi)*pSat
         )/psi;
     }
 
@@ -57,8 +57,8 @@
     p =
     (
         rho
-      - (1.0 - gamma)*rhol0
-      - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+      - gamma2*rhol0
+      - ((gamma*psiv + gamma2*psil) - psi)*pSat
     )/psi;
 
     p.correctBoundaryConditions();
diff --git a/applications/solvers/multiphase/cavitatingFoam/continuityErrs.H b/applications/solvers/multiphase/cavitatingFoam/continuityErrs.H
index c5e93c7efb3..ce618ee40f2 100644
--- a/applications/solvers/multiphase/cavitatingFoam/continuityErrs.H
+++ b/applications/solvers/multiphase/cavitatingFoam/continuityErrs.H
@@ -1,5 +1,5 @@
 {
-    volScalarField thermoRho = psi*p + (1.0 - gamma)*rhol0;
+    volScalarField thermoRho = psi*p + gamma2*rhol0;
 
     dimensionedScalar totalMass = fvc::domainIntegrate(rho);
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/createFields.H b/applications/solvers/multiphase/cavitatingFoam/createFields.H
index dbacf1dbd3f..c7241406504 100644
--- a/applications/solvers/multiphase/cavitatingFoam/createFields.H
+++ b/applications/solvers/multiphase/cavitatingFoam/createFields.H
@@ -49,6 +49,8 @@
     volScalarField& gamma(twoPhaseProperties.alpha1());
     gamma.oldTime();
 
+    volScalarField& gamma2(twoPhaseProperties.alpha2());
+
     Info<< "Creating compressibilityModel\n" << endl;
     autoPtr<barotropicCompressibilityModel> psiModel =
         barotropicCompressibilityModel::New
@@ -62,8 +64,8 @@
     rho == max
     (
         psi*p
-      + (1.0 - gamma)*rhol0
-      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
+      + gamma2*rhol0
+      + ((gamma*psiv + gamma2*psil) - psi)*pSat,
         rhoMin
     );
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/gammaPsi.H b/applications/solvers/multiphase/cavitatingFoam/gammaPsi.H
index b259ddd3222..4edda7336de 100644
--- a/applications/solvers/multiphase/cavitatingFoam/gammaPsi.H
+++ b/applications/solvers/multiphase/cavitatingFoam/gammaPsi.H
@@ -1,5 +1,6 @@
 {
     gamma = max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));
+    gamma2 = 1.0 - gamma;
 
     Info<< "max-min gamma: " << max(gamma).value()
         << " " << min(gamma).value() << endl;
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index ea35f79fa14..b2c7d953b53 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -4,8 +4,8 @@
         p =
         (
             rho
-          - (1.0 - gamma)*rhol0
-          - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+          - gamma2*rhol0
+          - ((gamma*psiv + gamma2*psil) - psi)*pSat
         )/psi;
     }
 
@@ -49,8 +49,8 @@
     rho == max
     (
         psi*p
-      + (1.0 - gamma)*rhol0
-      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
+      + gamma2*rhol0
+      + ((gamma*psiv + gamma2*psil) - psi)*pSat,
         rhoMin
     );
 
@@ -59,8 +59,8 @@
     p =
     (
         rho
-      - (1.0 - gamma)*rhol0
-      - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+      - gamma2*rhol0
+      - ((gamma*psiv + gamma2*psil) - psi)*pSat
     )/psi;
 
     p.correctBoundaryConditions();
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
index 5f5ac824b8a..faae1976703 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
@@ -18,7 +18,7 @@
             )
           + fvc::flux
             (
-                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+                -fvc::flux(-phir, alpha2, alpharScheme),
                 alpha1,
                 alpharScheme
             )
@@ -27,6 +27,7 @@
         MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
         //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
 
+        alpha2 = 1.0 - alpha1;
         rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
     }
 
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
index 6c82f940499..57c78027a47 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
@@ -23,4 +23,4 @@ else
     #include "alphaEqn.H"
 }
 
-rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;
+rho == alpha1*rho1 + alpha2*rho2;
diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H
index b96dcf898da..a2720e20eee 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqn.H
@@ -18,7 +18,7 @@
             )
           + fvc::flux
             (
-                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+                -fvc::flux(-phir, alpha2, alpharScheme),
                 alpha1,
                 alpharScheme
             )
@@ -26,6 +26,7 @@
 
         MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
 
+        alpha2 = 1.0 - alpha1;
         rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
     }
 
diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
index 6c82f940499..57c78027a47 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
@@ -23,4 +23,4 @@ else
     #include "alphaEqn.H"
 }
 
-rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;
+rho == alpha1*rho1 + alpha2*rho2;
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index df1edf04e1e..c607c685489 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -33,6 +33,7 @@
     twoPhaseMixture twoPhaseProperties(U, phi);
 
     volScalarField& alpha1(twoPhaseProperties.alpha1());
+    volScalarField& alpha2(twoPhaseProperties.alpha2());
 
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
@@ -48,7 +49,7 @@
             mesh,
             IOobject::READ_IF_PRESENT
         ),
-        alpha1*rho1 + (scalar(1) - alpha1)*rho2,
+        alpha1*rho1 + alpha2*rho2,
         alpha1.boundaryField().types()
     );
     rho.oldTime();
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H
index 139646e43c6..6bc5953a4a6 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H
@@ -12,7 +12,8 @@
 
     alpha1Eqn.solve();
 
+    alpha2 = 1.0 - alpha1;
     rhoPhi += alpha1Eqn.flux()*(rho1 - rho2);
 }
 
-rho = alpha1*rho1 + (scalar(1) - alpha1)*rho2;
+rho = alpha1*rho1 + alpha2*rho2;
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
index 6c82f940499..57c78027a47 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
@@ -23,4 +23,4 @@ else
     #include "alphaEqn.H"
 }
 
-rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;
+rho == alpha1*rho1 + alpha2*rho2;
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
index 0d01b9a9e5e..e98f06e74e5 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
@@ -32,6 +32,7 @@
     twoPhaseMixture twoPhaseProperties(U, phi);
 
     volScalarField& alpha1(twoPhaseProperties.alpha1());
+    volScalarField& alpha2(twoPhaseProperties.alpha2());
 
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
@@ -42,7 +43,7 @@
     dimensionedScalar alphatab(twoPhaseProperties.lookup("alphatab"));
 
     // Need to store rho for ddt(rho, U)
-    volScalarField rho("rho", alpha1*rho1 + (scalar(1) - alpha1)*rho2);
+    volScalarField rho("rho", alpha1*rho1 + alpha2*rho2);
     rho.oldTime();
 
 
diff --git a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.C b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.C
index 4ec3ce72f85..4b698138a77 100644
--- a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.C
+++ b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.C
@@ -54,7 +54,8 @@ Foam::twoPhaseMixture::twoPhaseMixture
 (
     const volVectorField& U,
     const surfaceScalarField& phi,
-    const word& alpha1Name
+    const word& alpha1Name,
+    const word& alpha2Name
 )
 :
     transportModel(U, phi),
@@ -102,6 +103,17 @@ Foam::twoPhaseMixture::twoPhaseMixture
         U_.mesh()
     ),
 
+    alpha2_
+    (
+        IOobject
+        (
+            found("phases") ? word("alpha" + phase2Name_) : alpha2Name,
+            U_.time().timeName(),
+            U_.db()
+        ),
+        1.0 - alpha1_
+    ),
+
     nu_
     (
         IOobject
diff --git a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H
index 025d95f7011..265f862bce2 100644
--- a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H
+++ b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H
@@ -70,6 +70,7 @@ protected:
         const surfaceScalarField& phi_;
 
         volScalarField alpha1_;
+        volScalarField alpha2_;
 
         volScalarField nu_;
 
@@ -89,7 +90,8 @@ public:
         (
             const volVectorField& U,
             const surfaceScalarField& phi,
-            const word& alpha1Name = "alpha1"
+            const word& alpha1Name = "alpha1",
+            const word& alpha2Name = "alpha2"
         );
 
 
@@ -122,6 +124,18 @@ public:
             return alpha1_;
         }
 
+        //- Return the phase-fraction of phase 2
+        const volScalarField& alpha2() const
+        {
+            return alpha2_;
+        }
+
+        //- Return the phase-fraction of phase 2
+        volScalarField& alpha2()
+        {
+            return alpha2_;
+        }
+
         //- Return const-access to phase1 viscosityModel
         const viscosityModel& nuModel1() const
         {
-- 
GitLab