compressibleTwoPhaseEulerFoam: Updated virtual-mass formulation

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H

@@ -28,12 +28,17 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
 
         U1Eqn =
         (
-            (scalar(1) + Cvm*rho2*alpha2/rho1)*
+            fvm::ddt(alpha1, U1)
+          + fvm::div(alphaPhi1, U1)
+          - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), U1)
+
+          + Cvm*rho2*alpha1*alpha2/rho1*
             (
-                fvm::ddt(alpha1, U1)
-              + fvm::div(alphaPhi1, U1)
-              - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), U1)
+                fvm::ddt(U1)
+              + fvm::div(phi1, U1)
+              - fvm::Sp(fvc::div(phi1), U1)
             )
+
           - fvm::laplacian(alpha1*nuEff1, U1)
           + fvc::div(alpha1*Rc1)
          ==
@@ -54,12 +59,17 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
 
         U2Eqn =
         (
-            (scalar(1) + Cvm*rho2*alpha1/rho2)*
+            fvm::ddt(alpha2, U2)
+          + fvm::div(alphaPhi2, U2)
+          - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), U2)
+
+          + Cvm*rho2*alpha1*alpha2/rho2*
             (
-                fvm::ddt(alpha2, U2)
-              + fvm::div(alphaPhi2, U2)
-              - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), U2)
+                fvm::ddt(U2)
+              + fvm::div(phi2, U2)
+              - fvm::Sp(fvc::div(phi2), U2)
             )
+
           - fvm::laplacian(alpha2*nuEff2, U2)
           + fvc::div(alpha2*Rc2)
          ==

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H

@@ -37,6 +37,9 @@
 
     surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
 
+    HbyA1 += (1.0/rho1)*rAU1*dragCoeff*U2;
+    HbyA2 += (1.0/rho2)*rAU2*dragCoeff*U1;
+
     surfaceScalarField Dp
     (
         "Dp",
@@ -110,10 +113,7 @@
             p.relax();
             mSfGradp = pEqnIncomp.flux()/Dp;
 
-            volVectorField U10 = U1;
-
             U1 = HbyA1
-               + (1.0/rho1)*rAU1*dragCoeff*U2
                + fvc::reconstruct
                  (
                      rAlphaAU1f*(g & mesh.Sf())
@@ -122,7 +122,6 @@
             U1.correctBoundaryConditions();
 
             U2 = HbyA2
-               + (1.0/rho2)*rAU2*dragCoeff*U10
                + fvc::reconstruct
                  (
                      rAlphaAU2f*(g & mesh.Sf())