diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createRDeltaT.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createRDeltaT.H
new file mode 100644
index 0000000000000000000000000000000000000000..b45dff137e4df5fb5005062e10f75c289812f4d1
--- /dev/null
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createRDeltaT.H
@@ -0,0 +1,47 @@
+bool LTS =
+    word(mesh.ddtScheme("default"))
+ == fv::localEulerDdtScheme<scalar>::typeName;
+
+tmp<volScalarField> trDeltaT;
+tmp<volScalarField> trSubDeltaT;
+
+if (LTS)
+{
+    scalar maxDeltaT
+    (
+        pimple.dict().lookupOrDefault<scalar>("maxDeltaT", GREAT)
+    );
+
+    trDeltaT = tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "rDeltaT",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    trSubDeltaT = tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "rSubDeltaT",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT)
+        )
+    );
+}
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
index 27826f12a1934a8d4c5745c09e0b0e52bc0b938c..b0434b45d67734e5617dc73f6c566ecb3d6f158b 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
@@ -38,6 +38,8 @@ Description
 #include "PhaseCompressibleTurbulenceModel.H"
 #include "fixedFluxPressureFvPatchScalarField.H"
 #include "pimpleControl.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,11 +52,16 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
+    #include "createRDeltaT.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
-    #include "readTimeControls.H"
-    #include "CourantNos.H"
-    #include "setInitialDeltaT.H"
+
+    if (!LTS)
+    {
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
 
     Switch faceMomentum
     (
@@ -78,8 +85,16 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "CourantNos.H"
-        #include "setDeltaT.H"
+
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "CourantNos.H"
+            #include "setDeltaT.H"
+        }
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/setRDeltaT.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/setRDeltaT.H
new file mode 100644
index 0000000000000000000000000000000000000000..99e08c3f565e9a3d151bcf60564bff905b615099
--- /dev/null
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/setRDeltaT.H
@@ -0,0 +1,35 @@
+{
+    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rSubDeltaT = trSubDeltaT();
+
+    scalar rDeltaTSmoothingCoeff
+    (
+        runTime.controlDict().lookupOrDefault<scalar>
+        (
+            "rDeltaTSmoothingCoeff",
+            0.02
+        )
+    );
+
+    // Set the reciprocal time-step from the local Courant number
+    rDeltaT.dimensionedInternalField() = max
+    (
+        1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+        fvc::surfaceSum(max(mag(phi1), mag(phi2)))().dimensionedInternalField()
+       /((2*maxCo)*mesh.V())
+    );
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+
+    const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+    label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+
+    rSubDeltaT = rDeltaT*nAlphaSubCycles;
+}