compressibleTwoPhaseEulerFoam/pEqn.H: Corrected handling of the kinetic theory...

compressibleTwoPhaseEulerFoam/pEqn.H: Corrected handling of the kinetic theory drag in the pressure equation

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H

@@ -22,27 +22,29 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
-    surfaceScalarField phiHbyA1
-    (
-        "phiHbyA1",
-        (fvc::interpolate(HbyA1) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
-      + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
-      + rAlphaAU1f*(g & mesh.Sf())
-    );
+    surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
 
     if (g0.value() > 0.0)
     {
-        phiHbyA1 -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
+        ppDrag -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
     }
 
     if (kineticTheory.on())
     {
-        phiHbyA1 -=
-            fvc::interpolate((1.0/rho1)*rAU1)
+        ppDrag -=
+            fvc::interpolate(1.0/rho1)*rAlphaAU1f
            *fvc::snGrad(kineticTheory.pa())*mesh.magSf();
     }
 
+    surfaceScalarField phiHbyA1
+    (
+        "phiHbyA1",
+        (fvc::interpolate(HbyA1) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
+      + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
+      + ppDrag
+      + rAlphaAU1f*(g & mesh.Sf())
+    );
     mrfZones.relativeFlux(phiHbyA1);
 
     surfaceScalarField phiHbyA2
@@ -135,7 +137,8 @@
             U1 = HbyA1
                + fvc::reconstruct
                  (
-                     rAlphaAU1f
+                     ppDrag
+                   + rAlphaAU1f
                     *(
                          (g & mesh.Sf())
                        + mSfGradp/fvc::interpolate(rho1)