diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
index 9bb0de54b1c1247593f74b30369fc9893aaddd49..28f6f81d9f0f2c1c8499d292d5d8473c0233bf84 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
@@ -10,6 +10,8 @@
       + sources(rho, U)
     );
 
+    UEqn.relax();
+
     sources.constrain(UEqn);
 
     pZones.addResistance(UEqn);
@@ -17,5 +19,6 @@
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+        K = 0.5*magSqr(U);
     }
 
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
index e333a14712a819d8b7828313924a3a3af0e779c8..60f6e1a55da8e3a958546f19f814fef0a558a510 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
@@ -35,6 +35,8 @@ if (solveSpecies)
               + sources(rho, Yi)
             );
 
+            YiEqn.relax();
+
             sources.constrain(YiEqn);
 
             YiEqn.solve(mesh.solver("Yi"));
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
index ab6b01f9cea220ce72b805249fae37dee8e68c03..a656d4ea8c524c5610748271c9af2300c3e04c77 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
@@ -98,3 +98,9 @@
         mesh,
         dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );
+
+    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/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
index 27a9f6c8b3d136db75fa0b9bdd033d8df09422d5..fc88ba785e4abfd0de285de86c70f93688711e47 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
@@ -1,56 +1,32 @@
 {
-    tmp<volScalarField> pWork
+    fvScalarMatrix hsEqn
     (
-        new volScalarField
-        (
-            IOobject
-            (
-                "pWork",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedScalar("zero", dimensionSet(1, -1, -3, 0, 0), 0.0)
-        )
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+      - (fvc::ddt(rho, K) + fvc::div(phi, K))
+      + parcels.Sh(hs)
+      + radiation->Shs(thermo)
+      + combustion->Sh()
+      + sources(rho, hs)
     );
 
-    if (pressureWork)
+    if (pressureWorkTimeDerivative)
     {
-        surfaceScalarField phiU("phiU", phi/fvc::interpolate(rho));
-
-        pWork() += fvc::div(phiU*fvc::interpolate(p)) - p*fvc::div(phiU);
-
-        if (pressureWorkTimeDerivative)
-        {
-            pWork() += fvc::ddt(p);
-        }
+        hsEqn -= dpdt;
     }
 
-    {
-        fvScalarMatrix hsEqn
-        (
-            fvm::ddt(rho, hs)
-          + mvConvection->fvmDiv(phi, hs)
-          - fvm::laplacian(turbulence->alphaEff(), hs)
-         ==
-            pWork()
-          + parcels.Sh(hs)
-          + radiation->Shs(thermo)
-          + combustion->Sh()
-          + sources(rho, hs)
-        );
+    hsEqn.relax();
 
-        sources.constrain(hsEqn);
+    sources.constrain(hsEqn);
 
-        hsEqn.solve();
+    hsEqn.solve();
 
-        thermo.correct();
+    thermo.correct();
 
-        radiation->correct();
+    radiation->correct();
 
-        Info<< "T gas min/max   = " << min(T).value() << ", "
-            << max(T).value() << endl;
-    }
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
 }
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
index 23645a1b7be45e76d6cb3fa04e9746d868c324d5..080faeccedd2c9cd2d4a85d762510d6ee796ea19 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
@@ -59,4 +59,7 @@
 
     U -= rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+    K = 0.5*magSqr(U);
+
+    dpdt = fvc::ddt(p);
 }
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H
index 8136af1ad3e9161f9a4769dfd99933fad95d1ab2..050c24fe7adf101cc7a235a5d8dc5c8c183caaf3 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H
@@ -1,9 +1,8 @@
 dictionary additional = mesh.solutionDict().subDict("additional");
 
 // pressure work term for enthalpy equation
-bool pressureWork = additional.lookupOrDefault("pressureWork", true);
 bool pressureWorkTimeDerivative =
