twoPhaseMixture: Add alpha2

applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H

@@ -4,8 +4,8 @@
         p =
         (
             rho
-          - (1.0 - gamma)*rhol0
+          - gamma2*rhol0
-          - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+          - ((gamma*psiv + gamma2*psil) - psi)*pSat
         )/psi;
     }
 
@@ -57,8 +57,8 @@
     p =
     (
         rho
-      - (1.0 - gamma)*rhol0
+      - gamma2*rhol0
-      - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+      - ((gamma*psiv + gamma2*psil) - psi)*pSat
     )/psi;
 
     p.correctBoundaryConditions();

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);
 

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
+      + gamma2*rhol0
-      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
+      + ((gamma*psiv + gamma2*psil) - psi)*pSat,
         rhoMin
     );
 

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;

applications/solvers/multiphase/cavitatingFoam/pEqn.H

@@ -4,8 +4,8 @@
         p =
         (
             rho
-          - (1.0 - gamma)*rhol0
+          - gamma2*rhol0
-          - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+          - ((gamma*psiv + gamma2*psil) - psi)*pSat
         )/psi;
     }
 
@@ -49,8 +49,8 @@
     rho == max
     (
         psi*p
-      + (1.0 - gamma)*rhol0
+      + gamma2*rhol0
-      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
+      + ((gamma*psiv + gamma2*psil) - psi)*pSat,
         rhoMin
     );
 
@@ -59,8 +59,8 @@
     p =
     (
         rho
-      - (1.0 - gamma)*rhol0
+      - gamma2*rhol0
-      - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+      - ((gamma*psiv + gamma2*psil) - psi)*pSat
     )/psi;
 
     p.correctBoundaryConditions();

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;
     }
 

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;

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;
     }
 

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;

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();

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;

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;

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();
 
 

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

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
         {