diff --git a/applications/solvers/combustion/reactingFoam/setRDeltaT.H b/applications/solvers/combustion/reactingFoam/setRDeltaT.H
index 73423173c7..7e806d0997 100644
--- a/applications/solvers/combustion/reactingFoam/setRDeltaT.H
+++ b/applications/solvers/combustion/reactingFoam/setRDeltaT.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,53 +28,43 @@ License
 
     const dictionary& pimpleDict = pimple.dict();
 
-    // Maximum flow Courant number
-    scalar maxCo(pimpleDict.lookup<scalar>("maxCo"));
-
-    // Maximum time scale
-    scalar maxDeltaT(pimpleDict.lookupOrDefault<scalar>("maxDeltaT", great));
-
-    // Smoothing parameter (0-1) when smoothing iterations > 0
-    scalar rDeltaTSmoothingCoeff
-    (
-        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
-    );
-
-    // Damping coefficient (1-0)
-    scalar rDeltaTDampingCoeff
-    (
-        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
-    );
-
-    // Maximum change in cell temperature per iteration
-    // (relative to previous value)
-    scalar alphaTemp(pimpleDict.lookupOrDefault("alphaTemp", 0.05));
-
-    // Maximum change in cell concentration per iteration
-    // (relative to reference value)
-    scalar alphaY(pimpleDict.lookupOrDefault("alphaY", 1.0));
-
     Info<< "Time scales min/max:" << endl;
 
     // Cache old reciprocal time scale field
-    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+    const volScalarField rDeltaT0("rDeltaT0", rDeltaT);
 
     // Flow time scale
     {
+        // Maximum flow Courant number
+        const scalar maxCo(pimpleDict.lookup<scalar>("maxCo"));
+
         rDeltaT.ref() =
         (
             fvc::surfaceSum(mag(phi))()()
            /((2*maxCo)*mesh.V()*rho())
         );
 
-        // Limit the largest time scale
-        rDeltaT.max(1/maxDeltaT);
+        // Limit the largest time step
+        if (pimpleDict.found("maxDeltaT"))
+        {
+            rDeltaT.max(1/pimpleDict.lookup<scalar>("maxDeltaT"));
+        }
+
+        // Limit the smallest time step
+        if (pimpleDict.found("minDeltaT"))
+        {
+            rDeltaT.min(1/pimpleDict.lookup<scalar>("minDeltaT"));
+        }
 
         Info<< "    Flow        = "
             << 1/gMax(rDeltaT.primitiveField()) << ", "
             << 1/gMin(rDeltaT.primitiveField()) << endl;
     }
 
+    // Maximum change in cell temperature per iteration
+    // (relative to previous value)
+    const scalar alphaTemp(pimpleDict.lookupOrDefault("alphaTemp", 0.05));
+
     // Heat release rate time scale
     if (alphaTemp < 1)
     {
@@ -91,9 +81,13 @@ License
     }
 
     // Reaction rate time scale
