isothermalFluid::correctBuoyantPressure: Relax the net force momentum equation source

rather than p_rgh which introduces an imbalance between the pressure and
buoyancy forces.  The relaxation factor for p_rgh specified in fvSolution is
used to relax the net force as the intent is to relax the pressure and this
provides convenient usage and backwards-compatibility.

Optional relaxation for the thermodynamic pressure p is also available for case
this provides convergence benefit for steady cases by relaxing the pressure work
term is the energy equation.

Note these changes only relate to the operation of the isothermalFluid solver
module for buoyant cases.

applications/solvers/modules/fluid/isothermalFluid/correctBuoyantPressure.C

@@ -102,8 +102,9 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure()
            *fvc::relative(phiHbyA, rho, U)
         );
 
-        const fvScalarMatrix divPhidp(fvm::div(phid, p));
-        phiHbyA -= divPhidp.flux();
+        // Subtract the compressible part
+        // The resulting flux will be zero for a perfect gas
+        phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
         if (pimple.consistent())
         {
@@ -120,7 +121,7 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure()
         fvScalarMatrix p_rghDDtEqn
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-          + fvc::div(phiHbyA) + divPhidp
+          + fvc::div(phiHbyA) + fvm::div(phid, p)
          ==
             fvModels().source(psi, p_rgh, rho.name())
         );
@@ -183,6 +184,21 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure()
 
     phi = phiHbyA + p_rghEqn.flux();
 
+    // Calculate and relax the net pressure-buoyancy force
+    netForce.ref().relax
+    (
+        fvc::reconstruct((ghGradRhof + p_rghEqn.flux()/rhorAAtUf)),
+        p_rgh.relaxationFactor()
+    );
+
+    // Correct the momentum source with the pressure gradient flux
+    // calculated from the relaxed pressure
+    U = HbyA + rAAtU*netForce();
+    U.correctBoundaryConditions();
+    fvConstraints().constrain(U);
+
+    K = 0.5*magSqr(U);
+
     if (!mesh.schemes().steady())
     {
         p = p_rgh + rho*gh + pRef;
@@ -202,19 +218,8 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure()
 
     continuityErrors();
 
-    // Explicitly relax pressure for momentum corrector
-    p_rgh.relax();
-
     p = p_rgh + rho*gh + pRef;
 
-    // Correct the momentum source with the pressure gradient flux
-    // calculated from the relaxed pressure
-    U = HbyA + rAAtU*fvc::reconstruct((ghGradRhof + p_rghEqn.flux()/rhorAAtUf));
-    U.correctBoundaryConditions();
-    fvConstraints().constrain(U);
-
-    K = 0.5*magSqr(U);
-
     if (mesh.schemes().steady())
     {
         if (fvConstraints().constrain(p))
@@ -234,6 +239,9 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure()
         p_rgh.correctBoundaryConditions();
     }
 
+    // Optionally relax pressure for the thermophysics
+    p.relax();
+
     if (mesh.schemes().steady() || pimple.simpleRho() || adjustMass)
     {
         rho = thermo.rho();

applications/solvers/modules/fluid/isothermalFluid/isothermalFluid.C

@@ -170,6 +170,21 @@ Foam::solvers::isothermalFluid::isothermalFluid
             thermo,
             pimple.dict()
         );
+
+        netForce = new volVectorField
+        (
+            IOobject
+            (
+                "netForce",
+                runTime.timeName(),
+                mesh
+            ),
+            fvc::reconstruct
+            (
+                (-buoyancy->ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh))
+               *mesh.magSf()
+            )
+        );
     }
 
     if (mesh.dynamic())

applications/solvers/modules/fluid/isothermalFluid/isothermalFluid.H

@@ -147,6 +147,10 @@ protected:
 
     // Cached temporary fields
 
+        //- Momentum equation net force source term
+        //  Used for buoyant simulations only
+        tmp<volVectorField> netForce;
+
         //- Pointer to the surface momentum field
         //  used to recreate the flux after mesh-change
         autoPtr<surfaceVectorField> rhoUf;

applications/solvers/modules/fluid/isothermalFluid/momentumPredictor.C

@@ -51,13 +51,7 @@ void Foam::solvers::isothermalFluid::momentumPredictor()
             (
                 UEqn
              ==
-                fvc::reconstruct
-                (
-                    (
-                      - buoyancy->ghf*fvc::snGrad(rho)
-                      - fvc::snGrad(p_rgh)
-                    )*mesh.magSf()
-                )
+                netForce()
             );
         }
         else