From 1eaa024a089175240c1f0685ba48ce580af466ed Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Thu, 21 Dec 2017 17:08:33 -0800
Subject: [PATCH] BUG: Change the calculation of dgdt term in the pEq plus
 reverting compressible pEqComp until further investigation on the
 consequences on dynamic mesh for compressibleInterDyMFoam. alphaSuSp.H has to
 be added in the solver folder in order to make it compatible with the alpha
 Eq.

---
 .../compressibleInterDyMFoam/alphaSuSp.H      | 43 +++++++++++++++++++
 .../compressibleInterDyMFoam/pEqn.H           | 39 +++++++----------
 2 files changed, 58 insertions(+), 24 deletions(-)
 create mode 100644 applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaSuSp.H

diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaSuSp.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaSuSp.H
new file mode 100644
index 0000000000..81309cd091
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaSuSp.H
@@ -0,0 +1,43 @@
+volScalarField::Internal Sp
+(
+    IOobject
+    (
+        "Sp",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("Sp", dgdt.dimensions(), 0)
+);
+
+volScalarField::Internal Su
+(
+    IOobject
+    (
+        "Su",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("Su", dgdt.dimensions(), 0)
+);
+
+forAll(dgdt, celli)
+{
+    if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
+    {
+        Sp[celli] -= dgdt[celli]*alpha1[celli];
+        Su[celli] += dgdt[celli]*alpha1[celli];
+    }
+    else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
+    {
+        Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
+    }
+}
+
+volScalarField::Internal divU
+(
+    mesh.moving()
+  ? fvc::div(phiCN() + mesh.phi())
+  : fvc::div(phiCN())
+);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index b1eace14f2..dd6c345b85 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -56,29 +56,13 @@
     }
     else
     {
-        #include "rhofs.H"
-
-        p_rghEqnComp1 =
-            pos(alpha1)
-           *(
-                (
-                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
-                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
-                )/rho1
-              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
-              + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh))
-            );
+         p_rghEqnComp1 =
+            fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
+          + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
 
         p_rghEqnComp2 =
-            pos(alpha2)
-           *(
-               (
-                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
-                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
-               )/rho2
-             - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
-             + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh))
-            );
+            fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
+          + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
     }
 
     // Cache p_rgh prior to solve for density update
@@ -94,7 +78,11 @@
 
         solve
         (
-            p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
+            (
+                (max(alpha1, scalar(0))/rho1)*p_rghEqnComp1()
+              + (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2()
+            )
+          + p_rghEqnIncomp,
             mesh.solver(p_rgh.select(pimple.finalInnerIter()))
         );
 
@@ -105,8 +93,8 @@
 
             dgdt =
             (
-                alpha1*(p_rghEqnComp2 & p_rgh)
-              - alpha2*(p_rghEqnComp1 & p_rgh)
+                pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
+              - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
             );
 
             phi = phiHbyA + p_rghEqnIncomp.flux();
@@ -131,8 +119,11 @@
     rho = alpha1*rho1 + alpha2*rho2;
 
     // Correct p_rgh for consistency with p and the updated densities
+    p = max(p_rgh + rho*gh, pMin);
     p_rgh = p - rho*gh;
     p_rgh.correctBoundaryConditions();
 
+
+
     K = 0.5*magSqr(U);
 }