rhoPimpleFoam: Updated transonic option to be consistent with sonicFoam

Improves stability on start-up and allows running at slightly larger time-steps.

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -1,3 +1,5 @@
+rho = thermo.rho();
+
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -10,7 +12,7 @@ if (pimple.nCorrPISO() <= 1)
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
-    fvc::flux(rho*HbyA)
+    fvc::interpolate(rho)*fvc::flux(HbyA)
   + rhorAUf*fvc::ddtCorr(rho, U, phi)
 );
 
@@ -26,7 +28,8 @@ if (pimple.transonic())
         "phid",
         (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
-    phiHbyA -= fvc::interpolate(p)*phid;
+
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -84,9 +87,11 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-pressureControl.limit(p);
-p.correctBoundaryConditions();
-rho = thermo.rho();
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
 
 if (!pimple.transonic())
 {

applications/solvers/compressible/rhoPimpleFoam/pcEqn.H

@@ -1,3 +1,5 @@
+rho = thermo.rho();
+
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -11,7 +13,7 @@ surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (
-        fvc::flux(rho*HbyA)
+        fvc::interpolate(rho)*fvc::flux(HbyA)
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
     )
 );
@@ -33,7 +35,7 @@ if (pimple.transonic())
 
     phiHbyA +=
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-      - fvc::interpolate(p)*phid;
+      - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
@@ -96,9 +98,11 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-pressureControl.limit(p);
-p.correctBoundaryConditions();
-rho = thermo.rho();
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
 
 if (!pimple.transonic())
 {

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -1,25 +1,25 @@
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+tUEqn.clear();
+
+bool closedVolume = false;
+
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(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())
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-    tUEqn.clear();
-
-    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 phid
     (
         "phid",
         (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
-        phiHbyA -= fvc::interpolate(p)*phid;
+
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     while (simple.correctNonOrthogonal())
     {
@@ -48,9 +48,9 @@
             phi = phiHbyA + pEqn.flux();
         }
     }
-    }
-    else
-    {
+}
+else
+{
     closedVolume = adjustPhi(phiHbyA, U, p);
 
     while (simple.correctNonOrthogonal())
@@ -76,34 +76,35 @@
             phi = phiHbyA + pEqn.flux();
         }
     }
-    }
+}
 
+#include "incompressible/continuityErrs.H"
 
-    #include "incompressible/continuityErrs.H"
+// Explicitly relax pressure for momentum corrector
+p.relax();
 
-    // Explicitly relax pressure for momentum corrector
-    p.relax();
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
 
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    fvOptions.correct(U);
+bool pLimited = pressureControl.limit(p);
 
-    pressureControl.limit(p);
-
-    // For closed-volume cases adjust the pressure and density levels
-    // to obey overall mass continuity
-    if (closedVolume)
-    {
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+if (closedVolume)
+{
     p += (initialMass - fvc::domainIntegrate(psi*p))
         /fvc::domainIntegrate(psi);
-    }
-
-    p.correctBoundaryConditions();
-
-    rho = thermo.rho();
-
-    if (!simple.transonic())
-    {
-        rho.relax();
-    }
+}
+
+if (pLimited || closedVolume)
+{
+    p.correctBoundaryConditions();
+}
+
+rho = thermo.rho();
+
+if (!simple.transonic())
+{
+    rho.relax();
 }

applications/solvers/compressible/rhoSimpleFoam/pcEqn.H

@@ -1,3 +1,5 @@
+rho = thermo.rho();
+
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -5,7 +7,7 @@ tUEqn.clear();
 
 bool closedVolume = false;
 
-surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
 volScalarField rhorAtU("rhorAtU", rho*rAtU);
@@ -23,7 +25,7 @@ if (simple.transonic())
 
     phiHbyA +=
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-      - fvc::interpolate(p)*phid;
+      - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
@@ -98,7 +100,7 @@ U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 
-pressureControl.limit(p);
+bool pLimited = pressureControl.limit(p);
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
@@ -108,9 +110,11 @@ if (closedVolume)
         /fvc::domainIntegrate(psi);
 }
 
-p.correctBoundaryConditions();
+if (pLimited || closedVolume)
+{
+    p.correctBoundaryConditions();
+}
 
-// Recalculate density from the relaxed pressure
 rho = thermo.rho();
 
 if (!simple.transonic())

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

@@ -207,14 +207,24 @@ Foam::pressureControl::pressureControl
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::pressureControl::limit(volScalarField& p) const
+bool Foam::pressureControl::limit(volScalarField& p) const
 {
-    Info<< "pressureControl: p max/min "
-        << max(p).value() << " "
-        << min(p).value() << endl;
+    scalar pMax = max(p).value();
+    scalar pMin = min(p).value();
+
+    if (pMax > pMax_.value() || pMin < pMin_.value())
+    {
+        Info<< "pressureControl: p max/min " << pMax << " " << pMin << endl;
 
         p = max(p, pMin_);
         p = min(p, pMax_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 

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

@@ -89,8 +89,8 @@ public:
         //- Return the pressure reference level
         inline scalar refValue() const;
 
-        //- Limit the pressure
-        void limit(volScalarField& p) const;
+        //- Limit the pressure if necessary and return true if so
+        bool limit(volScalarField& p) const;
 };