diff --git a/applications/solvers/combustion/PDRFoam/UEqn.H b/applications/solvers/combustion/PDRFoam/UEqn.H
index 70912c7cfae08052593d297d0bdf6d7e19f396d6..79d81663c4f24d5d849e8e07438fe9d626042fc3 100644
--- a/applications/solvers/combustion/PDRFoam/UEqn.H
+++ b/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);
     }
diff --git a/applications/solvers/combustion/PDRFoam/createFields.H b/applications/solvers/combustion/PDRFoam/createFields.H
index 15a0ecafea7dd4bca89b0f84d13e876ed5995e49..2d8451c9ae665d5a429a9187cbdbb16c94c5622a 100644
--- a/applications/solvers/combustion/PDRFoam/createFields.H
+++ b/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;
diff --git a/applications/solvers/combustion/PDRFoam/hEqn.H b/applications/solvers/combustion/PDRFoam/hEqn.H
index 7f5292d01a903f7931aa0015e81cb26c9250253a..4d8bb54588ae0fe2348d28862ce1aead46f18d7c 100644
--- a/applications/solvers/combustion/PDRFoam/hEqn.H
+++ b/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();
diff --git a/applications/solvers/combustion/PDRFoam/huEqn.H b/applications/solvers/combustion/PDRFoam/huEqn.H
index 3467bc6b751177800f70b927250e2b61926ff04b..f0fa7be5eabadca152c27d6ed24037d5d9a400ed 100644
--- a/applications/solvers/combustion/PDRFoam/huEqn.H
+++ b/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()
     );
 }
diff --git a/applications/solvers/combustion/PDRFoam/pEqn.H b/applications/solvers/combustion/PDRFoam/pEqn.H
index e2a2a471e28e1935c990af0c7f1f1f2959dc5d41..45f982ba4d7fd2c575ab517439c96f1fe62b562f 100644
--- a/applications/solvers/combustion/PDRFoam/pEqn.H
+++ b/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);
diff --git a/applications/solvers/combustion/XiFoam/UEqn.H b/applications/solvers/combustion/XiFoam/UEqn.H
index 1e626d75b85f487ec98f31131f0c0370f8d639c1..b9bc567aae9107cd8d687ec4666488586472f863 100644
--- a/applications/solvers/combustion/XiFoam/UEqn.H
+++ b/applications/solvers/combustion/XiFoam/UEqn.H
@@ -12,4 +12,5 @@
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+        K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/combustion/XiFoam/createFields.H b/applications/solvers/combustion/XiFoam/createFields.H
index 7a7e796b62c1df021b6d19b2439babc15c985bdb..1750c693e5cb986984a29a9a750689b99c39b9c4 100644
--- a/applications/solvers/combustion/XiFoam/createFields.H
+++ b/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
diff --git a/applications/solvers/combustion/XiFoam/hEqn.H b/applications/solvers/combustion/XiFoam/hEqn.H
index 513ae604419f28b96bbf455e8ed0bc1e75d258bb..f562466b200fd6fc4ba2f9b3bba63b2f8686887e 100644
--- a/applications/solvers/combustion/XiFoam/hEqn.H
+++ b/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();
diff --git a/applications/solvers/combustion/XiFoam/huEqn.H b/applications/solvers/combustion/XiFoam/huEqn.H
index 0b4068344bc8b74a46691ba7ed205c6cd4efc9ba..3059b9445e60dbec71fd7381a39d48b485edeb89 100644
--- a/applications/solvers/combustion/XiFoam/huEqn.H
+++ b/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()
     );
 }
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index 784e9ca6b71b90d114e13d880d7bf18b067c66f6..cb25e7029f6fe9e2b1a4b3863a416e5f23adfd62 100644
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/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);
diff --git a/applications/solvers/combustion/coldEngineFoam/createFields.H b/applications/solvers/combustion/coldEngineFoam/createFields.H
index 9705f04525e4221a3561c51c2f42b5cbef2b55b0..6286b1a5553fea4b0a5cd79f20ac3e0572e90adc 100644
--- a/applications/solvers/combustion/coldEngineFoam/createFields.H
+++ b/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));
diff --git a/applications/solvers/combustion/coldEngineFoam/hEqn.H b/applications/solvers/combustion/coldEngineFoam/hEqn.H
index ae60d3316ec77061804e629360ed13f6cd891f68..2e1a2703a8eee6fb85a3a2cfa173cdfd26a72e2b 100644
--- a/applications/solvers/combustion/coldEngineFoam/hEqn.H
+++ b/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();
diff --git a/applications/solvers/combustion/engineFoam/UEqn.H b/applications/solvers/combustion/engineFoam/UEqn.H
index eff9e1a5d288a5ffde71f0d800bba9dc80095784..f110051946816e9a0f3aab7160c3681af6317377 100644
--- a/applications/solvers/combustion/engineFoam/UEqn.H
+++ b/applications/solvers/combustion/engineFoam/UEqn.H
@@ -8,4 +8,5 @@
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+        K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index 580be2da56ef087f2c63ed7b30470234e4e40eeb..d1b7135649a270428d8280bf0d8e9a4a766fddcb 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/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);
diff --git a/applications/solvers/combustion/fireFoam/UEqn.H b/applications/solvers/combustion/fireFoam/UEqn.H
index a64e50a2d24d19b3f1dd25e4bc73acbcfc27a356..31d01e4327c7749a5f9b30ccb1a814bc5f5e8d7c 100644
--- a/applications/solvers/combustion/fireFoam/UEqn.H
+++ b/applications/solvers/combustion/fireFoam/UEqn.H
@@ -23,4 +23,6 @@
                 )*mesh.magSf()
             )
         );
