Performance update for free surface solver with immersed boundary

applications/solvers/immersedBoundary/interIbFoam/UEqn.H

@@ -11,7 +11,6 @@
       + fvm::div(rhoPhi, U)
       - fvm::laplacian(muEff, U)
       - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();
@@ -22,13 +21,7 @@
         (
             UEqn
          ==
-            fvc::reconstruct
-            (
-                (
-                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
-                ) * mesh.magSf()
-            )
+          - gh*fvc::grad(rho)
+          - fvc::grad(pd)
         );
     }

applications/solvers/immersedBoundary/interIbFoam/alphaEqn.H

@@ -2,9 +2,22 @@
     word alphaScheme("div(phi,alpha)");
     word alpharScheme("div(phirb,alpha)");
 
-    surfaceScalarField phic = mag(phi/mesh.magSf());
-    phic = min(interface.cAlpha()*phic, max(phic));
-    surfaceScalarField phir = faceIbMask*phic*interface.nHatf();
+    // New formulation of phir: Co condition
+    surfaceScalarField phir
+    (
+        "phir",
+        faceIbMask*interface.cAlpha()*interface.nHatf()*
+        min
+        (
+            scalar(0.5)/    // This is ref compression Co number for cAlpha = 1
+            (
+                mesh.surfaceInterpolation::deltaCoeffs()*
+                mesh.time().deltaT()
+            ),
+            mag(phi/mesh.magSf())
+          + dimensionedScalar("small", dimVelocity, 1e-3)
+        )
+    );
 
     fvScalarMatrix alpha1Eqn
     (
@@ -18,6 +31,7 @@
         )
     );
 
+    alpha1Eqn.relax();
     alpha1Eqn.solve();
 
     Info<< "Liquid phase volume fraction = "
@@ -26,12 +40,17 @@
         << "  Max(alpha1) = " << max(cellIbMask*alpha1).value()
         << endl;
 
+    rhoPhi = faceIbMask*(alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2);
+
     // Limit alpha to handle IB cells
-    alpha1.max(0);
-    alpha1.min(1);
+//     alpha1.max(0);
+//     alpha1.min(1);
 
     // Update of interface and density
     interface.correct();
+    twoPhaseProperties.correct();
 
-    rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;
+    // Calculate density using limited alpha1
+    rho = twoPhaseProperties.rho();
+    rho.correctBoundaryConditions();
 }

applications/solvers/immersedBoundary/interIbFoam/createFields.H

@@ -105,21 +105,8 @@
     scalar pdRefValue = 0.0;
     setRefCell(pd, pimple.dict(), pdRefCell, pdRefValue);
 
-    scalar pRefValue = 0.0;
-
-    if (pd.needReference())
-    {
-        pRefValue = readScalar(pimple.dict().lookup("pRefValue"));
-
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pdRefCell)
-        );
-    }
-
     mesh.schemesDict().setFluxRequired(pd.name());
+    mesh.schemesDict().setFluxRequired(alpha1.name());
 
     // Construct interface from alpha1 distribution
     interfaceProperties interface(alpha1, U, twoPhaseProperties);

applications/solvers/immersedBoundary/interIbFoam/interIbFoam.C

@@ -88,8 +88,6 @@ int main(int argc, char *argv[])
         {
             twoPhaseProperties.correct();
 
-#           include "alphaEqn.H"
-
 #           include "UEqn.H"
 
             // --- PISO loop
@@ -102,20 +100,10 @@ int main(int argc, char *argv[])
 
 #           include "limitU.H"
 
-            // Recalculate the mass fluxes
-            rhoPhi = phi*fvc::interpolate(rho);
+            p = pd + rho*gh;
+            p.correctBoundaryConditions();
 
-            p = pd + cellIbMask*rho*gh;
-
-            if (pd.needReference())
-            {
-                p += dimensionedScalar
-                (
-                    "p",
-                    p.dimensions(),
-                    pRefValue - getRefCellValue(p, pdRefCell)
-                );
-            }
+#           include "alphaEqn.H"
 
             turbulence->correct();
         }

applications/solvers/immersedBoundary/interIbFoam/pEqn.H

@@ -5,6 +5,7 @@
     U = rUA*UEqn.H();
     // Immersed boundary update
     U.correctBoundaryConditions();
+#   include "limitU.H"
 
     surfaceScalarField phiU
     (
@@ -12,10 +13,6 @@
         faceIbMask*(fvc::interpolate(U) & mesh.Sf())
     );
 
-    // Adjust immersed boundary fluxes
-    immersedBoundaryAdjustPhi(phiU, U);
-    adjustPhi(phiU, U, pd);
-
     phi = phiU
       + faceIbMask*
         (
@@ -23,12 +20,16 @@
           - ghf*fvc::snGrad(rho)
         )*rUAf*mesh.magSf();
 
+    // Adjust immersed boundary fluxes
+    immersedBoundaryAdjustPhi(phi, U);
+    adjustPhi(phi, U, pd);
+
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pdEqn
         (
-            fvm::laplacian(rUAf, pd) == fvc::div(phi)
+            fvm::laplacian(rUAf, pd, "laplacian(rAU,pd)") == fvc::div(phi)
         );
 
         pdEqn.setReference(pdRefCell, pdRefValue);
@@ -38,18 +39,14 @@
             mesh.solutionDict().solver(pd.select(pimple.finalInnerIter()))
         );
 
-
         if (pimple.finalNonOrthogonalIter())
         {
             phi -= pdEqn.flux();
         }
     }
 
-    // Explicitly relax pressure except for last corrector
-    if (!pimple.finalIter())
-    {
+    // Explicitly relax pressure
     pd.relax();
-    }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();