Added ddtPhiCorr support for moving mesh

applications/solvers/multiphase/interDyMFoam/interDyMFoam.C

@@ -74,7 +74,7 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // Make the fluxes absolute
-        if (mesh.changing() && correctPhi)
+        if (mesh.changing())
         {
             fvc::makeAbsolute(phi, U);
         }
@@ -99,8 +99,11 @@ int main(int argc, char *argv[])
             #include "correctPhi.H"
         }
 
+        // Keep the absolute fluxes for use in ddtPhiCorr
+        surfaceScalarField phiAbs("phiAbs", phi);
+
         // Make the fluxes relative to the mesh motion
-        if (mesh.changing() && correctPhi)
+        if (mesh.changing())
         {
             fvc::makeRelative(phi, U);
         }

applications/solvers/multiphase/interDyMFoam/pEqn.H

@@ -9,7 +9,7 @@
     (
         "phiU",
         (fvc::interpolate(U) & mesh.Sf())
-    //+ fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
+      + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
     );
 
     phi = phiU +

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

@@ -1,4 +1,4 @@
-/*---------------------------------------------------------------------------*\
+/*---------------------------------------------------------------------------* \
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
@@ -21,7 +21,7 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-    
+
 \*---------------------------------------------------------------------------*/
 
 #include "EulerDdtScheme.H"
@@ -279,7 +279,7 @@ EulerDdtScheme<Type>::fvmDdt
     scalar rDeltaT = 1.0/mesh().time().deltaT().value();
 
     fvm.diag() = rDeltaT*mesh().V();
-    
+
     if (mesh().moving())
     {
         fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
@@ -314,7 +314,7 @@ EulerDdtScheme<Type>::fvmDdt
     scalar rDeltaT = 1.0/mesh().time().deltaT().value();
 
     fvm.diag() = rDeltaT*rho.value()*mesh().V();
-    
+
     if (mesh().moving())
     {
         fvm.source() = rDeltaT
@@ -375,50 +375,30 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr
 (
     const volScalarField& rA,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
+    const fluxFieldType& phiAbs
 )
 {
     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
 
     IOobject ddtIOobject
     (
-        "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+        "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phiAbs.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
 
-    if (mesh().moving())
+    tmp<fluxFieldType> phiCorr =
-    {
+        phiAbs.oldTime() - (fvc::interpolate(U.oldTime()) & mesh().Sf());
-        return tmp<fluxFieldType>
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rDeltaT.dimensions()*rA.dimensions()*phi.dimensions(),
-                    pTraits<typename flux<Type>::type>::zero
-                )
-            )
-        );
-    }
-    else
-    {
-        tmp<fluxFieldType> phiCorr =
-            phi.oldTime() - (fvc::interpolate(U.oldTime()) & mesh().Sf());
 
-        return tmp<fluxFieldType>
+    return tmp<fluxFieldType>
+    (
+        new fluxFieldType
         (
-            new fluxFieldType
+            ddtIOobject,
-            (
+            fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime(), phiCorr())
-                ddtIOobject,
+           *fvc::interpolate(rDeltaT*rA)*phiCorr
-                fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr())
+        )
-               *fvc::interpolate(rDeltaT*rA)*phiCorr
+    );
-            )
-        );
-    }
 }
 
 
@@ -429,7 +409,7 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr
     const volScalarField& rA,
     const volScalarField& rho,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
