surfaceFilmModels::kinematicSingleLayer: Added support for density gr…

…adient effects on the hydrostatic pressure avoiding gradients in the extruded mesh geometry which causes serious errors at corners.
OpenFOAM · Dec 19, 2019 · f6688bb · f6688bb
1 parent 68a3840
commit f6688bb
Show file tree

Hide file tree

Showing 2 changed files with 30 additions and 7 deletions.
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
@@ -176,9 +176,23 @@ tmp<volScalarField> kinematicSingleLayer::pe()
 }
 
 
-tmp<volScalarField> kinematicSingleLayer::rhog() const
+tmp<surfaceScalarField> kinematicSingleLayer::rhog() const
 {
-    return -rho_*min(nHat() & g_, dimensionedScalar(g_.dimensions(), 0))*VbyA();
+    return
+        fvc::interpolate
+        (
+            max(nHat() & -g_, dimensionedScalar(g_.dimensions(), 0))*VbyA()
+        )*fvc::interpolate(rho_);
+}
+
+
+tmp<surfaceScalarField> kinematicSingleLayer::gGradRho() const
+{
+    return
+        fvc::interpolate
+        (
+            max(nHat() & -g_, dimensionedScalar(g_.dimensions(), 0))*VbyA()
+        )*fvc::snGrad(rho_);
 }
 
 
@@ -311,6 +325,8 @@ tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
 
     if (pimple_.momentumPredictor())
     {
+        const surfaceScalarField alphaf(fvc::interpolate(alpha_));
+
         solve
         (
             UEqn
@@ -319,11 +335,12 @@ tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
             (
                 constrainFilmField
                 (
-                    fvc::interpolate(alpha_)
+                    alphaf
                    *(
                         (
                             fvc::snGrad(pe + pc(), "snGrad(p)")
-                          + fvc::interpolate(rhog())*fvc::snGrad(alpha_)
+                          + gGradRho()*alphaf
+                          + rhog()*fvc::snGrad(alpha_)
                         )*regionMesh().magSf()
                       - fvc::interpolate(rho_)*(g_ & regionMesh().Sf())
                     ), 0
@@ -355,15 +372,18 @@ void kinematicSingleLayer::solveAlpha
     const surfaceScalarField alphaf(fvc::interpolate(alpha_));
     const surfaceScalarField rhof(fvc::interpolate(rho_));
     const surfaceScalarField alpharAUf(fvc::interpolate(alpha_*rAU));
-    const surfaceScalarField rhogf(fvc::interpolate(rhog()));
+    const surfaceScalarField rhogf(rhog());
 
     const surfaceScalarField phiu
     (
         "phiu",
         (
             constrainFilmField
             (
-                fvc::snGrad(pe + pc(), "snGrad(p)")*regionMesh().magSf()
+                (
+                    fvc::snGrad(pe + pc(), "snGrad(p)")
+                  + gGradRho()*alphaf
+                )*regionMesh().magSf()
               - rhof*(g_ & regionMesh().Sf()),
                 0
             )

diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H
@@ -232,7 +232,10 @@ protected:
         virtual void transferPrimaryRegionSourceFields();
 
         //- Hydrostatic pressure coefficient
-        tmp<volScalarField> rhog() const;
+        tmp<surfaceScalarField> rhog() const;
+
+        //- Hydrostatic pressure coefficient gradient
+        tmp<surfaceScalarField> gGradRho() const;
 
         //- Capillary pressure
         tmp<volScalarField> pc();