+
+        K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/combustion/fireFoam/YhsEqn.H b/applications/solvers/combustion/fireFoam/YhsEqn.H
index 299eb64ab7669e20893c3aa614d34e7a66af2653..b85396caa171cb21fc55efaad74359adb80c05a1 100644
--- a/applications/solvers/combustion/fireFoam/YhsEqn.H
+++ b/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)
diff --git a/applications/solvers/combustion/fireFoam/createFields.H b/applications/solvers/combustion/fireFoam/createFields.H
index afabd393c74464a69451c860f700aeef2222e759..e0c11335e3f691c3ad1ffe2250622372725e1e1f 100644
--- a/applications/solvers/combustion/fireFoam/createFields.H
+++ b/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;
diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index 65259391daf83589f2c939ebcf829f1c9a50fda4..872a0513a3d87e6311723dded401bbff4e10e998 100644
--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/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);
diff --git a/applications/solvers/combustion/reactingFoam/UEqn.H b/applications/solvers/combustion/reactingFoam/UEqn.H
index 1e626d75b85f487ec98f31131f0c0370f8d639c1..b9bc567aae9107cd8d687ec4666488586472f863 100644
--- a/applications/solvers/combustion/reactingFoam/UEqn.H
+++ b/applications/solvers/combustion/reactingFoam/UEqn.H
@@ -12,4 +12,5 @@
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+        K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/combustion/reactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/createFields.H
index 7cf94f754d9a205246efd838f58564f236aa5c0f..0b6e9ee2e1df3308148b832d222348f60206e531 100644
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/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;
 
diff --git a/applications/solvers/combustion/reactingFoam/hsEqn.H b/applications/solvers/combustion/reactingFoam/hsEqn.H
index de1a85fddf5365c2be897aaeff7ffa8a42cd35cc..642b7091bce54db6a81f5871fc2aebd5f5c2af9e 100644
--- a/applications/solvers/combustion/reactingFoam/hsEqn.H
+++ b/applications/solvers/combustion/reactingFoam/hsEqn.H
@@ -6,7 +6,8 @@
       - fvm::laplacian(turbulence->alphaEff(), hs)
 //      - fvm::laplacian(turbulence->muEff(), hs)  // unit lewis no.
      ==
