interPhaseChangeDyMFoam: Corrected correctPhi to include the phase-change dilatation effect

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

@@ -37,14 +37,17 @@ if (mesh.changing())
         pcorrTypes
     );
 
-    // dimensionedScalar rAUf("rAUf", dimTime/rho.dimensions(), 1.0);
     surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
+    #ifndef divUCorr
+    #define divUCorr
+    #endif
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi) divUCorr
         );
 
         pcorrEqn.setReference(pRefCell, pRefValue);
@@ -56,5 +59,7 @@ if (mesh.changing())
         }
     }
 
+    #undef divUCorr
+
     #include "continuityErrs.H"
 }

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

@@ -82,7 +82,16 @@
                 0
             );
 
-            tphiAlpha() += tphiAlphaCorr();
+            // Under-relax the correction for all but the 1st corrector
+            if (aCorr == 0)
+            {
+                tphiAlpha() += tphiAlphaCorr();
+            }
+            else
+            {
+                alpha1 = 0.5*alpha1 + 0.5*alpha10;
+                tphiAlpha() += 0.5*tphiAlphaCorr();
+            }
         }
         else
         {

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C

@@ -82,7 +82,6 @@ int main(int argc, char *argv[])
         dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0)
     );
 
-    #include "../interFoam/interDyMFoam/correctPhi.H"
     #include "createUf.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
@@ -106,6 +105,9 @@ int main(int argc, char *argv[])
         {
             if (pimple.firstIter() || moveMeshOuterCorrectors)
             {
+                // Store divU from the previous mesh for the correctPhi
+                volScalarField divU(fvc::div(fvc::absolute(phi, U)));
+
                 scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
 
                 mesh.update();
@@ -125,6 +127,7 @@ int main(int argc, char *argv[])
                     // Calculate absolute flux from the mapped surface velocity
                     phi = mesh.Sf() & Uf;
 
+                    #define divUCorr -divU
                     #include "../interFoam/interDyMFoam/correctPhi.H"
 
                     // Make the flux relative to the mesh motion

tutorials/multiphase/interPhaseChangeDyMFoam/propeller/system/fvSchemes

@@ -32,7 +32,6 @@ divSchemes
 
     div(phi,alpha)      Gauss vanLeer;
     div(phirb,alpha)    Gauss linear;
-    UD                  Gauss upwind;
     div(rhoPhi,U)       Gauss linearUpwind grad(U);
     div(phi,k)          Gauss upwind;
     div(phi,epsilon)    Gauss upwind;