rhoSimpleFoam: 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.

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -8,20 +8,23 @@
 
     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)
-           *fvc::flux(HbyA)
+            (fvc::interpolate(psi)/rhof)*phiHbyA
         );
-
-        MRF.makeRelative(fvc::interpolate(psi), phid);
+        phiHbyA -= fvc::interpolate(p)*phid;
 
         while (simple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
-                fvm::div(phid, p)
+                fvc::div(phiHbyA)
+              + fvm::div(phid, p)
               - fvm::laplacian(rhorAUf, p)
               ==
                 fvOptions(psi, p, rho.name())
@@ -40,7 +43,7 @@
 
             if (simple.finalNonOrthogonalIter())
             {
-                phi == pEqn.flux();
+                phi = phiHbyA + pEqn.flux();
             }
         }
     }

applications/solvers/compressible/rhoSimpleFoam/pcEqn.H

@@ -7,20 +7,19 @@ bool closedVolume = false;
 
 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)
-       *fvc::flux(HbyA)
+        (fvc::interpolate(psi)/rhof)*phiHbyA
     );
 
-    MRF.makeRelative(fvc::interpolate(psi), phid);
-
-    surfaceScalarField phic
-    (
-        "phic",
+    phiHbyA +=
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-    );
+      - fvc::interpolate(p)*phid;
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
@@ -31,7 +30,7 @@ if (simple.transonic())
         fvScalarMatrix pEqn
         (
             fvm::div(phid, p)
-          + fvc::div(phic)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rhorAtU, p)
          ==
             fvOptions(psi, p, rho.name())
@@ -50,7 +49,7 @@ if (simple.transonic())
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi == phic + pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }

tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution

@@ -40,7 +40,7 @@ SIMPLE
 {
     nNonOrthogonalCorrectors 0;
     pMinFactor      0.1;
-    pMaxFactor      1.5;
+    pMaxFactor      2;
     transonic       yes;
     consistent      yes;