interDyMFoam: Improved mesh-motion mapping and reinstated ddtPhiCorr

applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C

@@ -58,6 +58,9 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
+    surfaceScalarField phiAbs("phiAbs", phi);
+    fvc::makeAbsolute(phiAbs, U);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
     Info<< "\nStarting time loop\n" << endl;
 
@@ -67,9 +70,6 @@ int main(int argc, char *argv[])
         #include "alphaCourantNo.H"
         #include "CourantNo.H"
 
-        // Make the fluxes absolute
-        fvc::makeAbsolute(phi, U);
-
         #include "setDeltaT.H"
 
         runTime++;
@@ -78,8 +78,18 @@ int main(int argc, char *argv[])
 
         scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
 
-        // Do any mesh changes
+        {
-        mesh.update();
+            // Calculate the relative velocity used to map the relative flux phi
+            volVectorField Urel("Urel", U);
+
+            if (mesh.moving())
+            {
+                Urel -= fvc::reconstruct(fvc::meshPhi(U));
+            }
+
+            // Do any mesh changes
+            mesh.update();
+        }
 
         if (mesh.changing())
         {
@@ -96,9 +106,6 @@ int main(int argc, char *argv[])
             #include "correctPhi.H"
         }
 
-        // Make the fluxes relative to the mesh motion
-        fvc::makeRelative(phi, U);
-
         if (mesh.changing() && checkMeshCourantNo)
         {
             #include "meshCourantNo.H"

applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H

@@ -3,16 +3,19 @@
     surfaceScalarField rAUf(fvc::interpolate(rAU));
 
     U = rAU*UEqn.H();
-    surfaceScalarField phiU("phiU", (fvc::interpolate(U) & mesh.Sf()));
+
+    phiAbs =
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU, rho, U, phiAbs);
 
     if (p_rgh.needReference())
     {
-        fvc::makeRelative(phiU, U);
+        fvc::makeRelative(phiAbs, U);
-        adjustPhi(phiU, U, p_rgh);
+        adjustPhi(phiAbs, U, p_rgh);
-        fvc::makeAbsolute(phiU, U);
+        fvc::makeAbsolute(phiAbs, U);
     }
 
-    phi = phiU +
+    phi = phiAbs +
     (
         fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
       - ghf*fvc::snGrad(rho)
@@ -38,11 +41,13 @@
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+    U += rAU*fvc::reconstruct((phi - phiAbs)/rAUf);
     U.correctBoundaryConditions();
 
     #include "continuityErrs.H"
 
+    phiAbs = phi;
+
     // Make the fluxes relative to the mesh motion
     fvc::makeRelative(phi, U);