FreeFOAM The Cross-Platform CFD Toolkit
SpalartAllmaras.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 
26 #include "SpalartAllmaras.H"
28 
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 namespace compressible
36 {
37 namespace RASModels
38 {
39 
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 
42 defineTypeNameAndDebug(SpalartAllmaras, 0);
43 addToRunTimeSelectionTable(RASModel, SpalartAllmaras, dictionary);
44 
45 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
46 
47 tmp<volScalarField> SpalartAllmaras::chi() const
48 {
49  return rho_*nuTilda_/mu();
50 }
51 
52 
53 tmp<volScalarField> SpalartAllmaras::fv1(const volScalarField& chi) const
54 {
55  volScalarField chi3 = pow3(chi);
56  return chi3/(chi3 + pow3(Cv1_));
57 }
58 
59 
60 tmp<volScalarField> SpalartAllmaras::fv2
61 (
62  const volScalarField& chi,
63  const volScalarField& fv1
64 ) const
65 {
66  return 1.0/pow3(scalar(1) + chi/Cv2_);
67 }
68 
69 
70 tmp<volScalarField> SpalartAllmaras::fv3
71 (
72  const volScalarField& chi,
73  const volScalarField& fv1
74 ) const
75 {
76  volScalarField chiByCv2 = (1/Cv2_)*chi;
77 
78  return
79  (scalar(1) + chi*fv1)
80  *(1/Cv2_)
81  *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
82  /pow3(scalar(1) + chiByCv2);
83 }
84 
85 
86 tmp<volScalarField> SpalartAllmaras::fw(const volScalarField& Stilda) const
87 {
89  (
90  nuTilda_
91  /(
92  max(Stilda, dimensionedScalar("SMALL", Stilda.dimensions(), SMALL))
93  *sqr(kappa_*d_)
94  ),
95  scalar(10.0)
96  );
97  r.boundaryField() == 0.0;
98 
99  volScalarField g = r + Cw2_*(pow6(r) - r);
100 
101  return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
102 }
103 
104 
105 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
106 
107 SpalartAllmaras::SpalartAllmaras
108 (
109  const volScalarField& rho,
110  const volVectorField& U,
111  const surfaceScalarField& phi,
112  const basicThermo& thermophysicalModel
113 )
114 :
115  RASModel(typeName, rho, U, phi, thermophysicalModel),
116 
117  sigmaNut_
118  (
120  (
121  "sigmaNut",
122  coeffDict_,
123  0.66666
124  )
125  ),
126  kappa_
127  (
129  (
130  "kappa",
131  coeffDict_,
132  0.41
133  )
134  ),
135  Prt_
136  (
138  (
139  "Prt",
140  coeffDict_,
141  1.0
142  )
143  ),
144 
145  Cb1_
146  (
148  (
149  "Cb1",
150  coeffDict_,
151  0.1355
152  )
153  ),
154  Cb2_
155  (
157  (
158  "Cb2",
159  coeffDict_,
160  0.622
161  )
162  ),
163  Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
164  Cw2_
165  (
167  (
168  "Cw2",
169  coeffDict_,
170  0.3
171  )
172  ),
173  Cw3_
174  (
176  (
177  "Cw3",
178  coeffDict_,
179  2.0
180  )
181  ),
182  Cv1_
183  (
185  (
186  "Cv1",
187  coeffDict_,
188  7.1
189  )
190  ),
191  Cv2_
192  (
194  (
195  "Cv2",
196  coeffDict_,
197  5.0
198  )
199  ),
200 
201  nuTilda_
202  (
203  IOobject
204  (
205  "nuTilda",
206  runTime_.timeName(),
207  mesh_,
210  ),
211  mesh_
212  ),
213 
214  mut_
215  (
216  IOobject
217  (
218  "mut",
219  runTime_.timeName(),
220  mesh_,
223  ),
224  mesh_
225  ),
226 
227  alphat_
228  (
229  IOobject
230  (
231  "alphat",
232  runTime_.timeName(),
233  mesh_,
236  ),
237  autoCreateAlphat("alphat", mesh_)
238  ),
239 
240  d_(mesh_)
241 {
242  alphat_ = mut_/Prt_;
243  alphat_.correctBoundaryConditions();
244 
245  printCoeffs();
246 }
247 
248 
249 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
250 
251 tmp<volSymmTensorField> SpalartAllmaras::R() const
252 {
254  (
256  (
257  IOobject
258  (
259  "R",
260  runTime_.timeName(),
261  mesh_,
264  ),
265  ((2.0/3.0)*I)*k() - (mut_/rho_)*dev(twoSymm(fvc::grad(U_)))
266  )
267  );
268 }
269 
270 
271 tmp<volSymmTensorField> SpalartAllmaras::devRhoReff() const
272 {
274  (
276  (
277  IOobject
278  (
279  "devRhoReff",
280  runTime_.timeName(),
281  mesh_,
284  ),
285  -muEff()*dev(twoSymm(fvc::grad(U_)))
286  )
287  );
288 }
289 
290 
291 tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff(volVectorField& U) const
292 {
293  volScalarField muEff_ = muEff();
294 
295  return
296  (
297  - fvm::laplacian(muEff_, U)
298  - fvc::div(muEff_*dev2(fvc::grad(U)().T()))
299  );
300 }
301 
302 
303 bool SpalartAllmaras::read()
304 {
305  if (RASModel::read())
306  {
307  sigmaNut_.readIfPresent(coeffDict());
308  kappa_.readIfPresent(coeffDict());
309  Prt_.readIfPresent(coeffDict());
310 
311  Cb1_.readIfPresent(coeffDict());
312  Cb2_.readIfPresent(coeffDict());
313  Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
314  Cw2_.readIfPresent(coeffDict());
315  Cw3_.readIfPresent(coeffDict());
316  Cv1_.readIfPresent(coeffDict());
317  Cv2_.readIfPresent(coeffDict());
318 
319  return true;
320  }
321  else
322  {
323  return false;
324  }
325 }
326 
327 
329 {
330  if (!turbulence_)
331  {
332  // Re-calculate viscosity
333  mut_ = rho_*nuTilda_*fv1(chi());
334  mut_.correctBoundaryConditions();
335 
336  // Re-calculate thermal diffusivity
337  alphat_ = mut_/Prt_;
338  alphat_.correctBoundaryConditions();
339 
340  return;
341  }
342 
344 
345  if (mesh_.changing())
346  {
347  d_.correct();
348  }
349 
350  volScalarField chi = this->chi();
351  volScalarField fv1 = this->fv1(chi);
352 
353  volScalarField Stilda =
354  fv3(chi, fv1)*::sqrt(2.0)*mag(skew(fvc::grad(U_)))
355  + fv2(chi, fv1)*nuTilda_/sqr(kappa_*d_);
356 
357  tmp<fvScalarMatrix> nuTildaEqn
358  (
359  fvm::ddt(rho_, nuTilda_)
360  + fvm::div(phi_, nuTilda_)
361  - fvm::laplacian(DnuTildaEff(), nuTilda_)
362  - Cb2_/sigmaNut_*rho_*magSqr(fvc::grad(nuTilda_))
363  ==
364  Cb1_*rho_*Stilda*nuTilda_
365  - fvm::Sp(Cw1_*fw(Stilda)*nuTilda_*rho_/sqr(d_), nuTilda_)
366  );
367 
368  nuTildaEqn().relax();
369  solve(nuTildaEqn);
370  bound(nuTilda_, dimensionedScalar("0", nuTilda_.dimensions(), 0.0));
371  nuTilda_.correctBoundaryConditions();
372 
373  // Re-calculate viscosity
374  mut_.internalField() = fv1*nuTilda_.internalField()*rho_.internalField();
375  mut_.correctBoundaryConditions();
376 
377  // Re-calculate thermal diffusivity
378  alphat_ = mut_/Prt_;
379  alphat_.correctBoundaryConditions();
380 }
381 
382 
383 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
384 
385 } // End namespace RASModels
386 } // End namespace compressible
387 } // End namespace Foam
388 
389 // ************************ vim: set sw=4 sts=4 et: ************************ //