rhoPimpleFoam: Improved support for compressible liquids

See tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq for exapmle

pimpleControl: Added SIMPLErho option for running in SIMPLE mode

with large time-step/Courant number and relaxation.  With this option the
density is updated from thermodynamics rather than continuity after the pressure
equation which is better behaved if pressure is relaxed and/or solved to a
loose relative tolerance.  The need for this option is demonstrated in the
tutorials/compressible/rhoPimpleFoam/RAS/angledDuct tutorial which is unstable
without the option.

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -42,6 +42,8 @@ volVectorField U
 
 pressureControl pressureControl(p, rho, pimple.dict(), false);
 
+mesh.setFluxRequired(p.name());
+
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -54,8 +56,6 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-mesh.setFluxRequired(p.name());
-
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
 (

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -1,4 +1,11 @@
-rho = thermo.rho();
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -31,17 +38,17 @@ if (pimple.transonic())
 
     phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA) + fvm::div(phid, p)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     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())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
@@ -56,16 +63,17 @@ if (pimple.transonic())
 }
 else
 {
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -76,6 +84,9 @@ else
     }
 }
 
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -92,10 +103,9 @@ if (pressureControl.limit(p))
     p.correctBoundaryConditions();
     rho = thermo.rho();
 }
-
-if (!pimple.transonic())
+else if (pimple.SIMPLErho())
 {
-    rho.relax();
+    rho = thermo.rho();
 }
 
 if (thermo.dpdt())

applications/solvers/compressible/rhoPimpleFoam/pcEqn.H

@@ -1,4 +1,11 @@
-rho = thermo.rho();
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
@@ -39,17 +46,17 @@ if (pimple.transonic())
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA) + fvm::div(phid, p)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAtU, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
 
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
@@ -67,16 +74,17 @@ else
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAtU, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -87,6 +95,9 @@ else
     }
 }
 
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -103,10 +114,9 @@ if (pressureControl.limit(p))
     p.correctBoundaryConditions();
     rho = thermo.rho();
 }
-
-if (!pimple.transonic())
+else if (pimple.SIMPLErho())
 {
-    rho.relax();
+    rho = thermo.rho();
 }
 
 if (thermo.dpdt())

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

@@ -36,7 +36,7 @@ if (pimple.transonic())
     {
         fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, p)
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phiHbyA)
           + fvm::div(phid, p)
           - fvm::laplacian(rhorAUf, p)
@@ -62,7 +62,7 @@ else
         // Pressure corrector
         fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, p)
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phiHbyA)
           - fvm::laplacian(rhorAUf, p)
          ==
@@ -95,11 +95,6 @@ if (pressureControl.limit(p))
     rho = thermo.rho();
 }
 
-if (!pimple.transonic())
-{
-    rho.relax();
-}
-
 {
     rhoUf = fvc::interpolate(rho*U);
     surfaceVectorField n(mesh.Sf()/mesh.magSf());

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -119,6 +119,8 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"