From f42adf48abc51b660e07ce8e38cf853307fc871d Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Fri, 5 Feb 2016 17:48:36 +0000
Subject: [PATCH] Bugfix: regular formulation of enthalpy source terms

---
 .../steadyCompressibleMRFFoam/createFields.H      |  2 +-
 .../compressible/steadyCompressibleMRFFoam/hEqn.H |  7 ++++++-
 .../compressible/steadyCompressibleMRFFoam/iEqn.H | 15 +++------------
 .../compressible/steadyCompressibleMRFFoam/pEqn.H |  3 ++-
 4 files changed, 12 insertions(+), 15 deletions(-)

diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/createFields.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/createFields.H
index 00073568e..78c32d810 100644
--- a/applications/solvers/compressible/steadyCompressibleMRFFoam/createFields.H
+++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/createFields.H
@@ -93,4 +93,4 @@
         ),
         h
     );
-    i -= 0.5*magSqr(Urot);
+    i == h - 0.5*(magSqr(Urot) - magSqr(Urel));
diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H
index 08392027f..6bf1b5f52 100644
--- a/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H
+++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H
@@ -2,7 +2,12 @@
     // Solve the enthalpy equation
     T.storePrevIter();
 
-    surfaceScalarField faceU = phi/fvc::interpolate(rho);
+    // Calculate face velocity from flux
+    surfaceScalarField faceU
+    (
+        "faceU",
+        phi/fvc::interpolate(rho)
+    );
 
     fvScalarMatrix hEqn
     (
diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H
index ca500d7af..ef6c36c13 100644
--- a/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H
+++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H
@@ -26,9 +26,10 @@
         fvm::ddt(rho, i)
       + fvm::div(phi, i)
       - fvm::laplacian(turbulence->alphaEff(), i)
-        // u & gradP term (steady-state formulation)
-      + fvm::SuSp((fvc::div(faceU, p, "div(U,p)") - fvc::div(faceU)*p)/i, i)
      ==
+        // ddt(p) term removed: steady-state.  HJ, 27/Apr/2010
+        fvc::div(faceU, p, "div(U,p)")
+      - p*fvc::div(faceU)
         // Viscous heating: note sign (devRhoReff has a minus in it)
       - (turbulence->devRhoReff() && fvc::grad(Urel))
     );
@@ -42,16 +43,6 @@
     h = i + 0.5*magSqr(Urot);
     h.correctBoundaryConditions();
 
-    // Bound the enthalpy using TMin and TMax
-    volScalarField Cp = thermo.Cp();
-
-    h = Foam::min(h, TMax*Cp);
-    h = Foam::max(h, TMin*Cp);
-    h.correctBoundaryConditions();
-
-    // Re-initialise rothalpy based on limited enthalpy
-    i = h - 0.5*magSqr(Urot);
-
     thermo.correct();
     psis = thermo.psi()/thermo.Cp()*thermo.Cv();
 }
diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H
index 868f74040..f66193394 100644
--- a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H
+++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H
@@ -4,9 +4,10 @@
     surfaceScalarField psisf = fvc::interpolate(psis);
     surfaceScalarField rhof = fvc::interpolate(rho);
 
-    // Needs to be outside of loop since p is changing, but psi and rho are not.
+    // Needs to be outside of loop since p is changing, but psi and rho are not
     surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p);
 
+    // --- PISO loop
     for (int corr = 0; corr < nCorr; corr++)
     {
         U = rUA*UEqn.H();