FreeFOAM The Cross-Platform CFD Toolkit
nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.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) 1991-2010 OpenCFD Ltd.
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 
29 #include <finiteVolume/volFields.H>
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 namespace incompressible
37 {
38 namespace RASModels
39 {
40 
41 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
42 
43 tmp<scalarField>
44 nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcNut() const
45 {
46  const label patchI = patch().index();
47 
48  const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
49  const scalarField& y = rasModel.y()[patchI];
50  const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
51  const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
52 
53  // The flow velocity at the adjacent cell centre
54  const scalarField magUp = mag(Uw.patchInternalField() - Uw);
55 
56  tmp<scalarField> tyPlus = calcYPlus(magUp);
57  scalarField& yPlus = tyPlus();
58 
59  tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
60  scalarField& nutw = tnutw();
61 
62  forAll(yPlus, facei)
63  {
64  if (yPlus[facei] > yPlusLam_)
65  {
66  const scalar Re = magUp[facei]*y[facei]/nuw[facei];
67  nutw[facei] = nuw[facei]*(sqr(yPlus[facei])/Re - 1);
68  }
69  }
70 
71  return tnutw;
72 }
73 
74 
75 tmp<scalarField>
76 nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
77 (
78  const scalarField& magUp
79 ) const
80 {
81  const label patchI = patch().index();
82 
83  const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
84  const scalarField& y = rasModel.y()[patchI];
85  const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
86 
87  tmp<scalarField> tyPlus(new scalarField(patch().size(), 0.0));
88  scalarField& yPlus = tyPlus();
89 
90  if (roughnessHeight_ > 0.0)
91  {
92  // Rough Walls
93  const scalar c_1 = 1/(90 - 2.25) + roughnessConstant_;
94  static const scalar c_2 = 2.25/(90 - 2.25);
95  static const scalar c_3 = 2.0*atan(1.0)/log(90/2.25);
96  static const scalar c_4 = c_3*log(2.25);
97 
98  //if (KsPlusBasedOnYPlus_)
99  {
100  // If KsPlus is based on YPlus the extra term added to the law
101  // of the wall will depend on yPlus
102  forAll(yPlus, facei)
103  {
104  const scalar magUpara = magUp[facei];
105  const scalar Re = magUpara*y[facei]/nuw[facei];
106  const scalar kappaRe = kappa_*Re;
107 
108  scalar yp = yPlusLam_;
109  const scalar ryPlusLam = 1.0/yp;
110 
111  int iter = 0;
112  scalar yPlusLast = 0.0;
113  scalar dKsPlusdYPlus = roughnessHeight_/y[facei];
114 
115  // Enforce the roughnessHeight to be less than the distance to
116  // the first cell centre
117  if (dKsPlusdYPlus > 1)
118  {
119  dKsPlusdYPlus = 1;
120  }
121 
122  // Additional tuning parameter (fudge factor) - nominally = 1
123  dKsPlusdYPlus *= roughnessFudgeFactor_;
124 
125  do
126  {
127  yPlusLast = yp;
128 
129  // The non-dimensional roughness height
130  scalar KsPlus = yp*dKsPlusdYPlus;
131 
132  // The extra term in the law-of-the-wall
133  scalar G = 0.0;
134 
135  scalar yPlusGPrime = 0.0;
136 
137  if (KsPlus >= 90)
138  {
139  const scalar t_1 = 1 + roughnessConstant_*KsPlus;
140  G = log(t_1);
141  yPlusGPrime = roughnessConstant_*KsPlus/t_1;
142  }
143  else if (KsPlus > 2.25)
144  {
145  const scalar t_1 = c_1*KsPlus - c_2;
146  const scalar t_2 = c_3*log(KsPlus) - c_4;
147  const scalar sint_2 = sin(t_2);
148  const scalar logt_1 = log(t_1);
149  G = logt_1*sint_2;
150  yPlusGPrime =
151  (c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2));
152  }
153 
154  scalar denom = 1.0 + log(E_*yp) - G - yPlusGPrime;
155  if (mag(denom) > VSMALL)
156  {
157  yp = (kappaRe + yp*(1 - yPlusGPrime))/denom;
158  }
159  } while
160  (
161  mag(ryPlusLam*(yp - yPlusLast)) > 0.0001
162  && ++iter < 10
163  && yp > VSMALL
164  );
165 
166  yPlus[facei] = max(0.0, yp);
167  }
168  }
169  }
170  else
171  {
172  // Smooth Walls
173  forAll(yPlus, facei)
174  {
175  const scalar magUpara = magUp[facei];
176  const scalar Re = magUpara*y[facei]/nuw[facei];
177  const scalar kappaRe = kappa_*Re;
178 
179  scalar yp = yPlusLam_;
180  const scalar ryPlusLam = 1.0/yp;
181 
182  int iter = 0;
183  scalar yPlusLast = 0.0;
184 
185  do
186  {
187  yPlusLast = yp;
188  yp = (kappaRe + yp)/(1.0 + log(E_*yp));
189 
190  } while(mag(ryPlusLam*(yp - yPlusLast)) > 0.0001 && ++iter < 10);
191 
192  yPlus[facei] = max(0.0, yp);
193  }
194  }
195 
196  return tyPlus;
197 }
198 
199 
200 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
201 
204 (
205  const fvPatch& p,
207 )
208 :
210  roughnessHeight_(pTraits<scalar>::zero),
211  roughnessConstant_(pTraits<scalar>::zero),
212  roughnessFudgeFactor_(pTraits<scalar>::zero)
213 {}
214 
215 
218 (
220  const fvPatch& p,
222  const fvPatchFieldMapper& mapper
223 )
224 :
225  nutWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
226  roughnessHeight_(ptf.roughnessHeight_),
227  roughnessConstant_(ptf.roughnessConstant_),
228  roughnessFudgeFactor_(ptf.roughnessFudgeFactor_)
229 {}
230 
231 
234 (
235  const fvPatch& p,
237  const dictionary& dict
238 )
239 :
241  roughnessHeight_(readScalar(dict.lookup("roughnessHeight"))),
242  roughnessConstant_(readScalar(dict.lookup("roughnessConstant"))),
243  roughnessFudgeFactor_(readScalar(dict.lookup("roughnessFudgeFactor")))
244 {}
245 
246 
249 (
251 )
252 :
254  roughnessHeight_(rwfpsf.roughnessHeight_),
255  roughnessConstant_(rwfpsf.roughnessConstant_),
256  roughnessFudgeFactor_(rwfpsf.roughnessFudgeFactor_)
257 {}
258 
259 
262 (
265 )
266 :
268  roughnessHeight_(rwfpsf.roughnessHeight_),
269  roughnessConstant_(rwfpsf.roughnessConstant_),
270  roughnessFudgeFactor_(rwfpsf.roughnessFudgeFactor_)
271 {}
272 
273 
274 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
275 
278 {
279  const label patchI = patch().index();
280 
281  const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
282  const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
283  const scalarField magUp = mag(Uw.patchInternalField() - Uw);
284 
285  return calcYPlus(magUp);
286 }
287 
288 
290 (
291  Ostream& os
292 ) const
293 {
295  writeLocalEntries(os);
296  os.writeKeyword("roughnessHeight")
297  << roughnessHeight_ << token::END_STATEMENT << nl;
298  os.writeKeyword("roughnessConstant")
299  << roughnessConstant_ << token::END_STATEMENT << nl;
300  os.writeKeyword("roughnessFudgeFactor")
301  << roughnessFudgeFactor_ << token::END_STATEMENT << nl;
302  writeEntry("value", os);
303 }
304 
305 
306 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
307 
309 (
312 );
313 
314 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
315 
316 } // End namespace RASModels
317 } // End namespace incompressible
318 } // End namespace Foam
319 
320 // ************************ vim: set sw=4 sts=4 et: ************************ //