-        DpDt
+        dpdt
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
       + combustion->Sh()
     );
 
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index 784e9ca6b71b90d114e13d880d7bf18b067c66f6..cb25e7029f6fe9e2b1a4b3863a416e5f23adfd62 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/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);
diff --git a/applications/solvers/combustion/rhoReactingFoam/UEqn.H b/applications/solvers/combustion/rhoReactingFoam/UEqn.H
index 1e626d75b85f487ec98f31131f0c0370f8d639c1..b9bc567aae9107cd8d687ec4666488586472f863 100644
--- a/applications/solvers/combustion/rhoReactingFoam/UEqn.H
+++ b/applications/solvers/combustion/rhoReactingFoam/UEqn.H
@@ -12,4 +12,5 @@
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+        K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/combustion/rhoReactingFoam/createFields.H b/applications/solvers/combustion/rhoReactingFoam/createFields.H
index 65860a67eb26605f0501062cdbfbebf9da5949f8..690cdc4b6054bde13a5e0904c8604f894ea43f05 100644
--- a/applications/solvers/combustion/rhoReactingFoam/createFields.H
+++ b/applications/solvers/combustion/rhoReactingFoam/createFields.H
@@ -67,11 +67,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;
 
diff --git a/applications/solvers/combustion/rhoReactingFoam/hsEqn.H b/applications/solvers/combustion/rhoReactingFoam/hsEqn.H
index 01f85ac1da7432acf7a54a555dc8283af58979d3..ba95be845274b234f80da59e698851297724e4be 100644
--- a/applications/solvers/combustion/rhoReactingFoam/hsEqn.H
+++ b/applications/solvers/combustion/rhoReactingFoam/hsEqn.H
@@ -6,7 +6,8 @@
       - fvm::laplacian(turbulence->alphaEff(), hs)
 //    - fvm::laplacian(turbulence->muEff(), hs)  // unit lewis no.
      ==
-        DpDt
+        dpdt
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
       + combustion->Sh()
     );
 
diff --git a/applications/solvers/combustion/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
index b22fc890794a469f5fef975b198c8592fa3af963..4e94de243c6429a38740558e91144a5659d524f8 100644
--- a/applications/solvers/combustion/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
@@ -86,6 +86,7 @@
 
     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);
 }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
index e6c74938ee00d4cc77c232de32e537576c33c9f3..8ff402a2b1aa414ad9583bd6925138f0b298a01a 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
@@ -5,8 +5,7 @@
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
-      - (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
+      - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
     );
 
     hEqn.relax();
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H
index 6e607f89e87ce91aeef4a66c7a9d1ad8d71cce78..24ed135c08160d6d3d9cfa8bfb8f4396919bb35a 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H
@@ -5,8 +5,7 @@
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
-      - p*fvc::div(phi/fvc::interpolate(rho))
+      - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
     );
 
     hEqn.relax();
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H
index 57de44e4a0b39ecfe5c17b719df318866e3bbdab..2c422a5cab526e69b5f38b2450413c24390f7d50 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H
@@ -5,8 +5,7 @@
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
-      - p*fvc::div(phi/fvc::interpolate(rho))
+      - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
       + radiation->Sh(thermo)
     );
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H
index 03b258042529f3e90c6ca942fdd946f094d78214..4505f77ed0aaf0f0321a028108793d900d30f4c6 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H
@@ -5,8 +5,7 @@
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turb.alphaEff(), h)
      ==
-        fvc::div(phi/fvc::interpolate(rho), p, "div(U,p)")
-      - p*fvc::div(phi/fvc::interpolate(rho))
+      - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
       + rad.Sh(thermo)
     );
 
diff --git a/tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/fvSchemes b/tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/fvSchemes
index c68d0f8a2f547dc62665f4bb862ccf8e6b70ce0b..7de6b2228e3b9b1cc180cdabc902af973af11e64 100644
--- a/tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/fvSchemes
+++ b/tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/fvSchemes
@@ -31,7 +31,7 @@ divSchemes
     default         none;
 
     div(phi,U)      Gauss limitedLinearV 1;
-    div(phiU,p)     Gauss limitedLinear 1;
+    div(phi,K)      Gauss limitedLinear 1;
     div(phid,p)     Gauss limitedLinear 1;
     div(phi,k)      Gauss limitedLinear 1;
     div(phi,epsilon) Gauss limitedLinear 1;
diff --git a/tutorials/combustion/XiFoam/les/pitzDaily/system/fvSchemes b/tutorials/combustion/XiFoam/les/pitzDaily/system/fvSchemes
index 1aff419daf0198d111d34f424dc3f05a629f8618..303b193a018ea01e4a5157446431e609277091d3 100644
--- a/tutorials/combustion/XiFoam/les/pitzDaily/system/fvSchemes
+++ b/tutorials/combustion/XiFoam/les/pitzDaily/system/fvSchemes
@@ -29,7 +29,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss linear;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,k)      Gauss limitedLinear 0.1;
     div(phiXi,Xi)   Gauss limitedLinear01 0.1;
     div(phiXi,Su)   Gauss limitedLinear01 0.1;
