From c3854609a41b9833df9fd0ef06796715009fa827 Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Wed, 9 Dec 2015 15:11:16 -0800
Subject: [PATCH] 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

---
 .../buoyantPimpleFoam/createFields.H          | 24 +++++++++
 .../heatTransfer/buoyantPimpleFoam/pEqn.H     | 53 ++++++++++++++++---
 .../chtMultiRegionFoam/fluid/pEqn.H           | 27 +++++++---
 .../hotRoom/system/fvSolution                 |  2 +
 4 files changed, 92 insertions(+), 14 deletions(-)

diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
index 8c6402f8f4..af2b8926aa 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/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));
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index e2f99671a3..76b4925917 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/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
     (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-      + fvc::div(phiHbyA)
-      ==
-        fvOptions(psi, p_rgh, rho.name())
+        new fvScalarMatrix(p_rgh, dimMass/dimTime)
     );
 
+    if (compressible)
+    {
+        p_rghDDtEqn =
+        (
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+        ==
+            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);
     }
 
-    #include "rhoEqn.H"
+    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;
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index 002faddb59..cd42cc00fd 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/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
-        (
-            fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-          + fvc::div(phiHbyA)
-        );
+        if (compressible)
+        {
+            p_rghDDtEqn =
+            (
+                fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+             ==
+                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);
     }
 
-    // Solve continuity
-    #include "rhoEqn.H"
+    if (compressible)
+    {
+        // Solve continuity
+        #include "rhoEqn.H"
+    }
 
     // Update continuity errors
     #include "compressibleContinuityErrors.H"
diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution
index ef63b73506..78d0b84613 100644
--- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution
@@ -60,6 +60,8 @@ PIMPLE
     nOuterCorrectors 1;
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       1e5;
 }