Updated steadyCompressible*Foam solvers & steadyCompressibleFoam tutorials

Conflicts:
	applications/solvers/compressible/steadyCompressibleFoam/pEqn.H

applications/solvers/compressible/steadyCompressibleFoam/createFields.H

@@ -9,7 +9,8 @@
     volScalarField& p = thermo.p();
     volScalarField& h = thermo.h();
     const volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo.psi();
+    volScalarField psis("psi", thermo.psi()/thermo.Cp()*thermo.Cv());
+    psis.oldTime();
 
     volScalarField rho
     (

applications/solvers/compressible/steadyCompressibleFoam/hEqn.H

@@ -29,4 +29,5 @@
 
     // Bounding of enthalpy taken out
     thermo.correct();
+    psis = thermo.psi()/thermo.Cp()*thermo.Cv();
 }

applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H

@@ -11,13 +11,17 @@
     {
         U = rUA*UEqn.H();
 
-        phi = rhoReff*(fvc::interpolate(U) & mesh.Sf());
+        // Calculate phi for boundary conditions
+        phi = rhof*fvc::interpolate(U) & mesh.Sf();
 
-        surfaceScalarField phid("phid", psisf/rhoReff*phi);
+        surfaceScalarField phid2 = rhoReff/rhof*phi;
+
+        surfaceScalarField phid("phid", psisf/rhof*phi);
 
         // Make fluxes relative within the MRF zone
-        mrfZones.relativeFlux(psisf, phid);
         mrfZones.relativeFlux(rhoReff, phi);
+        mrfZones.relativeFlux(psisf, phid);
+        mrfZones.relativeFlux(rhoReff, phid2);
 
         p.storePrevIter();
 
@@ -27,7 +31,7 @@
             (
                 fvm::ddt(psis, p)
               + fvm::div(phid, p)
-              + fvc::div(phi)
+              + fvc::div(phid2)
               - fvm::laplacian(rho*rUA, p)
             );
 
@@ -42,13 +46,13 @@
             // Calculate the flux
             if (nonOrth == nNonOrthCorr)
             {
-                phi += pEqn.flux();
+                phi = phid2 + pEqn.flux();
             }
         }
 
 #       include "compressibleContinuityErrs.H"
 
-        // Explicitly relax the pressure for momentum corrector
+        // Relax the pressure
         p.relax();
 
         U -= rUA*fvc::grad(p);

applications/solvers/compressible/steadyCompressibleSRFFoam/createFields.H

@@ -9,7 +9,8 @@
     volScalarField& p = thermo.p();
     volScalarField& h = thermo.h();
     const volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo.psi();
+    volScalarField psis("psi", thermo.psi()/thermo.Cp()*thermo.Cv());
+    psis.oldTime();
 
     volScalarField rho
     (

applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H

@@ -1,5 +1,5 @@
 {
-    // Solve the enthalpy equation
+    // Solve the rothalpy equation
     T.storePrevIter();
 
     // Calculate face velocity from flux
@@ -38,4 +38,5 @@
 
     // Bounding of enthalpy taken out
     thermo.correct();
+    psis = thermo.psi()/thermo.Cp()*thermo.Cv();
 }

applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H

@@ -1,24 +1,22 @@
 {
     volScalarField rUrelA = 1.0/UrelEqn.A();
 
+    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.
+    surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p);
+
     for (int corr = 0; corr < nCorr; corr++)
     {
         Urel = rUrelA*UrelEqn.H();
 
-        // Execute ddtPhiCorr before recalculating flux
-        // HJ, 27/Apr/2010
-        surfaceScalarField phid
-        (
-            "phid",
-            fvc::interpolate(thermo.psi())*
-            (
-                (fvc::interpolate(Urel) & mesh.Sf())
-              + fvc::ddtPhiCorr(rUrelA, rho, Urel, phi)
-            )
-        );
-
         // Calculate phi for boundary conditions
-        phi = fvc::interpolate(rho*Urel) & mesh.Sf();
+        phi = rhof*fvc::interpolate(Urel) & mesh.Sf();
+
+        surfaceScalarField phid2 = rhoReff/rhof*phi;
+
+        surfaceScalarField phid("phid", psisf/rhof*phi);
 
         p.storePrevIter();
 
@@ -26,8 +24,9 @@
         {
             fvScalarMatrix pEqn
             (
-                fvm::ddt(psi, p)
+                fvm::ddt(psis, p)
               + fvm::div(phid, p)
+              + fvc::div(phid2)
               - fvm::laplacian(rho*rUrelA, p)
             );
 
@@ -42,14 +41,13 @@
             // Calculate the flux
             if (nonOrth == nNonOrthCorr)
             {
-                phi = pEqn.flux();
+                phi = phid2 + pEqn.flux();
             }
         }
 
 #       include "compressibleContinuityErrs.H"
 
-
+        // Relax the pressure
-        // Explicitly relax the pressure for momentum corrector
         p.relax();
 
         Urel -= rUrelA*fvc::grad(p);

tutorials/compressible/steadyCompressibleFoam/2bump/system/fvSolution

@@ -58,11 +58,9 @@ PIMPLE
 relaxationFactors
 {
     // Note: under-relaxation factors used in wave-transmissive schemes
-    U               0.5;
+    U               0.7;
-    p               0.2;
+    p               0.7;
-    h               0.5;
+    h               0.7;
-    rho             0.5;
-    T               1;
 }
 
 fieldBounds

tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSchemes

@@ -20,7 +20,7 @@ ddtSchemes
 
     ddt(rho,U)      steadyState;
     ddt(rho,h)      steadyState;
-    ddt(psi,p)      steadyInertial phi rho 0.5;
+    ddt(psi,p)      steadyInertial phi rho 0.95;
 
     U               steadyState;
     T               steadyState;

tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution

@@ -76,8 +76,6 @@ relaxationFactors
     U               0.1;
     p               0.25;
     h               0.1;
-    rho             0.5;
-    T               0.5;
 }
 
 fieldBounds