diff --git a/tutorials/combustion/XiFoam/les/pitzDaily3D/system/fvSchemes b/tutorials/combustion/XiFoam/les/pitzDaily3D/system/fvSchemes
index 1aff419daf0198d111d34f424dc3f05a629f8618..303b193a018ea01e4a5157446431e609277091d3 100644
--- a/tutorials/combustion/XiFoam/les/pitzDaily3D/system/fvSchemes
+++ b/tutorials/combustion/XiFoam/les/pitzDaily3D/system/fvSchemes
@@ -29,7 +29,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss linear;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,k)      Gauss limitedLinear 0.1;
     div(phiXi,Xi)   Gauss limitedLinear01 0.1;
     div(phiXi,Su)   Gauss limitedLinear01 0.1;
diff --git a/tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSchemes b/tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSchemes
index 059d89339cb116f47afdc7cb85881d21924ae085..39ac52720671675285c768704aa12ffb4b173853 100644
--- a/tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSchemes
+++ b/tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSchemes
@@ -30,7 +30,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss limitedLinearV 1;
     div(phid,p)     Gauss limitedLinear 1;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss limitedLinear 1;
     div(phi,k)      Gauss limitedLinear 1;
     div(phi,epsilon) Gauss limitedLinear 1;
     div(phi,R)      Gauss limitedLinear 1;
diff --git a/tutorials/combustion/engineFoam/kivaTest/system/fvSchemes b/tutorials/combustion/engineFoam/kivaTest/system/fvSchemes
index ffb8a769e5591152e774f4e5439147e0950a11b3..62b296f3ed42a7239133a9641726d128a9dea836 100644
--- a/tutorials/combustion/engineFoam/kivaTest/system/fvSchemes
+++ b/tutorials/combustion/engineFoam/kivaTest/system/fvSchemes
@@ -29,8 +29,8 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
+    div(phi,K)      Gauss upwind;
     div(phid,p)     Gauss upwind;
-    div(phiU,p)     Gauss linear;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes
index 8fec68c55f476eef6c2c64b66118b578a2a32908..03db056e8b26e67fb249df01a89d3715cba994bf 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes
@@ -40,7 +40,7 @@ divSchemes
         hs              limitedLinear 1;
     };
     div((muEff*dev2(T(grad(U))))) Gauss linear;
-    div(phiU,p)         Gauss linear;
+    div(phi,K)          Gauss limitedLinear 1;
     div(Ji,Ii_h)        Gauss upwind;
 }
 
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSchemes b/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSchemes
index 6a10661e2400d2bdf384758d48b3d7ff3b82cc41..775f4cb396ae2457033199179cfb61fd5bbdd16b 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSchemes
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSchemes
@@ -41,7 +41,7 @@ divSchemes
     };
     div((muEff*dev2(T(grad(U))))) Gauss linear;
     div(phi,omega)  Gauss limitedLinear 1;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss limitedLinear 1;
     div(U)          Gauss linear;
     div(Ji,Ii_h)    Gauss upwind;
 }
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/G b/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/G
index c0b1c57ad9bdb629fb091cc842723200b9c701d6..fe75242d76939285cab48e4fc610b174e845687a 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/G
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/G
@@ -24,7 +24,8 @@ boundaryField
     {
         type            MarshakRadiation;
         T               T;
-        emissivity      1;
+        emissivityMode  lookup;
+        emissivity      uniform 1.0;
         value           uniform 0;
     }
 }
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/polyMesh/boundary b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/polyMesh/boundary
index 54beb4a79dbb458432d8117589fd83a510953f2e..5f99cb1862ebf281c53bac07be9cfea3822f01f4 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/polyMesh/boundary
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/polyMesh/boundary
@@ -15,25 +15,31 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-3
+4
 (
     base
     {
         type            patch;
-        nFaces          3600;
+        nFaces          3456;
         startFace       637200;
     }
     outlet
     {
         type            patch;
         nFaces          3600;
-        startFace       640800;
+        startFace       640656;
     }
     sides
     {
         type            patch;
         nFaces          14400;
-        startFace       644400;
+        startFace       644256;
+    }
+    inlet
+    {
+        type            patch;
+        nFaces          144;
+        startFace       658656;
     }
 )
 
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSchemes b/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSchemes
index c7a3c9b507a71d312cb66743ecee61665cb92abb..0a76845c6f5327f30c0567fd50ebf8ea8bb4179b 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSchemes
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSchemes
@@ -29,6 +29,7 @@ divSchemes
 {
     default        none;
     div(phi,U)      Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,k)      Gauss limitedLinear 0.1;
     div(phi,Yi_hs) Gauss multivariateSelection
     {
@@ -40,7 +41,6 @@ divSchemes
         hs              limitedLinear 1;
     };
     div((muEff*dev2(T(grad(U))))) Gauss linear;
