diff --git a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
index 4688633cb37109e35085ea3045380e88a346319f..1e39e983e42a476d62c73a65a0ef531771f43991 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
@@ -14,4 +14,5 @@ volScalarField rAU(1.0/UEqn().A());
 if (pimple.momentumPredictor())
 {
     solve(UEqn() == -fvc::grad(p));
+    K = 0.5*magSqr(U);
 }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index 070f651a999296d7214d5f86e7fec2e021df2b58..e2c7b950de7ba519f7eee4eec39b89667d4ef8fd 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/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/compressible/rhoPimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPimpleFoam/hEqn.H
index 3125cc3ffa86ce120e7dbbf774c9b46941105418..b1d9e4e8b0237ca0b9f28f6d36c8f5162121c846 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/hEqn.H
@@ -5,7 +5,8 @@
       + fvm::div(phi, h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        DpDt
+        dpdt
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
     );
 
     hEqn.relax();
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 15aa4b9300fb6663531a259ba76ba80a8a525526..f7bb932fdc28349536810212f41932b3436f917f 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -83,5 +83,6 @@ Info<< "rho max/min : " << max(rho).value()
 
 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/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
index cb7ecaa29d7652134db213612e523a8803987ae8..831f139e817ee9f882b786a05a6a4ad0a67e088e 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
@@ -85,5 +85,6 @@ Info<< "rho max/min : " << max(rho).value()
 
 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/rhoPorousMRFSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/hEqn.H
index ff0b91bcb7aa71baafec4132306043b07b6ffd04..cb989b804f3a01533fabd085f4c1bf0e05b74b4a 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/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)")
     );
 
     pZones.addEnthalpySource(thermo, rho, hEqn);
diff --git a/applications/solvers/compressible/sonicFoam/eEqn.H b/applications/solvers/compressible/sonicFoam/eEqn.H
index 1d1d6aef0b2d53efa849347af2a0c4d6f9cbf912..6733f528854a92ebba609a2769535f4b28fe36a3 100644
--- a/applications/solvers/compressible/sonicFoam/eEqn.H
+++ b/applications/solvers/compressible/sonicFoam/eEqn.H
@@ -4,7 +4,7 @@
         fvm::ddt(rho, e)
       + fvm::div(phi, e)
       - fvm::laplacian(turbulence->alphaEff(), e)
-      ==
+     ==
       - p*fvc::div(phi/fvc::interpolate(rho))
     );
 
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
index 7dff741e663b150ccb55b177ed5b725dcdc2a0a2..9597f6ac08ddcac42bba70642fe61c013f1c87ed 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
@@ -23,4 +23,5 @@
                 )*mesh.magSf()
             )
         );
+        K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
index e39d0bb80333b1cb16731cb0d63356c100aaf093..1c2a2b94e49f1d1bc7f14bff9fce5d0684d6f727 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@@ -74,9 +74,8 @@
     // Force p_rgh to be consistent with p
     p_rgh = p - rho*gh;
 
-    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));
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H
index 3125cc3ffa86ce120e7dbbf774c9b46941105418..b1d9e4e8b0237ca0b9f28f6d36c8f5162121c846 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H
@@ -5,7 +5,8 @@
       + fvm::div(phi, h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        DpDt
+        dpdt
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
     );
 
     hEqn.relax();
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index 5c60d78b35acfc20d7ef5e7372494db4d10d9d90..397c1c8bf488d90f794fa8eb2b6084c17e2b173a 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -47,6 +47,7 @@
             // calculated from the relaxed pressure
             U += rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
+            K = 0.5*magSqr(U);
         }
     }
 
@@ -55,7 +56,7 @@
     // Second part of thermodynamic density update
     thermo.rho() += psi*p_rgh;
 
-    DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+    dpdt = fvc::ddt(p);
 
     #include "rhoEqn.H"
     #include "compressibleContinuityErrs.H"
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
index 3d8531e1bb6bb98ac532724d71f3a62337950744..14f51f4c507b6eecb9f6e802766432aa1ed8f235 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
@@ -1,7 +1,7 @@
     // Initialise fluid field pointer lists
     PtrList<basicRhoThermo> thermoFluid(fluidRegions.size());
     PtrList<volScalarField> rhoFluid(fluidRegions.size());
-    PtrList<volScalarField> KFluid(fluidRegions.size());
+    PtrList<volScalarField> kappaFluid(fluidRegions.size());
     PtrList<volVectorField> UFluid(fluidRegions.size());
     PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
     PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
@@ -50,15 +50,15 @@
             )
         );
 
