Thermodynamics: Changed h-eqn to conserve total energy

applications/solvers/combustion/PDRFoam/UEqn.H

@@ -13,4 +13,5 @@
     {
         U = invA & (UEqn.H() - betav*fvc::grad(p));
         U.correctBoundaryConditions();
+        K = 0.5*magSqr(U);
     }

applications/solvers/combustion/PDRFoam/createFields.H

@@ -44,7 +44,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::RASModel> turbulence
@@ -58,12 +58,12 @@
         )
     );
 
-    Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt
-    (
-        "DpDt",
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
-    );
+
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt("dpdt", fvc::ddt(p));
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));
 
 
     Info<< "Creating the unstrained laminar flame speed\n" << endl;

applications/solvers/combustion/PDRFoam/hEqn.H

@@ -5,7 +5,8 @@
       + mvConvection->fvmDiv(phi, h)
       - fvm::laplacian(Db, h)
      ==
-        betav*DpDt
+        betav*dpdt
+      - betav*(fvc::ddt(rho, K) + fvc::div(phi, K))
     );
 
     thermo.correct();

applications/solvers/combustion/PDRFoam/huEqn.H

@@ -13,6 +13,6 @@ if (ign.ignited())
     //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), hu)
 
      ==
-        betav*DpDt*rho/thermo.rhou()
+        betav*(dpdt - (fvc::ddt(rho, K) + fvc::div(phi, K)))*rho/thermo.rhou()
     );
 }

applications/solvers/combustion/PDRFoam/pEqn.H

@@ -64,5 +64,6 @@ else
 
 U -= invA & (betav*fvc::grad(p));
 U.correctBoundaryConditions();
+K = 0.5*magSqr(U);
 
-DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+dpdt = fvc::ddt(p);

applications/solvers/combustion/XiFoam/UEqn.H

@@ -12,4 +12,5 @@
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+        K = 0.5*magSqr(U);
     }

applications/solvers/combustion/XiFoam/createFields.H

@@ -60,12 +60,11 @@
         )
     );
 
-    Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt
-    (
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
-    );
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt("dpdt", fvc::ddt(p));
 
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));
 
     Info<< "Creating field Xi\n" << endl;
     volScalarField Xi

applications/solvers/combustion/XiFoam/hEqn.H

@@ -5,7 +5,8 @@
       + mvConvection->fvmDiv(phi, h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        DpDt
+        dpdt
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
     );
 
     hEqn.relax();

applications/solvers/combustion/XiFoam/huEqn.H

@@ -12,7 +12,7 @@ if (ign.ignited())
     //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), hu)
     //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), hu)
 
-     ==
-        DpDt*rho/thermo.rhou()
+    ==
+        (dpdt - (fvc::ddt(rho, K) + fvc::div(phi, K)))*rho/thermo.rhou()
     );
 }

applications/solvers/combustion/XiFoam/pEqn.H

@@ -64,5 +64,6 @@ else
 
 U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+K = 0.5*magSqr(U);
 
-DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+dpdt = fvc::ddt(p);

applications/solvers/combustion/coldEngineFoam/createFields.H

@@ -54,8 +54,8 @@
         )
     );
 
-    Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt
-    (
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
-    );
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt("dpdt", fvc::ddt(p));
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));

applications/solvers/combustion/coldEngineFoam/hEqn.H

@@ -5,7 +5,7 @@
       + fvm::div(phi, h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        DpDt
+      - fvc::div(phi, 0.5*magSqr(U))
     );
 
     thermo.correct();

applications/solvers/combustion/engineFoam/UEqn.H

@@ -8,4 +8,5 @@
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+        K = 0.5*magSqr(U);
     }

applications/solvers/combustion/engineFoam/pEqn.H

@@ -57,5 +57,6 @@ else
 
 U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+K = 0.5*magSqr(U);
 
-DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+dpdt = fvc::ddt(p);

applications/solvers/combustion/fireFoam/UEqn.H

@@ -23,4 +23,6 @@
                 )*mesh.magSf()
             )
         );
+
+        K = 0.5*magSqr(U);
     }

applications/solvers/combustion/fireFoam/YhsEqn.H

@@ -53,7 +53,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
       + mvConvection->fvmDiv(phi, hs)
       - fvm::laplacian(turbulence->alphaEff(), hs)
      ==
-        DpDt
+        dpdt
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
       + combustion->Sh()
       + radiation->Shs(thermo)
       + parcels.Sh(hs)

applications/solvers/combustion/fireFoam/createFields.H

@@ -83,12 +83,12 @@
         dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );
 
-    Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt
-    (
-        "DpDt",
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
-    );
+
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt("dpdt", fvc::ddt(p));
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));
 
 
     Info<< "Calculating field g.h\n" << endl;

applications/solvers/combustion/fireFoam/pEqn.H

@@ -43,5 +43,6 @@ p = p_rgh + rho*gh;
 
 U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
 U.correctBoundaryConditions();
+K = 0.5*magSqr(U);
 
-DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+dpdt = fvc::ddt(p);

applications/solvers/combustion/reactingFoam/UEqn.H

@@ -12,4 +12,5 @@
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+        K = 0.5*magSqr(U);
     }

applications/solvers/combustion/reactingFoam/createFields.H

@@ -65,11 +65,13 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 
-Info<< "Creating field DpDt\n" << endl;
-volScalarField DpDt
-(
-    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
-);
+
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt("dpdt", fvc::ddt(p));
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;