-    additional.lookupOrDefault("pressureWorkTimeDerivative", true);
+    additional.lookupOrDefault("pressureWorkTimeDerivative", false);
 
 // flag to activate solve transport for each specie (Y vector)
 bool solveSpecies = additional.lookupOrDefault("solveSpecies", true);
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
index cb77f67e760ebc3b93dce7907228ce56d7af76f6..857dd7637cb2adbc0a5c18925c55104dadda266e 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
@@ -645,20 +645,18 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
                     phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
                 )
                *(
-                    fvc::interpolate(rA*rho.oldTime())
+                    fvc::interpolate(rA)
                    *(
                        coefft0*phiAbs.oldTime()
-                      /fvc::interpolate(rho.oldTime())
                      - coefft00*phiAbs.oldTime().oldTime()
-                      /fvc::interpolate(rho.oldTime().oldTime())
                     )
                   - (
                         fvc::interpolate
                         (
-                            rA*rho.oldTime()*
+                            rA*
                             (
-                                coefft0*U.oldTime()
-                              - coefft00*U.oldTime().oldTime()
+                                coefft0*rho.oldTime()*U.oldTime()
+                              - coefft00*rho.oldTime().oldTime()*U.oldTime().oldTime()
                             )
                         ) & mesh().Sf()
                     )
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/polyMesh/boundary b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/polyMesh/boundary
index 2016cf671ba467dcf3f8f338c06f985c7260dae7..c98262911110db1c328760886467e026b70f6b29 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/polyMesh/boundary
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/polyMesh/boundary
@@ -21,40 +21,40 @@ FoamFile
     {
         type            wall;
         nFaces          172;
-        startFace       3334;
+        startFace       3294;
     }
     inlet
     {
         type            patch;
         nFaces          20;
-        startFace       3506;
+        startFace       3466;
     }
     outlet
     {
         type            patch;
         nFaces          20;
-        startFace       3526;
+        startFace       3486;
     }
     cycLeft_half0
     {
         type            cyclic;
-        nFaces          0;
-        startFace       3546;
+        nFaces          20;
+        startFace       3506;
         matchTolerance  0.0001;
         neighbourPatch  cycLeft_half1;
     }
     cycLeft_half1
     {
         type            cyclic;
-        nFaces          0;
-        startFace       3546;
+        nFaces          20;
+        startFace       3526;
         matchTolerance  0.0001;
         neighbourPatch  cycLeft_half0;
     }
     cycRight_half0
     {
         type            cyclic;
-        nFaces          0;
+        nFaces          20;
         startFace       3546;
         matchTolerance  0.0001;
         neighbourPatch  cycRight_half1;
@@ -62,8 +62,8 @@ FoamFile
     cycRight_half1
     {
         type            cyclic;
-        nFaces          0;
-        startFace       3546;
+        nFaces          20;
+        startFace       3566;
         matchTolerance  0.0001;
         neighbourPatch  cycRight_half0;
     }
@@ -71,7 +71,7 @@ FoamFile
     {
         type            empty;
         nFaces          3440;
-        startFace       3546;
+        startFace       3586;
     }
 )
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution
index 574d2ac6a419ce943ed08c9a7cb330f2c96ffc43..672f614cf5e5aff9be0f46a9ff3e0a9f880b1cbe 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution
@@ -85,8 +85,7 @@ PIMPLE
 
 additional
 {
-    pressureWork    true;
-    pressureWorkTimeDerivative true;
+    pressureWorkTimeDerivative false;
 }
 
 relaxationFactors
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes
index 9b46cd8f1cc706b1e5439be1f3c1a031eafeb690..a74c78c1fcd48af64736074f8268b1b22ec5e4d3 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(phi,K)      Gauss linear;
+    div(phi,K)      Gauss upwind;
     div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
index 2927bf0ea22c1315df025babf1c54bbd6fbafa64..ff633dd7d5e9cfb9cb0e3853802a1c86b97c7889 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
@@ -83,19 +83,18 @@ PIMPLE
 
 additional
 {
-    pressureWork    true;
-    pressureWorkTimeDerivative true;
+    pressureWorkTimeDerivative false;
 }
 
 relaxationFactors
 {
     fields
     {
-        ".*Final"       1;
+        ".*"       1;
     }
     equations
     {
-        ".*Final"       1;
+        ".*"       1;
     }
 }
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes
index 7a96d901447dcddc0c0030f4353b902498605e6b..15fb830b80142e26c4420acc5ef19a3f78d3f82c 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(phi,K)      Gauss linear;
+    div(phi,K)      Gauss upwind;
     div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
index f5ac901b11cb2e0dd831d58c915b37b30d19edc6..3aeefebadd78710314a2f7ea2f720cbbf482ec91 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
@@ -89,7 +89,6 @@ potentialFlow
 
 additional
 {
-    pressureWork    false;
     pressureWorkTimeDerivative false;
 }
 
@@ -97,11 +96,11 @@ relaxationFactors
 {
     fields
     {
-        ".*Final"       1;
+        ".*"       1;
     }
     equations
     {
-        ".*Final"       1;
+        ".*"       1;
     }
 }