reactingTwoPhaseEulerFoam: Added LTS support to the face-momentum algorithm

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H

@@ -29,37 +29,25 @@ for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
     }
 
     {
-        tmp<fvScalarMatrix> E2eqn(phase2.heEqn());
+        tmp<fvScalarMatrix> E2Eqn(phase2.heEqn());
 
-        if (E2eqn.valid())
+        if (E2Eqn.valid())
         {
-            E2eqn =
+            E2Eqn =
             (
-                E2eqn
+                E2Eqn
              ==
                *heatTransfer[phase2.name()]
               + alpha2*rho2*(U2&g)
               + fvOptions(alpha2, rho2, phase2.thermo().he())
             );
 
-            E2eqn->relax();
+            E2Eqn->relax();
-            fvOptions.constrain(E2eqn.ref());
+            fvOptions.constrain(E2Eqn.ref());
-            E2eqn->solve();
+            E2Eqn->solve();
         }
     }
 
     fluid.correctThermo();
     fluid.correct();
 }
-
-Info<< phase1.name() << " min/max T "
-    << min(phase1.thermo().T()).value()
-    << " - "
-    << max(phase1.thermo().T()).value()
-    << endl;
-
-Info<< phase2.name() << " min/max T "
-    << min(phase2.thermo().T()).value()
-    << " - "
-    << max(phase2.thermo().T()).value()
-    << endl;

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createRDeltaTf.H

@@ -0,0 +1,21 @@
+tmp<surfaceScalarField> trDeltaTf;
+
+if (LTS && faceMomentum)
+{
+    trDeltaTf = tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                fv::localEulerDdt::rDeltaTfName,
+                runTime.timeName(),
+                mesh,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("one", dimless/dimTime, 1)
+        )
+    );
+}

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

@@ -1,9 +1,2 @@
-    ddtPhi1 =
+ddtPhi1 = fvc::ddt(phi1);
-    (
+ddtPhi2 = fvc::ddt(phi2);
-        (phi1 - phi1.oldTime())/runTime.deltaT()
-    );
-
-    ddtPhi2 =
-    (
-        (phi2 - phi2.oldTime())/runTime.deltaT()
-    );

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

@@ -1,9 +1,2 @@
-    surfaceScalarField ddtPhi1
+surfaceScalarField ddtPhi1(fvc::ddt(phi1));
-    (
+surfaceScalarField ddtPhi2(fvc::ddt(phi2));
-        (phi1 - phi1.oldTime())/runTime.deltaT()
-    );
-
-    surfaceScalarField ddtPhi2
-    (
-        (phi2 - phi2.oldTime())/runTime.deltaT()
-    );

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

@@ -27,15 +27,17 @@ surfaceScalarField Kdf("Kdf", fluid.Kdf());
 // Virtual-mass coefficient
 surfaceScalarField Vmf("Vmf", fluid.Vmf());
 
+const surfaceScalarField& rDeltaTf = fv::localEulerDdt::localRDeltaTf(mesh);
+
 surfaceScalarField rAUf1
 (
     IOobject::groupName("rAUf", phase1.name()),
     1.0
    /(
-        (alphaRhof10 + Vmf)/runTime.deltaT()
+        rDeltaTf*(alphaRhof10 + Vmf)
       + fvc::interpolate(U1Eqn.A())
       + Kdf
-   )
+    )
 );
 
 surfaceScalarField rAUf2
@@ -43,10 +45,10 @@ surfaceScalarField rAUf2
     IOobject::groupName("rAUf", phase2.name()),
     1.0
    /(
-        (alphaRhof20 + Vmf)/runTime.deltaT()
+        rDeltaTf*(alphaRhof20 + Vmf)
       + fvc::interpolate(U2Eqn.A())
       + Kdf
-   )
+    )
 );
 
 
@@ -157,7 +159,7 @@ while (pimple.correct())
         rAUf1
        *(
             (alphaRhof10 + Vmf)
-           *MRF.absolute(phi1.oldTime())/runTime.deltaT()
+           *rDeltaTf*MRF.absolute(phi1.oldTime())
           + fvc::flux(U1Eqn.H())
           + Vmf*ddtPhi2
           + Kdf*MRF.absolute(phi2)
@@ -175,7 +177,7 @@ while (pimple.correct())
         rAUf2
        *(
             (alphaRhof20 + Vmf)
-           *MRF.absolute(phi2.oldTime())/runTime.deltaT()
+           *rDeltaTf*MRF.absolute(phi2.oldTime())
           + fvc::flux(U2Eqn.H())
           + Vmf*ddtPhi1
           + Kdf*MRF.absolute(phi1)

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C

@@ -73,6 +73,7 @@ int main(int argc, char *argv[])
         )
     );
 
+    #include "createRDeltaTf.H"
     #include "pUf/createDDtU.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -91,6 +92,10 @@ int main(int argc, char *argv[])
         if (LTS)
         {
             #include "setRDeltaT.H"
+            if (faceMomentum)
+            {
+                #include "setRDeltaTf.H"
+            }
         }
         else
         {

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/setRDeltaTf.H

@@ -0,0 +1 @@
+trDeltaTf.ref() = fvc::interpolate(fv::localEulerDdt::localRDeltaT(mesh));

src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C

@@ -317,6 +317,33 @@ EulerDdtScheme<Type>::fvcDdt
 }
 
 
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+EulerDdtScheme<Type>::fvcDdt
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sf
+)
+{
+    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddt("+sf.name()+')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    return tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            ddtIOobject,
+            rDeltaT*(sf - sf.oldTime())
+        )
+    );
+}
+
+
 template<class Type>
 tmp<fvMatrix<Type>>
 EulerDdtScheme<Type>::fvmDdt

