multiphase: Updated all multiphase solvers to use consistent naming o…

…f phases and phase variables

Improved handling of pressure BCs
Added run-time selectable equation of state to compressibleInterFoam
Added MRF support to twoPhaseEulerFoam and compressibleTwoPhaseEulerFoam

applications/solvers/multiphase/bubbleFoam/DDtU.H

@@ -1,11 +1,11 @@
 {
-    DDtUa =
-        fvc::ddt(Ua)
-      + fvc::div(phia, Ua)
-      - fvc::div(phia)*Ua;
+    DDtU1 =
+        fvc::ddt(U1)
+      + fvc::div(phi1, U1)
+      - fvc::div(phi1)*U1;
 
-    DDtUb =
-        fvc::ddt(Ub)
-      + fvc::div(phib, Ub)
-      - fvc::div(phib)*Ub;
+    DDtU2 =
+        fvc::ddt(U2)
+      + fvc::div(phi2, U2)
+      - fvc::div(phi2)*U2;
 }

applications/solvers/multiphase/bubbleFoam/UEqns.H

@@ -1,72 +1,74 @@
-fvVectorMatrix UaEqn(Ua, Ua.dimensions()*dimVol/dimTime);
-fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
+fvVectorMatrix U1Eqn(U1, U1.dimensions()*dimVol/dimTime);
+fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
 
 {
-    volTensorField Rca(-nuEffa*(T(fvc::grad(Ua))));
-    Rca = Rca + (2.0/3.0)*sqr(Ct)*I*k - (2.0/3.0)*I*tr(Rca);
+    volTensorField Rc1(-nuEff1*(T(fvc::grad(U1))));
+    Rc1 = Rc1 + (2.0/3.0)*sqr(Ct)*I*k - (2.0/3.0)*I*tr(Rc1);
 
-    surfaceScalarField phiRa
+    surfaceScalarField phiR1
     (
-      - fvc::interpolate(nuEffa)
-       *mesh.magSf()*fvc::snGrad(alpha)/fvc::interpolate(alpha + scalar(0.001))
+      - fvc::interpolate(nuEff1)
+       *mesh.magSf()*fvc::snGrad(alpha1)
+       /fvc::interpolate(alpha1 + scalar(0.001))
     );
 
-    UaEqn =
+    U1Eqn =
     (
-        (scalar(1) + Cvm*rhob*beta/rhoa)*
+        (scalar(1) + Cvm*rho2*alpha2/rho1)*
         (
-            fvm::ddt(Ua)
-          + fvm::div(phia, Ua, "div(phia,Ua)")
-          - fvm::Sp(fvc::div(phia), Ua)
+            fvm::ddt(U1)
+          + fvm::div(phi1, U1, "div(phi1,U1)")
+          - fvm::Sp(fvc::div(phi1), U1)
         )
 
-      - fvm::laplacian(nuEffa, Ua)
-      + fvc::div(Rca)
+      - fvm::laplacian(nuEff1, U1)
+      + fvc::div(Rc1)
 
-      + fvm::div(phiRa, Ua, "div(phia,Ua)")
-      - fvm::Sp(fvc::div(phiRa), Ua)
-      + (fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca)
+      + fvm::div(phiR1, U1, "div(phi1,U1)")
+      - fvm::Sp(fvc::div(phiR1), U1)
+      + (fvc::grad(alpha1)/(fvc::average(alpha1) + scalar(0.001)) & Rc1)
      ==
     //  g                        // Buoyancy term transfered to p-equation
-      - fvm::Sp(beta/rhoa*dragCoef, Ua)
-    //+ beta/rhoa*dragCoef*Ub    // Explicit drag transfered to p-equation
-      - beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb)
+      - fvm::Sp(alpha2/rho1*dragCoef, U1)
+    //+ alpha2/rho1*dragCoef*U2    // Explicit drag transfered to p-equation
+      - alpha2/rho1*(liftCoeff - Cvm*rho2*DDtU2)
     );
 
-    UaEqn.relax();
+    U1Eqn.relax();
 
 
-    volTensorField Rcb(-nuEffb*T(fvc::grad(Ub)));
-    Rcb = Rcb + (2.0/3.0)*I*k - (2.0/3.0)*I*tr(Rcb);
+    volTensorField Rc2(-nuEff2*T(fvc::grad(U2)));
+    Rc2 = Rc2 + (2.0/3.0)*I*k - (2.0/3.0)*I*tr(Rc2);
 
-    surfaceScalarField phiRb
+    surfaceScalarField phiR2
     (
-       - fvc::interpolate(nuEffb)
-        *mesh.magSf()*fvc::snGrad(beta)/fvc::interpolate(beta + scalar(0.001))
+       - fvc::interpolate(nuEff2)
+        *mesh.magSf()*fvc::snGrad(alpha2)
+        /fvc::interpolate(alpha2 + scalar(0.001))
     );
 