-        Info<< "    Adding to KFluid\n" << endl;
-        KFluid.set
+        Info<< "    Adding to kappaFluid\n" << endl;
+        kappaFluid.set
         (
             i,
             new volScalarField
             (
                 IOobject
                 (
-                    "K",
+                    "kappa",
                     runTime.timeName(),
                     fluidRegions[i],
                     IOobject::NO_READ,
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index 862581ef790689a49475e6254d39de6f238e8437..73b72af316a7c3bf3cddc9eb05af4c35bc326d26 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -71,5 +71,5 @@
         << max(rho).value() << endl;
 
     // Update thermal conductivity
-    K = thermo.Cp()*turb.alphaEff();
+    kappa = thermo.Cp()*turb.alphaEff();
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
index 1b6a22db0da022cf756921210b2c59d3abcbb00c..cbad65e1c2bc3918515866f7cf0e7bc1fb3e6487 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
@@ -2,7 +2,7 @@
 
     basicRhoThermo& thermo = thermoFluid[i];
     volScalarField& rho = rhoFluid[i];
-    volScalarField& K = KFluid[i];
+    volScalarField& kappa = kappaFluid[i];
     volVectorField& U = UFluid[i];
     surfaceScalarField& phi = phiFluid[i];
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/setRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/setRegionSolidFields.H
index 179480249d8853e0c831e173fd7d7407c1d28ef7..7ff98e6211ea01aac48e362c43228da061dda68f 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/setRegionSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/setRegionSolidFields.H
@@ -7,8 +7,8 @@
     tmp<volScalarField> tcp = thermo.Cp();
     const volScalarField& cp = tcp();
 
-    tmp<volScalarField> tK = thermo.K();
-    //tmp<volSymmTensorField> tK = thermo.directionalK();
-    const volScalarField& K = tK();
+    tmp<volScalarField> tkappa = thermo.K();
+    //tmp<volSymmTensorField> tkappa = thermo.directionalkappa();
+    const volScalarField& kappa = tkappa();
 
     volScalarField& T = thermo.T();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H
index 4e2a2c7b19b57d4e17a9ea3dd4e769333f72115f..0df54bea01e371bba3ed74edc91eb6473269debc 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H
@@ -3,7 +3,7 @@
     {
         fvScalarMatrix tEqn
         (
-            -fvm::laplacian(K, T)
+            -fvm::laplacian(kappa, T)
         );
         tEqn.relax();
         tEqn.solve();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
index 119eeb4fd7d314bbd353934519093edb327e39c4..69a23fa1d2f7dcf543110315c0ad5a862c205382 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
@@ -22,5 +22,6 @@
                 )*mesh.magSf()
             ),
             mesh.solver(U.select(finalIter))
-        );
+         );
+         K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index f6788d26eb3ecccb2bea1484949a781e7f3010a3..7090319d9402eb2ff0583e69a7b6cd6850896b36 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -1,7 +1,7 @@
     // Initialise fluid field pointer lists
     PtrList<basicRhoThermo> thermoFluid(fluidRegions.size());
     PtrList<volScalarField> rhoFluid(fluidRegions.size());
-    PtrList<volScalarField> KFluid(fluidRegions.size());
+    PtrList<volScalarField> kappaFluid(fluidRegions.size());
     PtrList<volVectorField> UFluid(fluidRegions.size());
     PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
     PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
@@ -10,7 +10,8 @@
     PtrList<volScalarField> ghFluid(fluidRegions.size());
     PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
     PtrList<radiation::radiationModel> radiation(fluidRegions.size());
-    PtrList<volScalarField> DpDtFluid(fluidRegions.size());
+    PtrList<volScalarField> KFluid(fluidRegions.size());
+    PtrList<volScalarField> dpdtFluid(fluidRegions.size());
 
     List<scalar> initialMassFluid(fluidRegions.size());
 
@@ -45,15 +46,15 @@
             )
         );
 
