From f6688bb9b31ee92c71a2dcdd83b99b55aa195d9d Mon Sep 17 00:00:00 2001
From: Henry Weller <http://openfoam.org>
Date: Thu, 19 Dec 2019 12:27:41 +0000
Subject: [PATCH] surfaceFilmModels::kinematicSingleLayer: Added support for
 density gradient effects on the hydrostatic pressure

avoiding gradients in the extruded mesh geometry which causes serious errors at
corners.
---
 .../kinematicSingleLayer.C                    | 32 +++++++++++++++----
 .../kinematicSingleLayer.H                    |  5 ++-
 2 files changed, 30 insertions(+), 7 deletions(-)

diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
index e24b0c86b5..e5d0ddd6d0 100644
--- 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,7 +372,7 @@ 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
     (
@@ -363,7 +380,10 @@ void kinematicSingleLayer::solveAlpha
         (
             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
index f064f174e2..cd0ce92bc3 100644
--- 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();