reactingTwoPhaseEulerFoam: Added experimental run-time selectable LTS support

Select LTS via the ddtScheme:

    ddtSchemes
    {
        default         localEuler rDeltaT;
    }

The LTS algorithm is currently controlled with the standard settings in
controlDict, e.g.:

    maxCo           0.5;
    maxDeltaT       2e-8;

with the addition of the optional rDeltaT smoothing coefficient:

    rDeltaTSmoothingCoeff 0.02;

which defaults to 0.02.

    ddtSchemes
    {
        default         localEuler rDeltaT;
    }

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)
+        )
+    );
+}

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"
+    if (!LTS)
-    #include "setInitialDeltaT.H"
+    {
+        #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;

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;
+}