From 6f2c609a1446d29d81e997b575c97c1e8aa31c67 Mon Sep 17 00:00:00 2001
From: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
Date: Tue, 6 Dec 2022 11:52:55 +0000
Subject: [PATCH] ENH: contactAngleForce: add optional contact-line movement
 force

---
 .../contactAngleForce/contactAngleForce.C     | 18 +++++++++++++++---
 .../contactAngleForce/contactAngleForce.H     | 19 ++++++++++++++++++-
 2 files changed, 33 insertions(+), 4 deletions(-)

diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C
index 32ff65e6f2..ce46bad3d8 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C
@@ -62,6 +62,7 @@ contactAngleForce::contactAngleForce
 :
     force(typeName, film, dict),
     Ccf_(coeffDict_.get<scalar>("Ccf")),
+    hCrit_(coeffDict_.getOrDefault<scalar>("hCrit", GREAT)),
     mask_
     (
         IOobject
@@ -113,15 +114,20 @@ tmp<faVectorMatrix> contactAngleForce::correct(areaVectorField& U)
     const labelUList& nbr = film().regionMesh().neighbour();
 
     const DimensionedField<scalar, areaMesh>& magSf = film().regionMesh().S();
+    const scalarField& magSff = magSf.field();
 
     tmp<areaScalarField> talpha = film().alpha();
     const areaScalarField& sigma = film().sigma();
 
+    const areaScalarField& mu = film().mu();
     const areaScalarField& rhof = film().rho();
 
     tmp<areaScalarField> ttheta = theta();
     const areaScalarField& theta = ttheta();
 
+    const areaVectorField& Uf = film().Uf();
+    const areaScalarField& hf = film().h();
+
     const areaVectorField gradAlpha(fac::grad(talpha()));
 
     forAll(nbr, edgei)
@@ -148,14 +154,20 @@ tmp<faVectorMatrix> contactAngleForce::correct(areaVectorField& U)
             );
             const scalar cosTheta = cos(degToRad(theta[facei]));
 
+            // (MHDX:Eq. 13)
             force[facei] +=
-                Ccf_*n*sigma[facei]*(1 - cosTheta)/invDx/rhof[facei];
+                Ccf_*n*sigma[facei]*(1 - cosTheta)
+               /invDx/rhof[facei]/magSff[facei];
+
+            // (NDPC:Eq. 11)
+            if (hf[facei] > hCrit_)
+            {
+                force[facei] -= mu[facei]*Uf[facei]/hCrit_;
+            }
         }
     }
 
 
-    force /= magSf.field();
-
     if (film().regionMesh().time().writeTime())
     {
         tForce().write();
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.H b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.H
index f99aab9ff2..8209c6d080 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.H
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,6 +32,20 @@ Description
     The effect of the contact angle force can be ignored over a specified
     distance from patches.
 
+    Reference:
+    \verbatim
+        Governing equations (tag:MHDX):
+            Meredith, K. V., Heather, A., De Vries, J., & Xin, Y. (2011).
+            A numerical model for partially-wetted flow of thin liquid films.
+            Computational Methods in Multiphase Flow VI, 70, 239.
+
+        Contact line movement (tag:NDPC):
+            Novák, M., Devaradja, R., Papper, J., & Černý, M. (2020).
+            Efficient CFD methods for assessment of water management.
+            In 20. Internationales Stuttgarter Symposium (pp. 151-170).
+            Springer Vieweg, Wiesbaden.
+    \endverbatim
+
 SourceFiles
     contactAngleForce.C
 
@@ -64,6 +78,9 @@ class contactAngleForce
         //- Coefficient applied to the contact angle force
         scalar Ccf_;
 
+        //- Critical film thickness [m]
+        scalar hCrit_;
+
         //- Mask over which force is applied
         areaScalarField mask_;