From 7d6845defaa6f0840e68ebf239d024fa0d739127 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 24 Feb 2017 16:20:06 +0000
Subject: [PATCH] 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.
---
 .../solvers/compressible/rhoSimpleFoam/pEqn.H | 15 +++++++++------
 .../compressible/rhoSimpleFoam/pcEqn.H        | 19 +++++++++----------
 .../squareBend/system/fvSolution              |  2 +-
 3 files changed, 19 insertions(+), 17 deletions(-)

diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 06fa890c05..5f092ada2f 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/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();
             }
         }
     }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
index 0dbde82960..af2005febd 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
+++ b/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();
         }
     }
 }
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution
index 2ffe1fe4be..03707f600d 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution
+++ b/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;