-        Info<< "    Adding to KFluid\n" << endl;
-        KFluid.set
+        Info<< "    Adding to kappaFluid\n" << endl;
+        kappaFluid.set
         (
             i,
             new volScalarField
             (
                 IOobject
                 (
-                    "K",
+                    "kappa",
                     runTime.timeName(),
                     fluidRegions[i],
                     IOobject::NO_READ,
@@ -175,22 +176,25 @@
 
         initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
 
-        Info<< "    Adding to DpDtFluid\n" << endl;
-        DpDtFluid.set
+        Info<< "    Adding to KFluid\n" << endl;
+        KFluid.set
         (
             i,
             new volScalarField
             (
-                "DpDt",
-                fvc::DDt
-                (
-                    surfaceScalarField
-                    (
-                        "phiU",
-                        phiFluid[i]/fvc::interpolate(rhoFluid[i])
-                    ),
-                    thermoFluid[i].p()
-                )
+                "K",
+                0.5*magSqr(UFluid[i])
+            )
+        );
+
+        Info<< "    Adding to dpdtFluid\n" << endl;
+        dpdtFluid.set
+        (
+            i,
+            new volScalarField
+            (
+                "dpdt",
+                fvc::ddt(thermoFluid[i].p())
             )
         );
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
index 0bb3c8b47c61f70aedd822071aab0a6753957055..293a2001c0e47f5c8ada847eb4eadf850b3c73ec 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
@@ -5,7 +5,8 @@
       + fvm::div(phi, h)
       - fvm::laplacian(turb.alphaEff(), h)
      ==
-        DpDt
+        dpdt
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
       + rad.Sh(thermo)
     );
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index 07ab9c92e81799f1f80653cb414c095829c0e927..1df38d15bb63c704da4cdd51ce6f4734191e087b 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -68,11 +68,13 @@
     // Correct velocity field
     U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
     U.correctBoundaryConditions();
+    K = 0.5*magSqr(U);
 
     p = p_rgh + rho*gh;
 
-    // Update pressure substantive derivative
-    DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+    // Update pressure time derivative
+    dpdt = fvc::ddt(p);
+
 
     // Solve continuity
     #include "rhoEqn.H"
@@ -91,5 +93,5 @@
     }
 
     // Update thermal conductivity
-    K = thermoFluid[i].Cp()*turb.alphaEff();
+    kappa = thermoFluid[i].Cp()*turb.alphaEff();
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 81c6d25bb0ccc68597a0b154e0b7c3baa98552ba..6529c1393144bfbf2b444c14c459d0187a3970a6 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -2,12 +2,13 @@
 
     basicRhoThermo& thermo = thermoFluid[i];
     volScalarField& rho = rhoFluid[i];
-    volScalarField& K = KFluid[i];
+    volScalarField& kappa = kappaFluid[i];
     volVectorField& U = UFluid[i];
     surfaceScalarField& phi = phiFluid[i];
 
     compressible::turbulenceModel& turb = turbulence[i];
-    volScalarField& DpDt = DpDtFluid[i];
+    volScalarField& K = KFluid[i];
+    volScalarField& dpdt = dpdtFluid[i];
 
     volScalarField& p = thermo.p();
     const volScalarField& psi = thermo.psi();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
index a843ed8bd75f562bbba0dfb157356b36aaaaf5fd..c013b733db95dab9f2f4d87a3b15bf878f7c170d 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
@@ -7,9 +7,9 @@
     tmp<volScalarField> tcp = thermo.Cp();
     const volScalarField& cp = tcp();
 
-    tmp<volScalarField> tK = thermo.K();
-    //tmp<volSymmTensorField> tK = thermo.directionalK();
-    const volScalarField& K = tK();
-    //const volSymmTensorField& K = tK();
+    tmp<volScalarField> tkappa = thermo.K();
+    //tmp<volSymmTensorField> tkappa = thermo.directionalkappa();
+    const volScalarField& kappa = tkappa();
+    //const volSymmTensorField& kappa = tkappa();
 
     volScalarField& T = thermo.T();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.C
index 3dbab74858368248c4abe6ec2d507f011e2ac0ff..437ea78327462d81682bd59c630399111835002f 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.C
@@ -31,7 +31,7 @@ Foam::scalar Foam::solidRegionDiffNo
     const fvMesh& mesh,
     const Time& runTime,
     const volScalarField& Cprho,
-    const volScalarField& K
+    const volScalarField& kappa
 )
 {
     scalar DiNum = 0.0;
@@ -39,16 +39,16 @@ Foam::scalar Foam::solidRegionDiffNo
 
     //- Take care: can have fluid domains with 0 cells so do not test for
     //  zero internal faces.
-    surfaceScalarField KrhoCpbyDelta
+    surfaceScalarField kapparhoCpbyDelta
     (
         mesh.surfaceInterpolation::deltaCoeffs()
-      * fvc::interpolate(K)
+      * fvc::interpolate(kappa)
       / fvc::interpolate(Cprho)
     );
 
-    DiNum = gMax(KrhoCpbyDelta.internalField())*runTime.deltaT().value();
+    DiNum = gMax(kapparhoCpbyDelta.internalField())*runTime.deltaT().value();
 
-    meanDiNum = (average(KrhoCpbyDelta)).value()*runTime.deltaT().value();
+    meanDiNum = (average(kapparhoCpbyDelta)).value()*runTime.deltaT().value();
 
     Info<< "Region: " << mesh.name() << " Diffusion Number mean: " << meanDiNum
         << " max: " << DiNum << endl;
@@ -62,26 +62,26 @@ Foam::scalar Foam::solidRegionDiffNo
     const fvMesh& mesh,
     const Time& runTime,
     const volScalarField& Cprho,
-    const volSymmTensorField& Kdirectional
+    const volSymmTensorField& kappadirectional
 )
 {
     scalar DiNum = 0.0;
     scalar meanDiNum = 0.0;
 
-    volScalarField K(mag(Kdirectional));
+    volScalarField kappa(mag(kappadirectional));
 
     //- Take care: can have fluid domains with 0 cells so do not test for
     //  zero internal faces.
-    surfaceScalarField KrhoCpbyDelta
+    surfaceScalarField kapparhoCpbyDelta
     (
         mesh.surfaceInterpolation::deltaCoeffs()
-      * fvc::interpolate(K)
+      * fvc::interpolate(kappa)
       / fvc::interpolate(Cprho)
     );
 
-    DiNum = gMax(KrhoCpbyDelta.internalField())*runTime.deltaT().value();
+    DiNum = gMax(kapparhoCpbyDelta.internalField())*runTime.deltaT().value();
 
-    meanDiNum = (average(KrhoCpbyDelta)).value()*runTime.deltaT().value();
+    meanDiNum = (average(kapparhoCpbyDelta)).value()*runTime.deltaT().value();
 
     Info<< "Region: " << mesh.name() << " Diffusion Number mean: " << meanDiNum
         << " max: " << DiNum << endl;
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.H
index 2fe2efa95e10628cc90c9af1c108a84bae337626..cb4dd30b6a3801b7674abfe4a17c27e516a819e6 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.H
@@ -39,7 +39,7 @@ namespace Foam
         const fvMesh& mesh,
         const Time& runTime,
         const volScalarField& Cprho,
-        const volScalarField& K
+        const volScalarField& kappa
     );
 
     scalar solidRegionDiffNo
