contactAngleForce.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2017 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 = readScalar(coeffDict_.lookup("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  ),
84  dimensionedScalar("y", dimLength, GREAT)
85  );
86  dist.correct(y);
87 
88  mask_ = pos0(y - dimensionedScalar("dLim", dimLength, dLim));
89  }
90 }
91 
92 
93 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
94 
95 contactAngleForce::contactAngleForce
96 (
97  const word& typeName,
98  surfaceFilmModel& film,
99  const dictionary& dict
100 )
101 :
102  force(typeName, film, dict),
103  Ccf_(readScalar(coeffDict_.lookup("Ccf"))),
104  mask_
105  (
106  IOobject
107  (
108  typeName + ":contactForceMask",
113  ),
115  dimensionedScalar("mask", dimless, 1.0)
116  )
117 {
118  initialise();
119 }
120 
121 
122 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
123 
125 {}
126 
127 
128 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
129 
131 {
132  tmp<volVectorField> tForce
133  (
134  new volVectorField
135  (
136  IOobject
137  (
138  typeName + ":contactForce",
143  ),
146  )
147  );
148 
149  vectorField& force = tForce.ref().primitiveFieldRef();
150 
151  const labelUList& own = filmModel_.regionMesh().owner();
152  const labelUList& nbr = filmModel_.regionMesh().neighbour();
153 
154  const scalarField& magSf = filmModel_.magSf();
155 
158 
159  const tmp<volScalarField> ttheta = theta();
160  const volScalarField& theta = ttheta();
161 
162  const volVectorField gradAlpha(fvc::grad(alpha));
163 
164  forAll(nbr, facei)
165  {
166  const label cellO = own[facei];
167  const label cellN = nbr[facei];
168 
169  label celli = -1;
170  if ((alpha[cellO] > 0.5) && (alpha[cellN] < 0.5))
171  {
172  celli = cellO;
173  }
174  else if ((alpha[cellO] < 0.5) && (alpha[cellN] > 0.5))
175  {
176  celli = cellN;
177  }
178 
179  if (celli != -1 && mask_[celli] > 0.5)
180  {
181  const scalar invDx = filmModel_.regionMesh().deltaCoeffs()[facei];
182  const vector n =
183  gradAlpha[celli]/(mag(gradAlpha[celli]) + ROOTVSMALL);
184  const scalar cosTheta = cos(degToRad(theta[celli]));
185  force[celli] += Ccf_*n*sigma[celli]*(1 - cosTheta)/invDx;
186  }
187  }
188 
189  forAll(alpha.boundaryField(), patchi)
190  {
191  if (!filmModel_.isCoupledPatch(patchi))
192  {
193  const fvPatchField<scalar>& alphaPf = alpha.boundaryField()[patchi];
194  const fvPatchField<scalar>& maskPf = mask_.boundaryField()[patchi];
195  const fvPatchField<scalar>& sigmaPf = sigma.boundaryField()[patchi];
196  const fvPatchField<scalar>& thetaPf = theta.boundaryField()[patchi];
197  const scalarField& invDx = alphaPf.patch().deltaCoeffs();
198  const labelUList& faceCells = alphaPf.patch().faceCells();
199 
200  forAll(alphaPf, facei)
201  {
202  if (maskPf[facei] > 0.5)
203  {
204  label cellO = faceCells[facei];
205 
206  if ((alpha[cellO] > 0.5) && (alphaPf[facei] < 0.5))
207  {
208  const vector n =
209  gradAlpha[cellO]
210  /(mag(gradAlpha[cellO]) + ROOTVSMALL);
211  const scalar cosTheta = cos(degToRad(thetaPf[facei]));
212  force[cellO] +=
213  Ccf_*n*sigmaPf[facei]*(1 - cosTheta)/invDx[facei];
214  }
215  }
216  }
217  }
218  }
219 
220  force /= magSf;
221 
223  {
224  tForce().write();
225  }
226 
228  (
230  );
231 
232  tfvm.ref() += tForce;
233 
234  return tfvm;
235 }
236 
237 
238 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
239 
240 } // End namespace surfaceFilmModels
241 } // End namespace regionModels
242 } // End namespace Foam
243 
244 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:424
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:428
U
Definition: pEqn.H:83
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:137
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.
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
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:644
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:243
surfaceFilmModel & filmModel_
Reference to the film surface film model.
const dictionary coeffDict_
Coefficients dictionary.
Definition: subModelBase.H:79
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:288
bool writeTime() const
Return true if this is a write time.
Definition: TimeStateI.H:65
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
virtual tmp< fvVectorMatrix > correct(volVectorField &U)
Correct.
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
scalar y
dimensionedScalar cos(const dimensionedScalar &ds)
Calculate the gradient of the given field.
virtual const volScalarField & alpha() const =0
Return the film coverage, 1 = covered, 0 = uncovered / [].
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:91
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
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:61
const fvPatch & patch() const
Return patch.
Definition: fvPatchField.H:340
virtual const volScalarField & sigma() const =0
Return the film surface tension [N/m].
dimensionedScalar pos0(const dimensionedScalar &ds)
const dimensionSet dimForce
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:282
virtual const volScalarField & magSf() const
Return the face area magnitudes / [m2].
label patchi
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
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
fvMatrix< vector > fvVectorMatrix
Definition: fvMatricesFwd.H:45
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
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
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
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.
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:576