Adding 'compressible' flag in pEq for chtMultiRegionFoam and buoyantPimpleFoam, rhoThermo and

transient based solvers to account for incompressible Eq of State laws. It avoids taking into account
the term ddt(rho) as mass contribution due to compressibility effects

applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H

@@ -90,3 +90,27 @@ volScalarField dpdt
 
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+
+if (p_rgh.needReference())
+{
+    setRefCell
+    (
+        p,
+        p_rgh,
+        pimple.dict(),
+        pRefCell,
+        pRefValue
+    );
+
+    p += dimensionedScalar
+    (
+        "p",
+        p.dimensions(),
+        pRefValue - getRefCellValue(p, pRefCell)
+    );
+}
+
+dimensionedScalar initialMass("initialMass", fvc::domainIntegrate(rho));

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -1,4 +1,9 @@
 {
+    bool closedVolume = p_rgh.needReference();
+
+    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    bool compressible = (compressibility.value() > SMALL);
+
     rho = thermo.rho();
 
     // Thermodynamic density needs to be updated by psi*d(p) after the
@@ -36,19 +41,27 @@
         )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
     );
 
-    fvScalarMatrix p_rghDDtEqn
+    tmp<fvScalarMatrix> p_rghDDtEqn
+    (
+        new fvScalarMatrix(p_rgh, dimMass/dimTime)
+    );
+
+    if (compressible)
+    {
+        p_rghDDtEqn =
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-      + fvc::div(phiHbyA)
         ==
             fvOptions(psi, p_rgh, rho.name())
         );
+    }
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            p_rghDDtEqn
+            p_rghDDtEqn()
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rAUf, p_rgh)
         );
 
@@ -81,6 +94,32 @@
         dpdt = fvc::ddt(p);
     }
 
+    if (compressible)
+    {
         #include "rhoEqn.H"
+    }
     #include "compressibleContinuityErrs.H"
+
+    if (closedVolume)
+    {
+        if (!compressible)
+        {
+            p += dimensionedScalar
+            (
+                "p",
+                p.dimensions(),
+                pRefValue - getRefCellValue(p, pRefCell)
+            );
+        }
+        else
+        {
+            p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+                /compressibility;
+            rho = thermo.rho();
+        }
+        p_rgh = p - rho*gh;
+    }
+
+    Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
+        << endl;
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

@@ -36,12 +36,21 @@
         )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
     );
 
+    tmp<fvScalarMatrix> p_rghDDtEqn
+    (
+        new fvScalarMatrix(p_rgh, dimMass/dimTime)
+    );
+
     {
-        fvScalarMatrix p_rghDDtEqn
+        if (compressible)
+        {
+            p_rghDDtEqn =
             (
                 fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-          + fvc::div(phiHbyA)
+             ==
+                fvOptions(psi, p_rgh, rho.name())
             );
+        }
 
         // Thermodynamic density needs to be updated by psi*d(p) after the
         // pressure solution - done in 2 parts. Part 1:
@@ -52,6 +61,7 @@
             fvScalarMatrix p_rghEqn
             (
                 p_rghDDtEqn
+              + fvc::div(phiHbyA)
               - fvm::laplacian(rhorAUf, p_rgh)
             );
 
@@ -93,8 +103,11 @@
         dpdt = fvc::ddt(p);
     }
 
+    if (compressible)
+    {
         // Solve continuity
         #include "rhoEqn.H"
+    }
 
     // Update continuity errors
     #include "compressibleContinuityErrors.H"

tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution

@@ -60,6 +60,8 @@ PIMPLE
     nOuterCorrectors 1;
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       1e5;
 }