All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-2020 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 "momentumTransportModel.H"
29 #include "fvMatrix.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
35 
36 // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
37 
39 {
40  if (master_ != -1)
41  {
42  return;
43  }
44 
45  const volScalarField& epsilon =
46  static_cast<const volScalarField&>(this->internalField());
47 
48  const volScalarField::Boundary& bf = epsilon.boundaryField();
49 
50  label master = -1;
51  forAll(bf, patchi)
52  {
53  if (isA<epsilonWallFunctionFvPatchScalarField>(bf[patchi]))
54  {
56 
57  if (master == -1)
58  {
59  master = patchi;
60  }
61 
62  epf.master() = master;
63  }
64  }
65 }
66 
67 
69 {
70  const volScalarField& epsilon =
71  static_cast<const volScalarField&>(this->internalField());
72 
73  const volScalarField::Boundary& bf = epsilon.boundaryField();
74 
75  const fvMesh& mesh = epsilon.mesh();
76 
77  if (initialised_ && !mesh.changing())
78  {
79  return;
80  }
81 
82  volScalarField weights
83  (
84  IOobject
85  (
86  "weights",
87  mesh.time().timeName(),
88  mesh,
91  false // do not register
92  ),
93  mesh,
95  );
96 
97  DynamicList<label> epsilonPatches(bf.size());
98  forAll(bf, patchi)
99  {
100  if (isA<epsilonWallFunctionFvPatchScalarField>(bf[patchi]))
101  {
102  epsilonPatches.append(patchi);
103 
104  const labelUList& faceCells = bf[patchi].patch().faceCells();
105  forAll(faceCells, i)
106  {
107  weights[faceCells[i]]++;
108  }
109  }
110  }
111 
112  cornerWeights_.setSize(bf.size());
113  forAll(epsilonPatches, i)
114  {
115  label patchi = epsilonPatches[i];
116  const fvPatchScalarField& wf = weights.boundaryField()[patchi];
118  }
119 
120  G_.setSize(internalField().size(), 0.0);
122 
123  initialised_ = true;
124 }
125 
126 
129 {
130  const volScalarField& epsilon =
131  static_cast<const volScalarField&>(this->internalField());
132 
133  const volScalarField::Boundary& bf = epsilon.boundaryField();
134 
136  refCast<const epsilonWallFunctionFvPatchScalarField>(bf[patchi]);
137 
138  return const_cast<epsilonWallFunctionFvPatchScalarField&>(epf);
139 }
140 
141 
143 (
144  const momentumTransportModel& turbulence,
145  scalarField& G0,
146  scalarField& epsilon0
147 )
148 {
149  // Accumulate all of the G and epsilon contributions
151  {
152  if (!cornerWeights_[patchi].empty())
153  {
155 
156  const List<scalar>& w = cornerWeights_[patchi];
157 
158  epf.calculate(turbulence, w, epf.patch(), G0, epsilon0);
159  }
160  }
161 
162  // Apply zero-gradient condition for epsilon
164  {
165  if (!cornerWeights_[patchi].empty())
166  {
168 
169  epf == scalarField(epsilon0, epf.patch().faceCells());
170  }
171  }
172 }
173 
174 
176 (
177  const momentumTransportModel& turbModel,
178  const List<scalar>& cornerWeights,
179  const fvPatch& patch,
180  scalarField& G0,
181  scalarField& epsilon0
182 )
183 {
184  const label patchi = patch.index();
185 
187  nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi);
188 
189  const scalarField& y = turbModel.y()[patchi];
190 
191  const tmp<scalarField> tnuw = turbModel.nu(patchi);
192  const scalarField& nuw = tnuw();
193 
194  const tmp<volScalarField> tk = turbModel.k();
195  const volScalarField& k = tk();
196 
197  const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
198 
199  const scalarField magGradUw(mag(Uw.snGrad()));
200 
201  const scalar Cmu25 = pow025(nutw.Cmu());
202  const scalar Cmu75 = pow(nutw.Cmu(), 0.75);
203 
204  // Set epsilon and G
205  forAll(nutw, facei)
206  {
207  const label celli = patch.faceCells()[facei];
208 
209  const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei];
210 
211  const scalar w = cornerWeights[facei];
212 
213  if (yPlus > nutw.yPlusLam())
214  {
215  epsilon0[celli] +=
216  w*Cmu75*pow(k[celli], 1.5)/(nutw.kappa()*y[facei]);
217 
218  G0[celli] +=
219  w
220  *(nutw[facei] + nuw[facei])
221  *magGradUw[facei]
222  *Cmu25*sqrt(k[celli])
223  /(nutw.kappa()*y[facei]);
224  }
225  else
226  {
227  epsilon0[celli] += w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
228  }
229  }
230 }
231 
232 
233 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
234 
237 (
238  const fvPatch& p,
240 )
241 :
243  G_(),
244  epsilon_(),
245  initialised_(false),
246  master_(-1),
248 {}
249 
250 
253 (
254  const fvPatch& p,
256  const dictionary& dict
257 )
258 :
260  G_(),
261  epsilon_(),
262  initialised_(false),
263  master_(-1),
265 {
266  // Apply zero-gradient condition on start-up
268 }
269 
270 
273 (
275  const fvPatch& p,
277  const fvPatchFieldMapper& mapper
278 )
279 :
280  fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
281  G_(),
282  epsilon_(),
283  initialised_(false),
284  master_(-1),
286 {}
287 
288 
291 (
293 )
294 :
296  G_(),
297  epsilon_(),
298  initialised_(false),
299  master_(-1),
301 {}
302 
303 
306 (
309 )
310 :
311  fixedValueFvPatchField<scalar>(ewfpsf, iF),
312  G_(),
313  epsilon_(),
314  initialised_(false),
315  master_(-1),
317 {}
318 
319 
320 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
321 
323 {
324  if (patch().index() == master_)
325  {
326  if (init)
327  {
328  G_ = 0.0;
329  }
330 
331  return G_;
332  }
333 
334  return epsilonPatch(master_).G();
335 }
336 
337 
339 (
340  bool init
341 )
342 {
343  if (patch().index() == master_)
344  {
345  if (init)
346  {
347  epsilon_ = 0.0;
348  }
349 
350  return epsilon_;
351  }
352 
353  return epsilonPatch(master_).epsilon(init);
354 }
355 
356 
358 {
359  if (updated())
360  {
361  return;
362  }
363 
364  const momentumTransportModel& turbModel =
366  (
368  (
369  momentumTransportModel::typeName,
370  internalField().group()
371  )
372  );
373 
374  setMaster();
375 
376  if (patch().index() == master_)
377  {
379  calculateTurbulenceFields(turbModel, G(true), epsilon(true));
380  }
381 
382  const scalarField& G0 = this->G();
383  const scalarField& epsilon0 = this->epsilon();
384 
385  typedef DimensionedField<scalar, volMesh> FieldType;
386 
387  FieldType& G =
388  const_cast<FieldType&>
389  (
390  db().lookupObject<FieldType>(turbModel.GName())
391  );
392 
393  FieldType& epsilon = const_cast<FieldType&>(internalField());
394 
395  forAll(*this, facei)
396  {
397  label celli = patch().faceCells()[facei];
398 
399  G[celli] = G0[celli];
400  epsilon[celli] = epsilon0[celli];
401  }
402 
404 }
405 
406 
408 (
409  const scalarField& weights
410 )
411 {
412  if (updated())
413  {
414  return;
415  }
416 
417  const momentumTransportModel& turbModel =
419  (
421  (
422  momentumTransportModel::typeName,
423  internalField().group()
424  )
425  );
426 
427  setMaster();
428 
429  if (patch().index() == master_)
430  {
432  calculateTurbulenceFields(turbModel, G(true), epsilon(true));
433  }
434 
435  const scalarField& G0 = this->G();
436  const scalarField& epsilon0 = this->epsilon();
437 
438  typedef DimensionedField<scalar, volMesh> FieldType;
439 
440  FieldType& G =
441  const_cast<FieldType&>
442  (
443  db().lookupObject<FieldType>(turbModel.GName())
444  );
445 
446  FieldType& epsilon = const_cast<FieldType&>(internalField());
447 
448  scalarField& epsilonf = *this;
449 
450  // Only set the values if the weights are > tolerance
451  forAll(weights, facei)
452  {
453  scalar w = weights[facei];
454 
455  if (w > tolerance_)
456  {
457  label celli = patch().faceCells()[facei];
458 
459  G[celli] = (1.0 - w)*G[celli] + w*G0[celli];
460  epsilon[celli] = (1.0 - w)*epsilon[celli] + w*epsilon0[celli];
461  epsilonf[facei] = epsilon[celli];
462  }
463  }
464 
466 }
467 
468 
470 (
471  fvMatrix<scalar>& matrix
472 )
473 {
474  if (manipulatedMatrix())
475  {
476  return;
477  }
478 
479  matrix.setValues(patch().faceCells(), patchInternalField());
480 
482 }
483 
484 
486 (
487  fvMatrix<scalar>& matrix,
488  const Field<scalar>& weights
489 )
490 {
491  if (manipulatedMatrix())
492  {
493  return;
494  }
495 
496  DynamicList<label> constraintCells(weights.size());
497  DynamicList<scalar> constraintEpsilon(weights.size());
498  const labelUList& faceCells = patch().faceCells();
499 
501  = internalField();
502 
503  label nConstrainedCells = 0;
504 
505 
506  forAll(weights, facei)
507  {
508  // Only set the values if the weights are > tolerance
509  if (weights[facei] > tolerance_)
510  {
511  nConstrainedCells++;
512 
513  label celli = faceCells[facei];
514 
515  constraintCells.append(celli);
516  constraintEpsilon.append(epsilon[celli]);
517  }
518  }
519 
520  if (debug)
521  {
522  Pout<< "Patch: " << patch().name()
523  << ": number of constrained cells = " << nConstrainedCells
524  << " out of " << patch().size()
525  << endl;
526  }
527 
528  matrix.setValues
529  (
530  constraintCells,
531  scalarField(constraintEpsilon)
532  );
533 
535 }
536 
537 
538 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
539 
540 namespace Foam
541 {
543  (
546  );
547 }
548 
549 
550 // ************************************************************************* //
virtual tmp< volScalarField > nu() const =0
Return the laminar viscosity.
const char *const group
Group name for atomic constants.
bool changing() const
Is mesh changing (topology changing and/or moving)
Definition: polyMesh.H:540
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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 tmp< Field< Type > > snGrad() const
Return patch-normal gradient.
Definition: fvPatchField.C:186
label index() const
Return the index of this patch in the fvBoundaryMesh.
Definition: fvPatch.H:173
virtual void calculate(const momentumTransportModel &turbulence, const List< scalar > &cornerWeights, const fvPatch &patch, scalarField &G, scalarField &epsilon)
Calculate the epsilon and G.
const dimensionedScalar & epsilon0
Electric constant: default SI units: [F/m].
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
bool updated() const
Return true if the boundary condition has already been updated.
Definition: fvPatchField.H:356
const Boundary & boundaryField() const
Return const-reference to the boundary field.
virtual void calculateTurbulenceFields(const momentumTransportModel &turbulence, scalarField &G0, scalarField &epsilon0)
Main driver to calculate the turbulence fields.
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:164
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
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
const nearWallDist & y() const
Return the near wall distances.
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:622
label k
Boltzmann constant.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
bool manipulatedMatrix() const
Return true if the matrix has already been manipulated.
Definition: fvPatchField.H:362
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type.
Macros for easy insertion into run-time selection tables.
virtual void operator==(const fvPatchField< Type > &)
Definition: fvPatchField.C:461
static const nutWallFunctionFvPatchScalarField & nutw(const momentumTransportModel &turbModel, const label patchi)
Return the nut patchField for the given wall patch.
scalar y
This boundary condition provides a turbulent kinematic viscosity condition when using wall functions...
dynamicFvMesh & mesh
virtual label & master()
Return non-const access to the master patch ID.
const dimensionedScalar & e
Elementary charge.
Definition: doubleScalar.H:105
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:484
Foam::fvPatchFieldMapper.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
scalarField epsilon_
Local copy of turbulence epsilon field.
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:262
This boundary condition provides a turbulence dissipation wall constraint for low- and high-Reynolds ...
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
const fvPatch & patch() const
Return patch.
Definition: fvPatchField.H:325
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:155
const dimensionedScalar & G0
Conductance quantum: default SI units: [S].
virtual epsilonWallFunctionFvPatchScalarField & epsilonPatch(const label patchi)
Helper function to return non-const access to an epsilon patch.
virtual tmp< Field< Type > > patchInternalField() const
Return internal field next to patch as patch field.
Definition: fvPatchField.C:194
const Mesh & mesh() const
Return mesh.
virtual tmp< volScalarField > k() const =0
Return the turbulence kinetic energy.
Abstract base class for turbulence models (RAS, LES and laminar).
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
void setSize(const label)
Reset size of List.
Definition: List.C:281
const volVectorField & U() const
Access function to velocity field.
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.
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:229
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:167
dimensioned< scalar > mag(const dimensioned< Type > &)
word GName() const
Helper function to return the name of the turbulence G field.
virtual const word & name() const
Return name.
Definition: fvPatch.H:143
virtual void setMaster()
Set the master patch - master is responsible for updating all.
static scalar yPlusLam(const scalar kappa, const scalar E)
Calculate the Y+ at the edge of the laminar sublayer.
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:171
Namespace for OpenFOAM.
const DimensionedField< Type, volMesh > & internalField() const
Return dimensioned internal field reference.
Definition: fvPatchField.H:332
scalarField & epsilon(bool init=false)
Return non-const access to the master&#39;s epsilon field.