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