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-2015 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"
30 #include "fvPatchField.H"
32 
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37 namespace regionModels
38 {
39 namespace surfaceFilmModels
40 {
41 
42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43 
44 defineTypeNameAndDebug(contactAngleForce, 0);
45 addToRunTimeSelectionTable(force, contactAngleForce, dictionary);
46 
47 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
48 
49 void contactAngleForce::initialise()
50 {
51  const wordReList zeroForcePatches(coeffDict_.lookup("zeroForcePatches"));
52 
53  if (zeroForcePatches.size())
54  {
55  const polyBoundaryMesh& pbm = owner_.regionMesh().boundaryMesh();
56  scalar dLim = readScalar(coeffDict_.lookup("zeroForceDistance"));
57 
58  Info<< " Assigning zero contact force within " << dLim
59  << " of patches:" << endl;
60 
61  labelHashSet patchIDs = pbm.patchSet(zeroForcePatches);
62 
63  forAllConstIter(labelHashSet, patchIDs, iter)
64  {
65  label patchI = iter.key();
66  Info<< " " << pbm[patchI].name() << endl;
67  }
68 
69  // Temporary implementation until run-time selection covers this case
70  patchDistMethods::meshWave dist(owner_.regionMesh(), patchIDs);
72  (
73  IOobject
74  (
75  "y",
80  false
81  ),
83  dimensionedScalar("y", dimLength, GREAT)
84  );
85  dist.correct(y);
86 
87  mask_ = pos(y - dimensionedScalar("dLim", dimLength, dLim));
88  }
89 }
90 
91 
92 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
93 
94 contactAngleForce::contactAngleForce
95 (
96  surfaceFilmModel& owner,
97  const dictionary& dict
98 )
99 :
100  force(typeName, owner, dict),
101  Ccf_(readScalar(coeffDict_.lookup("Ccf"))),
102  rndGen_(label(0), -1),
103  distribution_
104  (
106  (
107  coeffDict_.subDict("contactAngleDistribution"),
108  rndGen_
109  )
110  ),
111  mask_
112  (
113  IOobject
114  (
115  typeName + ":contactForceMask",
116  owner_.time().timeName(),
117  owner_.regionMesh(),
120  ),
121  owner_.regionMesh(),
122  dimensionedScalar("mask", dimless, 1.0)
123  )
124 {
125  initialise();
126 }
127 
128 
129 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
130 
132 {}
133 
134 
135 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
136 
138 {
139  tmp<volVectorField> tForce
140  (
141  new volVectorField
142  (
143  IOobject
144  (
145  typeName + ":contactForce",
146  owner_.time().timeName(),
147  owner_.regionMesh(),
150  ),
151  owner_.regionMesh(),
153  )
154  );
155 
156  vectorField& force = tForce().internalField();
157 
158  const labelUList& own = owner_.regionMesh().owner();
159  const labelUList& nbr = owner_.regionMesh().neighbour();
160 
161  const scalarField& magSf = owner_.magSf();
162 
163  const volScalarField& alpha = owner_.alpha();
164  const volScalarField& sigma = owner_.sigma();
165 
166  volVectorField gradAlpha(fvc::grad(alpha));
167 
168  forAll(nbr, faceI)
169  {
170  const label cellO = own[faceI];
171  const label cellN = nbr[faceI];
172 
173  label cellI = -1;
174  if ((alpha[cellO] > 0.5) && (alpha[cellN] < 0.5))
175  {
176  cellI = cellO;
177  }
178  else if ((alpha[cellO] < 0.5) && (alpha[cellN] > 0.5))
179  {
180  cellI = cellN;
181  }
182 
183  if (cellI != -1 && mask_[cellI] > 0.5)
184  {
185  const scalar invDx = owner_.regionMesh().deltaCoeffs()[faceI];
186  const vector n =
187  gradAlpha[cellI]/(mag(gradAlpha[cellI]) + ROOTVSMALL);
188  scalar theta = cos(degToRad(distribution_->sample()));
189  force[cellI] += Ccf_*n*sigma[cellI]*(1.0 - theta)/invDx;
190  }
191  }
192 
193  forAll(alpha.boundaryField(), patchI)
194  {
195  if (!owner().isCoupledPatch(patchI))
196  {
197  const fvPatchField<scalar>& alphaf = alpha.boundaryField()[patchI];
198  const fvPatchField<scalar>& maskf = mask_.boundaryField()[patchI];
199  const scalarField& invDx = alphaf.patch().deltaCoeffs();
200  const labelUList& faceCells = alphaf.patch().faceCells();
201 
202  forAll(alphaf, faceI)
203  {
204  if (maskf[faceI] > 0.5)
205  {
206  label cellO = faceCells[faceI];
207 
208  if ((alpha[cellO] > 0.5) && (alphaf[faceI] < 0.5))
209  {
210  const vector n =
211  gradAlpha[cellO]
212  /(mag(gradAlpha[cellO]) + ROOTVSMALL);
213  scalar theta = cos(degToRad(distribution_->sample()));
214  force[cellO] +=
215  Ccf_*n*sigma[cellO]*(1.0 - theta)/invDx[faceI];
216  }
217  }
218  }
219  }
220  }
221 
222  force /= magSf;
223 
224  if (owner_.regionMesh().time().outputTime())
225  {
226  tForce().write();
227  }
228 
231 
232  tfvm() += tForce;
233 
234  return tfvm;
235 }
236 
237 
238 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
239 
240 } // End namespace surfaceFilmModels
241 } // End namespace regionModels
242 } // End namespace Foam
243 
244 // ************************************************************************* //
surfaceFilmModel & owner_
Reference to the owner surface film model.
const dimensionSet dimForce
dimensioned< scalar > mag(const dimensioned< Type > &)
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
defineTypeNameAndDebug(kinematicSingleLayer, 0)
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
List< wordRe > wordReList
A List of wordRe (word or regular expression)
Definition: wordReList.H:50
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
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
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:741
const dictionary coeffDict_
Coefficients dictionary.
Definition: subModelBase.H:79
addToRunTimeSelectionTable(surfaceFilmModel, kinematicSingleLayer, mesh)
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:282
const fvMesh & regionMesh() const
Return the region mesh database.
Definition: regionModelI.H:61
const surfaceScalarField & deltaCoeffs() const
Return reference to cell-centre difference coefficients.
const Time & time() const
Return the reference to the time database.
Definition: regionModelI.H:37
messageStream Info
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 dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:638
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
Namespace for OpenFOAM.
fvMatrix< vector > fvVectorMatrix
Definition: fvMatricesFwd.H:45
virtual tmp< fvVectorMatrix > correct(volVectorField &U)
Correct.
Calculate the gradient of the given field.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
label n
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
const surfaceFilmModel & owner() const
Return const access to the owner surface film model.
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
virtual const volScalarField & alpha() const =0
Return the film coverage, 1 = covered, 0 = uncovered / [].
#define forAll(list, i)
Definition: UList.H:421
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
virtual const volScalarField & magSf() const
Return the face area magnitudes / [m2].
dimensionedScalar cos(const dimensionedScalar &ds)
Macros for easy insertion into run-time selection tables.
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Unit conversion functions.
static autoPtr< distributionModel > New(const dictionary &dict, cachedRandom &rndGen)
Selector.
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:65
virtual const volScalarField & sigma() const =0
Return the film surface tension [N/m].
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:452
scalar y
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
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
bool outputTime() const
Return true if this is an output time (primary or secondary)
Definition: TimeState.C:91
static const Vector zero
Definition: Vector.H:80
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
bool isCoupledPatch(const label regionPatchI) const
Return true if patchI on the local region is a coupled.
Definition: regionModelI.H:152
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
dimensionedScalar pos(const dimensionedScalar &ds)
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:288
const fvPatch & patch() const
Return patch.
Definition: fvPatchField.H:300
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
U
Definition: pEqn.H:82