Bugfix: regular formulation of enthalpy source terms

applications/solvers/compressible/steadyCompressibleMRFFoam/createFields.H

@@ -93,4 +93,4 @@
         ),
         h
     );
-    i -= 0.5*magSqr(Urot);
+    i == h - 0.5*(magSqr(Urot) - magSqr(Urel));

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
     (

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();
 }

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();