omegaWallFunctionFvPatchScalarField.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 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 scalar omegaWallFunctionFvPatchScalarField::tolerance_ = 1e-5;
40 
41 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
42 
44 {
45  if (master_ != -1)
46  {
47  return;
48  }
49 
50  const volScalarField& omega =
51  static_cast<const volScalarField&>(this->internalField());
52 
53  const volScalarField::Boundary& bf = omega.boundaryField();
54 
55  label master = -1;
56  forAll(bf, patchi)
57  {
58  if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchi]))
59  {
61 
62  if (master == -1)
63  {
64  master = patchi;
65  }
66 
67  opf.master() = master;
68  }
69  }
70 }
71 
72 
74 {
75  const volScalarField& omega =
76  static_cast<const volScalarField&>(this->internalField());
77 
78  const volScalarField::Boundary& bf = omega.boundaryField();
79 
80  const fvMesh& mesh = omega.mesh();
81 
82  if (initialised_ && !mesh.changing())
83  {
84  return;
85  }
86 
87  volScalarField weights
88  (
89  IOobject
90  (
91  "weights",
92  mesh.time().timeName(),
93  mesh,
96  false // do not register
97  ),
98  mesh,
100  );
101 
102  DynamicList<label> omegaPatches(bf.size());
103  forAll(bf, patchi)
104  {
105  if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchi]))
106  {
107  omegaPatches.append(patchi);
108 
109  const labelUList& faceCells = bf[patchi].patch().faceCells();
110  forAll(faceCells, i)
111  {
112  label celli = faceCells[i];
113  weights[celli]++;
114  }
115  }
116  }
117 
118  cornerWeights_.setSize(bf.size());
119  forAll(omegaPatches, i)
120  {
121  label patchi = omegaPatches[i];
122  const fvPatchScalarField& wf = weights.boundaryField()[patchi];
124  }
125 
126  G_.setSize(internalField().size(), 0.0);
127  omega_.setSize(internalField().size(), 0.0);
128 
129  initialised_ = true;
130 }
131 
132 
135 {
136  const volScalarField& omega =
137  static_cast<const volScalarField&>(this->internalField());
138 
139  const volScalarField::Boundary& bf = omega.boundaryField();
140 
142  refCast<const omegaWallFunctionFvPatchScalarField>(bf[patchi]);
143 
144  return const_cast<omegaWallFunctionFvPatchScalarField&>(opf);
145 }
146 
147 
149 (
150  const momentumTransportModel& turbModel,
151  scalarField& G0,
152  scalarField& omega0
153 )
154 {
155  // accumulate all of the G and omega contributions
157  {
158  if (!cornerWeights_[patchi].empty())
159  {
161 
162  const List<scalar>& w = cornerWeights_[patchi];
163 
164  opf.calculate(turbModel, w, opf.patch(), G0, omega0);
165  }
166  }
167 
168  // apply zero-gradient condition for omega
170  {
171  if (!cornerWeights_[patchi].empty())
172  {
174 
175  opf == scalarField(omega0, opf.patch().faceCells());
176  }
177  }
178 }
179 
180 
182 (
183  const momentumTransportModel& turbModel,
184  const List<scalar>& cornerWeights,
185  const fvPatch& patch,
186  scalarField& G0,
187  scalarField& omega0
188 )
189 {
190  const label patchi = patch.index();
191 
193  nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi);
194 
195  const scalarField& y = turbModel.y()[patchi];
196 
197  const tmp<volScalarField> tk = turbModel.k();
198  const volScalarField& k = tk();
199 
200  const tmp<scalarField> tnuw = turbModel.nu(patchi);
201  const scalarField& nuw = tnuw();
202 
203  const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
204 
205  const scalarField magGradUw(mag(Uw.snGrad()));
206 
207  const FieldType& G =
208  db().lookupObject<FieldType>(turbModel.GName());
209 
210  const scalar Cmu25 = pow025(nutw.Cmu());
211  const scalar Cmu5 = sqrt(nutw.Cmu());
212 
213  // Set omega and G
214  forAll(nutw, facei)
215  {
216  const label celli = patch.faceCells()[facei];
217  const scalar w = cornerWeights[facei];
218 
219  const scalar Rey = y[facei]*sqrt(k[celli])/nuw[facei];
220  const scalar yPlus = Cmu25*Rey;
221  const scalar uPlus = (1/nutw.kappa())*log(nutw.E()*yPlus);
222 
223  if (blended_)
224  {
225  const scalar lamFrac = exp(-Rey/11);
226  const scalar turbFrac = 1 - lamFrac;
227 
228  const scalar uStar = sqrt
229  (
230  lamFrac*nuw[facei]*magGradUw[facei] + turbFrac*Cmu5*k[celli]
231  );
232 
233  const scalar omegaVis = 6*nuw[facei]/(beta1_*sqr(y[facei]));
234  const scalar omegaLog = uStar/(Cmu5*nutw.kappa()*y[facei]);
235 
236  omega0[celli] += w*(lamFrac*omegaVis + turbFrac*omegaLog);
237 
238  G0[celli] +=
239  w
240  *(
241  lamFrac*G[celli]
242 
243  + turbFrac
244  *sqr(uStar*magGradUw[facei]*y[facei]/uPlus)
245  /(nuw[facei]*nutw.kappa()*yPlus)
246  );
247  }
248  else
249  {
250  if (yPlus < nutw.yPlusLam())
251  {
252  const scalar omegaVis = 6*nuw[facei]/(beta1_*sqr(y[facei]));
253 
254  omega0[celli] += w*omegaVis;
255 
256  G0[celli] += w*G[celli];
257  }
258  else
259  {
260  const scalar uStar = sqrt(Cmu5*k[celli]);
261  const scalar omegaLog = uStar/(Cmu5*nutw.kappa()*y[facei]);
262 
263  omega0[celli] += w*omegaLog;
264 
265  G0[celli] +=
266  w*
267  sqr(uStar*magGradUw[facei]*y[facei]/uPlus)
268  /(nuw[facei]*nutw.kappa()*yPlus);
269  }
270  }
271  }
272 }
273 
274 
275 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
276 
278 (
279  const fvPatch& p,
281 )
282 :
284  beta1_(0.075),
285  blended_(false),
286  G_(),
287  omega_(),
288  initialised_(false),
289  master_(-1),
291 {}
292 
293 
295 (
296  const fvPatch& p,
298  const dictionary& dict
299 )
300 :
302  beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075)),
303  blended_(dict.lookupOrDefault<Switch>("blended", false)),
304  G_(),
305  omega_(),
306  initialised_(false),
307  master_(-1),
309 {
310  // apply zero-gradient condition on start-up
312 }
313 
314 
316 (
318  const fvPatch& p,
320  const fvPatchFieldMapper& mapper
321 )
322 :
323  fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
324  beta1_(ptf.beta1_),
325  blended_(ptf.blended_),
326  G_(),
327  omega_(),
328  initialised_(false),
329  master_(-1),
331 {}
332 
333 
335 (
338 )
339 :
340  fixedValueFvPatchField<scalar>(owfpsf, iF),
341  beta1_(owfpsf.beta1_),
342  blended_(owfpsf.blended_),
343  G_(),
344  omega_(),
345  initialised_(false),
346  master_(-1),
348 {}
349 
350 
351 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
352 
354 {
355  if (patch().index() == master_)
356  {
357  if (init)
358  {
359  G_ = 0.0;
360  }
361 
362  return G_;
363  }
364 
365  return omegaPatch(master_).G();
366 }
367 
368 
370 {
371  if (patch().index() == master_)
372  {
373  if (init)
374  {
375  omega_ = 0.0;
376  }
377 
378  return omega_;
379  }
380 
381  return omegaPatch(master_).omega(init);
382 }
383 
384 
386 {
387  if (updated())
388  {
389  return;
390  }
391 
392  const momentumTransportModel& turbModel =
394  (
396  (
397  momentumTransportModel::typeName,
398  internalField().group()
399  )
400  );
401 
402  setMaster();
403 
404  if (patch().index() == master_)
405  {
407  calculateTurbulenceFields(turbModel, G(true), omega(true));
408  }
409 
410  const scalarField& G0 = this->G();
411  const scalarField& omega0 = this->omega();
412 
413  FieldType& G =
414  const_cast<FieldType&>
415  (
416  db().lookupObject<FieldType>(turbModel.GName())
417  );
418 
419  FieldType& omega = const_cast<FieldType&>(internalField());
420 
421  forAll(*this, facei)
422  {
423  label celli = patch().faceCells()[facei];
424 
425  G[celli] = G0[celli];
426  omega[celli] = omega0[celli];
427  }
428 
430 }
431 
432 
434 (
435  const scalarField& weights
436 )
437 {
438  if (updated())
439  {
440  return;
441  }
442 
443  const momentumTransportModel& turbModel =
445  (
447  (
448  momentumTransportModel::typeName,
449  internalField().group()
450  )
451  );
452 
453  setMaster();
454 
455  if (patch().index() == master_)
456  {
458  calculateTurbulenceFields(turbModel, G(true), omega(true));
459  }
460 
461  const scalarField& G0 = this->G();
462  const scalarField& omega0 = this->omega();
463 
464  FieldType& G =
465  const_cast<FieldType&>
466  (
467  db().lookupObject<FieldType>(turbModel.GName())
468  );
469 
470  FieldType& omega = const_cast<FieldType&>(internalField());
471 
472  scalarField& omegaf = *this;
473 
474  // only set the values if the weights are > tolerance
475  forAll(weights, facei)
476  {
477  scalar w = weights[facei];
478 
479  if (w > tolerance_)
480  {
481  label celli = patch().faceCells()[facei];
482 
483  G[celli] = (1.0 - w)*G[celli] + w*G0[celli];
484  omega[celli] = (1.0 - w)*omega[celli] + w*omega0[celli];
485  omegaf[facei] = omega[celli];
486  }
487  }
488 
490 }
491 
492 
494 (
495  fvMatrix<scalar>& matrix
496 )
497 {
498  if (manipulatedMatrix())
499  {
500  return;
501  }
502 
503  matrix.setValues(patch().faceCells(), patchInternalField());
504 
506 }
507 
508 
510 (
511  fvMatrix<scalar>& matrix,
512  const Field<scalar>& weights
513 )
514 {
515  if (manipulatedMatrix())
516  {
517  return;
518  }
519 
520  DynamicList<label> constraintCells(weights.size());
521  DynamicList<scalar> constraintomega(weights.size());
522  const labelUList& faceCells = patch().faceCells();
523 
525  = internalField();
526 
527  label nConstrainedCells = 0;
528 
529 
530  forAll(weights, facei)
531  {
532  // only set the values if the weights are > tolerance
533  if (weights[facei] > tolerance_)
534  {
535  nConstrainedCells++;
536 
537  label celli = faceCells[facei];
538 
539  constraintCells.append(celli);
540  constraintomega.append(omega[celli]);
541  }
542  }
543 
544  if (debug)
545  {
546  Pout<< "Patch: " << patch().name()
547  << ": number of constrained cells = " << nConstrainedCells
548  << " out of " << patch().size()
549  << endl;
550  }
551 
552  matrix.setValues
553  (
554  constraintCells,
555  scalarField(constraintomega)
556  );
557 
559 }
560 
561 
563 {
564  writeEntry(os, "beta1", beta1_);
565  writeEntry(os, "blended", blended_);
567 }
568 
569 
570 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
571 
573 (
576 );
577 
578 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
579 
580 } // End namespace Foam
581 
582 // ************************************************************************* //
virtual tmp< volScalarField > nu() const =0
Return the laminar viscosity.
const char *const group
Group name for atomic constants.
scalar uPlus
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:325
Switch blended_
Blending switch (defaults to false)
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 updateCoeffs()
Update the coefficients associated with the patch field.
virtual void write(Ostream &) const
Write.
virtual tmp< Field< Type > > snGrad() const
Return patch-normal gradient.
Definition: fvPatchField.C:166
scalar Rey
dimensionedScalar log(const dimensionedScalar &ds)
label index() const
Return the index of this patch in the fvBoundaryMesh.
Definition: fvPatch.H:174
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
omegaWallFunctionFvPatchScalarField(const fvPatch &, const DimensionedField< scalar, volMesh > &)
Construct from patch and internal field.
bool updated() const
Return true if the boundary condition has already been updated.
Definition: fvPatchField.H:355
const Boundary & boundaryField() const
Return const-reference to the boundary field.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
virtual void updateWeightedCoeffs(const scalarField &weights)
Update the coefficients associated with the patch field.
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:62
dimensionedScalar pow025(const dimensionedScalar &ds)
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none/any.
Definition: Switch.H:60
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:66
virtual void createAveragingWeights()
Create the averaging weights for cells which are bounded by.
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:636
const dimensionSet dimless
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:361
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:441
scalarField & G(bool init=false)
Return non-const access to the master&#39;s G field.
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 omegaWallFunctionFvPatchScalarField & omegaPatch(const label patchi)
Helper function to return non-const access to an omega patch.
dimensionedScalar exp(const dimensionedScalar &ds)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual const labelUList & faceCells() const
Return faceCells.
Definition: fvPatch.C:99
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
static scalar tolerance_
Tolerance used in weighted calculations.
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:242
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:337
virtual label size() const
Return size.
Definition: fvPatch.H:156
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
virtual tmp< Field< Type > > patchInternalField() const
Return internal field next to patch as patch field.
Definition: fvPatchField.C:174
const Mesh & mesh() const
Return mesh.
List< List< scalar > > cornerWeights_
List of averaging corner weights.
void writeEntry(Ostream &os, const HashTable< T, Key, Hash > &ht)
Definition: HashTableIO.C:96
virtual void calculateTurbulenceFields(const momentumTransportModel &turbulence, scalarField &G0, scalarField &omega0)
Main driver to calculate the turbulence fields.
scalarField G_
Local copy of turbulence G field.
virtual tmp< volScalarField > k() const =0
Return the turbulence kinetic energy.
Abstract base class for turbulence models (RAS, LES and laminar).
scalarField & omega(bool init=false)
Return non-const access to the master&#39;s omega field.
void setSize(const label)
Reset size of List.
Definition: List.C:281
const volVectorField & U() const
Access function to velocity field.
label patchi
This boundary condition provides a wall constraint on turbulnce specific dissipation, omega for both low and high Reynolds number turbulence models.
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const dimensionedScalar G0
Conductance quantum: default SI units: [S].
scalar yPlus
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
Definition: fvPatchField.C:209
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:147
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual void setMaster()
Set the master patch - master is responsible for updating all.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
word GName() const
Helper function to return the name of the turbulence G field.
virtual const word & name() const
Return name.
Definition: fvPatch.H:144
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 label & master()
Return non-const access to the master patch ID.
virtual void calculate(const momentumTransportModel &turbulence, const List< scalar > &cornerWeights, const fvPatch &patch, scalarField &G, scalarField &omega)
Calculate the omega and G.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
virtual void manipulateMatrix(fvMatrix< scalar > &matrix)
Manipulate matrix.
makePatchTypeField(fvPatchScalarField, atmBoundaryLayerInletEpsilonFvPatchScalarField)
label size() const
Return the number of elements in the UList.
Definition: ListI.H:171
scalarField omega_
Local copy of turbulence omega field.
Namespace for OpenFOAM.
const DimensionedField< Type, volMesh > & internalField() const
Return dimensioned internal field reference.
Definition: fvPatchField.H:343