rhoPimpleFoam: Added support for transonic flow of liquids and real gases

Both stardard SIMPLE and the SIMPLEC (using the 'consistent' option in
fvSolution) are now supported for both subsonic and transonic flow of all
fluid types.

rhoPimpleFoam now instantiates the lower-level fluidThermo which instantiates
either a psiThermo or rhoThermo according to the 'type' specification in
thermophysicalProperties, see also commit 655fc78748

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -2,11 +2,11 @@
 
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<psiThermo> pThermo
+autoPtr<fluidThermo> pThermo
 (
-    psiThermo::New(mesh)
+    fluidThermo::New(mesh)
 );
-psiThermo& thermo = pThermo();
+fluidThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
 volScalarField& p = thermo.p();
@@ -40,27 +40,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-dimensionedScalar rhoMax
+pressureControl pressureControl(p, rho, pimple.dict(), false);
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        pimple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        pimple.dict(),
-        dimDensity,
-        0
-    )
-);
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -1,8 +1,3 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -12,55 +7,54 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (
+        fvc::flux(rho*HbyA)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+    )
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
 if (pimple.transonic())
 {
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
-       *(
-            fvc::flux(HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
-        )
     );
-
+    phiHbyA -= fvc::interpolate(p)*phid;
-    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
           + fvm::div(phid, p)
           - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
 
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi == pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
 else
 {
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
-    );
-
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -87,19 +81,20 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-Info<< "rho max/min : " << max(rho).value()
-    << " " << min(rho).value() << endl;
-
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+pressureControl.limit(p);
+p.correctBoundaryConditions();
+rho = thermo.rho();
+
+if (!pimple.transonic())
+{
+    rho.relax();
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);

applications/solvers/compressible/rhoPimpleFoam/pcEqn.H

@@ -1,8 +1,3 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -12,72 +7,64 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (
+        fvc::flux(rho*HbyA)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+    )
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
+
 if (pimple.transonic())
 {
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
-       *(
-            fvc::flux(HbyA)
-          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
-           /fvc::interpolate(rho)
-        )
     );
 
-    MRF.makeRelative(fvc::interpolate(psi), phid);
+    phiHbyA +=
-
-    surfaceScalarField phic
-    (
-        "phic",
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-    );
+      - fvc::interpolate(p)*phid;
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField rhorAtU("rhorAtU", rho*rAtU);
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
           + fvm::div(phid, p)
-          + fvc::div(phic)
           - fvm::laplacian(rhorAtU, p)
          ==
             fvOptions(psi, p, rho.name())
         );
 
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi == phic + pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
 else
 {
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
-        )
-    );
-
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField rhorAtU("rhorAtU", rho*rAtU);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -109,19 +96,16 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-if (thermo.dpdt())
+pressureControl.limit(p);
-{
+p.correctBoundaryConditions();
-    dpdt = fvc::ddt(p);
-}
-
-// Recalculate density from the relaxed pressure
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
 
 if (!pimple.transonic())
 {
     rho.relax();
 }
 
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H

