diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/Make/options b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/Make/options
index 5216efaa9c603cfef603c058b70812fd3b0a650b..622878cc520de18f4f27471139afb41d1788a077 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/Make/options
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -ggdb3 \
     -ItwoPhaseSystem/lnInclude \
     -I../phaseSystems/lnInclude \
     -I../interfacialModels/lnInclude \
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
index 3b225202e917f8839dd4d821cc15c2cc47bbd6d2..f268beae620af98d438252d28848bf34f2e87746 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -7,8 +7,7 @@ volScalarField rAU1
     1.0
    /(
         U1Eqn.A()
-      + max(phase1.residualAlpha() - alpha1, scalar(0))
-       *rho1/runTime.deltaT()
+      + byDt(max(phase1.residualAlpha() - alpha1, scalar(0))*rho1)
     )
 );
 volScalarField rAU2
@@ -17,8 +16,7 @@ volScalarField rAU2
     1.0
    /(
         U2Eqn.A()
-      + max(phase2.residualAlpha() - alpha2, scalar(0))
-       *rho2/runTime.deltaT()
+      + byDt(max(phase2.residualAlpha() - alpha2, scalar(0))*rho2)
     )
 );
 
@@ -101,8 +99,8 @@ while (pimple.correct())
         rAU1
        *(
             U1Eqn.H()
-          + max(phase1.residualAlpha() - alpha1, scalar(0))
-           *rho1*U1.oldTime()/runTime.deltaT()
+          + byDt(max(phase1.residualAlpha() - alpha1, scalar(0))*rho1)
+           *U1.oldTime()
         );
 
     volVectorField HbyA2
@@ -114,8 +112,8 @@ while (pimple.correct())
         rAU2
        *(
             U2Eqn.H()
-         +  max(phase2.residualAlpha() - alpha2, scalar(0))
-           *rho2*U2.oldTime()/runTime.deltaT()
+         +  byDt(max(phase2.residualAlpha() - alpha2, scalar(0))*rho2)
+           *U2.oldTime()
         );
 
     surfaceScalarField ghSnGradRho
@@ -175,11 +173,9 @@ while (pimple.correct())
     (
         IOobject::groupName("phiHbyA", phase1.name()),
         fvc::flux(HbyA1)
-      + phiCorrCoeff1*fvc::interpolate(alpha1.oldTime()*rho1.oldTime()*rAU1)
-       *(
-            MRF.absolute(phi1.oldTime())
-          - fvc::flux(U1.oldTime())
-        )/runTime.deltaT()
+      + phiCorrCoeff1
+       *fvc::interpolate(byDt(alpha1.oldTime()*rho1.oldTime()*rAU1))
+       *(MRF.absolute(phi1.oldTime()) - fvc::flux(U1.oldTime()))
       - phiF1()
       - phig1
     );
@@ -189,11 +185,9 @@ while (pimple.correct())
     (
         IOobject::groupName("phiHbyA", phase2.name()),
         fvc::flux(HbyA2)
-      + phiCorrCoeff2*fvc::interpolate(alpha2.oldTime()*rho2.oldTime()*rAU2)
-       *(
-            MRF.absolute(phi2.oldTime())
-          - fvc::flux(U2.oldTime())
-        )/runTime.deltaT()
+      + phiCorrCoeff2
+       *fvc::interpolate(byDt(alpha2.oldTime()*rho2.oldTime()*rAU2))
+       *(MRF.absolute(phi2.oldTime()) - fvc::flux(U2.oldTime()))
       - phiF2()
       - phig2
     );
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/DDtU.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/DDtU.H
index c7bccf66eac5e0a56750e2837027e16394903c56..557eb02678c560ccd8d32ce9dca262b5f05ac4f0 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/DDtU.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/DDtU.H
@@ -1,2 +1,2 @@
-ddtPhi1 = fvc::ddt(phi1);
-ddtPhi2 = fvc::ddt(phi2);
+tddtPhi1 = fvc::ddt(phi1);
+tddtPhi2 = fvc::ddt(phi2);
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createDDtU.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createDDtU.H
index 14e66e47a833f9cbf06eba9d03904b3c675599ba..b4c6de5f1a5be5a5391d4382e8b9cff163233719 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createDDtU.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createDDtU.H
@@ -1,2 +1,7 @@
-surfaceScalarField ddtPhi1(fvc::ddt(phi1));
-surfaceScalarField ddtPhi2(fvc::ddt(phi2));
+tmp<surfaceScalarField> tddtPhi1;
+tmp<surfaceScalarField> tddtPhi2;
+
+if (faceMomentum)
+{
+    #include "DDtU.H"
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createRDeltaTf.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createRDeltaTf.H
similarity index 100%
rename from applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createRDeltaTf.H
rename to applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createRDeltaTf.H
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index 7950d76b94470df97aaeebb600a3898177dee675..05bca24c45cf9fe63926cc936d949e33f0edddc9 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -159,7 +159,7 @@ while (pimple.correct())
             (alphaRhof10 + Vmf)
            *byDt(MRF.absolute(phi1.oldTime()))
           + fvc::flux(U1Eqn.H())
-          + Vmf*ddtPhi2
+          + Vmf*tddtPhi2()
           + Kdf*MRF.absolute(phi2)
           - Ff1()
         );
@@ -177,7 +177,7 @@ while (pimple.correct())
             (alphaRhof20 + Vmf)
            *byDt(MRF.absolute(phi2.oldTime()))
           + fvc::flux(U2Eqn.H())
-          + Vmf*ddtPhi1
+          + Vmf*tddtPhi1()
           + Kdf*MRF.absolute(phi1)
           - Ff2()
        );
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
index bebe5620011e4b038e11949e315d5f332e1d0e98..47b6fe63894497e69b59ef31a029dd414b8f64b2 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
@@ -42,6 +42,18 @@ Description
 
 namespace Foam
 {
+    tmp<volScalarField> byDt(const volScalarField& vf)
+    {
+        if (fv::localEulerDdt::enabled(vf.mesh()))
+        {
+            return fv::localEulerDdt::localRDeltaT(vf.mesh())*vf;
+        }
+        else
+        {
+            return vf/vf.mesh().time().deltaT();
+        }
+    }
+
     tmp<surfaceScalarField> byDt(const surfaceScalarField& sf)
     {
         if (fv::localEulerDdt::enabled(sf.mesh()))
@@ -89,7 +101,7 @@ int main(int argc, char *argv[])
         )
     );
 
-    #include "createRDeltaTf.H"
+    #include "pUf/createRDeltaTf.H"
     #include "pUf/createDDtU.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //