curvatureSeparation.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 "curvatureSeparation.H"
27 #include "kinematicSingleLayer.H"
28 #include "fvcGrad.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 namespace regionModels
36 {
37 namespace surfaceFilmModels
38 {
39 
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 
42 defineTypeNameAndDebug(curvatureSeparation, 0);
44 (
45  ejectionModel,
46  curvatureSeparation,
47  dictionary
48 );
49 
50 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
51 
52 tmp<volScalarField> curvatureSeparation::calcInvR1
53 (
54  const volVectorField& U
55 ) const
56 {
57  const dimensionedScalar smallU("smallU", dimVelocity, rootVSmall);
58  const volVectorField UHat(U/(mag(U) + smallU));
59 
60  tmp<volScalarField> tinvR1
61  (
62  volScalarField::New("invR1", UHat & (UHat & gradNHat_))
63  );
64 
65  scalarField& invR1 = tinvR1.ref().primitiveFieldRef();
66 
67  // Apply defined patch radii
68  const scalar rMin = 1e-6;
69  const fvMesh& mesh = film().regionMesh();
70  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
71 
72  forAll(definedPatchRadii_, i)
73  {
74  const label patchi = definedPatchRadii_[i].first();
75  const scalar definedInvR1 =
76  1.0/max(rMin, definedPatchRadii_[i].second());
77  UIndirectList<scalar>(invR1, pbm[patchi].faceCells()) = definedInvR1;
78  }
79 
80  // Filter out large radii
81  const scalar rMax = 1e6;
82  forAll(invR1, i)
83  {
84  if (mag(invR1[i]) < 1/rMax)
85  {
86  invR1[i] = -1.0;
87  }
88  }
89 
90  if (debug && mesh.time().writeTime())
91  {
92  tinvR1().write();
93  }
94 
95  return tinvR1;
96 }
97 
98 
100 (
101  const surfaceScalarField& phi
102 ) const
103 {
104  const fvMesh& mesh = film().regionMesh();
105  const vectorField nf(mesh.Sf()/mesh.magSf());
106  const unallocLabelList& own = mesh.owner();
107  const unallocLabelList& nbr = mesh.neighbour();
108 
109  scalarField phiMax(mesh.nCells(), -great);
110  scalarField cosAngle(mesh.nCells(), 0);
111 
112  forAll(nbr, facei)
113  {
114  const label cellO = own[facei];
115  const label cellN = nbr[facei];
116 
117  if (phi[facei] > phiMax[cellO])
118  {
119  phiMax[cellO] = phi[facei];
120  cosAngle[cellO] = -gHat_ & nf[facei];
121  }
122  if (-phi[facei] > phiMax[cellN])
123  {
124  phiMax[cellN] = -phi[facei];
125  cosAngle[cellN] = -gHat_ & -nf[facei];
126  }
127  }
128 
129  forAll(phi.boundaryField(), patchi)
130  {
131  const fvsPatchScalarField& phip = phi.boundaryField()[patchi];
132  const fvPatch& pp = phip.patch();
133  const labelList& faceCells = pp.faceCells();
134  const vectorField nf(pp.nf());
135  forAll(phip, i)
136  {
137  const label celli = faceCells[i];
138  if (phip[i] > phiMax[celli])
139  {
140  phiMax[celli] = phip[i];
141  cosAngle[celli] = -gHat_ & nf[i];
142  }
143  }
144  }
145 
146  if (debug && mesh.time().writeTime())
147  {
148  volScalarField volCosAngle
149  (
150  IOobject
151  (
152  "cosAngle",
153  mesh.time().timeName(),
154  mesh,
156  ),
157  mesh,
159  zeroGradientFvPatchScalarField::typeName
160  );
161  volCosAngle.primitiveFieldRef() = cosAngle;
162  volCosAngle.correctBoundaryConditions();
163  volCosAngle.write();
164  }
165 
166  return max(min(cosAngle, scalar(1)), scalar(-1));
167 }
168 
169 
170 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
171 
173 (
175  const dictionary& dict
176 )
177 :
178  ejectionModel(type(), film, dict),
179  gradNHat_(fvc::grad(film.nHat())),
180  deltaByR1Min_
181  (
182  coeffDict_.lookupOrDefault<scalar>("deltaByR1Min", scalar(0))
183  ),
184  deltaStable_(coeffDict_.lookupOrDefault("deltaStable", scalar(0))),
185  definedPatchRadii_(),
186  magG_(mag(film.g().value())),
187  gHat_(Zero)
188 {
189  if (magG_ < rootVSmall)
190  {
192  << "Acceleration due to gravity must be non-zero"
193  << exit(FatalError);
194  }
195 
196  gHat_ = film.g().value()/magG_;
197 
198  const List<Tuple2<word, scalar>> prIn
199  (
200  coeffDict_.lookup("definedPatchRadii")
201  );
202  const wordList& allPatchNames = film.regionMesh().boundaryMesh().names();
203 
204  DynamicList<Tuple2<label, scalar>> prData(allPatchNames.size());
205 
206  labelHashSet uniquePatchIDs;
207 
208  forAllReverse(prIn, i)
209  {
210  labelList patchIDs = findStrings(prIn[i].first(), allPatchNames);
211  forAll(patchIDs, j)
212  {
213  const label patchi = patchIDs[j];
214 
215  if (!uniquePatchIDs.found(patchi))
216  {
217  const scalar radius = prIn[i].second();
218  prData.append(Tuple2<label, scalar>(patchi, radius));
219 
220  uniquePatchIDs.insert(patchi);
221  }
222  }
223  }
224 
225  definedPatchRadii_.transfer(prData);
226 }
227 
228 
229 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
230 
232 {}
233 
234 
235 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
236 
238 (
239  scalarField& availableMass,
240  scalarField& massToEject,
241  scalarField& diameterToEject
242 )
243 {
244  const kinematicSingleLayer& film =
245  refCast<const kinematicSingleLayer>(this->film());
246  const fvMesh& mesh = film.regionMesh();
247 
248  const volScalarField& delta = film.delta();
249  const volVectorField& U = film.U();
250  const surfaceScalarField& phi = film.phi();
251  const volScalarField& rho = film.rho();
252  const scalarField magSqrU(magSqr(film.U()));
253 
254  const tmp<volScalarField> tsigma = film.sigma();
255  const volScalarField::Internal& sigma = tsigma();
256 
257  const scalarField invR1(calcInvR1(U));
258  const scalarField cosAngle(calcCosAngle(phi));
259 
260  // Calculate force balance
261  const scalar Fthreshold = 1e-10;
262  scalarField Fnet(mesh.nCells(), 0);
263  scalarField separated(mesh.nCells(), 0);
264 
265  forAll(invR1, i)
266  {
267  if ((invR1[i] > 0) && (delta[i]*invR1[i] > deltaByR1Min_))
268  {
269  const scalar R1 = 1.0/(invR1[i] + rootVSmall);
270  const scalar R2 = R1 + delta[i];
271 
272  // Inertial force
273  const scalar Fi = -delta[i]*rho[i]*magSqrU[i]*72.0/60.0*invR1[i];
274 
275  // Body force
276  const scalar Fb =
277  - 0.5*rho[i]*magG_*invR1[i]*(sqr(R1) - sqr(R2))*cosAngle[i];
278 
279  // Surface force
280  const scalar Fs = sigma[i]/R2;
281 
282  Fnet[i] = Fi + Fb + Fs;
283 
284  if (Fnet[i] + Fthreshold < 0 && delta[i] > deltaStable_)
285  {
286  separated[i] = 1.0;
287  }
288  }
289  }
290 
291  // Eject all available mass
292  massToEject = separated*availableMass;
293  diameterToEject = separated*delta;
294  availableMass -= separated*availableMass;
295 
296  addToEjectedMass(sum(separated*availableMass));
297 
298  if (debug && mesh.time().writeTime())
299  {
300  volScalarField volFnet
301  (
302  IOobject
303  (
304  "Fnet",
305  mesh.time().timeName(),
306  mesh,
308  ),
309  mesh,
311  zeroGradientFvPatchScalarField::typeName
312  );
313  volFnet.primitiveFieldRef() = Fnet;
314  volFnet.correctBoundaryConditions();
315  volFnet.write();
316  }
317 
319 }
320 
321 
322 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
323 
324 } // End namespace surfaceFilmModels
325 } // End namespace regionModels
326 } // End namespace Foam
327 
328 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
scalar delta
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
const volVectorField & U() const
Return the film velocity [m/s].
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const surfaceVectorField & Sf() const
Return cell face area vectors.
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
tmp< scalarField > calcCosAngle(const surfaceScalarField &phi) const
Calculate the cosine of the angle between gravity vector and.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Kinematic form of single-cell layer surface film model.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
A 2-tuple for storing two objects of different types.
Definition: HashTable.H:65
scalar deltaStable_
Stable film thickness - drips only formed if thickness.
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
tmp< vectorField > nf() const
Return face normals.
Definition: fvPatch.C:130
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Base class for film ejection models, handling mass transfer from the film.
Definition: ejectionModel.H:56
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
label nCells() const
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< scalar >> &)
Return a temporary field constructed from name,.
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:62
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: UList.H:446
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 top-level database.
Definition: fvMesh.H:239
tmp< volScalarField > calcInvR1(const volVectorField &U) const
Calculate local (inverse) radius of curvature.
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
const surfaceFilmRegionModel & film() const
Return const access to the film surface film model.
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:284
bool findStrings(const wordReListMatcher &matcher, const std::string &str)
Return true if string matches one of the regular expressions.
Definition: stringListOps.H:52
bool writeTime() const
Return true if this is a write time.
Definition: TimeStateI.H:65
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
dynamicFvMesh & mesh
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
Calculate the gradient of the given field.
virtual const labelUList & faceCells() const
Return faceCells.
Definition: fvPatch.C:99
wordList names() const
Return a list of patch names.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
const Type & value() const
Return const reference to value.
const volVectorField & nHat() const
Return the patch normal vectors.
const surfaceScalarField & phi() const
Return the film flux [kg m/s].
const fvMesh & regionMesh() const
Return the region mesh database.
Definition: regionModelI.H:55
static const zero Zero
Definition: zero.H:97
const dimensionSet dimForce
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
Foam::polyBoundaryMesh.
const dimensionedVector & g() const
Return the acceleration due to gravity.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
void addToEjectedMass(const scalar dMass)
Add to ejected mass.
Definition: ejectionModel.C:44
const dimensionSet dimVelocity
const volScalarField & delta() const
Return const access to the film thickness [m].
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
tmp< volScalarField > sigma() const
Return the surface tension coefficient [kg/s^2].
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:278
const fvPatch & patch() const
Return patch.
label patchi
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
const volScalarField & rho() const
Return the film density [kg/m^3].
A class for managing temporary objects.
Definition: PtrList.H:53
An indexed form of CGAL::Triangulation_face_base_2<K> used to keep track of the vertices in the trian...
Definition: indexedFace.H:48
void transfer(List< T > &)
Transfer contents of the argument List into this.
Definition: DynamicListI.H:285
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
defineTypeNameAndDebug(kinematicSingleLayer, 0)
curvatureSeparation(surfaceFilmRegionModel &film, const dictionary &dict)
Construct from surface film model.
scalar deltaByR1Min_
Minimum gravity driven film thickness (non-dimensionalised delta/R1)
An abstract base class with a fat-interface to all derived classes covering all possible ways in whic...
Definition: fvsPatchField.H:65
addToRunTimeSelectionTable(surfaceFilmModel, noFilm, mesh)
Namespace for OpenFOAM.