+    const fluxFieldType& phiAbs
 )
 {
     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
@@ -437,111 +417,94 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr
     IOobject ddtIOobject
     (
         "ddtPhiCorr("
-      + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+      + rA.name() + ','
+      + rho.name() + ','
+      + U.name() + ','
+      + phiAbs.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
 
-    if (mesh().moving())
+    if
+    (
+        U.dimensions() == dimVelocity
+     && phiAbs.dimensions() == dimVelocity*dimArea
+    )
     {
         return tmp<fluxFieldType>
         (
             new fluxFieldType
             (
                 ddtIOobject,
-                mesh(),
+                rDeltaT
-                dimensioned<typename flux<Type>::type>
+               *fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime())
+               *(
+                   fvc::interpolate(rA*rho.oldTime())*phiAbs.oldTime()
+                 - (fvc::interpolate(rA*rho.oldTime()*U.oldTime())
+                  & mesh().Sf())
+                )
+            )
+        );
+    }
+    else if
+    (
+        U.dimensions() == dimVelocity
+     && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                rDeltaT
+               *fvcDdtPhiCoeff
                 (
-                    "0",
+                    U.oldTime(),
-                    rDeltaT.dimensions()*rA.dimensions()
+                    phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
-                   *rho.dimensions()*phi.dimensions(),
+                )
-                    pTraits<typename flux<Type>::type>::zero
+               *(
+                   fvc::interpolate(rA*rho.oldTime())
+                  *phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
+                 - (
+                       fvc::interpolate
+                       (
+                           rA*rho.oldTime()*U.oldTime()
+                       ) & mesh().Sf()
+                 )
+               )
+            )
+        );
+    }
+    else if
+    (
+        U.dimensions() == rho.dimensions()*dimVelocity
+     && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                rDeltaT
+               *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phiAbs.oldTime())
+               *(
+                   fvc::interpolate(rA)*phiAbs.oldTime()
+                 - (fvc::interpolate(rA*U.oldTime()) & mesh().Sf())
                 )
             )
         );
     }
     else
     {
-        if
+        FatalErrorIn
         (
-            U.dimensions() == dimVelocity
+            "EulerDdtScheme<Type>::fvcDdtPhiCorr"
-         && phi.dimensions() == dimVelocity*dimArea
+        )   << "dimensions of phiAbs are not correct"
-        )
+            << abort(FatalError);
-        {
-            return tmp<fluxFieldType>
-            (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
-                   *(
-                        fvc::interpolate(rA*rho.oldTime())*phi.oldTime()
-                      - (fvc::interpolate(rA*rho.oldTime()*U.oldTime())
-                      & mesh().Sf())
-                    )
-                )
-            );
-        }
-        else if 
-        (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
-            (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    rDeltaT
-                   *fvcDdtPhiCoeff
-                    (
-                        U.oldTime(),
-                        phi.oldTime()/fvc::interpolate(rho.oldTime())
-                    )
-                   *(
-                        fvc::interpolate(rA*rho.oldTime())
-                       *phi.oldTime()/fvc::interpolate(rho.oldTime())
-                      - (
-                            fvc::interpolate
-                            (
-                                rA*rho.oldTime()*U.oldTime()
-                            ) & mesh().Sf()
-                        )
-                    )
-                )
-            );
-        }
-        else if 
-        (
-            U.dimensions() == rho.dimensions()*dimVelocity
-         && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
-            (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    rDeltaT
-                   *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())
-                   *(
-                        fvc::interpolate(rA)*phi.oldTime()
-                      - (fvc::interpolate(rA*U.oldTime()) & mesh().Sf())
-                    )
-                )
-            );
-        }
-        else
-        {
-            FatalErrorIn
-            (
-                "EulerDdtScheme<Type>::fvcDdtPhiCorr"
-            )   << "dimensions of phi are not correct"
-                << abort(FatalError);
 
-            return fluxFieldType::null();
+        return fluxFieldType::null();
-        }
     }
 }
 

src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C

@@ -21,7 +21,7 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-    
+
 \*---------------------------------------------------------------------------*/
 
 #include "backwardDdtScheme.H"
@@ -531,43 +531,160 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
         mesh()
     );
 