@@ -1,8 +1,3 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -12,55 +7,53 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    fvc::flux(rho*HbyA)
+  + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
+);
+
+fvc::makeRelative(phiHbyA, rho, U);
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
 if (pimple.transonic())
 {
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
-       *(
-            fvc::flux(HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
-        )
     );
-
+    phiHbyA -= fvc::interpolate(p)*phid;
-    fvc::makeRelative(phid, psi, U);
-    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
           + fvm::div(phid, p)
           - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
 
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi == pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
 else
 {
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        fvc::flux(rho*HbyA)
-      + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
-    );
-
-    fvc::makeRelative(phiHbyA, rho, U);
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
     while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
@@ -88,19 +81,20 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-Info<< "rho max/min : " << max(rho).value()
-    << " " << min(rho).value() << endl;
-
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+pressureControl.limit(p);
+p.correctBoundaryConditions();
+rho = thermo.rho();
+
+if (!pimple.transonic())
+{
+    rho.relax();
+}
+
 {
     rhoUf = fvc::interpolate(rho*U);
     surfaceVectorField n(mesh.Sf()/mesh.magSf());

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,10 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    rhoPimpleFoam
+    rhoPimpleDyMFoam
-
-Group
-    grpCompressibleSolvers grpMovingMeshSolvers
 
 Description
     Transient solver for turbulent flow of compressible fluids for HVAC and
@@ -38,10 +35,11 @@ Description
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
-#include "psiThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "CorrectPhi.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,10 +34,11 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "psiThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"

applications/solvers/compressible/rhoSimpleFoam/createFields.H

@@ -7,6 +7,8 @@ autoPtr<fluidThermo> pThermo
 fluidThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
+volScalarField& p = thermo.p();
+
 volScalarField rho
 (
     IOobject
@@ -20,8 +22,6 @@ volScalarField rho
     thermo.rho()
 );
 
-volScalarField& p = thermo.p();
-
 Info<< "Reading field U\n" << endl;
 volVectorField U
 (

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -6,16 +6,18 @@
 
     bool closedVolume = false;
 
+    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
     if (simple.transonic())
     {
-        surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-        surfaceScalarField rhof(fvc::interpolate(rho));
-        MRF.makeRelative(rhof, phiHbyA);
-
         surfaceScalarField phid
         (
             "phid",
-            (fvc::interpolate(psi)/rhof)*phiHbyA
+            (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
         );
         phiHbyA -= fvc::interpolate(p)*phid;
 
@@ -49,14 +51,8 @@
     }
     else
     {
-        surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-        MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
         closedVolume = adjustPhi(phiHbyA, U, p);
 
-        // Update the pressure BCs to ensure flux consistency
-        constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
         while (simple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn

applications/solvers/compressible/rhoSimpleFoam/pcEqn.H

@@ -5,16 +5,20 @@ tUEqn.clear();
 
 bool closedVolume = false;
 
+surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
+
 if (simple.transonic())
 {
-    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-    surfaceScalarField rhof(fvc::interpolate(rho));
-    MRF.makeRelative(rhof, phiHbyA);
-
     surfaceScalarField phid
     (
         "phid",
-        (fvc::interpolate(psi)/rhof)*phiHbyA
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
 
     phiHbyA +=
@@ -23,14 +27,12 @@ if (simple.transonic())
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField rhorAtU("rhorAtU", rho*rAtU);
-
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::div(phid, p)
+            fvc::div(phiHbyA)
-          + fvc::div(phiHbyA)
+          + fvm::div(phid, p)
           - fvm::laplacian(rhorAtU, p)
          ==
             fvOptions(psi, p, rho.name())
@@ -55,19 +57,11 @@ if (simple.transonic())
 }
 else
 {
-    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
     closedVolume = adjustPhi(phiHbyA, U, p);
 
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField rhorAtU("rhorAtU", rho*rAtU);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
-
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H

@@ -1,59 +0,0 @@
-Info<< "Reading thermophysical properties\n" << endl;
-
-autoPtr<fluidThermo> pThermo
-(
-    fluidThermo::New(mesh)
-);
-fluidThermo& thermo = pThermo();
-thermo.validate(args.executable(), "h", "e");
-
-volScalarField rho
-(
-    IOobject
-    (
-        "rho",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    thermo.rho()
-);
-
-volScalarField& p = thermo.p();
-
-Info<< "Reading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "compressibleCreatePhi.H"
-
-pressureControl pressureControl(p, rho, simple.dict());
-
-mesh.setFluxRequired(p.name());
-
-Info<< "Creating turbulence model\n" << endl;
-autoPtr<compressible::turbulenceModel> turbulence
-(
-    compressible::turbulenceModel::New
-    (
-        rho,
-        U,
-        phi,
-        thermo
-    )
-);
-
-dimensionedScalar initialMass = fvc::domainIntegrate(rho);
-
-#include "createMRF.H"

src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C

@@ -32,105 +32,171 @@ Foam::pressureControl::pressureControl
 (
     const volScalarField& p,
     const volScalarField& rho,
-    const dictionary& dict
+    const dictionary& dict,
+    const bool pRefRequired
 )
 :
-    refCell_(0),
+    refCell_(-1),
     refValue_(0),
     pMax_("pMax", dimPressure, 0),
     pMin_("pMin", dimPressure, GREAT)
 {
-    if (setRefCell(p, dict, refCell_, refValue_))
+    bool pLimits = false;
+
+    // Set the reference cell and value for closed domain simulations
+    if (pRefRequired && setRefCell(p, dict, refCell_, refValue_))
     {
+        pLimits = true;
+
         pMax_.value() = refValue_;
         pMin_.value() = refValue_;
     }
 
-    const volScalarField::Boundary& pbf = p.boundaryField();
+    if (dict.found("pMax") && dict.found("pMin"))
-    const volScalarField::Boundary& rhobf = rho.boundaryField();
-
-    scalar rhoRefMax = -GREAT;
-    scalar rhoRefMin = GREAT;
-    bool rhoLimits = false;
-
-    forAll(pbf, patchi)
-    {
-        if (pbf[patchi].fixesValue())
-        {
-            rhoLimits = true;
-
-            pMax_.value() = max(pMax_.value(), max(pbf[patchi]));
-            pMin_.value() = min(pMin_.value(), min(pbf[patchi]));
-
-            rhoRefMax = max(rhoRefMax, max(rhobf[patchi]));
-            rhoRefMin = min(rhoRefMin, min(rhobf[patchi]));
-        }
-    }
-
-    if (dict.found("pMax"))
     {
         pMax_.value() = readScalar(dict.lookup("pMax"));
-    }
-    else if (dict.found("pMaxFactor"))
-    {
-        const scalar pMaxFactor(readScalar(dict.lookup("pMaxFactor")));
-        pMax_ *= pMaxFactor;
-    }
-    else if (dict.found("rhoMax"))
-    {
-        // For backward-compatibility infer the pMax from rhoMax
-
-        IOWarningInFunction(dict)
-            << "'rhoMax' specified rather than 'pMax' or 'pMaxFactor'" << nl
-            << "    This is supported for backward-compatibility but "
-               "'pMax' or 'pMaxFactor' are more reliable." << endl;
-
-        if (!rhoLimits)
-        {
-            FatalIOErrorInFunction(dict)
-                << "'rhoMax' specified rather than 'pMaxFactor'" << nl
-                << "    but the corresponding reference density cannot"
-                   " be evaluated from the boundary conditions." << nl
-                << "Please specify 'pMaxFactor' rather than 'rhoMax'"
-                << exit(FatalError);
-        }
-
-        dimensionedScalar rhoMax("rhoMax", dimDensity, dict);
-
-        pMax_ *= max(rhoMax.value()/rhoRefMax, 1);
-    }
-
-    if (dict.found("pMin"))
-    {
         pMin_.value() = readScalar(dict.lookup("pMin"));
     }
-    else if (dict.found("pMinFactor"))
+    else
     {
-        const scalar pMinFactor(readScalar(dict.lookup("pMinFactor")));
+        const volScalarField::Boundary& pbf = p.boundaryField();
-        pMin_ *= pMinFactor;
+        const volScalarField::Boundary& rhobf = rho.boundaryField();
-    }
-    else if (dict.found("rhoMin"))
-    {
-        // For backward-compatibility infer the pMin from rhoMin
 
-        IOWarningInFunction(dict)
+        scalar rhoRefMax = -GREAT;
-            << "'rhoMin' specified rather than 'pMin' or 'pMinFactor'" << nl
+        scalar rhoRefMin = GREAT;
-            << "    This is supported for backward-compatibility but"
+        bool rhoLimits = false;
-               "'pMin' or 'pMinFactor' are more reliable." << endl;
 
-        if (!rhoLimits)
+        forAll(pbf, patchi)
         {
-            FatalIOErrorInFunction(dict)
+            if (pbf[patchi].fixesValue())
-                << "'rhoMin' specified rather than 'pMinFactor'" << nl
+            {
-                << "    but the corresponding reference density cannot"
+                pLimits = true;
-                   " be evaluated from the boundary conditions." << nl
+                rhoLimits = true;
-                << "Please specify 'pMinFactor' rather than 'rhoMin'"
+
-                << exit(FatalError);
+                pMax_.value() = max(pMax_.value(), max(pbf[patchi]));
+                pMin_.value() = min(pMin_.value(), min(pbf[patchi]));
+
+                rhoRefMax = max(rhoRefMax, max(rhobf[patchi]));
+                rhoRefMin = min(rhoRefMin, min(rhobf[patchi]));
+            }
         }
 
-        dimensionedScalar rhoMin("rhoMin", dimDensity, dict);
+        reduce(rhoLimits, andOp<bool>());
+        if (rhoLimits)
+        {
+            reduce(pMax_.value(), maxOp<scalar>());
+            reduce(pMin_.value(), minOp<scalar>());
 
-        pMin_ *= min(rhoMin.value()/rhoRefMin, 1);
+            reduce(rhoRefMax, maxOp<scalar>());
+            reduce(rhoRefMin, minOp<scalar>());
+        }
+
+        if (dict.found("pMax"))
+        {
+            pMax_.value() = readScalar(dict.lookup("pMax"));
+        }
+        else if (dict.found("pMaxFactor"))
+        {
+            if (!pLimits)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "'pMaxFactor' specified rather than 'pMax'" << nl
+                    << "    but the corresponding reference pressure cannot"
+                       " be evaluated from the boundary conditions." << nl
+                    << "    Please specify 'pMax' rather than 'pMaxFactor'"
+                    << exit(FatalIOError);
+            }
+
+            const scalar pMaxFactor(readScalar(dict.lookup("pMaxFactor")));
+            pMax_ *= pMaxFactor;
+        }
+        else if (dict.found("rhoMax"))
+        {
+            // For backward-compatibility infer the pMax from rhoMax
+
+            IOWarningInFunction(dict)
+                 << "'rhoMax' specified rather than 'pMax' or 'pMaxFactor'"
+                 << nl
+                 << "    This is supported for backward-compatibility but "
+                    "'pMax' or 'pMaxFactor' are more reliable." << endl;
+
+            if (!pLimits)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "'rhoMax' specified rather than 'pMax'" << nl
+                    << "    but the corresponding reference pressure cannot"
+                       " be evaluated from the boundary conditions." << nl
+                    << "    Please specify 'pMax' rather than 'rhoMax'"
+                    << exit(FatalIOError);
+            }
+
+            if (!rhoLimits)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "'rhoMax' specified rather than 'pMaxFactor'" << nl
+                    << "    but the corresponding reference density cannot"
+                       " be evaluated from the boundary conditions." << nl
+                    << "    Please specify 'pMaxFactor' rather than 'rhoMax'"
+                    << exit(FatalIOError);
+            }
+
+            dimensionedScalar rhoMax("rhoMax", dimDensity, dict);
+
+            pMax_ *= max(rhoMax.value()/rhoRefMax, 1);
+        }
+
+        if (dict.found("pMin"))
+        {
+            pMin_.value() = readScalar(dict.lookup("pMin"));
+        }
+        else if (dict.found("pMinFactor"))
+        {
+            if (!pLimits)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "'pMinFactor' specified rather than 'pMin'" << nl
+                    << "    but the corresponding reference pressure cannot"
+                       " be evaluated from the boundary conditions." << nl
+                    << "    Please specify 'pMin' rather than 'pMinFactor'"
+                    << exit(FatalIOError);
+            }
+
+            const scalar pMinFactor(readScalar(dict.lookup("pMinFactor")));
+            pMin_ *= pMinFactor;
+        }
+        else if (dict.found("rhoMin"))
+        {
+            // For backward-compatibility infer the pMin from rhoMin
+
+            IOWarningInFunction(dict)
+                << "'rhoMin' specified rather than 'pMin' or 'pMinFactor'" << nl
+                << "    This is supported for backward-compatibility but"
+                   "'pMin' or 'pMinFactor' are more reliable." << endl;
+
+            if (!pLimits)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "'rhoMin' specified rather than 'pMin'" << nl
+                    << "    but the corresponding reference pressure cannot"
+                       " be evaluated from the boundary conditions." << nl
+                    << "    Please specify 'pMin' rather than 'rhoMin'"
+                    << exit(FatalIOError);
+            }
+
+            if (!rhoLimits)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "'rhoMin' specified rather than 'pMinFactor'" << nl
+                    << "    but the corresponding reference density cannot"
+                       " be evaluated from the boundary conditions." << nl
+                    << "    Please specify 'pMinFactor' rather than 'rhoMin'"
+                    << exit(FatalIOError);
+            }
+
+            dimensionedScalar rhoMin("rhoMin", dimDensity, dict);
+
+            pMin_ *= min(rhoMin.value()/rhoRefMin, 1);
+        }
     }
 
     Info<< "pressureControl" << nl

src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H

@@ -76,7 +76,8 @@ public:
         (
             const volScalarField& p,
             const volScalarField& rho,
-            const dictionary& dict
+            const dictionary& dict,
+            const bool pRefRequired = true
         );
 
 

tutorials/compressible/rhoPimpleDyMFoam/annularThermalMixer/system/fvSolution

@@ -63,8 +63,9 @@ PIMPLE
     nOuterCorrectors    3;
     nCorrectors         1;
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.5;
+
-    rhoMax          2.0;
+    pMaxFactor          1.2;
+    pMinFactor          0.8;
 }
 
 relaxationFactors

tutorials/compressible/rhoPimpleFoam/LES/pitzDaily/system/fvSolution

@@ -52,8 +52,9 @@ PIMPLE
     nOuterCorrectors 3;
     nCorrectors     1;
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.5;
+
-    rhoMax          2.0;
+    pMinFactor      0.5;
+    pMaxFactor      2.0;
 }
 
 relaxationFactors

tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/constant/thermophysicalProperties

@@ -23,7 +23,7 @@ thermoType
     thermo          hConst;
     equationOfState perfectGas;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture

tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSchemes

@@ -30,8 +30,9 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
+    div(phiv,p)     Gauss linear;
     div(phi,K)      Gauss linear;
-    div(phi,h)      Gauss upwind;
+    div(phi,e)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;

tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution

@@ -28,11 +28,10 @@ solvers
     pFinal
     {
         $p;
-        tolerance       1e-07;
         relTol          0;
     }
 
-    "(rho|U|h|k|epsilon|omega)"
+    "(rho|U|e|k|epsilon|omega)"
     {
         solver          smoothSolver;
         smoother        symGaussSeidel;
@@ -40,10 +39,9 @@ solvers
         relTol          0.1;
     }
 
-    "(rho|U|h|k|epsilon|omega)Final"
+    "(rho|U|e|k|epsilon|omega)Final"
     {
         $U;
-        tolerance       1e-06;
         relTol          0;
     }
 }
@@ -57,8 +55,8 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
     consistent          yes;
 
-    rhoMin          0.4;
+    pMaxFactor          1.5;
-    rhoMax          2.0;
+    pMinFactor          0.9;
 
     residualControl
     {
@@ -76,13 +74,13 @@ relaxationFactors
 {
     fields
     {
-        "p.*"           0.9;
+        "p.*"           1;
         "rho.*"         1;
     }
     equations
     {
         "U.*"           0.9;
-        "h.*"           0.7;
+        "e.*"           0.7;
         "(k|epsilon|omega).*" 0.8;
     }
 }

tutorials/compressible/rhoPimpleFoam/RAS/angledDuctLTS/system/fvSolution

@@ -56,8 +56,8 @@ PIMPLE
     nCorrectors       1;
     nNonOrthogonalCorrectors 0;
 
-    rhoMin            0.5;
+    pMaxFactor          1.5;
-    rhoMax            2.0;
+    pMinFactor          0.9;
 
     maxCo             0.2;
     rDeltaTSmoothingCoeff 0.1;

tutorials/compressible/rhoPimpleFoam/RAS/cavity/system/fvSolution

@@ -59,8 +59,9 @@ PIMPLE
     nOuterCorrectors 1;
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.5;
+
-    rhoMax          2.0;
+    pMax                1.2e5;
+    pMin                0.8e5;
 }
 
 

tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/fvSolution

@@ -84,14 +84,13 @@ solvers
 
 PIMPLE
 {
-    nOuterCorrectors 1;
+    nOuterCorrectors    1;
-    nCorrectors     2;
+    nCorrectors         2;
-    nNonOrthogonalCorrectors 0;
+    nNonOrthogonalCorrectors  0;
-    momentumPredictor yes;
+    momentumPredictor   yes;
-    rhoMin          0.5;
+
-    rhoMax          2.0;
+    pMax                1.2e5;
-    pRefCell        0;
+    pMin                0.8e5;
-    pRefValue       1e5;
 }
 
 relaxationFactors

tutorials/compressible/rhoPimpleFoam/laminar/helmholtzResonance/system/fvSolution

@@ -52,8 +52,9 @@ PIMPLE
     nOuterCorrectors 3;
     nCorrectors     1;
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.5;
+
-    rhoMax          2.0;
+    pMinFactor      0.5;
+    pMaxFactor      2.0;
 }
 
 relaxationFactors