@@ -47,7 +47,7 @@ namespace Foam
         const fvMesh& mesh,
         const Time& runTime,
         const volScalarField& Cprho,
-        const volSymmTensorField& K
+        const volSymmTensorField& kappa
     );
 
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
index 77dc6f04bf44ba748359c2f8bf3267307416b650..878780baf1df939beceb7efde5ff738d5da93d84 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
@@ -2,7 +2,7 @@
 
     forAll(solidRegions, i)
     {
-#       include "setRegionSolidFields.H"
+        #include "setRegionSolidFields.H"
 
         DiNum = max
         (
@@ -11,7 +11,7 @@
                 solidRegions[i],
                 runTime,
                 rho*cp,
-                K
+                kappa
             ),
             DiNum
         );
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
index d8aa03283b413d632634933e396f09ac4ab64e9b..d7781c4cd9e061568d236bf1133756774975a1c8 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
@@ -9,7 +9,7 @@ if (finalIter)
         tmp<fvScalarMatrix> TEqn
         (
             fvm::ddt(rho*cp, T)
-          - fvm::laplacian(K, T)
+          - fvm::laplacian(kappa, T)
         );
         TEqn().relax();
         TEqn().solve(mesh.solver(T.select(finalIter)));
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H
index f1f6173cd611633a7d7115f4f476f1dba583c2e1..264756c8331d46b3626d84bb8085811ef36da7c6 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H
@@ -124,20 +124,3 @@
         dimensionedScalar("one", dimless/dimTime, 1),
         zeroGradientFvPatchScalarField::typeName
     );
-
-    Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt
-    (
-        IOobject
-        (
-            "DpDt",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimPressure/dimTime, 0.0)
-    );
-
-    #include "setPressureWork.H"
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H
index 151e8e60bbbf98afdfc40487ab8799494c492e27..e8a4da764519b0eec8870a46cdc3d8f96525ba8d 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H
@@ -5,7 +5,7 @@
       + mvConvection->fvmDiv(phi, hs)
       - fvm::laplacian(turbulence->alphaEff(), hs)
      ==
-        DpDt
+      - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
       + parcels.Sh(hs)
       + radiation->Shs(thermo)
       + combustion->Sh()
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
index 27eb76e42fdf29a47ba4a68853ecc4e57b0ea9f6..141faa30bf2ea0fe7e738720fa5977a22d338bb2 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
@@ -65,7 +65,5 @@
     rho = max(rho, rhoMin);
     rho = min(rho, rhoMax);
 
-    #include "setPressureWork.H"
-
     Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
 }
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/setPressureWork.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/setPressureWork.H
deleted file mode 100644
index c5dd4b745152e5ba46b6c2db65ee98b3c861e0b8..0000000000000000000000000000000000000000
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/setPressureWork.H
+++ /dev/null
@@ -1,13 +0,0 @@
-DpDt == dimensionedScalar("zero", DpDt.dimensions(), 0.0);
-
-if (pressureWork)
-{
-    surfaceScalarField phiU("phiU", phi/fvc::interpolate(rho));
-
-    DpDt += fvc::div(phiU*fvc::interpolate(p)) - p*fvc::div(phiU);
-
-    if (pressureWorkTimeDerivative)
-    {
-        DpDt += fvc::ddt(p);
-    }
-}
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H
index 43ebced834180d302f19109d4f138155f5d552e0..24277a9f45acd53d2105b8d803b265d082e4c582 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H
@@ -74,10 +74,10 @@ Info<< "Time scales min/max:" << endl;
             runTime.deltaTValue()
            *mag
             (
-                DpDt
-              + parcels.hsTrans()/(mesh.V()*runTime.deltaT())
-//              + sources(rho, hs)
+                parcels.hsTrans()/(mesh.V()*runTime.deltaT())
+           // + sources(rho, hs)
               + combustion->Sh()()
+              - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")()
             )
            /rho
         );
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
index b7321872c34700aa13cc8f848aee24e58ec3453a..fe7b047d9a24e41d4baa7f54324f75e3cc08f160 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
@@ -17,4 +17,5 @@
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+        K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
index 13e8d94e72431a3242bdecacf0482d7f6f608cd9..8c786ae251beeeab4fd7b50fbc96d0c14d3a44da 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
@@ -115,12 +115,11 @@
     // Set the turbulence into the combustion model
     combustion->setTurbulence(turbulence());
 
