reactingTwoPhaseEulerFoam: Updated LTS support in the cell-based momentum algorithm

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/Make/options

-EXE_INC = \
+EXE_INC = -ggdb3 \
     -ItwoPhaseSystem/lnInclude \
     -I../phaseSystems/lnInclude \
     -I../interfacialModels/lnInclude \

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

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/DDtU.H

-ddtPhi1 = fvc::ddt(phi1);
-ddtPhi2 = fvc::ddt(phi2);
+tddtPhi1 = fvc::ddt(phi1);
+tddtPhi2 = fvc::ddt(phi2);

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createDDtU.H

-surfaceScalarField ddtPhi1(fvc::ddt(phi1));
-surfaceScalarField ddtPhi2(fvc::ddt(phi2));
+tmp<surfaceScalarField> tddtPhi1;
+tmp<surfaceScalarField> tddtPhi2;
+
+if (faceMomentum)
+{
+    #include "DDtU.H"
+}

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()
        );

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"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //