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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016, 2019 OpenFOAM Foundation
9 Copyright (C) 2017-2022 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
31#include "turbulenceModel.H"
32#include "fvMatrix.H"
34
35// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36
38
39// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
40
42{
43 if (master_ != -1)
44 {
45 return;
46 }
47
48 const auto& omega =
49 static_cast<const volScalarField&>(this->internalField());
50
51 const volScalarField::Boundary& bf = omega.boundaryField();
52
53 label master = -1;
54 forAll(bf, patchi)
55 {
56 if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchi]))
57 {
59
60 if (master == -1)
61 {
62 master = patchi;
63 }
64
65 opf.master() = master;
66 }
67 }
68}
69
70
72{
73 const auto& omega =
74 static_cast<const volScalarField&>(this->internalField());
75
76 const volScalarField::Boundary& bf = omega.boundaryField();
77
78 const fvMesh& mesh = omega.mesh();
79
80 if (initialised_ && !mesh.changing())
81 {
82 return;
83 }
84
85 volScalarField weights
86 (
88 (
89 "weights",
90 mesh.time().timeName(),
91 mesh,
94 false // do not register
95 ),
96 mesh,
98 );
99
100 DynamicList<label> omegaPatches(bf.size());
101 forAll(bf, patchi)
102 {
103 if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchi]))
104 {
105 omegaPatches.append(patchi);
106
107 const labelUList& faceCells = bf[patchi].patch().faceCells();
108 for (const auto& celli : faceCells)
109 {
110 ++weights[celli];
111 }
112 }
113 }
114
115 cornerWeights_.setSize(bf.size());
116 for (const auto& patchi : omegaPatches)
117 {
118 const fvPatchScalarField& wf = weights.boundaryField()[patchi];
119 cornerWeights_[patchi] = 1.0/wf.patchInternalField();
120 }
121
122 G_.setSize(internalField().size(), 0.0);
123 omega_.setSize(internalField().size(), 0.0);
124
125 initialised_ = true;
126}
127
128
131(
132 const label patchi
133)
134{
135 const auto& omega =
136 static_cast<const volScalarField&>(this->internalField());
137
138 const volScalarField::Boundary& bf = omega.boundaryField();
139
140 const auto& opf =
141 refCast<const omegaWallFunctionFvPatchScalarField>(bf[patchi]);
142
143 return const_cast<omegaWallFunctionFvPatchScalarField&>(opf);
144}
145
146
148(
149 const turbulenceModel& turbModel,
150 scalarField& G0,
151 scalarField& omega0
152)
153{
154 // accumulate all of the G and omega contributions
155 forAll(cornerWeights_, patchi)
156 {
157 if (!cornerWeights_[patchi].empty())
158 {
159 omegaWallFunctionFvPatchScalarField& opf = omegaPatch(patchi);
160
161 const List<scalar>& w = cornerWeights_[patchi];
162
163 opf.calculate(turbModel, w, opf.patch(), G0, omega0);
164 }
165 }
166
167 // apply zero-gradient condition for omega
168 forAll(cornerWeights_, patchi)
169 {
170 if (!cornerWeights_[patchi].empty())
171 {
172 omegaWallFunctionFvPatchScalarField& opf = omegaPatch(patchi);
173
174 opf == scalarField(omega0, opf.patch().faceCells());
175 }
176 }
177}
178
179
181(
182 const turbulenceModel& turbModel,
183 const List<scalar>& cornerWeights,
184 const fvPatch& patch,
185 scalarField& G0,
186 scalarField& omega0
187)
188{
189 const label patchi = patch.index();
190
191 const tmp<scalarField> tnutw = turbModel.nut(patchi);
192 const scalarField& nutw = tnutw();
193
194 const scalarField& y = turbModel.y()[patchi];
195
196 const tmp<scalarField> tnuw = turbModel.nu(patchi);
197 const scalarField& nuw = tnuw();
198
199 const tmp<volScalarField> tk = turbModel.k();
200 const volScalarField& k = tk();
201
202 const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
203
204 const scalarField magGradUw(mag(Uw.snGrad()));
205
206 const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
207 const scalar kappa = wallCoeffs_.kappa();
208 const scalar yPlusLam = wallCoeffs_.yPlusLam();
209
210 // Set omega and G
211 forAll(nutw, facei)
212 {
213 const label celli = patch.faceCells()[facei];
214 const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei];
215 const scalar w = cornerWeights[facei];
216
217 // Contribution from the viscous sublayer
218 const scalar omegaVis = 6.0*nuw[facei]/(beta1_*sqr(y[facei]));
219
220 // Contribution from the inertial sublayer
221 const scalar omegaLog = sqrt(k[celli])/(Cmu25*kappa*y[facei]);
222
223 switch (blender_)
224 {
225 case blenderType::STEPWISE:
226 {
227 if (yPlus > yPlusLam)
228 {
229 omega0[celli] += w*omegaLog;
230 }
231 else
232 {
233 omega0[celli] += w*omegaVis;
234 }
235 break;
236 }
237
238 case blenderType::BINOMIAL:
239 {
240 omega0[celli] +=
241 w*pow
242 (
243 pow(omegaVis, n_) + pow(omegaLog, n_),
244 scalar(1)/n_
245 );
246 break;
247 }
248
249 case blenderType::MAX:
250 {
251 // (PH:Eq. 27)
252 omega0[celli] += max(omegaVis, omegaLog);
253 break;
254 }
255
256 case blenderType::EXPONENTIAL:
257 {
258 // (PH:Eq. 31)
259 const scalar Gamma = 0.01*pow4(yPlus)/(1 + 5*yPlus);
260 const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL);
261
262 omega0[celli] +=
263 w*(omegaVis*exp(-Gamma) + omegaLog*exp(-invGamma));
264 break;
265 }
266
267 case blenderType::TANH:
268 {
269 // (KAS:Eqs. 33-34)
270 const scalar phiTanh = tanh(pow4(0.1*yPlus));
271 const scalar b1 = omegaVis + omegaLog;
272 const scalar b2 =
273 pow(pow(omegaVis, 1.2) + pow(omegaLog, 1.2), 1.0/1.2);
274
275 omega0[celli] += phiTanh*b1 + (1 - phiTanh)*b2;
276 break;
277 }
278 }
279
280 if (!(blender_ == blenderType::STEPWISE) || yPlus > yPlusLam)
281 {
282 G0[celli] +=
283 w
284 *(nutw[facei] + nuw[facei])
285 *magGradUw[facei]
286 *Cmu25*sqrt(k[celli])
287 /(kappa*y[facei]);
288 }
289 }
290}
291
292
294(
295 Ostream& os
296) const
297{
299 os.writeEntryIfDifferent<scalar>("beta1", 0.075, beta1_);
300 wallCoeffs_.writeEntries(os);
301}
302
303
304// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
305
307(
308 const fvPatch& p,
310)
311:
312 fixedValueFvPatchField<scalar>(p, iF),
314 initialised_(false),
315 master_(-1),
316 beta1_(0.075),
317 wallCoeffs_(),
318 G_(),
319 omega_(),
320 cornerWeights_()
321{}
322
323
325(
327 const fvPatch& p,
329 const fvPatchFieldMapper& mapper
330)
331:
332 fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
334 initialised_(false),
335 master_(-1),
336 beta1_(ptf.beta1_),
337 wallCoeffs_(ptf.wallCoeffs_),
338 G_(),
339 omega_(),
340 cornerWeights_()
341{}
342
343
345(
346 const fvPatch& p,
348 const dictionary& dict
349)
350:
351 fixedValueFvPatchField<scalar>(p, iF, dict),
352 wallFunctionBlenders(dict, blenderType::BINOMIAL, scalar(2)),
353 initialised_(false),
354 master_(-1),
355 beta1_(dict.getOrDefault<scalar>("beta1", 0.075)),
356 wallCoeffs_(dict),
357 G_(),
358 omega_(),
359 cornerWeights_()
360{
361 // apply zero-gradient condition on start-up
363}
364
365
367(
369)
370:
371 fixedValueFvPatchField<scalar>(owfpsf),
372 wallFunctionBlenders(owfpsf),
373 initialised_(false),
374 master_(-1),
375 beta1_(owfpsf.beta1_),
376 wallCoeffs_(owfpsf.wallCoeffs_),
377 G_(),
378 omega_(),
379 cornerWeights_()
380{}
381
382
384(
387)
388:
389 fixedValueFvPatchField<scalar>(owfpsf, iF),
390 wallFunctionBlenders(owfpsf),
391 initialised_(false),
392 master_(-1),
393 beta1_(owfpsf.beta1_),
394 wallCoeffs_(owfpsf.wallCoeffs_),
395 G_(),
396 omega_(),
397 cornerWeights_()
398{}
399
400
401// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
402
404(
405 bool init
406)
407{
408 if (patch().index() == master_)
409 {
410 if (init)
411 {
412 G_ = 0.0;
413 }
414
415 return G_;
416 }
417
418 return omegaPatch(master_).G();
419}
420
421
423(
424 bool init
425)
426{
427 if (patch().index() == master_)
428 {
429 if (init)
430 {
431 omega_ = 0.0;
432 }
433
434 return omega_;
435 }
436
437 return omegaPatch(master_).omega(init);
438}
439
440
442{
443 if (updated())
444 {
445 return;
446 }
447
448 const auto& turbModel = db().lookupObject<turbulenceModel>
449 (
451 (
453 internalField().group()
454 )
455 );
456
457 setMaster();
458
459 if (patch().index() == master_)
460 {
461 createAveragingWeights();
462 calculateTurbulenceFields(turbModel, G(true), omega(true));
463 }
464
465 const scalarField& G0 = this->G();
466 const scalarField& omega0 = this->omega();
467
468 typedef DimensionedField<scalar, volMesh> FieldType;
469
470 FieldType& G = db().lookupObjectRef<FieldType>(turbModel.GName());
471
472 FieldType& omega = const_cast<FieldType&>(internalField());
473
474 forAll(*this, facei)
475 {
476 const label celli = patch().faceCells()[facei];
477
478 G[celli] = G0[celli];
479 omega[celli] = omega0[celli];
480 }
481
483}
484
485
487(
488 const scalarField& weights
489)
490{
491 if (updated())
492 {
493 return;
494 }
495
496 const auto& turbModel = db().lookupObject<turbulenceModel>
497 (
499 (
501 internalField().group()
502 )
503 );
504
505 setMaster();
506
507 if (patch().index() == master_)
508 {
509 createAveragingWeights();
510 calculateTurbulenceFields(turbModel, G(true), omega(true));
511 }
512
513 const scalarField& G0 = this->G();
514 const scalarField& omega0 = this->omega();
515
516 typedef DimensionedField<scalar, volMesh> FieldType;
517
518 FieldType& G = db().lookupObjectRef<FieldType>(turbModel.GName());
519
520 FieldType& omega = const_cast<FieldType&>(internalField());
521
522 scalarField& omegaf = *this;
523
524 // only set the values if the weights are > tolerance
525 forAll(weights, facei)
526 {
527 const scalar w = weights[facei];
528
529 if (w > tolerance_)
530 {
531 const label celli = patch().faceCells()[facei];
532
533 G[celli] = (1.0 - w)*G[celli] + w*G0[celli];
534 omega[celli] = (1.0 - w)*omega[celli] + w*omega0[celli];
535 omegaf[facei] = omega[celli];
536 }
537 }
538
540}
541
542
544(
545 fvMatrix<scalar>& matrix
546)
547{
548 if (manipulatedMatrix())
549 {
550 return;
551 }
552
553 matrix.setValues(patch().faceCells(), patchInternalField());
554
556}
557
558
560(
561 fvMatrix<scalar>& matrix,
562 const Field<scalar>& weights
563)
564{
565 if (manipulatedMatrix())
566 {
567 return;
568 }
569
570 DynamicList<label> constraintCells(weights.size());
571 DynamicList<scalar> constraintValues(weights.size());
572 const labelUList& faceCells = patch().faceCells();
573
574 const DimensionedField<scalar, volMesh>& fld = internalField();
575
576 forAll(weights, facei)
577 {
578 // only set the values if the weights are > tolerance
579 if (tolerance_ < weights[facei])
580 {
581 const label celli = faceCells[facei];
582
583 constraintCells.append(celli);
584 constraintValues.append(fld[celli]);
585 }
586 }
587
588 if (debug)
589 {
590 Pout<< "Patch: " << patch().name()
591 << ": number of constrained cells = " << constraintCells.size()
592 << " out of " << patch().size()
593 << endl;
594 }
595
596 matrix.setValues(constraintCells, constraintValues);
597
599}
600
601
603(
604 Ostream& os
605) const
606{
608 writeLocalEntries(os);
609 writeEntry("value", os);
610}
611
612
613// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
614
615namespace Foam
616{
618 (
621 );
622}
623
624
625// ************************************************************************* //
scalar y
label k
Macros for easy insertion into run-time selection tables.
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
static word groupName(StringType base, const word &group)
Create dot-delimited name.group string.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
void setSize(const label n)
Alias for resize()
Definition: List.H:218
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
Ostream & writeEntryIfDifferent(const word &key, const T &value1, const T &value2)
Write a keyword/value entry only when the two values differ.
Definition: Ostream.H:251
static word timeName(const scalar t, const int precision=precision_)
Definition: Time.C:780
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Smooth ATC in cells next to a set of patches supplied by type.
Definition: faceCells.H:59
This boundary condition supplies a fixed value constraint, and is the base class for a number of othe...
virtual bool write()
Write the output fields.
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:121
void setValues(const labelUList &cellLabels, const Type &value)
Definition: fvMatrix.C:969
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
A FieldMapper for finite-volume patch fields.
virtual tmp< Field< Type > > patchInternalField() const
Return internal field next to patch as patch field.
Definition: fvPatchField.C:237
virtual tmp< Field< Type > > snGrad() const
Return patch-normal gradient.
Definition: fvPatchField.C:229
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
Definition: fvPatchField.C:325
virtual void manipulateMatrix(fvMatrix< Type > &matrix)
Manipulate matrix.
Definition: fvPatchField.C:358
const DimensionedField< Type, volMesh > & internalField() const
Return dimensioned internal field reference.
Definition: fvPatchField.H:368
const fvPatch & patch() const
Return patch.
Definition: fvPatchField.H:362
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:71
virtual const labelUList & faceCells() const
Return faceCells.
Definition: fvPatch.C:113
This boundary condition provides a wall function for the specific dissipation rate (i....
virtual void calculate(const turbulenceModel &turbulence, const List< scalar > &cornerWeights, const fvPatch &patch, scalarField &G, scalarField &omega)
Calculate the omega and G.
void writeLocalEntries(Ostream &) const
Write local wall function variables.
virtual void manipulateMatrix(fvMatrix< scalar > &matrix)
Manipulate matrix.
static scalar tolerance_
Tolerance used in weighted calculations.
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
virtual label & master()
Return non-const access to the master patch ID.
virtual void updateWeightedCoeffs(const scalarField &weights)
Update the coefficients associated with the patch field.
scalarField & omega(bool init=false)
Return non-const access to the master's omega field.
virtual void calculateTurbulenceFields(const turbulenceModel &turbulence, scalarField &G0, scalarField &omega0)
Main driver to calculate the turbulence fields.
virtual omegaWallFunctionFvPatchScalarField & omegaPatch(const label patchi)
Helper function to return non-const access to an omega patch.
bool changing() const noexcept
Is mesh changing (topology changing and/or moving)
Definition: polyMesh.H:562
const volScalarField & omega() const
Access functions.
A class for managing temporary objects.
Definition: tmp.H:65
Abstract base class for turbulence models (RAS, LES and laminar).
const volVectorField & U() const
Access function to velocity field.
static const word propertiesName
Default name of the turbulence properties dictionary.
virtual tmp< volScalarField > k() const =0
Return the turbulence kinetic energy.
virtual tmp< volScalarField > nu() const =0
Return the laminar viscosity.
virtual tmp< volScalarField > nut() const =0
Return the turbulence viscosity.
const nearWallDist & y() const
Return the near wall distances.
The class wallFunctionBlenders is a base class that hosts common entries for various derived wall-fun...
void writeEntries(Ostream &) const
Write wall-function blending data as dictionary entries.
blenderType
Options for the blending treatment of viscous and inertial sublayers.
volScalarField & p
dynamicFvMesh & mesh
scalar yPlus
OBJstream os(runTime.globalPath()/outputName)
#define makePatchTypeField(PatchTypeField, typePatchTypeField)
Definition: fvPatchField.H:676
const dimensionedScalar G
Newtonian constant of gravitation.
Namespace for OpenFOAM.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedScalar exp(const dimensionedScalar &ds)
const dimensionSet dimless
Dimensionless.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar tanh(const dimensionedScalar &ds)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
tmp< faMatrix< Type > > operator==(const faMatrix< Type > &, const faMatrix< Type > &)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
dimensionedScalar pow4(const dimensionedScalar &ds)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
dimensionedScalar pow025(const dimensionedScalar &ds)
dictionary dict
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333