-    if (mesh().moving())
+    scalar deltaT = deltaT_();
+    scalar deltaT0 = deltaT0_(U);
+
+    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    scalar coefft0  = coefft + coefft00;
+
+    return tmp<fluxFieldType>
+    (
+        new fluxFieldType
+        (
+            ddtIOobject,
+            rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
+           *(
+                fvc::interpolate(rA)
+               *(
+                   coefft0*phi.oldTime()
+                 - coefft00*phi.oldTime().oldTime()
+                )
+              - (
+                    fvc::interpolate
+                    (
+                        rA*
+                        (
+                            coefft0*U.oldTime()
+                          - coefft00*U.oldTime().oldTime()
+                        )
+                    ) & mesh().Sf()
+                )
+            )
+        )
+    );
+}
+
+
+template<class Type>
+tmp<typename backwardDdtScheme<Type>::fluxFieldType>
+backwardDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const volScalarField& rA,
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phiAbs
+)
+{
+    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddtPhiCorr("
+      + rA.name() + ','
+      + rho.name() + ','
+      + U.name() + ','
+      + phiAbs.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    scalar deltaT = deltaT_();
+    scalar deltaT0 = deltaT0_(U);
+
+    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    scalar coefft0  = coefft + coefft00;
+
+    if
+    (
+        U.dimensions() == dimVelocity
+     && phiAbs.dimensions() == dimVelocity*dimArea
+    )
     {
         return tmp<fluxFieldType>
         (
             new fluxFieldType
             (
                 ddtIOobject,
-                mesh(),
+                rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime())
-                dimensioned<typename flux<Type>::type>
+               *(
-                (
+                    coefft0*fvc::interpolate(rA*rho.oldTime())
-                    "0",
+                   *phiAbs.oldTime()
-                    rDeltaT.dimensions()*rA.dimensions()*phi.dimensions(),
+                  - coefft00*fvc::interpolate(rA*rho.oldTime().oldTime())
-                    pTraits<typename flux<Type>::type>::zero
+                   *phiAbs.oldTime().oldTime()
+                  - (
+                        fvc::interpolate
+                        (
+                            rA*
+                            (
+                                coefft0*rho.oldTime()*U.oldTime()
+                              - coefft00*rho.oldTime().oldTime()
+                               *U.oldTime().oldTime()
+                            )
+                        ) & mesh().Sf()
+                    )
                 )
             )
         );
     }
-    else
+    else if
+    (
+        U.dimensions() == dimVelocity
+     && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
     {
-        scalar deltaT = deltaT_();
-        scalar deltaT0 = deltaT0_(U);
-
-        scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-        scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-        scalar coefft0  = coefft + coefft00;
-
         return tmp<fluxFieldType>
         (
             new fluxFieldType
             (
                 ddtIOobject,
-                rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
+                rDeltaT
+               *fvcDdtPhiCoeff
+                (
+                    U.oldTime(),
+                    phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
+                )
+               *(
+                    fvc::interpolate(rA*rho.oldTime())
+                   *(
+                       coefft0*phiAbs.oldTime()
+                      /fvc::interpolate(rho.oldTime())
+                     - coefft00*phiAbs.oldTime().oldTime()
+                      /fvc::interpolate(rho.oldTime().oldTime())
+                    )
+                  - (
+                        fvc::interpolate
+                        (
+                            rA*rho.oldTime()*
+                            (
+                                coefft0*U.oldTime()
+                              - coefft00*U.oldTime().oldTime()
+                            )
+                        ) & mesh().Sf()
+                    )
+                )
+            )
+        );
+    }
+    else if
+    (
+        U.dimensions() == rho.dimensions()*dimVelocity
+     && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                rDeltaT
+               *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phiAbs.oldTime())
                *(
                     fvc::interpolate(rA)
                    *(
-                       coefft0*phi.oldTime()
+                       coefft0*phiAbs.oldTime()
-                     - coefft00*phi.oldTime().oldTime()
+                     - coefft00*phiAbs.oldTime().oldTime()
                     )
                   - (
                         fvc::interpolate
@@ -583,170 +700,15 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
             )
         );
     }
-}
-
-
-template<class Type>
-tmp<typename backwardDdtScheme<Type>::fluxFieldType>
-backwardDdtScheme<Type>::fvcDdtPhiCorr
-(
-    const volScalarField& rA,
-    const volScalarField& rho,
-    const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
-)
-{
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
-
-    IOobject ddtIOobject
-    (
-        "ddtPhiCorr("
-      + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
-        mesh().time().timeName(),
-        mesh()
-    );
-
-    if (mesh().moving())
-    {
-        return tmp<fluxFieldType>
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rDeltaT.dimensions()*rA.dimensions()
-                   *rho.dimensions()*phi.dimensions(),
-                    pTraits<typename flux<Type>::type>::zero
-                )
-            )
-        );
-    }
     else
     {
-        scalar deltaT = deltaT_();
+        FatalErrorIn
-        scalar deltaT0 = deltaT0_(U);
-
-        scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-        scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-        scalar coefft0  = coefft + coefft00;
-
-        if
         (
-            U.dimensions() == dimVelocity
+            "backwardDdtScheme<Type>::fvcDdtPhiCorr"
-         && phi.dimensions() == dimVelocity*dimArea
+        )   << "dimensions of phiAbs are not correct"
-        )
+            << abort(FatalError);
-        {
-            return tmp<fluxFieldType>
-            (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
-                   *(
-                        coefft0*fvc::interpolate(rA*rho.oldTime())
-                       *phi.oldTime()
-                      - coefft00*fvc::interpolate(rA*rho.oldTime().oldTime())
-                       *phi.oldTime().oldTime()
-                      - (
-                            fvc::interpolate
-                            (
-                                rA*
-                                (
-                                    coefft0*rho.oldTime()*U.oldTime()
-                                  - coefft00*rho.oldTime().oldTime()
-                                   *U.oldTime().oldTime()
-                                )
-                            ) & mesh().Sf()
-                        )
-                    )
-                )
-            );
-        }
-        else if 
-        (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
-            (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    rDeltaT
-                   *fvcDdtPhiCoeff
-                    (
-                        U.oldTime(),
-                        phi.oldTime()/fvc::interpolate(rho.oldTime())
-                    )
-                   *(
-                        fvc::interpolate(rA*rho.oldTime())
-                       *(
-                           coefft0*phi.oldTime()
-                          /fvc::interpolate(rho.oldTime())
-                         - coefft00*phi.oldTime().oldTime()
-                          /fvc::interpolate(rho.oldTime().oldTime())
-                        )
-                      - (
-                            fvc::interpolate
-                            (
-                                rA*rho.oldTime()*
-                                (
-                                    coefft0*U.oldTime()
-                                  - coefft00*U.oldTime().oldTime()
-                                )
-                            ) & mesh().Sf()
-                        )
-                    )
-                )
-            );
-        }
-        else if 
-        (
-            U.dimensions() == rho.dimensions()*dimVelocity
-         && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
-            (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    rDeltaT
-                   *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())
-                   *(
-                        fvc::interpolate(rA)
-                       *(
-                           coefft0*phi.oldTime()
-                         - coefft00*phi.oldTime().oldTime()
-                        )
-                      - (
-                            fvc::interpolate
-                            (
-                                rA*
-                                (
-                                    coefft0*U.oldTime()
-                                  - coefft00*U.oldTime().oldTime()
-                                )
-                            ) & mesh().Sf()
-                        )
-                    )
-                )
-            );
-        }
-        else
-        {
-            FatalErrorIn
-            (
-                "backwardDdtScheme<Type>::fvcDdtPhiCorr"
-            )   << "dimensions of phi are not correct"
-                << abort(FatalError);
 
-            return fluxFieldType::null();
+        return fluxFieldType::null();
-        }
     }
 }