-    div(phiU,p)     Gauss linear;
     div(Ji,Ii_h)    Gauss upwind;
 }
 
diff --git a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/system/fvSchemes b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/system/fvSchemes
index f9d7d6ecf00d0447cfd70c4c6d72ca8ae8a128fa..6d4a329ff2af85dc7ad157ae8ef6b066599a0767 100644
--- a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/system/fvSchemes
+++ b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
     div(phi,U)      Gauss limitedLinearV 1;
     div(phi,Yi_h)   Gauss limitedLinear01 1;
     div(phi,h)      Gauss limitedLinear 1;
-    div(phiU,p)     Gauss limitedLinear 1;
+    div(phi,K)      Gauss limitedLinear 1;
     div(phid,p)     Gauss limitedLinear 1;
     div(phi,epsilon) Gauss limitedLinear 1;
     div(phi,k) Gauss limitedLinear 1;
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSchemes b/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSchemes
index 872b5d50fe35f9c139ec0c1694b825c1e83526d1..c9c95b180c2eb1e4f92092bd77559607b35226c4 100644
--- a/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSchemes
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSchemes
@@ -29,14 +29,14 @@ divSchemes
 {
     default             none;
 
-    div(phi,U)          Gauss upwind; //limitedLinearV 1; //upwind;
+    div(phi,U)          Gauss upwind;
     div((muEff*dev2(T(grad(U)))))      Gauss linear;
-    div(phi,h)          Gauss upwind; //limitedLinear 1; //upwind;
+    div(phi,h)          Gauss upwind;
     div(phi,epsilon)    Gauss upwind;
     div(phi,k)          Gauss upwind;
 
     div(phid,p)         Gauss upwind;
-    div(U,p)            Gauss linear;
+    div(phi,K)          Gauss upwind;
 }
 
 laplacianSchemes
@@ -46,24 +46,18 @@ laplacianSchemes
 
 interpolationSchemes
 {
-    default                 none;
-    interpolate(rho)        linear;
-    div(U,p)                upwind phi;
-    interpolate((psi*U))    linear;
-    interpolate(U)          linear;
-    UD                      upwind phid;
-    interpolate(p)          linear;
-    interpolate(((rho|(A(U)-H(1)))-(rho|A(U)))) linear;
+    default         linear;
+    UD              upwind phid;
 }
 
 snGradSchemes
 {
-    default                 corrected;
+    default         corrected;
 }
 
 fluxRequired
 {
-    default                 no;
+    default         no;
     p;
     pCorr;
 }
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution b/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution
index 8b5f9ccd5ff8c138a1d5a1251807e82e0483a2db..df2018d930bd1479894f25acee3e227938bcebfc 100644
--- a/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution
@@ -17,14 +17,6 @@ FoamFile
 
 solvers
 {
-    p0
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-08;
-        relTol          0.01;
-    }
-
     p
     {
         solver          GAMG;
@@ -40,120 +32,7 @@ solvers
         mergeLevels     1;
     }
 
