omegaWallFunctionFvPatchScalarField.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-2015 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 
27 #include "turbulenceModel.H"
28 #include "fvPatchFieldMapper.H"
29 #include "fvMatrix.H"
30 #include "volFields.H"
31 #include "wallFvPatch.H"
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 
42 scalar omegaWallFunctionFvPatchScalarField::tolerance_ = 1e-5;
43 
44 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
45 
47 {
48  if (!isA<wallFvPatch>(patch()))
49  {
50  FatalErrorIn("omegaWallFunctionFvPatchScalarField::checkType()")
51  << "Invalid wall function specification" << nl
52  << " Patch type for patch " << patch().name()
53  << " must be wall" << nl
54  << " Current patch type is " << patch().type() << nl << endl
55  << abort(FatalError);
56  }
57 }
58 
59 
61 {
62  os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
63  os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
64  os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
65  os.writeKeyword("beta1") << beta1_ << token::END_STATEMENT << nl;
66 }
67 
68 
70 {
71  if (master_ != -1)
72  {
73  return;
74  }
75 
76  const volScalarField& omega =
77  static_cast<const volScalarField&>(this->dimensionedInternalField());
78 
79  const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
80 
81  label master = -1;
82  forAll(bf, patchi)
83  {
84  if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchi]))
85  {
87 
88  if (master == -1)
89  {
90  master = patchi;
91  }
92 
93  opf.master() = master;
94  }
95  }
96 }
97 
98 
100 {
101  const volScalarField& omega =
102  static_cast<const volScalarField&>(this->dimensionedInternalField());
103 
104  const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
105 
106  const fvMesh& mesh = omega.mesh();
107 
108  if (initialised_ && !mesh.changing())
109  {
110  return;
111  }
112 
113  volScalarField weights
114  (
115  IOobject
116  (
117  "weights",
118  mesh.time().timeName(),
119  mesh,
122  false // do not register
123  ),
124  mesh,
125  dimensionedScalar("zero", dimless, 0.0)
126  );
127 
128  DynamicList<label> omegaPatches(bf.size());
129  forAll(bf, patchi)
130  {
131  if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchi]))
132  {
133  omegaPatches.append(patchi);
134 
135  const labelUList& faceCells = bf[patchi].patch().faceCells();
136  forAll(faceCells, i)
137  {
138  label cellI = faceCells[i];
139  weights[cellI]++;
140  }
141  }
142  }
143 
144  cornerWeights_.setSize(bf.size());
145  forAll(omegaPatches, i)
146  {
147  label patchi = omegaPatches[i];
148  const fvPatchScalarField& wf = weights.boundaryField()[patchi];
150  }
151 
154 
155  initialised_ = true;
156 }
157 
158 
161 {
162  const volScalarField& omega =
163  static_cast<const volScalarField&>(this->dimensionedInternalField());
164 
165  const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
166 
168  refCast<const omegaWallFunctionFvPatchScalarField>(bf[patchi]);
169 
170  return const_cast<omegaWallFunctionFvPatchScalarField&>(opf);
171 }
172 
173 
175 (
176  const turbulenceModel& turbModel,
177  scalarField& G0,
178  scalarField& omega0
179 )
180 {
181  // accumulate all of the G and omega contributions
183  {
184  if (!cornerWeights_[patchi].empty())
185  {
187 
188  const List<scalar>& w = cornerWeights_[patchi];
189 
190  opf.calculate(turbModel, w, opf.patch(), G0, omega0);
191  }
192  }
193 
194  // apply zero-gradient condition for omega
196  {
197  if (!cornerWeights_[patchi].empty())
198  {
200 
201  opf == scalarField(omega0, opf.patch().faceCells());
202  }
203  }
204 }
205 
206 
208 (
209  const turbulenceModel& turbModel,
210  const List<scalar>& cornerWeights,
211  const fvPatch& patch,
212  scalarField& G,
213  scalarField& omega
214 )
215 {
216  const label patchi = patch.index();
217 
218  const scalarField& y = turbModel.y()[patchi];
219 
220  const scalar Cmu25 = pow025(Cmu_);
221 
222  const tmp<volScalarField> tk = turbModel.k();
223  const volScalarField& k = tk();
224 
225  const tmp<scalarField> tnuw = turbModel.nu(patchi);
226  const scalarField& nuw = tnuw();
227 
228  const tmp<scalarField> tnutw = turbModel.nut(patchi);
229  const scalarField& nutw = tnutw();
230 
231  const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
232 
233  const scalarField magGradUw(mag(Uw.snGrad()));
234 
235  // Set omega and G
236  forAll(nutw, faceI)
237  {
238  label cellI = patch.faceCells()[faceI];
239 
240  scalar w = cornerWeights[faceI];
241 
242  scalar omegaVis = 6.0*nuw[faceI]/(beta1_*sqr(y[faceI]));
243 
244  scalar omegaLog = sqrt(k[cellI])/(Cmu25*kappa_*y[faceI]);
245 
246  omega[cellI] += w*sqrt(sqr(omegaVis) + sqr(omegaLog));
247 
248  G[cellI] +=
249  w
250  *(nutw[faceI] + nuw[faceI])
251  *magGradUw[faceI]
252  *Cmu25*sqrt(k[cellI])
253  /(kappa_*y[faceI]);
254  }
255 }
256 
257 
258 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
259 
261 (
262  const fvPatch& p,
264 )
265 :
267  Cmu_(0.09),
268  kappa_(0.41),
269  E_(9.8),
270  beta1_(0.075),
272  G_(),
273  omega_(),
274  initialised_(false),
275  master_(-1),
277 {
278  checkType();
279 }
280 
281 
283 (
285  const fvPatch& p,
287  const fvPatchFieldMapper& mapper
288 )
289 :
290  fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
291  Cmu_(ptf.Cmu_),
292  kappa_(ptf.kappa_),
293  E_(ptf.E_),
294  beta1_(ptf.beta1_),
295  yPlusLam_(ptf.yPlusLam_),
296  G_(),
297  omega_(),
298  initialised_(false),
299  master_(-1),
301 {
302  checkType();
303 }
304 
305 
307 (
308  const fvPatch& p,
310  const dictionary& dict
311 )
312 :
314  Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
315  kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
316  E_(dict.lookupOrDefault<scalar>("E", 9.8)),
317  beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075)),
319  G_(),
320  omega_(),
321  initialised_(false),
322  master_(-1),
324 {
325  checkType();
326 
327  // apply zero-gradient condition on start-up
329 }
330 
331 
333 (
335 )
336 :
338  Cmu_(owfpsf.Cmu_),
339  kappa_(owfpsf.kappa_),
340  E_(owfpsf.E_),
341  beta1_(owfpsf.beta1_),
342  yPlusLam_(owfpsf.yPlusLam_),
343  G_(),
344  omega_(),
345  initialised_(false),
346  master_(-1),
348 {
349  checkType();
350 }
351 
352 
354 (
357 )
358 :
359  fixedValueFvPatchField<scalar>(owfpsf, iF),
360  Cmu_(owfpsf.Cmu_),
361  kappa_(owfpsf.kappa_),
362  E_(owfpsf.E_),
363  beta1_(owfpsf.beta1_),
364  yPlusLam_(owfpsf.yPlusLam_),
365  G_(),
366  omega_(),
367  initialised_(false),
368  master_(-1),
370 {
371  checkType();
372 }
373 
374 
375 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
376 
378 {
379  if (patch().index() == master_)
380  {
381  if (init)
382  {
383  G_ = 0.0;
384  }
385 
386  return G_;
387  }
388 
389  return omegaPatch(master_).G();
390 }
391 
392 
394 {
395  if (patch().index() == master_)
396  {
397  if (init)
398  {
399  omega_ = 0.0;
400  }
401 
402  return omega_;
403  }
404 
405  return omegaPatch(master_).omega(init);
406 }
407 
408 
410 {
411  if (updated())
412  {
413  return;
414  }
415 
416  const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
417  (
419  (
422  )
423  );
424 
425  setMaster();
426 
427  if (patch().index() == master_)
428  {
430  calculateTurbulenceFields(turbModel, G(true), omega(true));
431  }
432 
433  const scalarField& G0 = this->G();
434  const scalarField& omega0 = this->omega();
435 
436  typedef DimensionedField<scalar, volMesh> FieldType;
437 
438  FieldType& G =
439  const_cast<FieldType&>
440  (
441  db().lookupObject<FieldType>(turbModel.GName())
442  );
443 
444  FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
445 
446  forAll(*this, faceI)
447  {
448  label cellI = patch().faceCells()[faceI];
449 
450  G[cellI] = G0[cellI];
451  omega[cellI] = omega0[cellI];
452  }
453 
455 }
456 
457 
459 (
460  const scalarField& weights
461 )
462 {
463  if (updated())
464  {
465  return;
466  }
467 
468  const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
469  (
471  (
474  )
475  );
476 
477  setMaster();
478 
479  if (patch().index() == master_)
480  {
482  calculateTurbulenceFields(turbModel, G(true), omega(true));
483  }
484 
485  const scalarField& G0 = this->G();
486  const scalarField& omega0 = this->omega();
487 
488  typedef DimensionedField<scalar, volMesh> FieldType;
489 
490  FieldType& G =
491  const_cast<FieldType&>
492  (
493  db().lookupObject<FieldType>(turbModel.GName())
494  );
495 
496  FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
497 
498  scalarField& omegaf = *this;
499 
500  // only set the values if the weights are > tolerance
501  forAll(weights, faceI)
502  {
503  scalar w = weights[faceI];
504 
505  if (w > tolerance_)
506  {
507  label cellI = patch().faceCells()[faceI];
508 
509  G[cellI] = (1.0 - w)*G[cellI] + w*G0[cellI];
510  omega[cellI] = (1.0 - w)*omega[cellI] + w*omega0[cellI];
511  omegaf[faceI] = omega[cellI];
512  }
513  }
514 
516 }
517 
518 
520 (
521  fvMatrix<scalar>& matrix
522 )
523 {
524  if (manipulatedMatrix())
525  {
526  return;
527  }
528 
529  matrix.setValues(patch().faceCells(), patchInternalField());
530 
532 }
533 
534 
536 (
537  fvMatrix<scalar>& matrix,
538  const Field<scalar>& weights
539 )
540 {
541  if (manipulatedMatrix())
542  {
543  return;
544  }
545 
546  DynamicList<label> constraintCells(weights.size());
547  DynamicList<scalar> constraintomega(weights.size());
548  const labelUList& faceCells = patch().faceCells();
549 
552 
553  label nConstrainedCells = 0;
554 
555 
556  forAll(weights, faceI)
557  {
558  // only set the values if the weights are > tolerance
559  if (weights[faceI] > tolerance_)
560  {
561  nConstrainedCells++;
562 
563  label cellI = faceCells[faceI];
564 
565  constraintCells.append(cellI);
566  constraintomega.append(omega[cellI]);
567  }
568  }
569 
570  if (debug)
571  {
572  Pout<< "Patch: " << patch().name()
573  << ": number of constrained cells = " << nConstrainedCells
574  << " out of " << patch().size()
575  << endl;
576  }
577 
578  matrix.setValues
579  (
580  constraintCells,
581  scalarField(constraintomega.xfer())
582  );
583 
585 }
586 
587 
589 {
590  writeLocalEntries(os);
592 }
593 
594 
595 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
596 
598 (
601 );
602 
603 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
604 
605 } // End namespace Foam
606 
607 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
virtual void manipulateMatrix(fvMatrix< Type > &matrix)
Manipulate matrix.
Definition: fvPatchField.C:336
virtual tmp< Field< Type > > snGrad() const
Return patch-normal gradient.
Definition: fvPatchField.C:207
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
scalarField omega_
Local copy of turbulence omega field.
scalarField & G(bool init=false)
Return non-const access to the master&#39;s G field.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
const dimensionedScalar G0
Conductance quantum: default SI units: [S].
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
dimensioned< scalar > mag(const dimensioned< Type > &)
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
Definition: fvPatchField.C:301
bool empty() const
Return true if the UList is empty (ie, size() is zero).
const char *const group
Group name for atomic constants.
word GName() const
Helper function to return the name of the turbulence G field.
bool changing() const
Is mesh changing (topology changing and/or moving)
Definition: polyMesh.H:521
virtual omegaWallFunctionFvPatchScalarField & omegaPatch(const label patchi)
Helper function to return non-const access to an omega patch.
virtual tmp< volScalarField > nu() const =0
Return the laminar viscosity.
virtual void calculate(const turbulenceModel &turbulence, const List< scalar > &cornerWeights, const fvPatch &patch, scalarField &G, scalarField &omega)
Calculate the omega and G.
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
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:310
Foam::fvPatchFieldMapper.
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:68
virtual const labelUList & faceCells() const
Return faceCells.
Definition: fvPatch.C:93
const Mesh & mesh() const
Return mesh.
dynamicFvMesh & mesh
bool updated() const
Return true if the boundary condition has already been updated.
Definition: fvPatchField.H:348
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
virtual void checkType()
Check the type of the patch.
Namespace for OpenFOAM.
const word & name() const
Return name.
Definition: fvPatch.H:149
virtual void createAveragingWeights()
Create the averaging weights for cells which are bounded by.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static word groupName(Name name, const word &group)
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
This boundary condition provides a wall function constraint on turbulnce specific dissipation...
dictionary dict
static const char nl
Definition: Ostream.H:260
const double e
Elementary charge.
Definition: doubleFloat.H:78
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
virtual tmp< Field< Type > > patchInternalField() const
Return internal field next to patch as patch field.
Definition: fvPatchField.C:215
volScalarField & p
Definition: createFields.H:51
virtual void writeLocalEntries(Ostream &) const
Write local wall function variables.
static scalar tolerance_
Tolerance used in weighted calculations.
virtual void manipulateMatrix(fvMatrix< scalar > &matrix)
Manipulate matrix.
static scalar yPlusLam(const scalar kappa, const scalar E)
Calculate the Y+ at the edge of the laminar sublayer.
#define forAll(list, i)
Definition: UList.H:421
label k
Boltzmann constant.
label patchi
dimensionedScalar pow025(const dimensionedScalar &ds)
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.
const nearWallDist & y() const
Return the near wall distances.
virtual void setMaster()
Set the master patch - master is responsible for updating all.
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:65
errorManip< error > abort(error &err)
Definition: errorManip.H:131
label index() const
Return the index of this patch in the fvBoundaryMesh.
Definition: fvPatch.H:179
Ostream & writeKeyword(const keyType &)
Write the keyword followed by an appropriate indentation.
Definition: Ostream.C:59
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
scalar yPlusLam_
Y+ at the edge of the laminar sublayer.
Abstract base class for turbulence models (RAS, LES and laminar).
scalarField G_
Local copy of turbulence G field.
virtual void calculateTurbulenceFields(const turbulenceModel &turbulence, scalarField &G0, scalarField &omega0)
Main driver to calculate the turbulence fields.
error FatalError
scalar y
scalarField & omega(bool init=false)
Return non-const access to the master&#39;s omega field.
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
const objectRegistry & db() const
Return local objectRegistry.
Definition: fvPatchField.C:188
bool manipulatedMatrix() const
Return true if the matrix has already been manipulated.
Definition: fvPatchField.H:354
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
virtual void write(Ostream &) const
Write.
label size() const
Return the number of elements in the UList.
List< List< scalar > > cornerWeights_
List of averaging corner weights.
virtual void operator==(const fvPatchField< Type > &)
Definition: fvPatchField.C:559
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
static const word propertiesName
Default name of the turbulence properties dictionary.
omegaWallFunctionFvPatchScalarField(const fvPatch &, const DimensionedField< scalar, volMesh > &)
Construct from patch and internal field.
virtual tmp< volScalarField > k() const =0
Return the turbulence kinetic energy.
const DimensionedField< Type, volMesh > & dimensionedInternalField() const
Return dimensioned internal field reference.
Definition: fvPatchField.H:307
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
const volVectorField & U() const
Access function to velocity field.
makePatchTypeField(fvPatchVectorField, SRFFreestreamVelocityFvPatchVectorField)
virtual tmp< volScalarField > nut() const =0
Return the turbulence viscosity.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
virtual label size() const
Return size.
Definition: fvPatch.H:161
const fvPatch & patch() const
Return patch.
Definition: fvPatchField.H:300
virtual label & master()
Return non-const access to the master patch ID.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
void setValues(const labelUList &cells, const UList< Type > &values)
Set solution in given cells to the specified values.
Definition: fvMatrix.C:488
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53