ENH: contactAngleForce: add optional contact-line movement force

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();

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_;