curvatureSeparation.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-2013 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 "curvatureSeparation.H"
28 #include "fvMesh.H"
29 #include "Time.H"
30 #include "volFields.H"
31 #include "kinematicSingleLayer.H"
32 #include "surfaceInterpolate.H"
33 #include "fvcDiv.H"
34 #include "fvcGrad.H"
35 #include "stringListOps.H"
36 #include "cyclicPolyPatch.H"
37 
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42 namespace regionModels
43 {
44 namespace surfaceFilmModels
45 {
46 
47 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
48 
49 defineTypeNameAndDebug(curvatureSeparation, 0);
51 (
52  injectionModel,
53  curvatureSeparation,
54  dictionary
55 );
56 
57 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
58 
59 tmp<volScalarField> curvatureSeparation::calcInvR1
60 (
61  const volVectorField& U
62 ) const
63 {
64  // method 1
65 /*
66  tmp<volScalarField> tinvR1
67  (
68  new volScalarField("invR1", fvc::div(owner().nHat()))
69  );
70 */
71 
72  // method 2
73  dimensionedScalar smallU("smallU", dimVelocity, ROOTVSMALL);
74  volVectorField UHat(U/(mag(U) + smallU));
75  tmp<volScalarField> tinvR1
76  (
77  new volScalarField("invR1", UHat & (UHat & gradNHat_))
78  );
79 
80 
81  scalarField& invR1 = tinvR1().internalField();
82 
83  // apply defined patch radii
84  const scalar rMin = 1e-6;
85  const fvMesh& mesh = owner().regionMesh();
86  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
87  forAll(definedPatchRadii_, i)
88  {
89  label patchI = definedPatchRadii_[i].first();
90  scalar definedInvR1 = 1.0/max(rMin, definedPatchRadii_[i].second());
91  UIndirectList<scalar>(invR1, pbm[patchI].faceCells()) = definedInvR1;
92  }
93 
94  // filter out large radii
95  const scalar rMax = 1e6;
96  forAll(invR1, i)
97  {
98  if (mag(invR1[i]) < 1/rMax)
99  {
100  invR1[i] = -1.0;
101  }
102  }
103 
104  if (debug && mesh.time().outputTime())
105  {
106  tinvR1().write();
107  }
108 
109  return tinvR1;
110 }
111 
112 
114 (
115  const surfaceScalarField& phi
116 ) const
117 {
118  const fvMesh& mesh = owner().regionMesh();
119  const vectorField nf(mesh.Sf()/mesh.magSf());
120  const unallocLabelList& own = mesh.owner();
121  const unallocLabelList& nbr = mesh.neighbour();
122 
123  scalarField phiMax(mesh.nCells(), -GREAT);
124  scalarField cosAngle(mesh.nCells(), 0.0);
125  forAll(nbr, faceI)
126  {
127  label cellO = own[faceI];
128  label cellN = nbr[faceI];
129 
130  if (phi[faceI] > phiMax[cellO])
131  {
132  phiMax[cellO] = phi[faceI];
133  cosAngle[cellO] = -gHat_ & nf[faceI];
134  }
135  if (-phi[faceI] > phiMax[cellN])
136  {
137  phiMax[cellN] = -phi[faceI];
138  cosAngle[cellN] = -gHat_ & -nf[faceI];
139  }
140  }
141 
142  forAll(phi.boundaryField(), patchI)
143  {
144  const fvsPatchScalarField& phip = phi.boundaryField()[patchI];
145  const fvPatch& pp = phip.patch();
146  const labelList& faceCells = pp.faceCells();
147  const vectorField nf(pp.nf());
148  forAll(phip, i)
149  {
150  label cellI = faceCells[i];
151  if (phip[i] > phiMax[cellI])
152  {
153  phiMax[cellI] = phip[i];
154  cosAngle[cellI] = -gHat_ & nf[i];
155  }
156  }
157  }
158 /*
159  // correction for cyclics - use cyclic pairs' face normal instead of
160  // local face normal
161  const fvBoundaryMesh& pbm = mesh.boundary();
162  forAll(phi.boundaryField(), patchI)
163  {
164  if (isA<cyclicPolyPatch>(pbm[patchI]))
165  {
166  const scalarField& phip = phi.boundaryField()[patchI];
167  const vectorField nf(pbm[patchI].nf());
168  const labelList& faceCells = pbm[patchI].faceCells();
169  const label sizeBy2 = pbm[patchI].size()/2;
170 
171  for (label face0=0; face0<sizeBy2; face0++)
172  {
173  label face1 = face0 + sizeBy2;
174  label cell0 = faceCells[face0];
175  label cell1 = faceCells[face1];
176 
177  // flux leaving half 0, entering half 1
178  if (phip[face0] > phiMax[cell0])
179  {
180  phiMax[cell0] = phip[face0];
181  cosAngle[cell0] = -gHat_ & -nf[face1];
182  }
183 
184  // flux leaving half 1, entering half 0
185  if (-phip[face1] > phiMax[cell1])
186  {
187  phiMax[cell1] = -phip[face1];
188  cosAngle[cell1] = -gHat_ & nf[face0];
189  }
190  }
191  }
192  }
193 */
194  // checks
195  if (debug && mesh.time().outputTime())
196  {
197  volScalarField volCosAngle
198  (
199  IOobject
200  (
201  "cosAngle",
202  mesh.time().timeName(),
203  mesh,
205  ),
206  mesh,
207  dimensionedScalar("zero", dimless, 0.0),
208  zeroGradientFvPatchScalarField::typeName
209  );
210  volCosAngle.internalField() = cosAngle;
211  volCosAngle.correctBoundaryConditions();
212  volCosAngle.write();
213  }
214 
215  return max(min(cosAngle, scalar(1.0)), scalar(-1.0));
216 }
217 
218 
219 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
220 
221 curvatureSeparation::curvatureSeparation
222 (
223  surfaceFilmModel& owner,
224  const dictionary& dict
225 )
226 :
227  injectionModel(type(), owner, dict),
228  gradNHat_(fvc::grad(owner.nHat())),
229  deltaByR1Min_(coeffDict_.lookupOrDefault<scalar>("deltaByR1Min", 0.0)),
230  definedPatchRadii_(),
231  magG_(mag(owner.g().value())),
232  gHat_(vector::zero)
233 {
234  if (magG_ < ROOTVSMALL)
235  {
237  (
238  "curvatureSeparation::curvatureSeparation"
239  "("
240  "const surfaceFilmModel&, "
241  "const dictionary&"
242  ")"
243  )
244  << "Acceleration due to gravity must be non-zero"
245  << exit(FatalError);
246  }
247 
248  gHat_ = owner.g().value()/magG_;
249 
250  List<Tuple2<word, scalar> > prIn(coeffDict_.lookup("definedPatchRadii"));
251  const wordList& allPatchNames = owner.regionMesh().boundaryMesh().names();
252 
253  DynamicList<Tuple2<label, scalar> > prData(allPatchNames.size());
254 
255  labelHashSet uniquePatchIDs;
256 
257  forAllReverse(prIn, i)
258  {
259  labelList patchIDs = findStrings(prIn[i].first(), allPatchNames);
260  forAll(patchIDs, j)
261  {
262  const label patchI = patchIDs[j];
263 
264  if (!uniquePatchIDs.found(patchI))
265  {
266  const scalar radius = prIn[i].second();
267  prData.append(Tuple2<label, scalar>(patchI, radius));
268 
269  uniquePatchIDs.insert(patchI);
270  }
271  }
272  }
273 
274  definedPatchRadii_.transfer(prData);
275 }
276 
277 
278 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
279 
281 {}
282 
283 
284 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
285 
287 (
288  scalarField& availableMass,
289  scalarField& massToInject,
290  scalarField& diameterToInject
291 )
292 {
293  const kinematicSingleLayer& film =
294  refCast<const kinematicSingleLayer>(this->owner());
295  const fvMesh& mesh = film.regionMesh();
296 
297  const volScalarField& delta = film.delta();
298  const volVectorField& U = film.U();
299  const surfaceScalarField& phi = film.phi();
300  const volScalarField& rho = film.rho();
301  const scalarField magSqrU(magSqr(film.U()));
302  const volScalarField& sigma = film.sigma();
303 
304  const scalarField invR1(calcInvR1(U));
305  const scalarField cosAngle(calcCosAngle(phi));
306 
307  // calculate force balance
308  const scalar Fthreshold = 1e-10;
309  scalarField Fnet(mesh.nCells(), 0.0);
310  scalarField separated(mesh.nCells(), 0.0);
311  forAll(invR1, i)
312  {
313  if ((invR1[i] > 0) && (delta[i]*invR1[i] > deltaByR1Min_))
314  {
315  scalar R1 = 1.0/(invR1[i] + ROOTVSMALL);
316  scalar R2 = R1 + delta[i];
317 
318  // inertial force
319  scalar Fi = -delta[i]*rho[i]*magSqrU[i]*72.0/60.0*invR1[i];
320 
321  // body force
322  scalar Fb =
323  - 0.5*rho[i]*magG_*invR1[i]*(sqr(R1) - sqr(R2))*cosAngle[i];
324 
325  // surface force
326  scalar Fs = sigma[i]/R2;
327 
328  Fnet[i] = Fi + Fb + Fs;
329 
330  if (Fnet[i] + Fthreshold < 0)
331  {
332  separated[i] = 1.0;
333  }
334  }
335  }
336 
337  // inject all available mass
338  massToInject = separated*availableMass;
339  diameterToInject = separated*delta;
340  availableMass -= separated*availableMass;
341 
342  addToInjectedMass(sum(separated*availableMass));
343 
344  if (debug && mesh.time().outputTime())
345  {
346  volScalarField volFnet
347  (
348  IOobject
349  (
350  "Fnet",
351  mesh.time().timeName(),
352  mesh,
354  ),
355  mesh,
356  dimensionedScalar("zero", dimForce, 0.0),
357  zeroGradientFvPatchScalarField::typeName
358  );
359  volFnet.internalField() = Fnet;
360  volFnet.correctBoundaryConditions();
361  volFnet.write();
362  }
363 
365 }
366 
367 
368 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
369 
370 } // End namespace surfaceFilmModels
371 } // End namespace regionModels
372 } // End namespace Foam
373 
374 // ************************************************************************* //
void transfer(List< T > &)
Transfer contents of the argument List into this.
Definition: DynamicListI.H:277
A 2-tuple for storing two objects of different types.
Definition: Tuple2.H:47
const dimensionSet dimForce
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
void addToInjectedMass(const scalar dMass)
Add to injected mass.
dimensioned< scalar > mag(const dimensioned< Type > &)
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
const surfaceVectorField & Sf() const
Return cell face area vectors.
defineTypeNameAndDebug(kinematicSingleLayer, 0)
const fvPatch & patch() const
Return patch.
virtual const volScalarField & rho() const
Return the film density [kg/m3].
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
scalar deltaByR1Min_
Minimum gravity driven film thickness (non-dimensionalised delta/R1)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
const dimensionedVector & g() const
Return the accleration due to gravity.
Base class for film injection models, handling mass transfer from the film.
const volScalarField & delta() const
Return const access to the film thickness / [m].
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
addToRunTimeSelectionTable(surfaceFilmModel, kinematicSingleLayer, mesh)
InternalField & internalField()
Return internal field.
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:282
label nCells() const
Surface Interpolation.
tmp< volScalarField > calcInvR1(const volVectorField &U) const
Calculate local (inverse) radius of curvature.
const fvMesh & regionMesh() const
Return the region mesh database.
Definition: regionModelI.H:61
tmp< scalarField > calcCosAngle(const surfaceScalarField &phi) const
Calculate the cosine of the angle between gravity vector and.
virtual const labelUList & faceCells() const
Return faceCells.
Definition: fvPatch.C:93
virtual const volVectorField & nHat() const
Return the patch normal vectors.
dynamicFvMesh & mesh
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
Namespace for OpenFOAM.
Calculate the gradient of the given field.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
wordList names() const
Return a list of patch names.
const double e
Elementary charge.
Definition: doubleFloat.H:78
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
virtual const volVectorField & U() const
Return the film velocity [m/s].
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
A List with indirect addressing.
Definition: fvMatrix.H:106
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
Foam::polyBoundaryMesh.
bool findStrings(const wordReListMatcher &matcher, const std::string &str)
Return true if string matches one of the regular expressions.
Definition: stringListOps.H:52
const surfaceFilmModel & owner() const
Return const access to the owner surface film model.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
const volScalarField & sigma() const
Return const access to the surface tension / [m/s2].
scalar delta
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
fileName::Type type(const fileName &)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:589
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
virtual const surfaceScalarField & phi() const
Return the film flux [kg.m/s].
tmp< vectorField > nf() const
Return face normals.
Definition: fvPatch.C:124
#define forAllReverse(list, i)
Definition: UList.H:424
Calculate the divergence of the given field.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
error FatalError
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
static const Vector zero
Definition: Vector.H:80
bool outputTime() const
Return true if this is an output time (primary or secondary)
Definition: TimeState.C:91
Operations on lists of strings.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:61
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:97
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
An indexed form of CGAL::Triangulation_face_base_2<K> used to keep track of the vertices in the trian...
Definition: indexedFace.H:48
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:288
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const dimensionSet dimVelocity
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
An abstract base class with a fat-interface to all derived classes covering all possible ways in whic...
Definition: fvsPatchField.H:65
const Type & value() const
Return const reference to value.