ENH: compressibleInterDyMFoam enhancements for mesh motion and sphere drop test case tutorial

applications/solvers/multiphase/compressibleInterFoam/TEqn.H

@@ -5,7 +5,7 @@
       + fvm::div(rhoPhi, T)
       - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
       + (
-            fvc::div(fvc::absolute(phi, U), p)
+            divU*p
           + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
         )
        *(

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C

@@ -82,12 +82,12 @@ int main(int argc, char *argv[])
     {
         #include "readControls.H"
 
-        {
-            // Store divU from the previous mesh so that it can be mapped
-            // and used in correctPhi to ensure the corrected phi has the
-            // same divergence
-            volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
+        // Store divU from the previous mesh so that it can be mapped
+        // and used in correctPhi to ensure the corrected phi has the
+        // same divergence
+        volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
 
+        {
             #include "CourantNo.H"
             #include "setDeltaT.H"
 
@@ -110,8 +110,12 @@ int main(int argc, char *argv[])
                 ghf = (g & mesh.Cf()) - ghRef;
             }
 
-            if (mesh.changing() && correctPhi)
+            if ((correctPhi && mesh.changing()) || mesh.topoChanging())
             {
+                // Calculate absolute flux from the mapped surface velocity
+                // SAF: temporary fix until mapped Uf is assessed
+                Uf = fvc::interpolate(U);
+
                 // Calculate absolute flux from the mapped surface velocity
                 phi = mesh.Sf() & Uf;
 
@@ -119,6 +123,8 @@ int main(int argc, char *argv[])
 
                 // Make the fluxes relative to the mesh motion
                 fvc::makeRelative(phi, U);
+
+                mesh.topoChanging(false);
             }
         }
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H

@@ -87,15 +87,6 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
-            p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
-
-            dgdt =
-            (
-                pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
-              - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
-            );
-
             phi = phiHbyA + p_rghEqnIncomp.flux();
 
             U = HbyA
@@ -116,10 +107,16 @@
 
     rho = alpha1*rho1 + alpha2*rho2;
 
-    // Correct p_rgh for consistency with p and the updated densities
+    p = max(p_rgh + rho*gh, pMin);
     p_rgh = p - rho*gh;
     p_rgh.correctBoundaryConditions();
 
+    dgdt =
+    (
+        pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
+      - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
+    );
+
     K = 0.5*magSqr(U);
 
     Info<< "max(U) " << max(mag(U)).value() << endl;

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -98,6 +98,7 @@ int main(int argc, char *argv[])
             solve(fvm::ddt(rho) + fvc::div(rhoPhi));
 
             #include "UEqn.H"
+            volScalarField divU(fvc::div(fvc::absolute(phi, U)));
             #include "TEqn.H"
 
             // --- Pressure corrector loop