ENH: Adding dynamic mesh functionality to icoReactingMultiphaseInterFoam

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/icoReactingMultiphaseInterFoam.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,16 +42,14 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "dynamicFvMesh.H"
 #include "subCycle.H"
 #include "multiphaseSystem.H"
 #include "turbulentFluidThermoModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 #include "radiationModel.H"
-#include "HashPtrTable.H"
-#include "fvcDDt.H"
-#include "zeroField.H"
+#include "CorrectPhi.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -70,14 +68,18 @@ int main(int argc, char *argv[])
 
     #include "setRootCaseLists.H"
     #include "createTime.H"
-    #include "createMesh.H"
+    #include "createDynamicFvMesh.H"
 
-    pimpleControl pimple(mesh);
+    #include "initContinuityErrs.H"
+    #include "createDyMControls.H"
 
     #include "createFields.H"
     #include "createFieldRefs.H"
+
+    #include "initCorrectPhi.H"
+    #include "createUfIfPresent.H"
+
     #include "createFvOptions.H"
-    #include "createTimeControls.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
@@ -89,7 +91,8 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readTimeControls.H"
+        #include "readDyMControls.H"
+
         #include "CourantNo.H"
         #include "alphaCourantNo.H"
         #include "setDeltaT.H"
@@ -105,6 +108,28 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    gh = (g & mesh.C()) - ghRef;
+                    ghf = (g & mesh.Cf()) - ghRef;
+
+                    if (correctPhi)
+                    {
+                        // Calculate absolute flux
+                        // from the mapped surface velocity
+                        phi = mesh.Sf() & Uf();
+
+                        #include "correctPhi.H"
+
+                        // Make the flux relative to the mesh motion
+                        fvc::makeRelative(phi, U);
+                    }
+                }
+            }
             #include "UEqn.H"
             #include "YEqns.H"
             #include "TEqn.H"

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H

@@ -1,17 +1,31 @@
 {
-    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+    if (correctPhi)
+    {
+        rAU.ref() = 1.0/UEqn.A();
+    }
+    else
+    {
+        rAU = 1.0/UEqn.A();
+    }
+
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
 
     volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn.H();
+    HbyA = rAU()*UEqn.H();
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
+      + fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi)
     );
-    adjustPhi(phiHbyA, U, p_rgh);
+
+    if (p_rgh.needReference())
+    {
+        fvc::makeRelative(phiHbyA, U);
+        adjustPhi(phiHbyA, U, p_rgh);
+        fvc::makeAbsolute(phiHbyA, U);
+    }
 
     surfaceScalarField phig
     (
@@ -50,7 +64,7 @@
 
             p_rgh.relax();
 
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU()*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
 
             U.correctBoundaryConditions();
             fvOptions.correct(U);
@@ -58,6 +72,14 @@
         }
     }
 
+    #include "continuityErrs.H"
+
+    // Correct Uf if the mesh is moving
+    fvc::correctUf(Uf, U, phi);
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+
     p == p_rgh + rho*gh;
 
     if (p_rgh.needReference())
@@ -70,4 +92,9 @@
         );
         p_rgh = p - rho*gh;
     }
+
+    if (!correctPhi)
+    {
+        rAU.clear();
+    }
 }