-    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<< "\nConstructing sources" << endl;
     IObasicSourceList sources(mesh);
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
index 6df66134c432a442faa7c758a8873c863508d9b4..6a28b325600595dcd42f3962cae54fe48da59ebe 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
@@ -5,7 +5,8 @@
       + mvConvection->fvmDiv(phi, hs)
       - fvm::laplacian(turbulence->alphaEff(), hs)
      ==
-        DpDt
+        dpdt
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
       + combustion->Sh()
       + coalParcels.Sh(hs)
       + limestoneParcels.Sh(hs)
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index 539521ca9b48ea9eca5dc7bd285e97d54a5a2443..56435d31ce12c0c67ef9b8730b49022048c824ef 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -74,5 +74,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/lagrangian/reactingParcelFilmFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H
index a64e50a2d24d19b3f1dd25e4bc73acbcfc27a356..81c587e2438ce31da3b78672923bcf326c7a3a19 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H
@@ -23,4 +23,5 @@
                 )*mesh.magSf()
             )
         );
+        K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
index 999515f4309aa2ca8526ee39a3daddaf8ed9d71b..345d795a1822fab0b35f2a3cc896c7f45943a3f5 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
@@ -69,13 +69,11 @@
     // Set the turbulence into the combustion model
     combustion->setTurbulence(turbulence());
 
-    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;
     volScalarField gh("gh", g & mesh.C());
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H
index 79c5e1c5b3cc5a239437601f3bd8e3b57c6224e3..a6124813a913d88e00a137575db4e7e777db84fc 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H
@@ -5,7 +5,8 @@
       + mvConvection->fvmDiv(phi, hs)
       - fvm::laplacian(turbulence->alphaEff(), hs)
      ==
-        DpDt
+        dpdt
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
       + parcels.Sh(hs)
       + surfaceFilm.Sh()
       + radiation->Shs(thermo)
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
index 65259391daf83589f2c939ebcf829f1c9a50fda4..872a0513a3d87e6311723dded401bbff4e10e998 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/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/lagrangian/reactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H
index 31e00d8f35a0fedd2ad576ab4c795e0951ef4c7e..eac49bfc331f66e87e4cf7f62c235d17172911ff 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H
@@ -13,4 +13,5 @@
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+        K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
index c318498d94c3efcd27e0c257e306795081ce99fb..6e764c0d6d5af3a793aee5d2541d7db8be6aa7fc 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
@@ -76,13 +76,11 @@
     // Set the turbulence into the combustion model
     combustion->setTurbulence(turbulence());
 
-    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));
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
     forAll(Y, i)
@@ -103,4 +101,4 @@
         ),
         mesh,
         dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
-    );
\ No newline at end of file
+    );
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H
index f7499813c7f63fcfdbe026d28cfc4485c3cff0df..12438380aa7e7fb97793202e6781674b61e45489 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H
@@ -5,10 +5,11 @@
       + mvConvection->fvmDiv(phi, hs)
       - fvm::laplacian(turbulence->alphaEff(), hs)
      ==
-        DpDt
-     +  parcels.Sh(hs)
-     +  radiation->Shs(thermo)
-     +  combustion->Sh()
+        dpdt
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
+      + parcels.Sh(hs)
+      + radiation->Shs(thermo)
+      + combustion->Sh()
     );
 
     hEqn.relax();
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index 8031af20d576e8ecbf8362316668f59adb649e86..5ee551a0ea5a50ba0fca611a4c4ece7843840c55 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -68,6 +68,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/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H
index 26015141d27b79c84bd2d5cb4706ccd79101dbc2..0af1e179462850f275779cacfe6f82bf85cb160a 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H
@@ -11,7 +11,7 @@
      ==
         heatTransferCoeff*T2/Cp1/rho1
       - fvm::Sp(heatTransferCoeff/Cp1/rho1, T1)
-      + alpha1*Dp1Dt/Cp1/rho1
+      + alpha1*(dpdt - (fvc::ddt(rho1, K1) + fvc::div(phi1, K1)))/Cp1/rho1
     );
 
     fvScalarMatrix T2Eqn
@@ -23,7 +23,7 @@
      ==
         heatTransferCoeff*T1/Cp2/rho2
       - fvm::Sp(heatTransferCoeff/Cp2/rho2, T2)
-      + alpha2*Dp2Dt/Cp2/rho2
+      + alpha2*(dpdt - (fvc::ddt(rho2, K2) + fvc::div(phi2, K2)))/Cp2/rho2
     );
 
     T1Eqn.relax();
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
index cd2bc8b6417a4405842f2c64f5db8b114f0d5423..b3eca666c269703747b35d39747aee0782961dbd 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
@@ -357,6 +357,10 @@
     volScalarField dgdt =
         pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001));
 
