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-2019 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();
150 
151  const tmp<volScalarField> ttheta = theta();
152  const volScalarField& theta = ttheta();
153 
154  const volVectorField gradCoverage(fvc::grad(coverage));
155 
156  forAll(nbr, facei)
157  {
158  const label cellO = own[facei];
159  const label cellN = nbr[facei];
160 
161  label celli = -1;
162  if ((coverage[cellO] > 0.5) && (coverage[cellN] < 0.5))
163  {
164  celli = cellO;
165  }
166  else if ((coverage[cellO] < 0.5) && (coverage[cellN] > 0.5))
167  {
168  celli = cellN;
169  }
170 
171  if (celli != -1 && mask_[celli] > 0.5)
172  {
173  const scalar invDx = filmModel_.regionMesh().deltaCoeffs()[facei];
174  const vector n =
175  gradCoverage[celli]/(mag(gradCoverage[celli]) + rootVSmall);
176  const scalar cosTheta = cos(degToRad(theta[celli]));
177  force[celli] += Ccf_*n*sigma[celli]*(1 - cosTheta)/invDx;
178  }
179  }
180 
181  forAll(coverage.boundaryField(), patchi)
182  {
183  if (!filmModel_.isCoupledPatch(patchi))
184  {
185  const fvPatchField<scalar>& coveragePf =
186  coverage.boundaryField()[patchi];
187  const fvPatchField<scalar>& maskPf = mask_.boundaryField()[patchi];
188  const fvPatchField<scalar>& sigmaPf = sigma.boundaryField()[patchi];
189  const fvPatchField<scalar>& thetaPf = theta.boundaryField()[patchi];
190  const scalarField& invDx = coveragePf.patch().deltaCoeffs();
191  const labelUList& faceCells = coveragePf.patch().faceCells();
192 
193  forAll(coveragePf, facei)
194  {
195  if (maskPf[facei] > 0.5)
196  {
197  const label cellO = faceCells[facei];
198 
199  if ((coverage[cellO] > 0.5) && (coveragePf[facei] < 0.5))
200  {
201  const vector n =
202  gradCoverage[cellO]
203  /(mag(gradCoverage[cellO]) + rootVSmall);
204 
205  const scalar cosTheta = cos(degToRad(thetaPf[facei]));
206  force[cellO] +=
207  Ccf_*n*sigmaPf[facei]*(1 - cosTheta)/invDx[facei];
208  }
209  }
210  }
211  }
212  }
213 
214  force /= V;
215 
217  {
218  tForce().write();
219  }
220 
222  (
223  new fvVectorMatrix(U, dimForce)
224  );
225 
226  tfvm.ref() += tForce;
227 
228  return tfvm;
229 }
230 
231 
232 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
233 
234 } // End namespace surfaceFilmModels
235 } // End namespace regionModels
236 } // End namespace Foam
237 
238 // ************************************************************************* //
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:158
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:174
Unit conversion functions.
const dimensionedScalar & sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
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:622
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 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
virtual tmp< fvVectorMatrix > correct(volVectorField &U)
Correct.
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
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:138
const fvMesh & regionMesh() const
Return the region mesh database.
Definition: regionModelI.H:61
static const zero Zero
Definition: zero.H:97
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:325
dimensionedScalar pos0(const dimensionedScalar &ds)
const dimensionSet dimForce
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:278
label patchi
U
Definition: pEqn.H:72
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
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
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
virtual const volScalarField & sigma() const =0
Return the film surface tension [N/m].
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
defineTypeNameAndDebug(kinematicSingleLayer, 0)
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:812