-    if (alphaY < 1)
+    if (pimpleDict.found("alphaY"))
     {
-        dictionary Yref(pimpleDict.subDict("Yref"));
+        // Maximum change in cell concentration per iteration
+        // (relative to reference value)
+        const scalar alphaY(pimpleDict.lookup<scalar>("alphaY"));
+
+        const dictionary Yref(pimpleDict.subDict("Yref"));
 
         volScalarField::Internal rDeltaTY
         (
@@ -152,6 +146,12 @@ License
     // Update the boundary values of the reciprocal time-step
     rDeltaT.correctBoundaryConditions();
 
+    // Smoothing parameter (0-1) when smoothing iterations > 0
+    const scalar rDeltaTSmoothingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
     // Spatially smooth the time scale field
     if (rDeltaTSmoothingCoeff < 1)
     {
@@ -163,10 +163,16 @@ License
     // - only increase at a fraction of old time scale
     if
     (
-        rDeltaTDampingCoeff < 1
+        pimpleDict.found("rDeltaTDampingCoeff")
      && runTime.timeIndex() > runTime.startTimeIndex() + 1
     )
     {
+        // Damping coefficient (1-0)
+        const scalar rDeltaTDampingCoeff
+        (
+            pimpleDict.lookup<scalar>("rDeltaTDampingCoeff")
+        );
+
         rDeltaT = max
         (
             rDeltaT,
diff --git a/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H b/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
index 15ccd343e5..44ef88e34b 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
@@ -3,34 +3,39 @@
 
     const dictionary& pimpleDict = pimple.dict();
 
-    scalar maxCo
+    const scalar maxCo
     (
         pimpleDict.lookupOrDefault<scalar>("maxCo", 0.8)
     );
 
-    scalar rDeltaTSmoothingCoeff
+    const scalar rDeltaTSmoothingCoeff
     (
         pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.02)
     );
 
-    scalar rDeltaTDampingCoeff
-    (
-        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
-    );
-
-    scalar maxDeltaT
+    const scalar maxDeltaT
     (
         pimpleDict.lookupOrDefault<scalar>("maxDeltaT", great)
     );
 
-    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+    const scalar minDeltaT
+    (
+        pimpleDict.lookupOrDefault<scalar>("minDeltaT", small)
+    );
+
+    const volScalarField rDeltaT0("rDeltaT0", rDeltaT);
 
     // Set the reciprocal time-step from the local Courant number
-    rDeltaT.ref() = max
+    // and maximum and minimum time-steps
+    rDeltaT.ref() = min
     (
-        1/dimensionedScalar(dimTime, maxDeltaT),
-        fvc::surfaceSum(mag(phi))()()
-       /((2*maxCo)*mesh.V()*rho())
+        1/dimensionedScalar(dimTime, minDeltaT),
+        max
+        (
+            1/dimensionedScalar(dimTime, maxDeltaT),
+            fvc::surfaceSum(mag(phi))()()
+           /((2*maxCo)*mesh.V()*rho())
+        )
     );
 
     if (pimple.transonic())
@@ -70,10 +75,16 @@
     // - only increase at a fraction of old time scale
     if
     (
-        rDeltaTDampingCoeff < 1.0
+        pimpleDict.found("rDeltaTDampingCoeff")
      && runTime.timeIndex() > runTime.startTimeIndex() + 1
     )
     {
+        // Damping coefficient (1-0)
+        const scalar rDeltaTDampingCoeff
+        (
+            pimpleDict.lookup<scalar>("rDeltaTDampingCoeff")
+        );
+
         rDeltaT =
             rDeltaT0
            *max(rDeltaT/rDeltaT0, scalar(1) - rDeltaTDampingCoeff);
diff --git a/applications/solvers/incompressible/pimpleFoam/setRDeltaT.H b/applications/solvers/incompressible/pimpleFoam/setRDeltaT.H
index 263012242c..3d40b4bbee 100644
--- a/applications/solvers/incompressible/pimpleFoam/setRDeltaT.H
+++ b/applications/solvers/incompressible/pimpleFoam/setRDeltaT.H
@@ -3,34 +3,44 @@
 
     const dictionary& pimpleDict = pimple.dict();
 
-    scalar maxCo
+    const scalar maxCo
     (
         pimpleDict.lookupOrDefault<scalar>("maxCo", 0.8)
     );
 
-    scalar rDeltaTSmoothingCoeff
+    const scalar rDeltaTSmoothingCoeff
     (
         pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.02)
     );
 
-    scalar rDeltaTDampingCoeff
+    const scalar rDeltaTDampingCoeff
     (
         pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
     );
 
-    scalar maxDeltaT
+    const scalar maxDeltaT
     (
         pimpleDict.lookupOrDefault<scalar>("maxDeltaT", great)
     );
 
-    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+    const scalar minDeltaT
+    (
+        pimpleDict.lookupOrDefault<scalar>("minDeltaT", small)
+    );
+
+    const volScalarField rDeltaT0("rDeltaT0", rDeltaT);
 
     // Set the reciprocal time-step from the local Courant number
-    rDeltaT.ref() = max
+    // and maximum and minimum time-steps
+    rDeltaT.ref() = min
     (
-        1/dimensionedScalar(dimTime, maxDeltaT),
-        fvc::surfaceSum(mag(phi))()()
-       /((2*maxCo)*mesh.V())
+        1/dimensionedScalar(dimTime, minDeltaT),
+        max
+        (
+            1/dimensionedScalar(dimTime, maxDeltaT),
+            fvc::surfaceSum(mag(phi))()()
+           /((2*maxCo)*mesh.V())
+        )
     );
 
     // Update the boundary values of the reciprocal time-step
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/setRDeltaT.H b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/setRDeltaT.H
index a297c72e4e..2cd961a153 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/setRDeltaT.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/setRDeltaT.H
@@ -3,17 +3,22 @@
 
     const dictionary& pimpleDict = pimple.dict();
 
-    scalar maxCo
+    const scalar maxCo
     (
         pimpleDict.lookupOrDefault<scalar>("maxCo", 0.2)
     );
 
-    scalar maxDeltaT
+    const scalar maxDeltaT
     (
         pimpleDict.lookupOrDefault<scalar>("maxDeltaT", great)
     );
 
-    scalar rDeltaTSmoothingCoeff
+    const scalar minDeltaT
+    (
+        pimpleDict.lookupOrDefault<scalar>("minDeltaT", small)
+    );
+
+    const scalar rDeltaTSmoothingCoeff
     (
         pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.02)
     );
@@ -26,11 +31,16 @@
     }
 
     // Set the reciprocal time-step from the local Courant number
-    rDeltaT.ref() = max
+    // and maximum and minimum time-steps
+    rDeltaT.ref() = min
     (
-        1/dimensionedScalar(dimTime, maxDeltaT),
-        fvc::surfaceSum(maxPhi)()()
-       /((2*maxCo)*mesh.V())
+        1/dimensionedScalar(dimTime, minDeltaT),
+        max
+        (
+            1/dimensionedScalar(dimTime, maxDeltaT),
+            fvc::surfaceSum(maxPhi)()()
+           /((2*maxCo)*mesh.V())
+        )
     );
 
     // Update the boundary values of the reciprocal time-step
diff --git a/src/twoPhaseModels/twoPhaseMixture/VoF/setRDeltaT.H b/src/twoPhaseModels/twoPhaseMixture/VoF/setRDeltaT.H
index 5748069779..29dfca5700 100644
--- a/src/twoPhaseModels/twoPhaseMixture/VoF/setRDeltaT.H
+++ b/src/twoPhaseModels/twoPhaseMixture/VoF/setRDeltaT.H
@@ -3,64 +3,69 @@
 
     const dictionary& pimpleDict = pimple.dict();
 
-    scalar maxCo
+    const scalar maxCo
     (
         pimpleDict.lookupOrDefault<scalar>("maxCo", 0.9)
     );
 
-    scalar maxAlphaCo
+    const scalar maxAlphaCo
     (
         pimpleDict.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
     );
 
-    scalar rDeltaTSmoothingCoeff
+    const scalar rDeltaTSmoothingCoeff
     (
         pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
     );
 
-    label nAlphaSpreadIter
+    const label nAlphaSpreadIter
     (
         pimpleDict.lookupOrDefault<label>("nAlphaSpreadIter", 1)
     );
 
-    scalar alphaSpreadDiff
+    const scalar alphaSpreadDiff
     (
         pimpleDict.lookupOrDefault<scalar>("alphaSpreadDiff", 0.2)
     );
 
-    scalar alphaSpreadMax
+    const scalar alphaSpreadMax
     (
         pimpleDict.lookupOrDefault<scalar>("alphaSpreadMax", 0.99)
     );
 
-    scalar alphaSpreadMin
+    const scalar alphaSpreadMin
     (
         pimpleDict.lookupOrDefault<scalar>("alphaSpreadMin", 0.01)
     );
 
-    label nAlphaSweepIter
+    const label nAlphaSweepIter
     (
         pimpleDict.lookupOrDefault<label>("nAlphaSweepIter", 5)
     );
 
-    scalar rDeltaTDampingCoeff
-    (
-        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
-    );
-
-    scalar maxDeltaT
+    const scalar maxDeltaT
     (
         pimpleDict.lookupOrDefault<scalar>("maxDeltaT", great)
     );
 
-    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+    const scalar minDeltaT
+    (
+        pimpleDict.lookupOrDefault<scalar>("minDeltaT", small)
+    );
+
+    const volScalarField rDeltaT0("rDeltaT0", rDeltaT);
 
     // Set the reciprocal time-step from the local Courant number
-    rDeltaT.ref() = max
+    // and maximum and minimum time-steps
+    rDeltaT.ref() = min
     (
-        1/dimensionedScalar(dimTime, maxDeltaT),
-        fvc::surfaceSum(mag(phi))()()
-       /((2*maxCo)*mesh.V())
+        1/dimensionedScalar(dimTime, minDeltaT),
+        max
+        (
+            1/dimensionedScalar(dimTime, maxDeltaT),
+            fvc::surfaceSum(mag(phi))()()
+           /((2*maxCo)*mesh.V())
+        )
     );
 
     if (maxAlphaCo < maxCo)
@@ -119,10 +124,16 @@
     // - only increase at a fraction of old time scale
     if
     (
-        rDeltaTDampingCoeff < 1.0
+        pimpleDict.found("rDeltaTDampingCoeff")
      && runTime.timeIndex() > runTime.startTimeIndex() + 1
     )
     {
+        // Damping coefficient (1-0)
+        const scalar rDeltaTDampingCoeff
+        (
+            pimpleDict.lookup<scalar>("rDeltaTDampingCoeff")
+        );
+
         rDeltaT = max
         (
             rDeltaT,