nutURoughWallFunctionFvPatchScalarField.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-2024 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 "momentumTransportModel.H"
28 #include "fieldMapper.H"
29 #include "volFields.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
38 
40 {
41  const label patchi = patch().index();
42 
43  const momentumTransportModel& turbModel =
44  db().lookupType<momentumTransportModel>(internalField().group());
45 
46  const scalarField& y = turbModel.yb()[patchi];
47  const tmp<scalarField> tnuw = turbModel.nu(patchi);
48  const scalarField& nuw = tnuw();
49 
50  const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
51  const scalarField magUp(mag(Uw.patchInternalField() - Uw));
52 
53  tmp<scalarField> tyPlus = yPlus(magUp);
54  scalarField& yPlus = tyPlus.ref();
55 
56  tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
57  scalarField& nutw = tnutw.ref();
58 
59  forAll(yPlus, facei)
60  {
61  const scalar Re = magUp[facei]*y[facei]/nuw[facei];
62 
63  if (sqr(yPlus[facei]) > Re)
64  {
65  nutw[facei] = nuw[facei]*(sqr(yPlus[facei])/Re - 1);
66  }
67  else
68  {
69  nutw[facei] = 0;
70  }
71  }
72 
73  return tnutw;
74 }
75 
76 
78 (
79  const scalarField& magUp
80 ) const
81 {
82  const label patchi = patch().index();
83 
84  const momentumTransportModel& turbModel =
85  db().lookupType<momentumTransportModel>(internalField().group());
86 
87  const scalarField& y = turbModel.yb()[patchi];
88  const tmp<scalarField> tnuw = turbModel.nu(patchi);
89  const scalarField& nuw = tnuw();
90 
91  tmp<scalarField> tyPlus(new scalarField(patch().size(), 0.0));
92  scalarField& yPlus = tyPlus.ref();
93 
94  static const scalar c2 = 2.25/(90 - 2.25);
95  static const scalar c3 = 2.0*atan(1.0)/log(90/2.25);
96  static const scalar c4 = c3*log(2.25);
97 
98  // If KsPlus is based on YPlus the extra term added to the law
99  // of the wall will depend on yPlus
100  forAll(yPlus, facei)
101  {
102  if (Ks_[facei] > 0)
103  {
104  // Rough Walls
105 
106  const scalar c1 = 1/(90 - 2.25) + Cs_[facei];
107 
108  const scalar magUpara = magUp[facei];
109  const scalar Re = magUpara*y[facei]/nuw[facei];
110  const scalar kappaRe = kappa_*Re;
111 
112  scalar yp = yPlusLam_;
113  const scalar ryPlusLam = 1/yp;
114 
115  int iter = 0;
116  scalar yPlusLast = 0.0;
117 
118  const scalar dKsPlusdYPlus = Ks_[facei]/y[facei];
119 
120  do
121  {
122  yPlusLast = yp;
123 
124  // The non-dimensional roughness height
125  scalar KsPlus = yp*dKsPlusdYPlus;
126 
127  // The extra term in the law-of-the-wall
128  scalar yPlusGPrime = 0;
129  scalar E = E_;
130 
131  if (KsPlus >= 90)
132  {
133  const scalar t1 = 1 + Cs_[facei]*KsPlus;
134  E = E_/t1;
135  yPlusGPrime = Cs_[facei]*KsPlus/t1;
136  }
137  else if (KsPlus > 2.25)
138  {
139  const scalar t1 = c1*KsPlus - c2;
140  const scalar t2 = c3*log(KsPlus) - c4;
141  const scalar sint2 = sin(t2);
142  const scalar logt1 = log(t1);
143  E = E_/pow(t1, sint2);
144  yPlusGPrime =
145  (c1*sint2*KsPlus/t1) + (c3*logt1*cos(t2));
146  }
147 
148  const scalar yPlusMin = constant::mathematical::e/E;
149 
150  if (kappa_*yPlusMin > 1)
151  {
152  yp = max
153  (
154  (kappaRe + yp*(1 - yPlusGPrime))
155  /(1 + log(E*yp) - yPlusGPrime),
156  sqrt(Re)
157  );
158  }
159  else
160  {
161  if (log(E*yp) < kappa_*yp)
162  {
163  yp = max
164  (
165  (kappaRe + yp*(1 - yPlusGPrime))
166  /(1 + log(E*yp) - yPlusGPrime),
167  yPlusMin
168  );
169  }
170  else
171  {
172  yp = sqrt(Re);
173  }
174  }
175  } while(mag(ryPlusLam*(yp - yPlusLast)) > 0.0001 && ++iter < 20);
176 
177  yPlus[facei] = yp;
178  }
179  else
180  {
181  // Smooth Walls
182  const scalar Re = magUp[facei]*y[facei]/nuw[facei];
183  const scalar ryPlusLam = 1/yPlusLam_;
184 
185  int iter = 0;
186  scalar yp = yPlusLam_;
187  scalar yPlusLast = yp;
188 
189  do
190  {
191  yPlusLast = yp;
192  if (yp > yPlusLam_)
193  {
194  yp = (kappa_*Re + yp)/(1 + log(E_*yp));
195  }
196  else
197  {
198  yp = sqrt(Re);
199  }
200  } while(mag(ryPlusLam*(yp - yPlusLast)) > 0.0001 && ++iter < 20);
201 
202  yPlus[facei] = yp;
203  }
204  }
205 
206  return tyPlus;
207 }
208 
209 
210 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
211 
213 (
214  const fvPatch& p,
216  const dictionary& dict
217 )
218 :
220  Ks_("Ks", dimLength, dict, p.size()),
221  Cs_("Cs", dimless, dict, p.size())
222 {}
223 
224 
226 (
228  const fvPatch& p,
230  const fieldMapper& mapper
231 )
232 :
233  nutUWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
234  Ks_(mapper(ptf.Ks_)),
235  Cs_(mapper(ptf.Cs_))
236 {}
237 
238 
240 (
243 )
244 :
246  Ks_(rwfpsf.Ks_),
247  Cs_(rwfpsf.Cs_)
248 {}
249 
250 
251 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
252 
254 (
255  const fvPatchScalarField& ptf,
256  const fieldMapper& mapper
257 )
258 {
260 
262  refCast<const nutURoughWallFunctionFvPatchScalarField>(ptf);
263 
264  mapper(Ks_, nrwfpsf.Ks_);
265  mapper(Cs_, nrwfpsf.Cs_);
266 }
267 
268 
270 (
271  const fvPatchScalarField& ptf
272 )
273 {
274  nutUWallFunctionFvPatchScalarField::reset(ptf);
275 
277  refCast<const nutURoughWallFunctionFvPatchScalarField>(ptf);
278 
279  Ks_.reset(nrwfpsf.Ks_);
280  Cs_.reset(nrwfpsf.Cs_);
281 }
282 
283 
285 {
287  writeLocalEntries(os);
288  writeEntry(os, "Cs", Cs_);
289  writeEntry(os, "Ks", Ks_);
290  writeEntry(os, "value", *this);
291 }
292 
293 
294 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
295 
297 (
300 );
301 
302 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
303 
304 } // End namespace Foam
305 
306 // ************************************************************************* //
scalar y
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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...
void reset(const Field< Type > &)
Reset the field values to the given field.
Definition: Field.C:456
const Boundary & boundaryField() const
Return const-reference to the boundary field.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:57
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
Abstract base class for field mapping.
Definition: fieldMapper.H:48
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:88
virtual void write(Ostream &) const
Write.
Definition: fvPatchField.C:229
virtual tmp< Field< Type > > patchInternalField() const
Return internal field next to patch as patch field.
Definition: fvPatchField.C:170
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:64
Abstract base class for turbulence models (RAS, LES and laminar).
const volScalarField::Boundary & yb() const
Return the near wall distance.
const volVectorField & U() const
Access function to velocity field.
virtual tmp< volScalarField > nu() const
Return the laminar viscosity.
This boundary condition provides a turbulent kinematic viscosity condition when using wall functions ...
virtual tmp< scalarField > nut() const
Calculate the turbulence viscosity.
nutURoughWallFunctionFvPatchScalarField(const fvPatch &, const DimensionedField< scalar, volMesh > &, const dictionary &)
Construct from patch, internal field and dictionary.
virtual void reset(const fvPatchScalarField &)
Reset the fvPatchField to the given fvPatchField.
virtual void map(const fvPatchScalarField &, const fieldMapper &)
Map the given fvPatchField onto this fvPatchField.
This boundary condition provides a turbulent kinematic viscosity condition when using wall functions,...
virtual tmp< scalarField > yPlus() const
Calculate and return the yPlus at the boundary.
scalar yPlusLam_
Y+ at the edge of the laminar sublayer.
virtual void writeLocalEntries(Ostream &) const
Write local wall function variables.
virtual void map(const fvPatchScalarField &, const fieldMapper &)
Map the given fvPatchField onto this fvPatchField.
static const nutWallFunctionFvPatchScalarField & nutw(const momentumTransportModel &turbModel, const label patchi)
Return the nut patchField for the given wall patch.
A class for managing temporary objects.
Definition: tmp.H:55
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
const scalar magUp
label patchi
const dimensionedScalar c1
First radiation constant: default SI units: [W/m^2].
const dimensionedScalar c2
Second radiation constant: default SI units: [m K].
Namespace for OpenFOAM.
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
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimless
dimensionedScalar sin(const dimensionedScalar &ds)
const dimensionSet dimLength
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar log(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
void writeEntry(Ostream &os, const HashTable< T, Key, Hash > &ht)
Definition: HashTableIO.C:96
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
makePatchTypeField(fvPatchScalarField, atmBoundaryLayerInletEpsilonFvPatchScalarField)
dimensionedScalar atan(const dimensionedScalar &ds)
dimensionedScalar cos(const dimensionedScalar &ds)
scalarField Re(const UList< complex > &cf)
Definition: complexFields.C:97
dictionary dict
volScalarField & p