epsilonWallFunctionFvPatchScalarField.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-2018 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 
28 #include "turbulenceModel.H"
29 #include "fvPatchFieldMapper.H"
30 #include "fvMatrix.H"
31 #include "volFields.H"
32 #include "wallFvPatch.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
38 
39 // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
40 
42 {
43  if (!isA<wallFvPatch>(patch()))
44  {
46  << "Invalid wall function specification" << nl
47  << " Patch type for patch " << patch().name()
48  << " must be wall" << nl
49  << " Current patch type is " << patch().type() << nl << endl
50  << abort(FatalError);
51  }
52 }
53 
54 
56 (
57  Ostream& os
58 ) const
59 {
60  os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
61  os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
62  os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
63 }
64 
65 
67 {
68  if (master_ != -1)
69  {
70  return;
71  }
72 
73  const volScalarField& epsilon =
74  static_cast<const volScalarField&>(this->internalField());
75 
76  const volScalarField::Boundary& bf = epsilon.boundaryField();
77 
78  label master = -1;
79  forAll(bf, patchi)
80  {
81  if (isA<epsilonWallFunctionFvPatchScalarField>(bf[patchi]))
82  {
84 
85  if (master == -1)
86  {
87  master = patchi;
88  }
89 
90  epf.master() = master;
91  }
92  }
93 }
94 
95 
97 {
98  const volScalarField& epsilon =
99  static_cast<const volScalarField&>(this->internalField());
100 
101  const volScalarField::Boundary& bf = epsilon.boundaryField();
102 
103  const fvMesh& mesh = epsilon.mesh();
104 
105  if (initialised_ && !mesh.changing())
106  {
107  return;
108  }
109 
110  volScalarField weights
111  (
112  IOobject
113  (
114  "weights",
115  mesh.time().timeName(),
116  mesh,
119  false // do not register
120  ),
121  mesh,
122  dimensionedScalar("zero", dimless, 0.0)
123  );
124 
125  DynamicList<label> epsilonPatches(bf.size());
126  forAll(bf, patchi)
127  {
128  if (isA<epsilonWallFunctionFvPatchScalarField>(bf[patchi]))
129  {
130  epsilonPatches.append(patchi);
131 
132  const labelUList& faceCells = bf[patchi].patch().faceCells();
133  forAll(faceCells, i)
134  {
135  weights[faceCells[i]]++;
136  }
137  }
138  }
139 
140  cornerWeights_.setSize(bf.size());
141  forAll(epsilonPatches, i)
142  {
143  label patchi = epsilonPatches[i];
144  const fvPatchScalarField& wf = weights.boundaryField()[patchi];
146  }
147 
148  G_.setSize(internalField().size(), 0.0);
150 
151  initialised_ = true;
152 }
153 
154 
157 {
158  const volScalarField& epsilon =
159  static_cast<const volScalarField&>(this->internalField());
160 
161  const volScalarField::Boundary& bf = epsilon.boundaryField();
162 
164  refCast<const epsilonWallFunctionFvPatchScalarField>(bf[patchi]);
165 
166  return const_cast<epsilonWallFunctionFvPatchScalarField&>(epf);
167 }
168 
169 
171 (
172  const turbulenceModel& turbulence,
173  scalarField& G0,
174  scalarField& epsilon0
175 )
176 {
177  // Accumulate all of the G and epsilon contributions
179  {
180  if (!cornerWeights_[patchi].empty())
181  {
183 
184  const List<scalar>& w = cornerWeights_[patchi];
185 
186  epf.calculate(turbulence, w, epf.patch(), G0, epsilon0);
187  }
188  }
189 
190  // Apply zero-gradient condition for epsilon
192  {
193  if (!cornerWeights_[patchi].empty())
194  {
196 
197  epf == scalarField(epsilon0, epf.patch().faceCells());
198  }
199  }
200 }
201 
202 
204 (
205  const turbulenceModel& turbModel,
206  const List<scalar>& cornerWeights,
207  const fvPatch& patch,
208  scalarField& G0,
209  scalarField& epsilon0
210 )
211 {
212  const label patchi = patch.index();
213 
214  const scalarField& y = turbModel.y()[patchi];
215 
216  const scalar Cmu25 = pow025(Cmu_);
217  const scalar Cmu75 = pow(Cmu_, 0.75);
218 
219  const tmp<volScalarField> tk = turbModel.k();
220  const volScalarField& k = tk();
221 
222  const tmp<scalarField> tnuw = turbModel.nu(patchi);
223  const scalarField& nuw = tnuw();
224 
225  const tmp<scalarField> tnutw = turbModel.nut(patchi);
226  const scalarField& nutw = tnutw();
227 
228  const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
229 
230  const scalarField magGradUw(mag(Uw.snGrad()));
231 
232  // Set epsilon and G
233  forAll(nutw, facei)
234  {
235  const label celli = patch.faceCells()[facei];
236 
237  const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei];
238 
239  const scalar w = cornerWeights[facei];
240 
241  if (yPlus > yPlusLam_)
242  {
243  epsilon0[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);
244 
245  G0[celli] +=
246  w
247  *(nutw[facei] + nuw[facei])
248  *magGradUw[facei]
249  *Cmu25*sqrt(k[celli])
250  /(kappa_*y[facei]);
251  }
252  else
253  {
254  epsilon0[celli] += w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
255  }
256  }
257 }
258 
259 
260 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
261 
264 (
265  const fvPatch& p,
267 )
268 :
270  Cmu_(0.09),
271  kappa_(0.41),
272  E_(9.8),
274  G_(),
275  epsilon_(),
276  initialised_(false),
277  master_(-1),
279 {
280  checkType();
281 }
282 
283 
286 (
287  const fvPatch& p,
289  const dictionary& dict
290 )
291 :
293  Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
294  kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
295  E_(dict.lookupOrDefault<scalar>("E", 9.8)),
297  G_(),
298  epsilon_(),
299  initialised_(false),
300  master_(-1),
302 {
303  checkType();
304 
305  // Apply zero-gradient condition on start-up
307 }
308 
309 
312 (
314  const fvPatch& p,
316  const fvPatchFieldMapper& mapper
317 )
318 :
319  fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
320  Cmu_(ptf.Cmu_),
321  kappa_(ptf.kappa_),
322  E_(ptf.E_),
323  yPlusLam_(ptf.yPlusLam_),
324  G_(),
325  epsilon_(),
326  initialised_(false),
327  master_(-1),
329 {
330  checkType();
331 }
332 
333 
336 (
338 )
339 :
341  Cmu_(ewfpsf.Cmu_),
342  kappa_(ewfpsf.kappa_),
343  E_(ewfpsf.E_),
344  yPlusLam_(ewfpsf.yPlusLam_),
345  G_(),
346  epsilon_(),
347  initialised_(false),
348  master_(-1),
350 {
351  checkType();
352 }
353 
354 
357 (
360 )
361 :
362  fixedValueFvPatchField<scalar>(ewfpsf, iF),
363  Cmu_(ewfpsf.Cmu_),
364  kappa_(ewfpsf.kappa_),
365  E_(ewfpsf.E_),
366  yPlusLam_(ewfpsf.yPlusLam_),
367  G_(),
368  epsilon_(),
369  initialised_(false),
370  master_(-1),
372 {
373  checkType();
374 }
375 
376 
377 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
378 
380 {
381  if (patch().index() == master_)
382  {
383  if (init)
384  {
385  G_ = 0.0;
386  }
387 
388  return G_;
389  }
390 
391  return epsilonPatch(master_).G();
392 }
393 
394 
396 (
397  bool init
398 )
399 {
400  if (patch().index() == master_)
401  {
402  if (init)
403  {
404  epsilon_ = 0.0;
405  }
406 
407  return epsilon_;
408  }
409 
410  return epsilonPatch(master_).epsilon(init);
411 }
412 
413 
415 {
416  if (updated())
417  {
418  return;
419  }
420 
421  const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
422  (
424  (
426  internalField().group()
427  )
428  );
429 
430  setMaster();
431 
432  if (patch().index() == master_)
433  {
435  calculateTurbulenceFields(turbModel, G(true), epsilon(true));
436  }
437 
438  const scalarField& G0 = this->G();
439  const scalarField& epsilon0 = this->epsilon();
440 
441  typedef DimensionedField<scalar, volMesh> FieldType;
442 
443  FieldType& G =
444  const_cast<FieldType&>
445  (
446  db().lookupObject<FieldType>(turbModel.GName())
447  );
448 
449  FieldType& epsilon = const_cast<FieldType&>(internalField());
450 
451  forAll(*this, facei)
452  {
453  label celli = patch().faceCells()[facei];
454 
455  G[celli] = G0[celli];
456  epsilon[celli] = epsilon0[celli];
457  }
458 
460 }
461 
462 
464 (
465  const scalarField& weights
466 )
467 {
468  if (updated())
469  {
470  return;
471  }
472 
473  const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
474  (
476  (
478  internalField().group()
479  )
480  );
481 
482  setMaster();
483 
484  if (patch().index() == master_)
485  {
487  calculateTurbulenceFields(turbModel, G(true), epsilon(true));
488  }
489 
490  const scalarField& G0 = this->G();
491  const scalarField& epsilon0 = this->epsilon();
492 
493  typedef DimensionedField<scalar, volMesh> FieldType;
494 
495  FieldType& G =
496  const_cast<FieldType&>
497  (
498  db().lookupObject<FieldType>(turbModel.GName())
499  );
500 
501  FieldType& epsilon = const_cast<FieldType&>(internalField());
502 
503  scalarField& epsilonf = *this;
504 
505  // Only set the values if the weights are > tolerance
506  forAll(weights, facei)
507  {
508  scalar w = weights[facei];
509 
510  if (w > tolerance_)
511  {
512  label celli = patch().faceCells()[facei];
513 
514  G[celli] = (1.0 - w)*G[celli] + w*G0[celli];
515  epsilon[celli] = (1.0 - w)*epsilon[celli] + w*epsilon0[celli];
516  epsilonf[facei] = epsilon[celli];
517  }
518  }
519 
521 }
522 
523 
525 (
526  fvMatrix<scalar>& matrix
527 )
528 {
529  if (manipulatedMatrix())
530  {
531  return;
532  }
533 
534  matrix.setValues(patch().faceCells(), patchInternalField());
535 
537 }
538 
539 
541 (
542  fvMatrix<scalar>& matrix,
543  const Field<scalar>& weights
544 )
545 {
546  if (manipulatedMatrix())
547  {
548  return;
549  }
550 
551  DynamicList<label> constraintCells(weights.size());
552  DynamicList<scalar> constraintEpsilon(weights.size());
553  const labelUList& faceCells = patch().faceCells();
554 
556  = internalField();
557 
558  label nConstrainedCells = 0;
559 
560 
561  forAll(weights, facei)
562  {
563  // Only set the values if the weights are > tolerance
564  if (weights[facei] > tolerance_)
565  {
566  nConstrainedCells++;
567 
568  label celli = faceCells[facei];
569 
570  constraintCells.append(celli);
571  constraintEpsilon.append(epsilon[celli]);
572  }
573  }
574 
575  if (debug)
576  {
577  Pout<< "Patch: " << patch().name()
578  << ": number of constrained cells = " << nConstrainedCells
579  << " out of " << patch().size()
580  << endl;
581  }
582 
583  matrix.setValues
584  (
585  constraintCells,
586  scalarField(constraintEpsilon.xfer())
587  );
588 
590 }
591 
592 
594 {
595  writeLocalEntries(os);
597 }
598 
599 
600 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
601 
602 namespace Foam
603 {
605  (
608  );
609 }
610 
611 
612 // ************************************************************************* //
const char *const group
Group name for atomic constants.
bool changing() const
Is mesh changing (topology changing and/or moving)
Definition: polyMesh.H:530
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:313
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 void write(Ostream &) const
Write.
virtual void writeLocalEntries(Ostream &) const
Write local wall function variables.
virtual tmp< Field< Type > > snGrad() const
Return patch-normal gradient.
Definition: fvPatchField.C:216
const dimensionedScalar G0
Conductance quantum: default SI units: [S].
scalar yPlusLam() const
Return the Y+ at the edge of the laminar sublayer.
virtual void calculate(const turbulenceModel &turbulence, const List< scalar > &cornerWeights, const fvPatch &patch, scalarField &G, scalarField &epsilon)
Calculate the epsilon and G.
const volVectorField & U() const
Access function to velocity field.
label index() const
Return the index of this patch in the fvBoundaryMesh.
Definition: fvPatch.H:179
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
bool updated() const
Return true if the boundary condition has already been updated.
Definition: fvPatchField.H:371
const Boundary & boundaryField() const
Return const-reference to the boundary field.
scalarField G_
Local copy of turbulence G field.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
virtual tmp< volScalarField > nu() const =0
Return the laminar viscosity.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:61
dimensionedScalar pow025(const dimensionedScalar &ds)
virtual void manipulateMatrix(fvMatrix< scalar > &matrix)
Manipulate matrix.
static scalar tolerance_
Tolerance used in weighted calculations.
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:626
label k
Boltzmann constant.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
bool manipulatedMatrix() const
Return true if the matrix has already been manipulated.
Definition: fvPatchField.H:377
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type.
word GName() const
Helper function to return the name of the turbulence G field.
const dimensionedScalar epsilon0
Electric constant: default SI units: [F/m].
Abstract base class for turbulence models (RAS, LES and laminar).
const word & name() const
Return name.
Definition: fvPatch.H:149
Macros for easy insertion into run-time selection tables.
virtual void operator==(const fvPatchField< Type > &)
Definition: fvPatchField.C:561
virtual tmp< volScalarField > k() const =0
Return the turbulence kinetic energy.
scalar y
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:98
dynamicFvMesh & mesh
virtual label & master()
Return non-const access to the master patch ID.
static const word propertiesName
Default name of the turbulence properties dictionary.
virtual void calculateTurbulenceFields(const turbulenceModel &turbulence, scalarField &G0, scalarField &epsilon0)
Main driver to calculate the turbulence fields.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual const labelUList & faceCells() const
Return faceCells.
Definition: fvPatch.C:93
static word groupName(Name name, const word &group)
void setValues(const labelUList &cells, const UList< Type > &values)
Set solution in given cells to the specified values.
Definition: fvMatrix.C:486
scalar yPlusLam_
y+ at the edge of the laminar sublayer
Foam::fvPatchFieldMapper.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:292
scalarField epsilon_
Local copy of turbulence epsilon field.
const nearWallDist & y() const
Return the near wall distances.
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:72
virtual void manipulateMatrix(fvMatrix< Type > &matrix)
Manipulate matrix.
Definition: fvPatchField.C:344
This boundary condition provides a turbulence dissipation wall constraint for low- and high-Reynolds ...
errorManip< error > abort(error &err)
Definition: errorManip.H:131
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
scalarField & G(bool init=false)
Return non-const access to the master&#39;s G field.
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
virtual label size() const
Return size.
Definition: fvPatch.H:161
virtual epsilonWallFunctionFvPatchScalarField & epsilonPatch(const label patchi)
Helper function to return non-const access to an epsilon patch.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
static const char nl
Definition: Ostream.H:265
virtual tmp< Field< Type > > patchInternalField() const
Return internal field next to patch as patch field.
Definition: fvPatchField.C:224
const Mesh & mesh() const
Return mesh.
Ostream & writeKeyword(const keyType &)
Write the keyword followed by an appropriate indentation.
Definition: Ostream.C:54
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
void setSize(const label)
Reset size of List.
Definition: List.C:281
label patchi
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
List< List< scalar > > cornerWeights_
List of averaging corner weights.
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
virtual void updateWeightedCoeffs(const scalarField &weights)
Update the coefficients associated with the patch field.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
scalar yPlus
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
Definition: fvPatchField.C:311
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
const objectRegistry & db() const
Return local objectRegistry.
Definition: fvPatchField.C:197
virtual void checkType()
Check the type of the patch.
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual void setMaster()
Set the master patch - master is responsible for updating all.
virtual tmp< volScalarField > nut() const =0
Return the turbulence viscosity.
volScalarField & p
A class for managing temporary objects.
Definition: PtrList.H:53
virtual void createAveragingWeights()
Create the averaging weights for cells which are bounded by.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
epsilonWallFunctionFvPatchScalarField(const fvPatch &, const DimensionedField< scalar, volMesh > &)
Construct from patch and internal field.
makePatchTypeField(fvPatchScalarField, atmBoundaryLayerInletEpsilonFvPatchScalarField)
label size() const
Return the number of elements in the UList.
Definition: ListI.H:170
Namespace for OpenFOAM.
const DimensionedField< Type, volMesh > & internalField() const
Return dimensioned internal field reference.
Definition: fvPatchField.H:347
scalarField & epsilon(bool init=false)
Return non-const access to the master&#39;s epsilon field.