Modification on rhoPimpleFoam pEq's for handling rho thermo and incompressible EoS. Adding rho limiters if p is limited.

This is important when LTS stepping or large Co number is used.

Updating rhoBuoyantPimpleFoam to handle closed domain for rho thermo and incompressible Eos.
Consolidating chtMultiRegionSimpleFoam and chtMultiRegionFoam pEqs to use the same formulation as rhoBuoyantPimpleFoam and
rhoBuoyantSimpleFoam

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -73,3 +73,26 @@ Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 
 #include "createMRF.H"
+
+
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        pimple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        pimple.dict(),
+        dimDensity,
+        0
+    )
+);

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -84,9 +84,6 @@ else
     }
 }
 
-// Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
-
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -101,12 +98,20 @@ K = 0.5*magSqr(U);
 if (pressureControl.limit(p))
 {
     p.correctBoundaryConditions();
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
 else if (pimple.SIMPLErho())
 {
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
+else
+{
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
+}
+
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
 
 if (thermo.dpdt())
 {

applications/solvers/compressible/rhoPimpleFoam/pcEqn.H

@@ -95,9 +95,6 @@ else
     }
 }
 
-// Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
-
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -112,12 +109,21 @@ K = 0.5*magSqr(U);
 if (pressureControl.limit(p))
 {
     p.correctBoundaryConditions();
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
+
 else if (pimple.SIMPLErho())
 {
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
+else
+{
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
+}
+
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
 
 if (thermo.dpdt())
 {

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -1,3 +1,7 @@
+
+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
@@ -41,6 +45,12 @@ while (pimple.correctNonOrthogonal())
       - fvm::laplacian(rhorAUf, p_rgh)
     );
 
+    p_rghEqn.setReference
+    (
+        pRefCell,
+        compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+    );
+
     p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
@@ -63,12 +73,40 @@ while (pimple.correctNonOrthogonal())
 p = p_rgh + rho*gh;
 
 // Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
+//thermo.correctRho(psi*p - psip0);
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+if (p_rgh.needReference())
+{
+    if (!compressible)
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }
+    else
+    {
+        p += (initialMass - fvc::domainIntegrate(psi*p))
+            /compressibility;
+        thermo.correctRho(psi*p - psip0);
+        rho = thermo.rho();
+
+    }
+    p_rgh = p - rho*gh;
+}
+else
+{
+    thermo.correctRho(psi*p - psip0);
+}
+
+rho = thermo.rho();
 
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
 }
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"

applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H

@@ -88,3 +88,26 @@ dimensionedScalar totalVolume = sum(mesh.V());
 
 #include "createMRF.H"
 #include "createRadiationModel.H"
+
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        simple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        simple.dict(),
+        dimDensity,
+        0
+    )
+);
+

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -1,7 +1,4 @@
 {
-    rho = thermo.rho();
-    rho.relax();
-
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
     volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
@@ -24,6 +21,8 @@
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
+    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    bool compressible = (compressibility.value() > SMALL);
 
     while (simple.correctNonOrthogonal())
     {
@@ -32,7 +31,12 @@
             fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.setReference
+        (
+            pRefCell,
+            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+        );
+
         p_rghEqn.solve();
 
         if (simple.finalNonOrthogonalIter())
@@ -51,12 +55,9 @@
         }
     }
 
-    #include "continuityErrs.H"
-
     p = p_rgh + rho*gh;
 
-    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
-    bool compressible = (compressibility.value() > SMALL);
+    #include "continuityErrs.H"
 
     // For closed-volume cases adjust the pressure level
     // to obey overall mass continuity

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

@@ -1,8 +1,10 @@
 {
+    /*
     rho = thermo.rho();
     rho = max(rho, rhoMin[i]);
     rho = min(rho, rhoMax[i]);
     rho.relax();
+    */
 
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -67,10 +69,22 @@
 
     // For closed-volume cases adjust the pressure level
     // to obey overall mass continuity
-    if (closedVolume && compressible)
+    if (closedVolume)
     {
-        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
-            /compressibility;
+        if (!compressible)
+        {
+            p += dimensionedScalar
+            (
+                "p",
+                p.dimensions(),
+                pRefValue - getRefCellValue(p, pRefCell)
+            );
+        }
+        else
+        {
+            p += (initialMass - fvc::domainIntegrate(psi*p))
+                /compressibility;
+        }
         p_rgh = p - rho*gh;
     }
 

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

@@ -19,6 +19,9 @@ List<bool> frozenFlowFluid(fluidRegions.size(), false);
 PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
 PtrList<fv::options> fluidFvOptions(fluidRegions.size());
 
+List<label> refCellFluid(fluidRegions.size());
+List<scalar> refValueFluid(fluidRegions.size());
+
 // Populate fluid field pointer lists
 forAll(fluidRegions, i)
 {
@@ -248,4 +251,20 @@ forAll(fluidRegions, i)
     );
 
     turbulence[i].validate();
+
+    refCellFluid[i] = 0;
+    refValueFluid[i] = 0.0;
+
+    if (p_rghFluid[i].needReference())
+    {
+        setRefCell
+        (
+            thermoFluid[i].p(),
+            p_rghFluid[i],
+            pimpleDict,
+            refCellFluid[i],
+            refValueFluid[i]
+        );
+    }
+
 }

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

