Thermodynamics: Changed all eEqn to EEqn and reformulated to conserve E in sonic solvers

To support these changes the need for "Sp" corrections on div-terms has been
eliminated by introducing a "bounded" convection scheme which subtracts the Sp
term from the selected scheme.  The equivalent will be needed for the ddt term.

A warning message is generated for steady-state solvers in which the "bounded"
scheme is not selected for the convection terms.

applications/solvers/compressible/rhoSimpleFoam/EEqn.H

+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::div(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+    );
+
+    EEqn.relax();
+    EEqn.solve();
+
+    thermo.correct();
+}

applications/solvers/compressible/rhoSimpleFoam/createFields.H

@@ -5,6 +5,7 @@
         psiThermo::New(mesh)
     );
     psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
 
     volScalarField rho
     (
@@ -20,7 +21,6 @@
     );
 
     volScalarField& p = thermo.p();
-    volScalarField& e = thermo.he();
     const volScalarField& psi = thermo.psi();
 
     Info<< "Reading field U\n" << endl;

applications/solvers/compressible/rhoSimpleFoam/eEqn.H

-{
-    // Kinetic + pressure energy
-    volScalarField Ekp("Ekp", 0.5*magSqr(U) + p/rho);
-
-    fvScalarMatrix eEqn
-    (
-        fvm::div(phi, e)
-      - fvm::Sp(fvc::div(phi), e)
-      - fvm::laplacian(turbulence->alphaEff(), e)
-     ==
-        fvc::div(phi)*Ekp - fvc::div(phi, Ekp)
-    );
-
-    eEqn.relax();
-    eEqn.solve();
-
-    thermo.correct();
-}

applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H

+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::div(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+    );
+
+    pZones.addEnergySource(thermo, rho, EEqn);
+
+    EEqn.relax();
+    EEqn.solve();
+
+    thermo.correct();
+}

applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H

@@ -5,6 +5,7 @@
         rhoThermo::New(mesh)
     );
     rhoThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
 
     volScalarField rho
     (
@@ -20,7 +21,6 @@
     );
 
     volScalarField& p = thermo.p();
-    volScalarField& e = thermo.he();
 
     Info<< "Reading field U\n" << endl;
     volVectorField U

applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H

-{
-    // Kinetic + pressure energy
-    volScalarField Ekp("Ekp", 0.5*magSqr(U) + p/rho);
-
-    fvScalarMatrix eEqn
-    (
-        fvm::div(phi, e)
-      - fvm::Sp(fvc::div(phi), e)
-      - fvm::laplacian(turbulence->alphaEff(), e)
-     ==
-        fvc::div(phi)*Ekp - fvc::div(phi, Ekp)
-    );
-
-    pZones.addEnergySource(thermo, rho, eEqn);
-
-    eEqn.relax();
-    eEqn.solve();
-
-    thermo.correct();
-}

applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C

@@ -63,7 +63,7 @@ int main(int argc, char *argv[])
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
-            #include "eEqn.H"
+            #include "EEqn.H"
             #include "pEqn.H"
         }
 

applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C

@@ -59,7 +59,7 @@ int main(int argc, char *argv[])
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
-            #include "eEqn.H"
+            #include "EEqn.H"
             #include "pEqn.H"
         }
 

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C

@@ -61,7 +61,7 @@ int main(int argc, char *argv[])
         // Velocity-pressure-enthalpy SIMPLEC corrector
         {
             #include "UEqn.H"
-            #include "eEqn.H"
+            #include "EEqn.H"
             #include "pEqn.H"
         }
 

applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H → applications/solvers/compressible/sonicFoam/EEqn.H

 {
     solve
     (
-        fvm::ddt(rho, e)
-      + fvm::div(phi, e)
+        fvm::ddt(rho, e) + fvm::div(phi, e)
+      + fvc::ddt(rho, K) + fvc::div(phi, volScalarField("Ekp", K + p/rho))
       - fvm::laplacian(turbulence->alphaEff(), e)
-     ==
-      - p*fvc::div(phi/fvc::interpolate(rho) + mesh.phi())
     );
 
     thermo.correct();

applications/solvers/compressible/sonicFoam/createFields.H

@@ -5,6 +5,7 @@
         psiThermo::New(mesh)
     );
     psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "e");
 
     volScalarField& p = thermo.p();
     volScalarField& e = thermo.he();

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -17,7 +17,8 @@ surfaceScalarField phid
 
 volScalarField Dp("Dp", rho*rAU);
 
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+// Non-orthogonal pressure corrector loop
+while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix pEqn
     (
@@ -28,7 +29,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
     pEqn.solve();
 
-    if (nonOrth == nNonOrthCorr)
+    if (pimple.finalNonOrthogonalIter())
     {
         phi = pEqn.flux();
     }

applications/solvers/compressible/sonicFoam/eEqn.H → applications/solvers/compressible/sonicFoam/sonicDyMFoam/EEqn.H

 {
     solve
     (
-        fvm::ddt(rho, e)
-      + fvm::div(phi, e)
+        fvm::ddt(rho, e) + fvm::div(phi, e)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + fvc::div(phi/fvc::interpolate(rho) + mesh.phi(), p, "div(phiv,p)")
       - fvm::laplacian(turbulence->alphaEff(), e)
-     ==
-      - (fvc::ddt(rho, K) + fvc::div(phi, volScalarField("Ekp", K + p/rho)))
     );
 
     thermo.correct();

applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C

@@ -34,6 +34,7 @@ Description
 #include "psiThermo.H"
 #include "turbulenceModel.H"
 #include "motionSolver.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +46,8 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -62,20 +65,24 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        #include "UEqn.H"
-
-        #include "eEqn.H"
-
-
-        // --- PISO loop
-
-        for (int corr=0; corr<nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
-            #include "pEqn.H"
+            #include "UEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
-        turbulence->correct();
-
         rho = thermo.rho();
 
         runTime.write();

applications/solvers/compressible/sonicFoam/sonicFoam.C

@@ -33,6 +33,7 @@ Description
 #include "fvCFD.H"
 #include "psiThermo.H"
 #include "turbulenceModel.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,8 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -52,21 +55,29 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
+        #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
 
         #include "rhoEqn.H"
-        #include "UEqn.H"
 
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
-            #include "eEqn.H"
-            #include "pEqn.H"
+            #include "UEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
-        turbulence->correct();
-
         rho = thermo.rho();
 
         runTime.write();

applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C

@@ -31,6 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,8 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -52,57 +55,59 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
+        #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
 
         #include "rhoEqn.H"
 
-        fvVectorMatrix UEqn
-        (
-            fvm::ddt(rho, U)
-          + fvm::div(phi, U)
-          - fvm::laplacian(mu, U)
-        );
-
-        solve(UEqn == -fvc::grad(p));
-
-
-        // --- PISO loop
-
-        for (int corr=0; corr<nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
-            volScalarField rAU(1.0/UEqn.A());
-            U = rAU*UEqn.H();
-
-            surfaceScalarField phid
+            fvVectorMatrix UEqn
             (
-                "phid",
-                psi
-               *(
-                    (fvc::interpolate(U) & mesh.Sf())
-                  + fvc::ddtPhiCorr(rAU, rho, U, phi)
-                )
+                fvm::ddt(rho, U)
+              + fvm::div(phi, U)
+              - fvm::laplacian(mu, U)
             );
 
-            phi = (rhoO/psi)*phid;
-            volScalarField Dp("Dp", rho*rAU);
+            solve(UEqn == -fvc::grad(p));
 
-            fvScalarMatrix pEqn
-            (
-                fvm::ddt(psi, p)
-              + fvc::div(phi)
-              + fvm::div(phid, p)
-              - fvm::laplacian(Dp, p)
-            );
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                volScalarField rAU(1.0/UEqn.A());
+                U = rAU*UEqn.H();
+
+                surfaceScalarField phid
+                (
+                    "phid",
+                    psi
+                   *(
+                       (fvc::interpolate(U) & mesh.Sf())
+                     + fvc::ddtPhiCorr(rAU, rho, U, phi)
+                   )
+                );
+
+                phi = (rhoO/psi)*phid;
+                volScalarField Dp("Dp", rho*rAU);
+
+                fvScalarMatrix pEqn
+                (
+                    fvm::ddt(psi, p)
+                  + fvc::div(phi)
+                  + fvm::div(phid, p)
+                  - fvm::laplacian(Dp, p)
+                );
 
-            pEqn.solve();
+                pEqn.solve();
 
-            phi += pEqn.flux();
+                phi += pEqn.flux();
 
-            #include "compressibleContinuityErrs.H"
+                #include "compressibleContinuityErrs.H"
 
-            U -= rAU*fvc::grad(p);
-            U.correctBoundaryConditions();
+                U -= rAU*fvc::grad(p);
+                U.correctBoundaryConditions();
+            }
         }
 
         rho = rhoO + psi*p;

applications/solvers/heatTransfer/buoyantPimpleFoam/EEqn.H

+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + fvm::div(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ep", p/rho))
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+    );
+
+    EEqn.relax();
+    EEqn.solve();
+
+    thermo.correct();
+}

applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C

@@ -75,7 +75,7 @@ int main(int argc, char *argv[])
         while (pimple.loop())
         {
             #include "UEqn.H"
-            #include "hEqn.H"
+            #include "EEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())

applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H

@@ -5,6 +5,7 @@
         rhoThermo::New(mesh)
     );
     rhoThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
 
     volScalarField rho
     (
@@ -20,7 +21,6 @@
     );
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.he();
     const volScalarField& psi = thermo.psi();
 
 

applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H

-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + fvm::div(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        dpdt
-      - (fvc::ddt(rho, K) + fvc::div(phi, K))
-    );
-
-    hEqn.relax();
-    hEqn.solve();
-
-    thermo.correct();
-}