-    Info<< "Creating field DpDt\n" << endl;
-    volScalarField Dp1Dt(fvc::DDt(phi1, p));
-    volScalarField Dp2Dt(fvc::DDt(phi2, p));
+
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt("dpdt", fvc::ddt(p));
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K1("K1", 0.5*magSqr(U1));
+    volScalarField K2("K2", 0.5*magSqr(U2));
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
index d6af50988faa90a997668c0563df9dbf3caba1e1..6e2fc44b0837b76de74363b2b973959416919d7b 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
@@ -127,6 +127,8 @@
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
-    Dp1Dt = fvc::DDt(phi1, p);
-    Dp2Dt = fvc::DDt(phi2, p);
+    K1 = 0.5*magSqr(U1);
+    K2 = 0.5*magSqr(U1);
+
+    dpdt = fvc::ddt(p);
 }
diff --git a/tutorials/compressible/rhoPimpleFoam/les/pitzDaily/system/fvSchemes b/tutorials/compressible/rhoPimpleFoam/les/pitzDaily/system/fvSchemes
index a27a7cb00d86138ff2e11f80e3d9dc1d4c9f182e..56aebd7a94fc31d1ae3a7d20453629f08aac63c2 100644
--- a/tutorials/compressible/rhoPimpleFoam/les/pitzDaily/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/les/pitzDaily/system/fvSchemes
@@ -30,7 +30,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss filteredLinear2V 0.2 0;
     div(phi,h)      Gauss filteredLinear2 0.2 0;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,k)      Gauss limitedLinear 1;
     div(phi,B)      Gauss limitedLinear 1;
     div(phi,muTilda) Gauss limitedLinear 1;
diff --git a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSchemes b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSchemes
index f8198db5e4be0e5868b72b72a91e46c0424498c0..1772c2a3fd2de453b29a4d86681ed374d8c0286d 100644
--- a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSchemes
@@ -31,7 +31,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/compressible/rhoPimpleFoam/ras/cavity/system/fvSchemes b/tutorials/compressible/rhoPimpleFoam/ras/cavity/system/fvSchemes
index 459736fe584d053660ad0afba527bb532a791887..3049da2bed4e6013c657ea02af19f736cffba999 100644
--- a/tutorials/compressible/rhoPimpleFoam/ras/cavity/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/ras/cavity/system/fvSchemes
@@ -31,7 +31,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 linear;
     div(phi,h)      Gauss limitedLinear 1;
     div(phi,k)      Gauss limitedLinear 1;
     div(phi,epsilon) Gauss limitedLinear 1;
diff --git a/tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/system/fvSchemes b/tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/system/fvSchemes
index a5ed5fdc7548934a80607d75c7a7fae3287465f2..00981f61fcb6a70666b768c6177a089b00cfc7c5 100644
--- a/tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/system/fvSchemes
@@ -31,7 +31,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes
index 45e2412652bb70ae63e7b4e412e7e19ab53e723f..dbca5c1607ad223ed5c440bd179c7708125c49b0 100644
--- a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
     div(phi,k)      Gauss limitedLinear 1;
     div(phi,epsilon) Gauss limitedLinear 1;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
-    div(phiU,p) Gauss linear;
+    div(phi,K)      Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes
index 4f47f1ef0d858ce76b8cb94b19dc27353580e270..d237bd73617ff269aa6435fc68b5269a8bca0fec 100644
--- a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes
@@ -32,7 +32,7 @@ divSchemes
     div(phi,h)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,k)      Gauss upwind;
-    div(U,p)        Gauss upwind;
+    div(phi,K)      Gauss upwind;
 }
 
 laplacianSchemes
diff --git a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
index 0abd1608aba0dcb6aa66c9488133a3c4b51c7588..a664893343e7216491ca9913ee0171ee31980aba 100644
--- a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
+++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
@@ -8,7 +8,7 @@
 FoamFile
 {
     version     2.0;
-    format      ascii;
+    format      binary;
     class       polyBoundaryMesh;
     location    "constant/polyMesh";
     object      boundary;
diff --git a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/system/fvSchemes b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/system/fvSchemes
index 491ee75580653f5ec676a6ee095d0e31ca1d1240..711e4e3233a2e700ada115f8ffeb80db86b6950d 100644
--- a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/system/fvSchemes
@@ -32,7 +32,7 @@ divSchemes
     div(phi,h)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,k)      Gauss upwind;
-    div(U,p)        Gauss upwind;
+    div(phi,K)      Gauss upwind;
 }
 
 laplacianSchemes
diff --git a/tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSchemes b/tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSchemes
index 6c0fad8e3e2116feeaf6a627673eb2d5b639330a..b4cf8d2e38d798122b9909b0a8e91b78f1c87e23 100644
--- a/tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSchemes
+++ b/tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss limitedLinear 1;
     div(phi,e)      Gauss limitedLinear 1;
-    div(phiU,p)     Gauss limitedLinear 1;
+    div(phi,K)      Gauss limitedLinear 1;
     div((muEff*dev2(T(grad(U))))) Gauss linear 1;
 }
 
