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 8a8edd0f8a2bf1978541ccc3338216fecbf5bdd1..5dc0e4b8d10353bab81ff6e4ca3326c4db7b8360 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C
@@ -44,13 +44,6 @@ namespace areaSurfaceFilmModels
defineTypeNameAndDebug(contactAngleForce, 0);
-// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
-
-void contactAngleForce::initialise()
-{
-}
-
-
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
contactAngleForce::contactAngleForce
@@ -62,28 +55,7 @@ contactAngleForce::contactAngleForce
:
force(typeName, film, dict),
Ccf_(coeffDict_.get<scalar>("Ccf")),
- hCrit_(coeffDict_.getOrDefault<scalar>("hCrit", GREAT)),
- mask_
- (
- IOobject
- (
- typeName + ":fContactForceMask",
- film.primaryMesh().time().timeName(),
- film.primaryMesh(),
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- film.regionMesh(),
- dimensionedScalar("mask", dimless, 1.0)
- )
-{
- initialise();
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-contactAngleForce::~contactAngleForce()
+ hCrit_(coeffDict_.getOrDefault<scalar>("hCrit", GREAT))
{}
@@ -91,73 +63,69 @@ contactAngleForce::~contactAngleForce()
tmp<faVectorMatrix> contactAngleForce::correct(areaVectorField& U)
{
- tmp<areaVectorField> tForce
+ auto tforce = tmp<areaVectorField>::New
(
- new areaVectorField
+ IOobject
(
- IOobject
- (
- typeName + ":contactForce",
- film().primaryMesh().time().timeName(),
- film().primaryMesh(),
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- film().regionMesh(),
- dimensionedVector(dimForce/dimDensity/dimArea, Zero)
- )
+ IOobject::scopedName(typeName, "contactForce"),
+ film().primaryMesh().time().timeName(),
+ film().primaryMesh(),
+ IOobject::NO_READ,
+ IOobject::NO_WRITE,
+ IOobject::NO_REGISTER
+ ),
+ film().regionMesh(),
+ dimensionedVector(dimForce/dimDensity/dimArea, Zero)
);
-
- vectorField& force = tForce.ref().primitiveFieldRef();
+ vectorField& force = tforce.ref().primitiveFieldRef();
const labelUList& own = film().regionMesh().owner();
- const labelUList& nbr = film().regionMesh().neighbour();
+ const labelUList& nei = 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 edgeScalarField& invDx = film().regionMesh().deltaCoeffs();
+ const areaScalarField& sigma = film().sigma();
const areaScalarField& mu = film().mu();
const areaScalarField& rhof = film().rho();
+ const areaVectorField& Uf = film().Uf();
+ const areaScalarField& hf = film().h();
+
+ tmp<areaScalarField> talpha = film().alpha();
+ const areaScalarField& alpha = talpha();
tmp<areaScalarField> ttheta = theta();
const areaScalarField& theta = ttheta();
- const areaVectorField& Uf = film().Uf();
- const areaScalarField& hf = film().h();
+ tmp<areaVectorField> tgradAlpha = fac::grad(alpha);
+ const areaVectorField& gradAlpha = tgradAlpha();
- const areaVectorField gradAlpha(fac::grad(talpha()));
-
- forAll(nbr, edgei)
+ forAll(nei, edgei)
{
const label faceO = own[edgei];
- const label faceN = nbr[edgei];
+ const label faceN = nei[edgei];
label facei = -1;
- if ((talpha()[faceO] > 0.5) && (talpha()[faceN] < 0.5))
+ if ((alpha[faceO] > 0.5) && (alpha[faceN] < 0.5))
{
facei = faceO;
}
- else if ((talpha()[faceO] < 0.5) && (talpha()[faceN] > 0.5))
+ else if ((alpha[faceO] < 0.5) && (alpha[faceN] > 0.5))
{
facei = faceN;
}
- if (facei != -1 && mask_[facei] > 0.5)
+ if (facei != -1)
{
- const scalar invDx = film().regionMesh().deltaCoeffs()[edgei];
- const vector n
- (
- gradAlpha[facei]/(mag(gradAlpha[facei]) + ROOTVSMALL)
- );
+ const vector n(normalised(gradAlpha[facei]));
const scalar cosTheta = cos(degToRad(theta[facei]));
// (MHDX:Eq. 13)
force[facei] +=
Ccf_*n*sigma[facei]*(1 - cosTheta)
- /invDx/rhof[facei]/magSff[facei];
+ /invDx[edgei]/rhof[facei]/magSff[facei];
// (NDPC:Eq. 11)
if (hf[facei] > hCrit_)
@@ -185,16 +153,13 @@ tmp<faVectorMatrix> contactAngleForce::correct(areaVectorField& U)
if (film().regionMesh().time().writeTime())
{
- tForce().write();
+ tforce().write();
gradAlpha.write();
}
- tmp<faVectorMatrix> tfvm
- (
- new faVectorMatrix(U, dimForce/dimDensity)
- );
+ auto tfvm = tmp<faVectorMatrix>::New(U, dimForce/dimDensity);
- tfvm.ref() += tForce;
+ tfvm.ref() += tforce;
return tfvm;
}
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 8209c6d08085ea72f0a1b15858afcb5e2b56163c..d763d37e98d339b430a90520a069e65d1f06f900 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.H
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.H
@@ -32,6 +32,36 @@ Description
The effect of the contact angle force can be ignored over a specified
distance from patches.
+ The contact-angle force is implemented as follows:
+ \f[
+ \vec{f}_\theta =
+ \beta \frac{\sigma (1 - cos(\theta))}{\Delta_{cl}} \vec{n}_{cl}
+ \frac{1}{\rho_f|\vec{S}_f|}
+ - \vec{f}_{cl}
+ \f]
+
+ with if \f$ h > h_{crit} \f$:
+
+ \f[
+ \vec{f}_{cl} = \nu_f \frac{\vec{U}_f}{h_{crit}}
+ \f]
+
+ where
+ \vartable
+ \vec{f}_\theta | Contact-angle force [N/m^2/rho]
+ \beta | Empirical constant [-]
+ \sigma | Liquid-air interfacial surface tension [kg/s^2]
+ \theta | Contact angle [rad]
+ \Delta_{cl} | Inverse width in direction normal to contact line [1/m]
+ \vec{n}_{cl} | Direction normal to contact line [-]
+ \rho_f | Film density [kg/m^3]
+ \vec{S}_f | Face-area vector [m^2]
+ \vec{f}_{cl} | Contact-line movement force [N/m^2/rho]
+ \nu_f | Film kinematic viscosity [m^2/s]
+ \vec{U}_f | Film velocity [m/s]
+ h_{crit} | Critical film height [m]
+ \endvartable
+
Reference:
\verbatim
Governing equations (tag:MHDX):
@@ -46,13 +76,38 @@ Description
Springer Vieweg, Wiesbaden.
\endverbatim
+Usage
+ Minimal example:
+ \verbatim
+ {
+ // Mandatory entries
+ Ccf <scalar>;
+
+ // Optional entries
+ hCrit <scalar>;
+
+ // Inherited entries
+ ...
+ }
+ \endverbatim
+
+ where the entries mean:
+ \table
+ Property | Description | Type | Reqd | Deflt
+ Ccf | Empirical coefficient | scalar | yes | -
+ hCrit | Critical film height [m] | scalar | no | GREAT
+ \endtable
+
+ The inherited entries are elaborated in:
+ - \link force.H \endlink
+
SourceFiles
contactAngleForce.C
\*---------------------------------------------------------------------------*/
-#ifndef contactAngleForce_H
-#define contactAngleForce_H
+#ifndef areaSurfaceFilmModels_contactAngleForce_H
+#define areaSurfaceFilmModels_contactAngleForce_H
#include "force.H"
@@ -81,15 +136,9 @@ class contactAngleForce
//- Critical film thickness [m]
scalar hCrit_;
- //- Mask over which force is applied
- areaScalarField mask_;
-
// Private Member Functions
- //- Initialise
- void initialise();
-
//- No copy construct
contactAngleForce(const contactAngleForce&) = delete;
@@ -121,7 +170,7 @@ public:
//- Destructor
- virtual ~contactAngleForce();
+ virtual ~contactAngleForce() = default;
// Member Functions