src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H

@@ -124,6 +124,11 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& psi
         );
 
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> fvcDdt
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+        );
+
         tmp<fvMatrix<Type>> fvmDdt
         (
             const GeometricField<Type, fvPatchField, volMesh>&

src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C

@@ -132,6 +132,22 @@ tmp<fvMatrix<Type>> ddtScheme<Type>::fvmDdt
 }
 
 
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> ddtScheme<Type>::fvcDdt
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sf
+)
+{
+    NotImplemented;
+
+    return tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+    (
+        GeometricField<Type, fvsPatchField, surfaceMesh>::null()
+    );
+}
+
+
+
 template<class Type>
 tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
 (

src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H

@@ -167,6 +167,11 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>&
         ) = 0;
 
+        virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> fvcDdt
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+        );
+
         virtual tmp<fvMatrix<Type>> fvmDdt
         (
             const GeometricField<Type, fvPatchField, volMesh>&

src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdt.C

@@ -29,6 +29,7 @@ License
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 Foam::word Foam::fv::localEulerDdt::rDeltaTName("rDeltaT");
+Foam::word Foam::fv::localEulerDdt::rDeltaTfName("rDeltaTf");
 Foam::word Foam::fv::localEulerDdt::rSubDeltaTName("rSubDeltaT");
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -53,6 +54,18 @@ const Foam::volScalarField& Foam::fv::localEulerDdt::localRDeltaT
 }
 
 
+const Foam::surfaceScalarField& Foam::fv::localEulerDdt::localRDeltaTf
+(
+    const fvMesh& mesh
+)
+{
+    return mesh.objectRegistry::lookupObject<surfaceScalarField>
+    (
+        rDeltaTfName
+    );
+}
+
+
 Foam::tmp<Foam::volScalarField> Foam::fv::localEulerDdt::localRSubDeltaT
 (
     const fvMesh& mesh,

src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C

@@ -46,6 +46,13 @@ const volScalarField& localEulerDdtScheme<Type>::localRDeltaT() const
 }
 
 
+template<class Type>
+const surfaceScalarField& localEulerDdtScheme<Type>::localRDeltaTf() const
+{
+    return localEulerDdt::localRDeltaTf(mesh());
+}
+
+
 template<class Type>
 tmp<GeometricField<Type, fvPatchField, volMesh>>
 localEulerDdtScheme<Type>::fvcDdt
@@ -333,6 +340,33 @@ localEulerDdtScheme<Type>::fvcDdt
 }
 
 
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+localEulerDdtScheme<Type>::fvcDdt
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sf
+)
+{
+    const surfaceScalarField& rDeltaT = localRDeltaTf();
+
+    IOobject ddtIOobject
+    (
+        "ddt("+sf.name()+')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    return tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            ddtIOobject,
+            rDeltaT*(sf - sf.oldTime())
+        )
+    );
+}
+
+
 template<class Type>
 tmp<fvMatrix<Type>>
 localEulerDdtScheme<Type>::fvmDdt

src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H

@@ -69,6 +69,9 @@ public:
     //- Name of the reciprocal local time-step field
     static word rDeltaTName;
 
+    //- Name of the reciprocal local face time-step field
+    static word rDeltaTfName;
+
     //- Name of the reciprocal local sub-cycling time-step field
     static word rSubDeltaTName;
 
@@ -88,6 +91,10 @@ public:
         //  looked-up from the objectRegistry
         static const volScalarField& localRDeltaT(const fvMesh& mesh);
 
+        //- Return the reciprocal of the local face time-step
+        //  looked-up from the objectRegistry
+        static const surfaceScalarField& localRDeltaTf(const fvMesh& mesh);
+
         //- Calculate and return the reciprocal of the local sub-cycling
         //  time-step
         static tmp<volScalarField> localRSubDeltaT
@@ -119,6 +126,9 @@ class localEulerDdtScheme
         //- Return the reciprocal of the local time-step
         const volScalarField& localRDeltaT() const;
 
+        //- Return the reciprocal of the local face time-step
+        const surfaceScalarField& localRDeltaTf() const;
+
 
 public:
 
@@ -178,6 +188,11 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& psi
         );
 
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> fvcDdt
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+        );
+
         tmp<fvMatrix<Type>> fvmDdt
         (
             const GeometricField<Type, fvPatchField, volMesh>&

src/finiteVolume/finiteVolume/fvc/fvcDdt.C

@@ -125,6 +125,21 @@ ddt
 }
 
 
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+ddt
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sf
+)
+{
+    return fv::ddtScheme<Type>::New
+    (
+        sf.mesh(),
+        sf.mesh().ddtScheme("ddt(" + sf.name() + ')')
+    ).ref().fvcDdt(sf);
+}
+
+
 template<class Type>
 tmp<GeometricField<Type, fvPatchField, volMesh>>
 ddt

src/finiteVolume/finiteVolume/fvc/fvcDdt.H

@@ -88,6 +88,12 @@ namespace fvc
         const GeometricField<Type, fvPatchField, volMesh>&
     );
 
+    template<class Type>
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> ddt
+    (
+        const GeometricField<Type, fvsPatchField, surfaceMesh>&
+    );
+
     template<class Type>
     tmp<GeometricField<Type, fvPatchField, volMesh>> ddt
     (