diff --git a/tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSchemes b/tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSchemes
index 6c0fad8e3e2116feeaf6a627673eb2d5b639330a..b4cf8d2e38d798122b9909b0a8e91b78f1c87e23 100644
--- a/tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSchemes
+++ b/tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss limitedLinear 1;
     div(phi,e)      Gauss limitedLinear 1;
-    div(phiU,p)     Gauss limitedLinear 1;
+    div(phi,K)      Gauss limitedLinear 1;
     div((muEff*dev2(T(grad(U))))) Gauss linear 1;
 }
 
diff --git a/tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSchemes b/tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSchemes
index a5fc7790daa6338fb9b5b0db793fd79ba6a26685..b12dc45dcf6d0ba6aaad7e83fbbbdcac2bffd7c9 100644
--- a/tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSchemes
+++ b/tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(phid,p)     Gauss limitedLinear 1;
-    div(phiU,p)     Gauss limitedLinear 1;
+    div(phi,K)      Gauss limitedLinear 1;
     div(phi,e)      Gauss limitedLinear 1;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
diff --git a/tutorials/compressible/sonicFoam/ras/prism/system/fvSchemes b/tutorials/compressible/sonicFoam/ras/prism/system/fvSchemes
index 11831249de6f36cde0898690b2978c37697924c0..e7b5a41a604f719f59b61610b77d2b2de68a5e4c 100644
--- a/tutorials/compressible/sonicFoam/ras/prism/system/fvSchemes
+++ b/tutorials/compressible/sonicFoam/ras/prism/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(phid,p)     Gauss limitedLinear 1;
-    div(phiU,p)     Gauss limitedLinear 1;
+    div(phi,K)      Gauss limitedLinear 1;
     div(phi,e)      Gauss limitedLinear 1;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes
index 5d814f89ed7bade3e906b6c3993cfbbb44b00fec..3a7b6154e54e9b0d16611a9b96e1fbb61d7c16e1 100644
--- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(R)          Gauss linear;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSchemes
index 211e990ac15d8747b0b030fa3c6f4b88f96a461b..e4be18c50c4cc9746d7d7568943cb8d899c6aa07 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSchemes
@@ -28,7 +28,7 @@ divSchemes
 {
     default        none;
     div(phi,U)      Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSchemes
index 7d87f15326e34e25d7051b755e0571ef57c64d22..865dfd9f993d8329406c98c469be9cf3e9d7695c 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSchemes
@@ -28,7 +28,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/bottomWater/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/bottomWater/fvSchemes
index 7d87f15326e34e25d7051b755e0571ef57c64d22..865dfd9f993d8329406c98c469be9cf3e9d7695c 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/bottomWater/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/bottomWater/fvSchemes
@@ -28,7 +28,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/topAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/topAir/fvSchemes
index 7d87f15326e34e25d7051b755e0571ef57c64d22..865dfd9f993d8329406c98c469be9cf3e9d7695c 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/topAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/topAir/fvSchemes
@@ -28,7 +28,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSchemes
index 7d87f15326e34e25d7051b755e0571ef57c64d22..865dfd9f993d8329406c98c469be9cf3e9d7695c 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSchemes
@@ -28,7 +28,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSchemes
index 7d87f15326e34e25d7051b755e0571ef57c64d22..865dfd9f993d8329406c98c469be9cf3e9d7695c 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSchemes
@@ -28,7 +28,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/system/fvSchemes b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/system/fvSchemes
index a5fc7790daa6338fb9b5b0db793fd79ba6a26685..b12dc45dcf6d0ba6aaad7e83fbbbdcac2bffd7c9 100644
--- a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/system/fvSchemes
+++ b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(phid,p)     Gauss limitedLinear 1;
-    div(phiU,p)     Gauss limitedLinear 1;
+    div(phi,K)      Gauss limitedLinear 1;
     div(phi,e)      Gauss limitedLinear 1;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
diff --git a/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvSchemes b/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvSchemes
index ab163d0856e6d306a27a6749561852108ca7c74f..1718307504e15d7634313c07da4ee3c0126bb3f8 100644
--- a/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvSchemes
+++ b/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phi,Yi_h)   Gauss upwind;
     div(phi,h)      Gauss upwind;
-    div(phiU,p)     Gauss upwind;
+    div(phi,K)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,k)      Gauss upwind;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
diff --git a/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvSolution b/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvSolution
index f653a66ba66821ed8fdab3f235ff2239ef321c95..e2d4fcfdeeaf346b9552e3baba6bc5e786762657 100644
--- a/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvSolution
+++ b/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvSolution
@@ -92,8 +92,6 @@ PIMPLE
 
 additional
 {
-    pressureWork    true;
-    pressureWorkTimeDerivative true;
     solveSpecies    true;
 }
 
diff --git a/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvSchemes b/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvSchemes
index 98dab9606b44d5b925f3b6cb622090e1c316928b..f0c610dcbf46eb29f0c8bc63b23f67e7587ef501 100644
--- a/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvSchemes
+++ b/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvSchemes
@@ -31,7 +31,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvSolution b/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvSolution
index eb09deb87a1631268cd31d0409bf960ed48101b0..82179b044e5f0078065cf042211e49eb3a01fa49 100644
--- a/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvSolution
+++ b/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvSolution
@@ -79,8 +79,6 @@ PIMPLE
 
 additional
 {
-    pressureWork    true;
-    pressureWorkTimeDerivative false; // true;
     solveSpecies    true;
 }
 
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSchemes b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSchemes
index 11223b58184cc90f40d46885c5e6f58e9c3eff2a..f2abd9112f30100f1b0157e750f21aaeb4ab075b 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSchemes
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSchemes
@@ -31,7 +31,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSchemes b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSchemes
index 11223b58184cc90f40d46885c5e6f58e9c3eff2a..f2abd9112f30100f1b0157e750f21aaeb4ab075b 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSchemes
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSchemes
@@ -31,7 +31,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes
index 68df90c908981dacd37fc017844a69b2098fec08..9b46cd8f1cc706b1e5439be1f3c1a031eafeb690 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes
@@ -31,7 +31,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes
index b65a4948da955a9834f6c3f5897d138669a5af16..7a96d901447dcddc0c0030f4353b902498605e6b 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/system/fvSchemes b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/system/fvSchemes
index 91117027586772e7972f0961630c0ddbd40b3da2..57fb0a2e123b417fe11b2e860de95b38eea1cc9d 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/system/fvSchemes
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/system/fvSchemes
@@ -30,7 +30,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
-    div(phiU,p)     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/lagrangian/reactingParcelFilmFoam/hotBoxes/system/fvSchemes b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/fvSchemes
index e7f90a15fbcb2685f5bb4d942162e21d85855f43..729ed886a08b2ef6ebab37497f8d32b322833d34 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/fvSchemes
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/fvSchemes
@@ -30,7 +30,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
-    div(phiU,p)     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/lagrangian/reactingParcelFilmFoam/rivuletPanel/system/fvSchemes b/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/system/fvSchemes
index e7f90a15fbcb2685f5bb4d942162e21d85855f43..729ed886a08b2ef6ebab37497f8d32b322833d34 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/system/fvSchemes
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/system/fvSchemes
@@ -30,7 +30,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
-    div(phiU,p)     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/lagrangian/reactingParcelFilmFoam/splashPanel/system/fvSchemes b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/system/fvSchemes
index e7f90a15fbcb2685f5bb4d942162e21d85855f43..729ed886a08b2ef6ebab37497f8d32b322833d34 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/system/fvSchemes
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/system/fvSchemes
@@ -30,7 +30,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
-    div(phiU,p)     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/lagrangian/reactingParcelFoam/evaporationTest/system/fvSchemes b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSchemes
index 11223b58184cc90f40d46885c5e6f58e9c3eff2a..f2abd9112f30100f1b0157e750f21aaeb4ab075b 100644
--- a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSchemes
+++ b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSchemes
@@ -31,7 +31,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/lagrangian/sprayFoam/aachenBomb/system/fvSchemes b/tutorials/lagrangian/sprayFoam/aachenBomb/system/fvSchemes
index 11223b58184cc90f40d46885c5e6f58e9c3eff2a..f2abd9112f30100f1b0157e750f21aaeb4ab075b 100644
--- a/tutorials/lagrangian/sprayFoam/aachenBomb/system/fvSchemes
+++ b/tutorials/lagrangian/sprayFoam/aachenBomb/system/fvSchemes
@@ -31,7 +31,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(phi,K)      Gauss linear;
     div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary
index 1bf12fb1ac720f9fb1d2d4ace54957b82decc7f5..56e0a545c15839416fcd3beca616ce0f609c5879 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary
@@ -1,4 +1,4 @@
-/*--------------------------------*- C++ -*----------------------------------* \
+/*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 |  \\    /   O peration     | Version:  dev                                   |
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
index ab8afe3d4d8c6030affe02d29997f2c6068a0db5..93b115565da213bde2ddf9929ba567fb6081ec6e 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
@@ -44,8 +44,8 @@ divSchemes
     div(phi,Theta)      Gauss limitedLinear 1;
     div(phid1,p)        Gauss upwind;
     div(phid2,p)        Gauss upwind;
-    div(phi1,p)         Gauss linear;
-    div(phi2,p)         Gauss linear;
+    div(phi1,K1)        Gauss limitedLinearV 1;
+    div(phi2,K2)        Gauss limitedLinearV 1;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSchemes b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSchemes
index 99b15f62aa3ec943aa0cf534e565881a02b84598..27ba67b86d244f6cea2d70740fc81486a5705f18 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSchemes
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSchemes
@@ -44,8 +44,8 @@ divSchemes
     div((alpha2*Rc2))        Gauss linear;
     div(phid1,p)        Gauss upwind;
     div(phid2,p)        Gauss upwind;
-    div(phi1,p)         Gauss linear;
-    div(phi2,p)         Gauss linear;
+    div(phi1,K1)        Gauss limitedLinearV 1;
+    div(phi2,K2)        Gauss limitedLinearV 1;
 }
 
 laplacianSchemes