From f9c325021f35df55c1b2f0c165e4574b386f5d56 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 11 Oct 2018 16:47:58 +0100
Subject: [PATCH] Simplification of the pimple IB solver: mesh update and flux
 correction

---
 .../pimpleDyMIbFoam/correctMeshMotion.H       | 45 +++++++++++++++++++
 .../pimpleDyMIbFoam/correctPhi.H              |  9 ++--
 .../pimpleDyMIbFoam/pimpleDyMIbFoam.C         | 24 +---------
 3 files changed, 53 insertions(+), 25 deletions(-)
 create mode 100644 applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctMeshMotion.H

diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctMeshMotion.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctMeshMotion.H
new file mode 100644
index 000000000..8534224c8
--- /dev/null
+++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctMeshMotion.H
@@ -0,0 +1,45 @@
+{
+    // Make the fluxes absolute (using the ddt(U) scheme)
+    fvc::makeAbsolute(phi, U);
+
+    // Motion inside PIMPLE loop: simple IB update
+    bool meshChanged = mesh.update();
+
+#   include "updateIbPatchFields.H"
+#   include "updateIbMasks.H"
+#   include "volContinuity.H"
+
+    if (runTime.outputTime())
+    {
+        volScalarField divMeshPhi
+        (
+            "divMeshPhi",
+            -fvc::div(mesh.phi())
+        );
+
+        divMeshPhi.internalField() +=
+            (1.0 - mesh.V0()/mesh.V())/runTime.deltaT().value();
+
+        divMeshPhi = mag(divMeshPhi);
+        divMeshPhi.write();
+
+        U.write();
+        mesh.phi().write();
+    }
+
+    // Fluxes will be corrected to absolute velocity
+    // HJ, 6/Feb/2009
+#   include "correctPhi.H"
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+
+    if (checkMeshCourantNo)
+    {
+#       include "immersedBoundaryCourantNo.H"
+    }
+
+    // Correct velocity boundary conditions after mesh motion, with relative
+    // fluxes
+    U.correctBoundaryConditions();
+}
diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H
index eda2e9eda..814e31803 100644
--- a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H
+++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H
@@ -1,5 +1,6 @@
 {
     volScalarField pcorr("pcorr", p);
+    pcorr *= 0;
 
     // Initialise flux with interpolated velocity
     phi = fvc::interpolate(U) & mesh.Sf();
@@ -8,11 +9,13 @@
 
     mesh.schemesDict().setFluxRequired(pcorr.name());
 
+    dimensionedScalar rUAf("(1|A(U))", dimTime, 1.0);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(1/aU, pcorr) == fvc::div(phi)
+            fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
         );
 
         pcorrEqn.setReference(pRefCell, pRefValue);
@@ -26,6 +29,6 @@
         // Fluxes are corrected to absolute velocity and further corrected
         // later.  HJ, 6/Feb/2009
     }
-#   include "continuityErrs.H"
-}
 
+    #   include "continuityErrs.H"
+}
diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C
index a1a4464d0..a4fc1f0ca 100644
--- a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C
+++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C
@@ -72,34 +72,14 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
 #       include "readControls.H"
-#       include "CourantNo.H"
+#       include "immersedBoundaryCourantNo.H"
 #       include "setDeltaT.H"
 
-        // Make the fluxes absolute
-        fvc::makeAbsolute(phi, U);
-
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        bool meshChanged = mesh.update();
-        reduce(meshChanged, orOp<bool>());
-
-#       include "updateIbPatchFields.H"
-#       include "updateIbMasks.H"
-#       include "volContinuity.H"
-
-        if (checkMeshCourantNo)
-        {
-#           include "meshCourantNo.H"
-        }
-
-        // Fluxes will be corrected to absolute velocity
-        // HJ, 6/Feb/2009
-#       include "correctPhi.H"
-
-        // Make the fluxes relative to the mesh motion
-        fvc::makeRelative(phi, U);
+#       include "correctMeshMotion.H"
 
         // --- PIMPLE loop
         while (pimple.loop())