From 8504e436579f891c4f6ee6fd0b3ddd2a12bc898d 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 06fa890c0..5f092ada2 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 0dbde8296..af2005feb 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 2ffe1fe4b..03707f600 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;