diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H
index 96f8473e9..2a90edea6 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqn.H
@@ -82,7 +82,7 @@
               ? fv::localEulerDdtScheme<scalar>
                 (
                     mesh,
-                    trSubDeltaT().name()
+                    nAlphaSubCycles > 1 ? "rSubDeltaT" : "rDeltaT"
                 ).fvmDdt(alpha1)
               : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
             )
diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
index 5db9f1654..71de07723 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
@@ -13,6 +13,16 @@ if (nAlphaSubCycles > 1)
         dimensionedScalar("0", rhoPhi.dimensions(), 0)
     );
 
+    tmp<volScalarField> trSubDeltaT;
+
+    if (LTS)
+    {
+        trSubDeltaT = tmp<volScalarField>
+        (
+            new volScalarField("rSubDeltaT", trDeltaT()*nAlphaSubCycles)
+        );
+    }
+
     for
     (
         subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
diff --git a/applications/solvers/multiphase/interFoam/createRDeltaT.H b/applications/solvers/multiphase/interFoam/createRDeltaT.H
deleted file mode 100644
index b45dff137..000000000
--- a/applications/solvers/multiphase/interFoam/createRDeltaT.H
+++ /dev/null
@@ -1,47 +0,0 @@
-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/interFoam/setRDeltaT.H b/applications/solvers/multiphase/interFoam/setRDeltaT.H
index 0752857a9..01bb8aa41 100644
--- a/applications/solvers/multiphase/interFoam/setRDeltaT.H
+++ b/applications/solvers/multiphase/interFoam/setRDeltaT.H
@@ -1,6 +1,5 @@
 {
     volScalarField& rDeltaT = trDeltaT();
-    volScalarField& rSubDeltaT = trSubDeltaT();
 
     const dictionary& pimpleDict = pimple.dict();
 
@@ -134,9 +133,4 @@
             << 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;
 }
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C
index 48b9703eb..3ce13f2c4 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C
@@ -175,6 +175,10 @@ void Foam::twoPhaseSystem::solve()
     label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
     label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
 
+    bool LTS =
+        word(mesh.ddtScheme("default"))
+     == fv::localEulerDdtScheme<scalar>::typeName;
+
     word alphaScheme("div(phi," + alpha1.name() + ')');
     word alpharScheme("div(phir," + alpha1.name() + ')');
 
@@ -316,6 +320,22 @@ void Foam::twoPhaseSystem::solve()
 
         if (nAlphaSubCycles > 1)
         {
+            tmp<volScalarField> trSubDeltaT;
+
+            if (LTS)
+            {
+                const volScalarField& rDeltaT =
+                    mesh.objectRegistry::lookupObject<volScalarField>
+                    (
+                        "rDeltaT"
+                    );
+
+                trSubDeltaT = tmp<volScalarField>
+                (
+                    new volScalarField("rSubDeltaT", rDeltaT*nAlphaSubCycles)
+                );
+            }
+
             for
             (
                 subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
index ecd7671fe..75180f09d 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
@@ -96,7 +96,10 @@ void Foam::MULES::correct
     if (LTS)
     {
         const volScalarField& rDeltaT =
-            mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
+            mesh.objectRegistry::lookupObject<volScalarField>
+            (
+                mesh.time().subCycling() ? "rSubDeltaT" : "rDeltaT"
+            );
 
         limitCorr
         (
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index 83eebdb15..b0205becc 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -118,7 +118,10 @@ void Foam::MULES::explicitSolve
     if (LTS)
     {
         const volScalarField& rDeltaT =
-            mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
+            mesh.objectRegistry::lookupObject<volScalarField>
+            (
+                mesh.time().subCycling() ? "rSubDeltaT" : "rDeltaT"
+            );
 
         limit
         (