sonicFoam: change internal energy equation formulation to conserve total energy

applications/solvers/compressible/sonicFoam/UEqn.H

@@ -6,3 +6,4 @@ fvVectorMatrix UEqn
 );
 
 solve(UEqn == -fvc::grad(p));
+K = 0.5*magSqr(U);

applications/solvers/compressible/sonicFoam/createFields.H

@@ -49,3 +49,6 @@
             thermo
         )
     );
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));

applications/solvers/compressible/sonicFoam/eEqn.H

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

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -39,3 +39,4 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+K = 0.5*magSqr(U);

applications/solvers/compressible/sonicFoam/sonicFoam.C

@@ -56,16 +56,12 @@ int main(int argc, char *argv[])
         #include "compressibleCourantNo.H"
 
         #include "rhoEqn.H"
-
         #include "UEqn.H"
 
-        #include "eEqn.H"
-
-
         // --- PISO loop
-
         for (int corr=0; corr<nCorr; corr++)
         {
+            #include "eEqn.H"
             #include "pEqn.H"
         }