-    U0
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-08;
-        relTol          0.1;
-    }
-
-    U1
-    {
-        solver          smoothSolver;
-        smoother        GaussSeidel;
-        nSweeps         1;
-        tolerance       1e-08;
-        relTol          0.1;
-    }
-
-    U
-    {
-        solver          GAMG;
-        tolerance       1e-08;
-        relTol          0.1;
-        smoother        GaussSeidel;
-        nPreSweeps      0;
-        nPostSweeps     2;
-        nFinestSweeps   2;
-        cacheAgglomeration true;
-        nCellsInCoarsestLevel 20;
-        agglomerator    faceAreaPair;
-        mergeLevels     1;
-    }
-
-    h0
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-08;
-        relTol          0.1;
-    }
-
-    h1
-    {
-        solver          smoothSolver;
-        smoother        GaussSeidel;
-        nSweeps         1;
-        tolerance       1e-08;
-        relTol          0.1;
-    }
-
-    h
-    {
-        solver          GAMG;
-        tolerance       1e-08;
-        relTol          0.1;
-        smoother        GaussSeidel;
-        nPreSweeps      0;
-        nPostSweeps     2;
-        nFinestSweeps   2;
-        cacheAgglomeration true;
-        nCellsInCoarsestLevel 20;
-        agglomerator    faceAreaPair;
-        mergeLevels     1;
-    }
-
-    k0
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-08;
-        relTol          0.1;
-    }
-
-    k1
-    {
-        solver          smoothSolver;
-        smoother        GaussSeidel;
-        nSweeps         1;
-        tolerance       1e-08;
-        relTol          0.1;
-    }
-
-    k
-    {
-        solver          GAMG;
-        tolerance       1e-08;
-        relTol          0.1;
-        smoother        GaussSeidel;
-        nPreSweeps      0;
-        nPostSweeps     2;
-        nFinestSweeps   2;
-        cacheAgglomeration true;
-        nCellsInCoarsestLevel 20;
-        agglomerator    faceAreaPair;
-        mergeLevels     1;
-    }
-
-    epsilon0
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-08;
-        relTol          0.1;
-    }
-
-    epsilon1
-    {
-        solver          smoothSolver;
-        smoother        GaussSeidel;
-        nSweeps         1;
-        tolerance       1e-08;
-        relTol          0.1;
-    }
-
-    epsilon
+    "(U|h|k|epsilon)"
     {
         solver          GAMG;
         tolerance       1e-08;
@@ -182,31 +61,15 @@ relaxationFactors
     fields
     {
         p               1;
-        rho             1; //0.1;
+        rho             1;
     }
     equations
     {
         U               0.9;
-        h               0.95;
+        h               0.8;
         k               0.9;
         epsilon         0.9;
     }
 }
 
-relaxationFactors0
-{
-    fields
-    {
-        p               0.3;
-        rho             0.1;
-    }
-    equations
-    {
-        U               0.7;
-        h               0.7;
-        k               0.7;
-        epsilon         0.7;
-    }
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/fvSchemes
index d19197862d6add401dff443b4af0733dd1519b74..02edc9816b94775af20c1716ecad09bf6ab7f8b9 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/fvSchemes
@@ -28,6 +28,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss limitedLinear 0.2;
+    div(phi,K)      Gauss limitedLinear 0.2;
     div(phi,h)      Gauss limitedLinear 0.2;
     div(phi,k)      Gauss limitedLinear 0.2;
     div(phi,epsilon) Gauss limitedLinear 0.2;
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/fvSchemes
index d19197862d6add401dff443b4af0733dd1519b74..02edc9816b94775af20c1716ecad09bf6ab7f8b9 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/fvSchemes
@@ -28,6 +28,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss limitedLinear 0.2;
+    div(phi,K)      Gauss limitedLinear 0.2;
     div(phi,h)      Gauss limitedLinear 0.2;
     div(phi,k)      Gauss limitedLinear 0.2;
     div(phi,epsilon) Gauss limitedLinear 0.2;
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes
index 593744fbe739b75781abd103b56da3704953ba48..f0019f66909410c050c47880c88b649476777274 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phi,h)      Gauss upwind;
+    div(phi,K)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes
index f35af85372788bc3119d1713d5c4d0e186c3adb0..572a6207d2280e9135103318008b254328051db3 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes
@@ -29,6 +29,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
+    div(phi,K)      Gauss upwind;
     div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes
index ad314dec6803cf961da2d20cafab2a50c3e3bec5..3cbd11840882a2c5ceeca7760d7932ee1f771620 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes
@@ -29,6 +29,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
+    div(phi,K)      Gauss upwind;
     div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;