nutURoughWallFunctionFvPatchScalarField.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 OpenFOAM Foundation
9 Copyright (C) 2019-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
30#include "turbulenceModel.H"
31#include "fvPatchFieldMapper.H"
32#include "volFields.H"
34
35
36// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
37
39Foam::nutURoughWallFunctionFvPatchScalarField::calcNut() const
40{
41 const label patchi = patch().index();
42
43 const auto& turbModel = db().lookupObject<turbulenceModel>
44 (
46 (
48 internalField().group()
49 )
50 );
51
52 const scalarField& y = turbModel.y()[patchi];
53
54 const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
55
56 const tmp<scalarField> tnuw = turbModel.nu(patchi);
57 const scalarField& nuw = tnuw();
58
59 // The flow velocity at the adjacent cell centre
60 const scalarField magUp(mag(Uw.patchInternalField() - Uw));
61
62 const scalar yPlusLam = wallCoeffs_.yPlusLam();
63
64 tmp<scalarField> tyPlus = calcYPlus(magUp);
65 scalarField& yPlus = tyPlus.ref();
66
67 auto tnutw = tmp<scalarField>::New(patch().size(), Zero);
68 auto& nutw = tnutw.ref();
69
70 forAll(yPlus, facei)
71 {
72 if (yPlusLam < yPlus[facei])
73 {
74 const scalar Re = magUp[facei]*y[facei]/nuw[facei] + ROOTVSMALL;
75 nutw[facei] = nuw[facei]*(sqr(yPlus[facei])/Re - 1);
76 }
77 }
78
79 return tnutw;
80}
81
82
84Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
85(
86 const scalarField& magUp
87) const
88{
89 const label patchi = patch().index();
90
91 const auto& turbModel = db().lookupObject<turbulenceModel>
92 (
94 (
96 internalField().group()
97 )
98 );
99
100 const scalarField& y = turbModel.y()[patchi];
101
102 const tmp<scalarField> tnuw = turbModel.nu(patchi);
103 const scalarField& nuw = tnuw();
104
105 const scalar kappa = wallCoeffs_.kappa();
106 const scalar E = wallCoeffs_.E();
107 const scalar yPlusLam = wallCoeffs_.yPlusLam();
108
109 auto tyPlus = tmp<scalarField>::New(patch().size(), Zero);
110 auto& yPlus = tyPlus.ref();
111
112 if (roughnessHeight_ > 0.0)
113 {
114 // Rough Walls
115 const scalar c_1 = 1.0/(90.0 - 2.25) + roughnessConstant_;
116 static const scalar c_2 = 2.25/(90.0 - 2.25);
117 static const scalar c_3 = 2.0*atan(1.0)/log(90.0/2.25);
118 static const scalar c_4 = c_3*log(2.25);
119
120 {
121 // If KsPlus is based on YPlus the extra term added to the law
122 // of the wall will depend on yPlus
123 forAll(yPlus, facei)
124 {
125 const scalar magUpara = magUp[facei];
126 const scalar Re = magUpara*y[facei]/nuw[facei];
127 const scalar kappaRe = kappa*Re;
128
129 scalar yp = yPlusLam;
130 const scalar ryPlusLam = 1.0/yp;
131
132 int iter = 0;
133 scalar yPlusLast = 0.0;
134 scalar dKsPlusdYPlus = roughnessHeight_/y[facei];
135
136 // Additional tuning parameter - nominally = 1
137 dKsPlusdYPlus *= roughnessFactor_;
138
139 do
140 {
141 yPlusLast = yp;
142
143 // The non-dimensional roughness height
144 const scalar KsPlus = yp*dKsPlusdYPlus;
145
146 // The extra term in the law-of-the-wall
147 scalar G = 0;
148
149 scalar yPlusGPrime = 0;
150
151 if (KsPlus >= 90.0)
152 {
153 const scalar t_1 = 1 + roughnessConstant_*KsPlus;
154 G = log(t_1);
155 yPlusGPrime = roughnessConstant_*KsPlus/t_1;
156 }
157 else if (KsPlus > 2.25)
158 {
159 const scalar t_1 = c_1*KsPlus - c_2;
160 const scalar t_2 = c_3*log(KsPlus) - c_4;
161 const scalar sint_2 = sin(t_2);
162 const scalar logt_1 = log(t_1);
163 G = logt_1*sint_2;
164 yPlusGPrime =
165 (c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2));
166 }
167
168 const scalar denom = 1.0 + log(E*yp) - G - yPlusGPrime;
169 if (mag(denom) > VSMALL)
170 {
171 yp = (kappaRe + yp*(1 - yPlusGPrime))/denom;
172 }
173 } while
174 (
175 mag(ryPlusLam*(yp - yPlusLast)) > tolerance_
176 && ++iter < maxIter_
177 && yp > VSMALL
178 );
179
180 yPlus[facei] = max(scalar(0), yp);
181 }
182 }
183 }
184 else
185 {
186 // Smooth Walls
187 forAll(yPlus, facei)
188 {
189 const scalar magUpara = magUp[facei];
190 const scalar Re = magUpara*y[facei]/nuw[facei];
191 const scalar kappaRe = kappa*Re;
192
193 scalar yp = yPlusLam;
194 const scalar ryPlusLam = 1.0/yp;
195
196 int iter = 0;
197 scalar yPlusLast = 0;
198
199 do
200 {
201 yPlusLast = yp;
202 yp = (kappaRe + yp)/(1.0 + log(E*yp));
203
204 }
205 while
206 (
207 mag(ryPlusLam*(yp - yPlusLast)) > tolerance_
208 && ++iter < maxIter_
209 );
210
211 yPlus[facei] = max(scalar(0), yp);
212 }
213 }
214
215 return tyPlus;
216}
217
218
219void Foam::nutURoughWallFunctionFvPatchScalarField::writeLocalEntries
220(
221 Ostream& os
222) const
223{
224 os.writeEntry("roughnessHeight", roughnessHeight_);
225 os.writeEntry("roughnessConstant", roughnessConstant_);
226 os.writeEntry("roughnessFactor", roughnessFactor_);
227 os.writeEntryIfDifferent<label>("maxIter", 10, maxIter_);
228 os.writeEntryIfDifferent<scalar>("tolerance", 0.0001, tolerance_);
229}
230
231
232// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
233
234Foam::nutURoughWallFunctionFvPatchScalarField::
235nutURoughWallFunctionFvPatchScalarField
236(
237 const fvPatch& p,
239)
240:
242 roughnessHeight_(Zero),
243 roughnessConstant_(Zero),
244 roughnessFactor_(Zero),
245 maxIter_(10),
246 tolerance_(0.0001)
247{}
248
249
250Foam::nutURoughWallFunctionFvPatchScalarField::
251nutURoughWallFunctionFvPatchScalarField
252(
254 const fvPatch& p,
256 const fvPatchFieldMapper& mapper
257)
258:
259 nutWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
260 roughnessHeight_(ptf.roughnessHeight_),
261 roughnessConstant_(ptf.roughnessConstant_),
262 roughnessFactor_(ptf.roughnessFactor_),
263 maxIter_(ptf.maxIter_),
264 tolerance_(ptf.tolerance_)
265{}
266
267
268Foam::nutURoughWallFunctionFvPatchScalarField::
269nutURoughWallFunctionFvPatchScalarField
270(
271 const fvPatch& p,
273 const dictionary& dict
274)
275:
277 roughnessHeight_(dict.get<scalar>("roughnessHeight")),
278 roughnessConstant_(dict.get<scalar>("roughnessConstant")),
279 roughnessFactor_(dict.get<scalar>("roughnessFactor")),
280 maxIter_(dict.getOrDefault<label>("maxIter", 10)),
281 tolerance_(dict.getOrDefault<scalar>("tolerance", 0.0001))
282{}
283
284
285Foam::nutURoughWallFunctionFvPatchScalarField::
286nutURoughWallFunctionFvPatchScalarField
287(
289)
290:
292 roughnessHeight_(rwfpsf.roughnessHeight_),
293 roughnessConstant_(rwfpsf.roughnessConstant_),
294 roughnessFactor_(rwfpsf.roughnessFactor_),
295 maxIter_(rwfpsf.maxIter_),
296 tolerance_(rwfpsf.tolerance_)
297{}
298
299
300Foam::nutURoughWallFunctionFvPatchScalarField::
301nutURoughWallFunctionFvPatchScalarField
302(
305)
306:
308 roughnessHeight_(rwfpsf.roughnessHeight_),
309 roughnessConstant_(rwfpsf.roughnessConstant_),
310 roughnessFactor_(rwfpsf.roughnessFactor_),
311 maxIter_(rwfpsf.maxIter_),
312 tolerance_(rwfpsf.tolerance_)
313{}
314
315
316// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
317
320{
321 const label patchi = patch().index();
322
323 const auto& turbModel = db().lookupObject<turbulenceModel>
324 (
326 (
328 internalField().group()
329 )
330 );
331
332 const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
333 tmp<scalarField> magUp = mag(Uw.patchInternalField() - Uw);
334
335 return calcYPlus(magUp());
336}
337
338
340(
341 Ostream& os
342) const
343{
345 writeLocalEntries(os);
346 writeEntry("value", os);
347}
348
349
350// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
351
352namespace Foam
353{
355 (
358 );
359}
360
361// ************************************************************************* //
scalar y
Macros for easy insertion into run-time selection tables.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
static word groupName(StringType base, const word &group)
Create dot-delimited name.group string.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
Definition: Ostream.H:239
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 autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
virtual bool write()
Write the output fields.
A FieldMapper for finite-volume patch fields.
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:71
This boundary condition provides a wall function on the turbulent viscosity (i.e. nut) based on veloc...
virtual tmp< scalarField > yPlus() const
Calculate and return the yPlus at the boundary.
The class nutWallFunction is an abstract base class that hosts calculation methods and common functi...
static const nutWallFunctionFvPatchScalarField & nutw(const turbulenceModel &turbModel, const label patchi)
Return the nut patchField for the given wall patch.
wallFunctionCoefficients wallCoeffs_
Wall-function coefficients.
A class for managing temporary objects.
Definition: tmp.H:65
Abstract base class for turbulence models (RAS, LES and laminar).
static const word propertiesName
Default name of the turbulence properties dictionary.
scalar yPlusLam() const noexcept
Return the object: yPlusLam.
U
Definition: pEqn.H:72
volScalarField & p
scalar yPlus
scalar magUp
OBJstream os(runTime.globalPath()/outputName)
#define makePatchTypeField(PatchTypeField, typePatchTypeField)
Definition: fvPatchField.H:676
ThermalDiffusivity< CompressibleTurbulenceModel< fluidThermo > > turbulenceModel
constexpr const char *const group
Group name for atomic constants.
const dimensionedScalar kappa
Coulomb constant: default SI units: [N.m2/C2].
const dimensionedScalar G
Newtonian constant of gravitation.
const std::string patch
OpenFOAM patch number as a std::string.
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
scalarField Re(const UList< complex > &cf)
Extract real component.
Definition: complexField.C:159
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sin(const dimensionedScalar &ds)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar log(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
dimensionedScalar atan(const dimensionedScalar &ds)
fvPatchField< vector > fvPatchVectorField
dimensionedScalar cos(const dimensionedScalar &ds)
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333