Consistency fix

Ratio of old volumes to current volumes needs to be taken into account when
creating time-step consistent face flux in ddt schemes.

src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C

@@ -660,7 +660,35 @@ CoEulerDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 )
 {
-    return (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
+    tmp<fluxFieldType> toldTimeFlux =
+        (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        toldTimeFlux() *= fvc::interpolate(V0ByV);
+    }
+
+    return toldTimeFlux;
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C

@@ -1202,15 +1202,38 @@ CrankNicolsonDdtScheme<Type>::fvcDdtConsistentPhiCorr
         );
     }
     
-    return
+    // Calculate old time flux
-        rAUf*
+    fluxFieldType oldTimeFlux =
+        rAUf*rDtCoef_(faceUDdt0)*(mesh().Sf() & faceU.oldTime());
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
         (
-            mesh().Sf()
+            IOobject
-          & (
+            (
-                rDtCoef_(faceUDdt0)*faceU.oldTime()
+                "V0ByV",
-              + offCentre_(faceUDdt0())
+                mesh().time().timeName(),
-            )
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
         );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        oldTimeFlux *= fvc::interpolate(V0ByV);
+    }
+
+    return
+        oldTimeFlux
+      + rAUf*rDtCoef_(faceUDdt0)*(mesh().Sf() & offCentre_(faceUDdt0()));
 }
 
 

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

@@ -522,7 +522,35 @@ EulerDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 )
 {
-    return (mesh().Sf() & faceU.oldTime())*rAUf/mesh().time().deltaT();
+    tmp<fluxFieldType> toldTimeFlux =
+        (mesh().Sf() & faceU.oldTime())*rAUf/mesh().time().deltaT();
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        toldTimeFlux() *= fvc::interpolate(V0ByV);
+    }
+
+    return toldTimeFlux;
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C

@@ -665,7 +665,35 @@ SLTSDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 )
 {
-    return (mesh().Sf() & faceU.oldTime())*rAUf*fvc::interpolate(SLrDeltaT());
+    tmp<fluxFieldType> toldTimeFlux =
+        (mesh().Sf() & faceU.oldTime())*rAUf*fvc::interpolate(SLrDeltaT());
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        toldTimeFlux() *= fvc::interpolate(V0ByV);
+    }
+
+    return toldTimeFlux;
 }
 
 

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

@@ -732,18 +732,66 @@ backwardDdtScheme<Type>::fvcDdtConsistentPhiCorr
 
     const scalar rDeltaT = 1.0/deltaT;
 
+    // Note: minus sign in gamma coefficient so we can simply add the fluxes
+    // together at the end
     const dimensionedScalar beta("beta", dimless/dimTime, coefft0*rDeltaT);
     const dimensionedScalar gamma("gamma", dimless/dimTime, -coefft00*rDeltaT);
 
-    return
+    // Calculate old and old-old flux contributions
-        rAUf*
+    fluxFieldType oldTimeFlux =
+        beta*rAUf*(mesh().Sf() & faceU.oldTime());
+    fluxFieldType oldOldTimeFlux =
+        gamma*rAUf*(mesh().Sf() & faceU.oldTime().oldTime());
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes for old flux contribution
+        volScalarField V0ByV
         (
-            mesh().Sf()
+            IOobject
-          & (
+            (
-                beta*faceU.oldTime()
+                "V0ByV",
-              + gamma*faceU.oldTime().oldTime()
+                mesh().time().timeName(),
-            )
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
         );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct old time flux contribution
+        oldTimeFlux *= fvc::interpolate(V0ByV);
+
+
+        // Also need to take into account the ratio between old-old and current
+        // cell volumes for old-old time flux contribution
+        volScalarField V00ByV
+        (
+            IOobject
+            (
+                "V00ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V00ByV.internalField() = mesh().V00()/mesh().V();
+        V00ByV.correctBoundaryConditions();
+
+        // Correct old-old time flux contribution
+        oldOldTimeFlux *= fvc::interpolate(V00ByV);
+    }
+
+    return oldTimeFlux + oldOldTimeFlux;
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/steadyInertialDdtScheme/steadyInertialDdtScheme.C

@@ -667,7 +667,35 @@ steadyInertialDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 )
 {
-    return (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
+    tmp<fluxFieldType> toldTimeFlux =
+        (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        toldTimeFlux() *= fvc::interpolate(V0ByV);
+    }
+
+    return toldTimeFlux;
 }