FreeFOAM The Cross-Platform CFD Toolkit
alphaSgsJayatillekeWallFunctionFvPatchScalarField.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) 2008-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>
32 
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37 namespace compressible
38 {
39 namespace LESModels
40 {
41 
42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43 
44 scalar alphaSgsJayatillekeWallFunctionFvPatchScalarField::maxExp_ = 50.0;
45 scalar alphaSgsJayatillekeWallFunctionFvPatchScalarField::tolerance_ = 0.01;
46 label alphaSgsJayatillekeWallFunctionFvPatchScalarField::maxIters_ = 10;
47 
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 
50 void alphaSgsJayatillekeWallFunctionFvPatchScalarField::checkType()
51 {
52  if (!isA<wallFvPatch>(patch()))
53  {
55  (
56  "alphaSgsJayatillekeWallFunctionFvPatchScalarField::checkType()"
57  )
58  << "Patch type for patch " << patch().name() << " must be wall\n"
59  << "Current patch type is " << patch().type() << nl
60  << exit(FatalError);
61  }
62 }
63 
64 
65 scalar alphaSgsJayatillekeWallFunctionFvPatchScalarField::Psmooth
66 (
67  const scalar Prat
68 ) const
69 {
70  return 9.24*(pow(Prat, 0.75) - 1.0)*(1.0 + 0.28*exp(-0.007*Prat));
71 }
72 
73 
74 scalar alphaSgsJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
75 (
76  const scalar P,
77  const scalar Prat
78 ) const
79 {
80  scalar ypt = 11.0;
81 
82  for (int i=0; i<maxIters_; i++)
83  {
84  scalar f = ypt - (log(E_*ypt)/kappa_ + P)/Prat;
85  scalar df = 1.0 - 1.0/(ypt*kappa_*Prat);
86  scalar yptNew = ypt - f/df;
87 
88  if (yptNew < VSMALL)
89  {
90  return 0;
91  }
92  else if (mag(yptNew - ypt) < tolerance_)
93  {
94  return yptNew;
95  }
96  else
97  {
98  ypt = yptNew;
99  }
100  }
101 
102  return ypt;
103 }
104 
105 
106 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
107 
108 alphaSgsJayatillekeWallFunctionFvPatchScalarField::
109 alphaSgsJayatillekeWallFunctionFvPatchScalarField
110 (
111  const fvPatch& p,
113 )
114 :
115  fixedValueFvPatchScalarField(p, iF),
116  Prt_(0.85),
117  kappa_(0.41),
118  E_(9.8),
119  hsName_("hs")
120 {
121  checkType();
122 }
123 
124 
125 alphaSgsJayatillekeWallFunctionFvPatchScalarField::
126 alphaSgsJayatillekeWallFunctionFvPatchScalarField
127 (
129  const fvPatch& p,
131  const fvPatchFieldMapper& mapper
132 )
133 :
134  fixedValueFvPatchScalarField(ptf, p, iF, mapper),
135  Prt_(ptf.Prt_),
136  kappa_(ptf.kappa_),
137  E_(ptf.E_),
138  hsName_(ptf.hsName_)
139 {}
140 
141 
142 alphaSgsJayatillekeWallFunctionFvPatchScalarField::
143 alphaSgsJayatillekeWallFunctionFvPatchScalarField
144 (
145  const fvPatch& p,
147  const dictionary& dict
148 )
149 :
150  fixedValueFvPatchScalarField(p, iF, dict),
151  Prt_(dict.lookupOrDefault<scalar>("Prt", 0.85)),
152  kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
153  E_(dict.lookupOrDefault<scalar>("E", 9.8)),
154  hsName_(dict.lookupOrDefault<word>("hs", "hs"))
155 {
156  checkType();
157 }
158 
159 
160 alphaSgsJayatillekeWallFunctionFvPatchScalarField::
161 alphaSgsJayatillekeWallFunctionFvPatchScalarField
162 (
164 )
165 :
166  fixedValueFvPatchScalarField(awfpsf),
167  Prt_(awfpsf.Prt_),
168  kappa_(awfpsf.kappa_),
169  E_(awfpsf.E_),
170  hsName_(awfpsf.hsName_)
171 {
172  checkType();
173 }
174 
175 
176 alphaSgsJayatillekeWallFunctionFvPatchScalarField::
177 alphaSgsJayatillekeWallFunctionFvPatchScalarField
178 (
181 )
182 :
183  fixedValueFvPatchScalarField(awfpsf, iF),
184  Prt_(awfpsf.Prt_),
185  kappa_(awfpsf.kappa_),
186  E_(awfpsf.E_),
187  hsName_(awfpsf.hsName_)
188 {
189  checkType();
190 }
191 
192 
193 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
194 
195 void alphaSgsJayatillekeWallFunctionFvPatchScalarField::evaluate
196 (
197  const Pstream::commsTypes
198 )
199 {
200  const LESModel& lesModel = db().lookupObject<LESModel>("LESProperties");
201 
202  // Field data
203  const label patchI = patch().index();
204 
205  const scalarField& muw = lesModel.mu().boundaryField()[patchI];
206  const scalarField muSgsw = lesModel.muSgs()().boundaryField()[patchI];
207 
208  const scalarField& alphaw = lesModel.alpha().boundaryField()[patchI];
209  scalarField& alphaSgsw = *this;
210 
211  const fvPatchVectorField& Uw = lesModel.U().boundaryField()[patchI];
212  const scalarField magUp = mag(Uw.patchInternalField() - Uw);
213  const scalarField magGradUw = mag(Uw.snGrad());
214 
215  const scalarField& rhow = lesModel.rho().boundaryField()[patchI];
216  const fvPatchScalarField& hw =
217  patch().lookupPatchField<volScalarField, scalar>(hsName_);
218 
219  const scalarField& ry = patch().deltaCoeffs();
220 
221  // Heat flux [W/m2] - lagging alphaSgsw
222  const scalarField qDot = (alphaw + alphaSgsw)*hw.snGrad();
223 
224  // Populate boundary values
225  forAll(alphaSgsw, faceI)
226  {
227  // Calculate uTau using Newton-Raphson iteration
228  scalar uTau =
229  sqrt((muSgsw[faceI] + muw[faceI])/rhow[faceI]*magGradUw[faceI]);
230 
231  if (uTau > ROOTVSMALL)
232  {
233  label iter = 0;
234  scalar err = GREAT;
235 
236  do
237  {
238  scalar kUu = min(kappa_*magUp[faceI]/uTau, maxExp_);
239  scalar fkUu = exp(kUu) - 1.0 - kUu*(1.0 + 0.5*kUu);
240 
241  scalar f =
242  - uTau/(ry[faceI]*muw[faceI]/rhow[faceI])
243  + magUp[faceI]/uTau
244  + 1.0/E_*(fkUu - 1.0/6.0*kUu*sqr(kUu));
245 
246  scalar df =
247  - 1.0/(ry[faceI]*muw[faceI]/rhow[faceI])
248  - magUp[faceI]/sqr(uTau)
249  - 1.0/E_*kUu*fkUu/uTau;
250 
251  scalar uTauNew = uTau - f/df;
252  err = mag((uTau - uTauNew)/uTau);
253  uTau = uTauNew;
254 
255  } while (uTau>VSMALL && err>tolerance_ && ++iter<maxIters_);
256 
257  scalar yPlus = uTau/ry[faceI]/(muw[faceI]/rhow[faceI]);
258 
259  // Molecular Prandtl number
260  scalar Pr = muw[faceI]/alphaw[faceI];
261 
262  // Molecular-to-turbulenbt Prandtl number ratio
263  scalar Prat = Pr/Prt_;
264 
265  // Thermal sublayer thickness
266  scalar P = Psmooth(Prat);
267  scalar yPlusTherm = this->yPlusTherm(P, Prat);
268 
269  // Evaluate new effective thermal diffusivity
270  scalar alphaEff = 0.0;
271  if (yPlus < yPlusTherm)
272  {
273  scalar A = qDot[faceI]*rhow[faceI]*uTau/ry[faceI];
274  scalar B = qDot[faceI]*Pr*yPlus;
275  scalar C = Pr*0.5*rhow[faceI]*uTau*sqr(magUp[faceI]);
276  alphaEff = A/(B + C + VSMALL);
277  }
278  else
279  {
280  scalar A = qDot[faceI]*rhow[faceI]*uTau/ry[faceI];
281  scalar B = qDot[faceI]*Prt_*(1.0/kappa_*log(E_*yPlus) + P);
282  scalar magUc = uTau/kappa_*log(E_*yPlusTherm) - mag(Uw[faceI]);
283  scalar C =
284  0.5*rhow[faceI]*uTau
285  *(Prt_*sqr(magUp[faceI]) + (Pr - Prt_)*sqr(magUc));
286  alphaEff = A/(B + C + VSMALL);
287  }
288 
289  // Update turbulent thermal diffusivity
290  alphaSgsw[faceI] = max(0.0, alphaEff - alphaw[faceI]);
291 
292  if (debug)
293  {
294  Info<< " uTau = " << uTau << nl
295  << " Pr = " << Pr << nl
296  << " Prt = " << Prt_ << nl
297  << " qDot = " << qDot[faceI] << nl
298  << " yPlus = " << yPlus << nl
299  << " yPlusTherm = " << yPlusTherm << nl
300  << " alphaEff = " << alphaEff << nl
301  << " alphaw = " << alphaw[faceI] << nl
302  << " alphaSgsw = " << alphaSgsw[faceI] << nl
303  << endl;
304  }
305  }
306  else
307  {
308  alphaSgsw[faceI] = 0.0;
309  }
310  }
311 }
312 
313 
314 void alphaSgsJayatillekeWallFunctionFvPatchScalarField::write(Ostream& os) const
315 {
317  os.writeKeyword("Prt") << Prt_ << token::END_STATEMENT << nl;
318  os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
319  os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
320  os.writeKeyword("hs") << hsName_ << token::END_STATEMENT << nl;
321  writeEntry("value", os);
322 }
323 
324 
325 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
326 
328 (
331 );
332 
333 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
334 
335 } // End namespace LESModels
336 } // End namespace compressible
337 } // End namespace Foam
338 
339 // ************************ vim: set sw=4 sts=4 et: ************************ //