contactAngleForce.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "contactAngleForce.H"
28 #include "fvcGrad.H"
29 #include "unitConversion.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 namespace regionModels
37 {
38 namespace surfaceFilmModels
39 {
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
43 defineTypeNameAndDebug(contactAngleForce, 0);
44 
45 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
46 
47 void contactAngleForce::initialise()
48 {
49  const wordReList zeroForcePatches
50  (
51  coeffDict_.lookupOrDefault<wordReList>("zeroForcePatches", wordReList())
52  );
53 
54  if (zeroForcePatches.size())
55  {
56  const polyBoundaryMesh& pbm = filmModel_.regionMesh().boundaryMesh();
57  scalar dLim = coeffDict_.lookup<scalar>("zeroForceDistance");
58 
59  Info<< " Assigning zero contact force within " << dLim
60  << " of patches:" << endl;
61 
62  labelHashSet patchIDs = pbm.patchSet(zeroForcePatches);
63 
64  forAllConstIter(labelHashSet, patchIDs, iter)
65  {
66  label patchi = iter.key();
67  Info<< " " << pbm[patchi].name() << endl;
68  }
69 
70  // Temporary implementation until run-time selection covers this case
71  patchDistMethods::meshWave dist(filmModel_.regionMesh(), patchIDs);
73  (
74  IOobject
75  (
76  "y",
81  false
82  ),
85  );
86  dist.correct(y);
87 
88  mask_ = pos0(y - dimensionedScalar(dimLength, dLim));
89  }
90 }
91 
92 
93 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
94 
96 (
97  const word& typeName,
99  const dictionary& dict
100 )
101 :
102  force(typeName, film, dict),
103  Ccf_(coeffDict_.lookup<scalar>("Ccf")),
104  mask_
105  (
106  IOobject
107  (
108  IOobject::modelName("contactForceMask", typeName),
113  ),
116  )
117 {
118  initialise();
119 }
120 
121 
122 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
123 
125 {}
126 
127 
128 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
129 
131 {
132  tmp<volVectorField> tForce
133  (
135  (
136  IOobject::modelName("contactForce", typeName),
139  )
140  );
141 
142  vectorField& force = tForce.ref().primitiveFieldRef();
143 
144  const labelUList& own = filmModel_.regionMesh().owner();
145  const labelUList& nbr = filmModel_.regionMesh().neighbour();
146  const scalarField& V = filmModel_.regionMesh().V();
147 
148  const volScalarField& coverage = filmModel_.coverage();
149 
150  const tmp<volScalarField> tsigma = filmModel_.sigma();
151  const volScalarField& sigma = tsigma();
152 
153  const tmp<volScalarField> ttheta = theta();
154  const volScalarField& theta = ttheta();
155 
156  const volVectorField gradCoverage(fvc::grad(coverage));
157 
158  forAll(nbr, facei)
159  {
160  const label cellO = own[facei];
161  const label cellN = nbr[facei];
162 
163  label celli = -1;
164  if ((coverage[cellO] > 0.5) && (coverage[cellN] < 0.5))
165  {
166  celli = cellO;
167  }
168  else if ((coverage[cellO] < 0.5) && (coverage[cellN] > 0.5))
169  {
170  celli = cellN;
171  }
172 
173  if (celli != -1 && mask_[celli] > 0.5)
174  {
175  const scalar invDx = filmModel_.regionMesh().deltaCoeffs()[facei];
176  const vector n =
177  gradCoverage[celli]/(mag(gradCoverage[celli]) + rootVSmall);
178  const scalar cosTheta = cos(degToRad(theta[celli]));
179  force[celli] += Ccf_*n*sigma[celli]*(1 - cosTheta)/invDx;
180  }
181  }
182 
183  forAll(coverage.boundaryField(), patchi)
184  {
185  if (!filmModel_.isCoupledPatch(patchi))
186  {
187  const fvPatchField<scalar>& coveragePf =
188  coverage.boundaryField()[patchi];
189  const fvPatchField<scalar>& maskPf = mask_.boundaryField()[patchi];
190  const fvPatchField<scalar>& sigmaPf = sigma.boundaryField()[patchi];
191  const fvPatchField<scalar>& thetaPf = theta.boundaryField()[patchi];
192  const scalarField& invDx = coveragePf.patch().deltaCoeffs();
193  const labelUList& faceCells = coveragePf.patch().faceCells();
194 
195  forAll(coveragePf, facei)
196  {
197  if (maskPf[facei] > 0.5)
198  {
199  const label cellO = faceCells[facei];
200 
201  if ((coverage[cellO] > 0.5) && (coveragePf[facei] < 0.5))
202  {
203  const vector n =
204  gradCoverage[cellO]
205  /(mag(gradCoverage[cellO]) + rootVSmall);
206 
207  const scalar cosTheta = cos(degToRad(thetaPf[facei]));
208  force[cellO] +=
209  Ccf_*n*sigmaPf[facei]*(1 - cosTheta)/invDx[facei];
210  }
211  }
212  }
213  }
214  }
215 
216  force /= V;
217 
219  {
220  tForce().write();
221  }
222 
224  (
225  new fvVectorMatrix(U, dimForce)
226  );
227 
228  tfvm.ref() += tForce;
229 
230  return tfvm;
231 }
232 
233 
234 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
235 
236 } // End namespace surfaceFilmModels
237 } // End namespace regionModels
238 } // End namespace Foam
239 
240 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
virtual tmp< volScalarField > theta() const =0
Return the contact angle field.
Base class for film (stress-based) force models.
Definition: force.H:55
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
const Boundary & boundaryField() const
Return const-reference to the boundary field.
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
Unit conversion functions.
static tmp< GeometricField< vector, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< vector >> &)
Return a temporary field constructed from name,.
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:66
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:636
const dimensionSet dimless
const Time & time() const
Return the reference to the time database.
Definition: regionModelI.H:37
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
const dictionary coeffDict_
Coefficients dictionary.
Definition: subModelBase.H:74
Macros for easy insertion into run-time selection tables.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
const dimensionSet dimLength
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:284
bool writeTime() const
Return true if this is a write time.
Definition: TimeStateI.H:65
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:211
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
virtual tmp< fvVectorMatrix > correct(volVectorField &U)
Correct.
scalar y
contactAngleForce(const word &typeName, surfaceFilmRegionModel &film, const dictionary &dict)
Construct from surface film model.
dimensionedScalar cos(const dimensionedScalar &ds)
Calculate the gradient of the given field.
A class for handling words, derived from string.
Definition: word.H:59
bool isCoupledPatch(const label regionPatchi) const
Return true if patchi on the local region is a coupled.
Definition: regionModelI.H:132
const fvMesh & regionMesh() const
Return the region mesh database.
Definition: regionModelI.H:55
static const zero Zero
Definition: zero.H:97
const dimensionSet dimForce
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
const fvPatch & patch() const
Return patch.
Definition: fvPatchField.H:337
dimensionedScalar pos0(const dimensionedScalar &ds)
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:278
label patchi
U
Definition: pEqn.H:72
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
force(surfaceFilmRegionModel &film)
Construct null.
Definition: force.C:44
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
fvMatrix< vector > fvVectorMatrix
Definition: fvMatricesFwd.H:45
const dimensionSet dimVolume
messageStream Info
surfaceFilmRegionModel & filmModel_
Reference to the film surface film model.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
virtual const volScalarField & coverage() const =0
Return the film coverage, 1 = covered, 0 = uncovered / [].
static word modelName(Name name, const word &model)
Return the name of the object within the given model.
List< wordRe > wordReList
A List of wordRe (word or regular expression)
Definition: wordReList.H:50
A class for managing temporary objects.
Definition: PtrList.H:53
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
defineTypeNameAndDebug(kinematicSingleLayer, 0)
virtual tmp< volScalarField > sigma() const =0
Return the film surface tension [N/m].
const surfaceScalarField & deltaCoeffs() const
Return reference to cell-centre difference coefficients.
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844