compressibleTwoPhaseEulerFoam: Improved handling of compressibility effect on continuity and boundedness

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H

@@ -21,6 +21,10 @@
     (
         fvm::ddt(alpha1, he1)
       + fvm::div(alphaPhi1, he1)
+
+        // Compressibity correction
+      - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), he1)
+
       - fvm::laplacian(k1, he1)
      ==
         heatTransferCoeff*(thermo2.T() - thermo1.T())/rho1
@@ -33,6 +37,10 @@
     (
         fvm::ddt(alpha2, he2)
       + fvm::div(alphaPhi2, he2)
+
+        // Compressibity correction
+      - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), he2)
+
       - fvm::laplacian(k2, he2)
      ==
         heatTransferCoeff*(thermo1.T() - thermo2.T())/rho2

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H

@@ -31,6 +31,9 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
             fvm::ddt(alpha1, U1)
           + fvm::div(alphaPhi1, U1)
 
+            // Compressibity correction
+          - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), U1)
+
           + Cvm*rho2*alpha1*alpha2/rho1*
             (
                 fvm::ddt(U1)
@@ -61,6 +64,9 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
             fvm::ddt(alpha2, U2)
           + fvm::div(alphaPhi2, U2)
 
+            // Compressibity correction
+          - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), U2)
+
           + Cvm*rho2*alpha1*alpha2/rho2*
             (
                 fvm::ddt(U2)

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H

@@ -16,6 +16,10 @@ if (turbulence)
     (
         fvm::ddt(alpha2, epsilon)
       + fvm::div(alphaPhi2, epsilon)
+
+        // Compressibity correction
+      - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), epsilon)
+
       - fvm::laplacian
         (
             alpha1Eps*nuEff2, epsilon,
@@ -41,6 +45,10 @@ if (turbulence)
     (
         fvm::ddt(alpha2, k)
       + fvm::div(alphaPhi2, k)
+
+        // Compressibity correction
+      - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), k)
+
       - fvm::laplacian
         (
             alpha1k*nuEff2, k,