-    UbEqn =
+    U2Eqn =
     (
-        (scalar(1) + Cvm*rhob*alpha/rhob)*
+        (scalar(1) + Cvm*rho2*alpha1/rho2)*
         (
-            fvm::ddt(Ub)
-          + fvm::div(phib, Ub, "div(phib,Ub)")
-          - fvm::Sp(fvc::div(phib), Ub)
+            fvm::ddt(U2)
+          + fvm::div(phi2, U2, "div(phi2,U2)")
+          - fvm::Sp(fvc::div(phi2), U2)
         )
 
-      - fvm::laplacian(nuEffb, Ub)
-      + fvc::div(Rcb)
+      - fvm::laplacian(nuEff2, U2)
+      + fvc::div(Rc2)
 
-      + fvm::div(phiRb, Ub, "div(phib,Ub)")
-      - fvm::Sp(fvc::div(phiRb), Ub)
+      + fvm::div(phiR2, U2, "div(phi2,U2)")
+      - fvm::Sp(fvc::div(phiR2), U2)
 
-      + (fvc::grad(beta)/(fvc::average(beta) + scalar(0.001)) & Rcb)
+      + (fvc::grad(alpha2)/(fvc::average(alpha2) + scalar(0.001)) & Rc2)
      ==
     //  g                        // Buoyancy term transfered to p-equation
-      - fvm::Sp(alpha/rhob*dragCoef, Ub)
-    //+ alpha/rhob*dragCoef*Ua   // Explicit drag transfered to p-equation
-      + alpha/rhob*(liftCoeff + Cvm*rhob*DDtUa)
+      - fvm::Sp(alpha1/rho2*dragCoef, U2)
+    //+ alpha1/rho2*dragCoef*U1   // Explicit drag transfered to p-equation
+      + alpha1/rho2*(liftCoeff + Cvm*rho2*DDtU1)
     );
 
-    UbEqn.relax();
+    U2Eqn.relax();
 }

applications/solvers/multiphase/bubbleFoam/alphaEqn.H

@@ -1,7 +1,7 @@
 {
-    word scheme("div(phi,alpha)");
+    word scheme("div(phi,alpha1)");
 
-    surfaceScalarField phir(phia - phib);
+    surfaceScalarField phir(phi1 - phi2);
 
     Info<< "Max Ur Courant Number = "
         << (
@@ -15,42 +15,47 @@
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
     {
-        fvScalarMatrix alphaEqn
+        fvScalarMatrix alpha1Eqn
         (
-             fvm::ddt(alpha)
-           + fvm::div(phi, alpha, scheme)
-           + fvm::div(-fvc::flux(-phir, beta, scheme), alpha, scheme)
+             fvm::ddt(alpha1)
+           + fvm::div(phi, alpha1, scheme)
+           + fvm::div(-fvc::flux(-phir, alpha2, scheme), alpha1, scheme)
         );
-        alphaEqn.relax();
-        alphaEqn.solve();
+        alpha1Eqn.relax();
+        alpha1Eqn.solve();
 
         /*
-        fvScalarMatrix betaEqn
+        fvScalarMatrix alpha2Eqn
         (
-              fvm::ddt(beta)
-            + fvm::div(phi, beta, scheme)
-            + fvm::div(-fvc::flux(phir, scalar(1) - beta, scheme), beta, scheme)
+            fvm::ddt(alpha2)
+          + fvm::div(phi, alpha2, scheme)
+          + fvm::div
+            (
+                -fvc::flux(phir, scalar(1) - alpha2, scheme),
+                alpha2,
+                scheme
+            )
         );
-        betaEqn.relax();
-        betaEqn.solve();
+        alpha2Eqn.relax();
+        alpha2Eqn.solve();
 
-        alpha =
+        alpha1 =
             0.5
            *(
                  scalar(1)
-               + sqr(scalar(1) - beta)
-               - sqr(scalar(1) - alpha)
+               + sqr(scalar(1) - alpha2)
+               - sqr(scalar(1) - alpha1)
             );
         */
 
-        beta = scalar(1) - alpha;
+        alpha2 = scalar(1) - alpha1;
     }
 
     Info<< "Dispersed phase volume fraction = "
-        << alpha.weightedAverage(mesh.V()).value()
-        << "  Min(alpha) = " << min(alpha).value()
-        << "  Max(alpha) = " << max(alpha).value()
+        << alpha1.weightedAverage(mesh.V()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
         << endl;
 }
 
-rho = alpha*rhoa + beta*rhob;
+rho = alpha1*rho1 + alpha2*rho2;

applications/solvers/multiphase/bubbleFoam/bubbleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -90,6 +90,7 @@ int main(int argc, char *argv[])
             if (pimple.turbCorr())
             {
                 #include "kEpsilon.H"
+                nuEff1 = sqr(Ct)*nut2 + nu1;
             }
         }