@@ -4,6 +4,10 @@ bool compressible = (compressibility.value() > SMALL);
 
 rho = thermo.rho();
 
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
 volScalarField rAU("rAU", 1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
@@ -32,10 +36,6 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
       + fvc::div(phiHbyA)
     );
 
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution
-    const volScalarField psip0(psi*p);
-
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix p_rghEqn
@@ -44,6 +44,12 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
           - fvm::laplacian(rhorAUf, p_rgh)
         );
 
+        p_rghEqn.setReference
+        (
+            pRefCell,
+            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+        );
+
         p_rghEqn.solve
         (
             mesh.solver
@@ -62,6 +68,9 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
         if (nonOrth == nNonOrthCorr)
         {
             phi = phiHbyA + p_rghEqn.flux();
+
+            p_rgh.relax();
+
             U = HbyA
               + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
@@ -73,15 +82,10 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
     p = p_rgh + rho*gh;
 
     // Thermodynamic density update
-    thermo.correctRho(psi*p - psip0);
+    //thermo.correctRho(psi*p - psip0);
 }
 
 
-// Update pressure time derivative if needed
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p);
-}
 
 // Solve continuity
 #include "rhoEqn.H"
@@ -91,10 +95,35 @@ if (thermo.dpdt())
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
-if (closedVolume && compressible)
+if (closedVolume)
 {
-    p += (initialMass - fvc::domainIntegrate(thermo.rho()))
-        /compressibility;
-    rho = thermo.rho();
+    if (!compressible)
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }
+    else
+    {
+        p += (initialMass - fvc::domainIntegrate(psi*p))
+            /compressibility;
+        thermo.correctRho(psi*p - psip0);
+        rho = thermo.rho();
+    }
     p_rgh = p - rho*gh;
 }
+else
+{
+    thermo.correctRho(psi*p - psip0);
+}
+
+rho = thermo.rho();
+
+// Update pressure time derivative if needed
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}

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

@@ -33,3 +33,6 @@
     );
 
     const bool frozenFlow = frozenFlowFluid[i];
+
+    const label pRefCell = refCellFluid[i];
+    const scalar pRefValue = refValueFluid[i];

src/thermophysicalModels/basic/fluidThermo/fluidThermo.H

@@ -105,7 +105,12 @@ public:
 
             //- Add the given density correction to the density field.
             //  Used to update the density field following pressure solution
-            virtual void correctRho(const volScalarField& deltaRho) = 0;
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho,
+                const dimensionedScalar& rhoMin,
+                const dimensionedScalar& rhoMax
+            ) = 0;
 
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const = 0;

src/thermophysicalModels/basic/psiThermo/psiThermo.C

@@ -102,10 +102,14 @@ Foam::tmp<Foam::scalarField> Foam::psiThermo::rho(const label patchi) const
 }
 
 
-void Foam::psiThermo::correctRho(const Foam::volScalarField& deltaRho)
+void Foam::psiThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho,
+    const dimensionedScalar& rhoMin,
+    const dimensionedScalar& rhoMax
+)
 {}
 
-
 const Foam::volScalarField& Foam::psiThermo::psi() const
 {
     return psi_;

src/thermophysicalModels/basic/psiThermo/psiThermo.H

@@ -118,7 +118,12 @@ public:
             //- Add the given density correction to the density field.
             //  Used to update the density field following pressure solution.
             //  For psiThermo does nothing.
-            virtual void correctRho(const volScalarField& deltaRho);
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho,
+                const dimensionedScalar& rhoMin,
+                const dimensionedScalar& rhoMax
+            );
 
             //- Density [kg/m^3] - uses current value of pressure
             virtual tmp<volScalarField> rho() const;

src/thermophysicalModels/basic/rhoThermo/rhoThermo.C

@@ -173,12 +173,18 @@ Foam::volScalarField& Foam::rhoThermo::rho()
 }
 
 
-void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)
+void Foam::rhoThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho,
+    const dimensionedScalar& rhoMin,
+    const dimensionedScalar& rhoMax
+)
 {
     rho_ += deltaRho;
+    rho_ = max(rho_, rhoMin);
+    rho_ = min(rho_, rhoMax);
 }
 
-
 const Foam::volScalarField& Foam::rhoThermo::psi() const
 {
     return psi_;

src/thermophysicalModels/basic/rhoThermo/rhoThermo.H

@@ -138,7 +138,13 @@ public:
 
             //- Add the given density correction to the density field.
             //  Used to update the density field following pressure solution
-            virtual void correctRho(const volScalarField& deltaRho);
+            //  Limit thermo rho between rhoMin and rhoMax
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho,
+                const dimensionedScalar& rhoMin,
+                const dimensionedScalar& rhoMax
+            );
 
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const;