Multiphase solvers: Update p_rgh following density changes

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C

@@ -146,6 +146,10 @@ int main(int argc, char *argv[])
 
         rho = alpha1*rho1 + alpha2*rho2;
 
+        // Correct p_rgh for consistency with p and the updated densities
+        p_rgh = p - rho*gh;
+        p_rgh.correctBoundaryConditions();
+
         runTime.write();
 
         Info<< "ExecutionTime = "

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H

@@ -116,6 +116,10 @@
 
     rho = alpha1*rho1 + alpha2*rho2;
 
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
+    p_rgh.correctBoundaryConditions();
+
     K = 0.5*magSqr(U);
 
     Info<< "max(U) " << max(mag(U)).value() << endl;

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -101,14 +101,16 @@
         }
     }
 
-    // p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
-
     // Update densities from change in p_rgh
     rho1 += psi1*(p_rgh - p_rgh_0);
     rho2 += psi2*(p_rgh - p_rgh_0);
 
     rho = alpha1*rho1 + alpha2*rho2;
 
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
+    p_rgh.correctBoundaryConditions();
+
     K = 0.5*magSqr(U);
 
     Info<< "max(U) " << max(mag(U)).value() << endl;

applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H

@@ -95,9 +95,6 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            // p = max(p_rgh + mixture.rho()*gh, pMin);
-            // p_rgh = p - mixture.rho()*gh;
-
             phasei = 0;
             forAllIter
             (
@@ -125,6 +122,10 @@
     mixture.correctRho(p_rgh - p_rgh_0);
     rho = mixture.rho();
 
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
+    p_rgh.correctBoundaryConditions();
+
     K = 0.5*magSqr(U);
 
     Info<< "max(U) " << max(mag(U)).value() << endl;

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H

@@ -92,6 +92,11 @@ while (pimple.correct())
     // Update continuity errors due to temperature changes
     fluid.correct();
 
+    volScalarField rho("rho", fluid.rho());
+
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
+
     PtrList<volVectorField> HbyAs(phases.size());
 
     forAll(phases, phasei)
@@ -121,9 +126,6 @@ while (pimple.correct())
             );
     }
 
-    // Mean density for buoyancy force and p_rgh -> p
-    volScalarField rho("rho", fluid.rho());
-
     surfaceScalarField ghSnGradRho
     (
         "ghSnGradRho",

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H

@@ -85,6 +85,7 @@ while (pimple.correct())
 {
     // Update continuity errors due to temperature changes
     fluid.correct();
+
     volScalarField rho("rho", fluid.rho());
 
     // Correct p_rgh for consistency with p and the updated densities

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H

@@ -95,6 +95,7 @@ while (pimple.correct())
 {
     // Update continuity errors due to temperature changes
     fluid.correct();
+
     volScalarField rho("rho", fluid.rho());
 
     // Correct p_rgh for consistency with p and the updated densities

applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H

@@ -83,6 +83,11 @@ while (pimple.correct())
     // Update continuity errors due to temperature changes
     #include "correctContErrs.H"
 
+    volScalarField rho("rho", fluid.rho());
+
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
+
     // Correct fixed-flux BCs to be consistent with the velocity BCs
     MRF.correctBoundaryFlux(U1, phi1);
     MRF.correctBoundaryFlux(U2, phi2);
@@ -113,9 +118,6 @@ while (pimple.correct())
            *rho2*U2.oldTime()/runTime.deltaT()
         );
 
-    // Mean density for buoyancy force and p_rgh -> p
-    volScalarField rho("rho", fluid.rho());
-
     surfaceScalarField ghSnGradRho
     (
         "ghSnGradRho",

applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H

@@ -95,6 +95,11 @@ while (pimple.correct())
     // Update continuity errors due to temperature changes
     #include "correctContErrs.H"
 
+    volScalarField rho("rho", fluid.rho());
+
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
+
     surfaceScalarField rhof1(fvc::interpolate(rho1));
     surfaceScalarField rhof2(fvc::interpolate(rho2));
 
@@ -114,7 +119,6 @@ while (pimple.correct())
         max(alphaf2, phase2.residualAlpha())*rAUf2
     );
 
-    volScalarField rho("rho", fluid.rho());
     surfaceScalarField ghSnGradRho
     (
         "ghSnGradRho",