ENH: Changes in handling topological changes in VOF solvers

1) Using divU instead of fvc::absolute(phi,U) in TEqn as the latter uses latest time meshPhi which is inconsistent
2) Adding fvc::interpolate(U) when topo changes
3) in pEq for compressible dgdt is updated using the latest rho1 and rho2 after compressible effects are considered

applications/solvers/multiphase/compressibleInterFoam/TEqn.H

@@ -5,8 +5,8 @@
       + 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)
+         + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
         )
        *(
            alpha1/mixture.thermo1().Cv()
@@ -20,4 +20,4 @@
     mixture.correct();
 
     Info<< "min(T) " << min(T).value() << endl;
-}
+}

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

@@ -70,6 +70,21 @@ int main(int argc, char *argv[])
 
     turbulence->validate();
 
+
+    volScalarField rAU
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rAU", dimTime/rho.dimensions(), 1.0)
+    );
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
     Info<< "\nStarting time loop\n" << endl;
 
@@ -77,45 +92,17 @@ int main(int argc, char *argv[])
     {
         #include "readControls.H"
 
-        {
+        // Store divU from the previous mesh so that it can be mapped
-            // Store divU from the previous mesh so that it can be mapped
+        // and used in correctPhi to ensure the corrected phi has the
-            // and used in correctPhi to ensure the corrected phi has the
+        // same divergence
-            // same divergence
+        volScalarField divU("divU0",fvc::div(fvc::absolute(phi, U)));
-            volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
 
-            #include "CourantNo.H"
+        #include "CourantNo.H"
-            #include "setDeltaT.H"
+        #include "setDeltaT.H"
 
-            runTime++;
+        runTime++;
 
-            Info<< "Time = " << runTime.timeName() << nl << endl;
+        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-            scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
-
-            // Do any mesh changes
-            mesh.update();
-
-            if (mesh.changing())
-            {
-                Info<< "Execution time for mesh.update() = "
-                    << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
-                    << " s" << endl;
-
-                gh = (g & mesh.C()) - ghRef;
-                ghf = (g & mesh.Cf()) - ghRef;
-            }
-
-            if (mesh.changing() && correctPhi)
-            {
-                // Calculate absolute flux from the mapped surface velocity
-                phi = mesh.Sf() & Uf;
-
-                #include "correctPhi.H"
-
-                // Make the fluxes relative to the mesh motion
-                fvc::makeRelative(phi, U);
-            }
-        }
 
         if (mesh.changing() && checkMeshCourantNo)
         {
@@ -127,6 +114,46 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            if (pimple.firstIter())
+            {
+                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    Info<< "Execution time for mesh.update() = "
+                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                        << " s" << endl;
+
+                    gh = (g & mesh.C()) - ghRef;
+                    ghf = (g & mesh.Cf()) - ghRef;
+                }
+
+                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);
+
+                    phi = mesh.Sf() & Uf;
+
+                    #include "correctPhi.H"
+
+                    fvc::makeRelative(phi, U);
+
+                    mixture.correct();
+
+                    mesh.topoChanging(false);
+                }
+
+                if (mesh.changing() && checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
+                }
+            }
+
+
             #include "alphaEqnsSubCycle.H"
 
             // correct interface on first PIMPLE corrector
@@ -145,6 +172,11 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         rho = alpha1*rho1 + alpha2*rho2;

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,8 +107,19 @@
 
     rho = alpha1*rho1 + alpha2*rho2;
 
+    p = max(p_rgh + rho*gh, pMin);
+    p_rgh = p - rho*gh;
+
+    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;
     Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
 }
+
+

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -93,6 +93,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

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -84,15 +84,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
@@ -101,14 +92,21 @@
         }
     }
 
-    // p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
-
     // Update densities from change in p_rgh
     rho1 += psi1*(p_rgh - p_rgh_0);
     rho2 += psi2*(p_rgh - p_rgh_0);
 
     rho = alpha1*rho1 + alpha2*rho2;
 
+    p = max(p_rgh + rho*gh, pMin);
+    p_rgh = p - rho*gh;
+
+    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/interFoam/interDyMFoam/interDyMFoam.C

@@ -133,9 +133,12 @@ int main(int argc, char *argv[])
                     ghf = (g & mesh.Cf()) - ghRef;
                 }
 
-                if (mesh.changing() && correctPhi)
+                if ((mesh.changing() && correctPhi) || mesh.topoChanging())
                 {
                     // Calculate absolute flux from the mapped surface velocity
+                    // SAF: temporary fix until mapped Uf is assessed
+                    Uf = fvc::interpolate(U);
+
                     phi = mesh.Sf() & Uf;
 
                     #include "correctPhi.H"
@@ -144,6 +147,8 @@ int main(int argc, char *argv[])
                     fvc::makeRelative(phi, U);
 
                     mixture.correct();
+
+                    mesh.topoChanging(false);
                 }
 
                 if (mesh.changing() && checkMeshCourantNo)