SpalartAllmarasDES.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-2021 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 
26 #include "SpalartAllmarasDES.H"
27 #include "fvModels.H"
28 #include "fvConstraints.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 namespace LESModels
35 {
36 
37 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
38 
39 template<class BasicMomentumTransportModel>
41 {
42  return volScalarField::New
43  (
44  modelName("chi"),
45  nuTilda_/this->nu()
46  );
47 }
48 
49 
50 template<class BasicMomentumTransportModel>
52 (
53  const volScalarField& chi
54 ) const
55 {
56  const volScalarField chi3("chi3", pow3(chi));
57  return volScalarField::New
58  (
59  modelName("fv1"),
60  chi3/(chi3 + pow3(Cv1_))
61  );
62 }
63 
64 
65 template<class BasicMomentumTransportModel>
68 (
69  const volScalarField::Internal& chi,
70  const volScalarField::Internal& fv1
71 ) const
72 {
74  (
75  modelName("fv2"),
76  1.0 - chi/(1.0 + chi*fv1)
77  );
78 }
79 
80 
81 template<class BasicMomentumTransportModel>
84 (
85  const volTensorField::Internal& gradU
86 ) const
87 {
89  (
90  modelName("Omega"),
91  sqrt(2.0)*mag(skew(gradU))
92  );
93 }
94 
95 
96 template<class BasicMomentumTransportModel>
99 (
100  const volScalarField::Internal& chi,
101  const volScalarField::Internal& fv1,
102  const volScalarField::Internal& Omega,
103  const volScalarField::Internal& dTilda
104 ) const
105 {
107  (
108  modelName("Stilda"),
109  max
110  (
111  Omega
112  + fv2(chi, fv1)*nuTilda_()/sqr(kappa_*dTilda),
113  Cs_*Omega
114  )
115  );
116 }
117 
118 
119 template<class BasicMomentumTransportModel>
121 (
122  const volScalarField::Internal& nur,
123  const volScalarField::Internal& Omega,
124  const volScalarField::Internal& dTilda
125 ) const
126 {
128  (
129  modelName("r"),
130  min
131  (
132  nur
133  /(
134  max
135  (
136  Omega,
137  dimensionedScalar(Omega.dimensions(), small)
138  )
139  *sqr(kappa_*dTilda)
140  ),
141  scalar(10)
142  )
143  );
144 }
145 
146 
147 template<class BasicMomentumTransportModel>
150 (
151  const volScalarField::Internal& Omega,
152  const volScalarField::Internal& dTilda
153 ) const
154 {
155  const volScalarField::Internal r(this->r(nuTilda_, Omega, dTilda));
156  const volScalarField::Internal g(modelName("g"), r + Cw2_*(pow6(r) - r));
157 
159  (
160  modelName("fw"),
161  g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0)
162  );
163 }
164 
165 
166 template<class BasicMomentumTransportModel>
169 (
170  const volScalarField::Internal& chi,
171  const volScalarField::Internal& fv1,
172  const volTensorField::Internal& gradU
173 ) const
174 {
176  (
177  modelName("dTilda"),
178  min(CDES_*this->delta()(), y_)
179  );
180 }
181 
182 
183 template<class BasicMomentumTransportModel>
185 (
186  const volScalarField::Internal& dTilda
187 ) const
188 {
189  if (this->mesh_.cacheTemporaryObject(modelName("LESRegion")))
190  {
192  (
193  modelName("LESRegion"),
194  neg(dTilda - y_())
195  );
196  }
197 }
198 
199 
200 template<class BasicMomentumTransportModel>
202 (
203  const volScalarField& fv1
204 )
205 {
206  this->nut_ = nuTilda_*fv1;
207  this->nut_.correctBoundaryConditions();
208  fvConstraints::New(this->mesh_).constrain(this->nut_);
209 }
210 
211 
212 template<class BasicMomentumTransportModel>
214 {
215  correctNut(fv1(this->chi()));
216 }
217 
218 
219 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
220 
221 template<class BasicMomentumTransportModel>
223 (
224  const alphaField& alpha,
225  const rhoField& rho,
226  const volVectorField& U,
227  const surfaceScalarField& alphaRhoPhi,
228  const surfaceScalarField& phi,
229  const transportModel& transport,
230  const word& type
231 )
232 :
234  (
235  type,
236  alpha,
237  rho,
238  U,
239  alphaRhoPhi,
240  phi,
241  transport
242  ),
243 
244  sigmaNut_
245  (
247  (
248  "sigmaNut",
249  this->coeffDict_,
250  0.66666
251  )
252  ),
253  kappa_
254  (
256  (
257  "kappa",
258  this->coeffDict_,
259  0.41
260  )
261  ),
262  Cb1_
263  (
265  (
266  "Cb1",
267  this->coeffDict_,
268  0.1355
269  )
270  ),
271  Cb2_
272  (
274  (
275  "Cb2",
276  this->coeffDict_,
277  0.622
278  )
279  ),
280  Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
281  Cw2_
282  (
284  (
285  "Cw2",
286  this->coeffDict_,
287  0.3
288  )
289  ),
290  Cw3_
291  (
293  (
294  "Cw3",
295  this->coeffDict_,
296  2.0
297  )
298  ),
299  Cv1_
300  (
302  (
303  "Cv1",
304  this->coeffDict_,
305  7.1
306  )
307  ),
308  Cs_
309  (
311  (
312  "Cs",
313  this->coeffDict_,
314  0.3
315  )
316  ),
317  CDES_
318  (
320  (
321  "CDES",
322  this->coeffDict_,
323  0.65
324  )
325  ),
326  ck_
327  (
329  (
330  "ck",
331  this->coeffDict_,
332  0.07
333  )
334  ),
335 
336  nuTilda_
337  (
338  IOobject
339  (
340  "nuTilda",
341  this->runTime_.timeName(),
342  this->mesh_,
345  ),
346  this->mesh_
347  ),
348 
349  y_(wallDist::New(this->mesh_).y())
350 {
351  if (type == typeName)
352  {
353  this->printCoeffs(type);
354  }
355 }
356 
357 
358 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
359 
360 template<class BasicMomentumTransportModel>
362 {
364  {
365  sigmaNut_.readIfPresent(this->coeffDict());
366  kappa_.readIfPresent(*this);
367 
368  Cb1_.readIfPresent(this->coeffDict());
369  Cb2_.readIfPresent(this->coeffDict());
370  Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
371  Cw2_.readIfPresent(this->coeffDict());
372  Cw3_.readIfPresent(this->coeffDict());
373  Cv1_.readIfPresent(this->coeffDict());
374  Cs_.readIfPresent(this->coeffDict());
375 
376  CDES_.readIfPresent(this->coeffDict());
377  ck_.readIfPresent(this->coeffDict());
378 
379  return true;
380  }
381  else
382  {
383  return false;
384  }
385 }
386 
387 
388 template<class BasicMomentumTransportModel>
390 DnuTildaEff() const
391 {
392  return volScalarField::New
393  (
394  "DnuTildaEff",
395  (nuTilda_ + this->nu())/sigmaNut_
396  );
397 }
398 
399 
400 template<class BasicMomentumTransportModel>
402 {
403  const volScalarField chi(this->chi());
404  const volScalarField fv1(this->fv1(chi));
405 
406  volScalarField dTildaExtrapolated
407  (
408  IOobject
409  (
410  "dTildaExtrapolated",
411  this->mesh_.time().timeName(),
412  this->mesh_
413  ),
414  this->mesh_,
415  dimLength,
416  zeroGradientFvPatchScalarField::typeName
417  );
418  dTildaExtrapolated.ref() = dTilda(chi, fv1, fvc::grad(this->U_));
419  dTildaExtrapolated.correctBoundaryConditions();
420 
422  (
424  (
425  modelName("k"),
426  sqr(this->nut()/ck_/dTildaExtrapolated)
427  )
428  );
429 
430  const fvPatchList& patches = this->mesh_.boundary();
431  volScalarField::Boundary& kBf = tk.ref().boundaryFieldRef();
432 
433  forAll(patches, patchi)
434  {
435  if (isA<wallFvPatch>(patches[patchi]))
436  {
437  kBf[patchi] = 0;
438  }
439  }
440 
441  return tk;
442 }
443 
444 
445 template<class BasicMomentumTransportModel>
447 {
448  if (!this->turbulence_)
449  {
450  return;
451  }
452 
453  // Local references
454  const alphaField& alpha = this->alpha_;
455  const rhoField& rho = this->rho_;
456  const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
457  const volVectorField& U = this->U_;
458  const Foam::fvModels& fvModels(Foam::fvModels::New(this->mesh_));
460  (
461  Foam::fvConstraints::New(this->mesh_)
462  );
463 
465 
466  const volScalarField chi(this->chi());
467  const volScalarField fv1(this->fv1(chi));
468 
469  tmp<volTensorField> tgradU = fvc::grad(U);
470  const volScalarField::Internal Omega(this->Omega(tgradU()));
471  const volScalarField::Internal dTilda(this->dTilda(chi, fv1, tgradU()));
472  const volScalarField::Internal Stilda
473  (
474  this->Stilda(chi, fv1, Omega, dTilda)
475  );
476  tgradU.clear();
477 
478  tmp<fvScalarMatrix> nuTildaEqn
479  (
480  fvm::ddt(alpha, rho, nuTilda_)
481  + fvm::div(alphaRhoPhi, nuTilda_)
482  - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
483  - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
484  ==
485  Cb1_*alpha()*rho()*Stilda*nuTilda_()
486  - fvm::Sp
487  (
488  Cw1_*alpha()*rho()*fw(Stilda, dTilda)*nuTilda_()/sqr(dTilda),
489  nuTilda_
490  )
491  + fvModels.source(alpha, rho, nuTilda_)
492  );
493 
494  nuTildaEqn.ref().relax();
495  fvConstraints.constrain(nuTildaEqn.ref());
496  solve(nuTildaEqn);
497  fvConstraints.constrain(nuTilda_);
498  bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), 0));
499  nuTilda_.correctBoundaryConditions();
500 
501  correctNut();
502 
503  // Optionally cache the LESRegion field
504  cacheLESRegion(dTilda);
505 }
506 
507 
508 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
509 
510 } // End namespace LESModels
511 } // End namespace Foam
512 
513 // ************************************************************************* //
Eddy viscosity LES SGS model base class.
scalar delta
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
BasicMomentumTransportModel::rhoField rhoField
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:237
static tmp< DimensionedField< Type, GeoMesh > > New(const word &name, const Mesh &mesh, const dimensionSet &)
Return a temporary field constructed from name, mesh.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedTensor skew(const dimensionedTensor &dt)
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
dimensionedSymmTensor sqr(const dimensionedVector &dv)
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< scalar >> &)
Return a temporary field constructed from name,.
dimensionedScalar sqrt(const dimensionedScalar &ds)
bool readIfPresent(const dictionary &)
Update the value of dimensioned<Type> if found in the dictionary.
virtual bool read()
Read model coefficients if they have changed.
Generic dimensioned Type class.
patches[0]
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
DimensionedField< scalar, volMesh > Internal
Type of the internal field from which this GeometricField is derived.
dimensionedScalar neg(const dimensionedScalar &ds)
const dimensionSet dimLength
virtual void correct()=0
Solve the turbulence equations and correct the turbulence viscosity.
scalar y
Foam::fvConstraints & fvConstraints
A class for handling words, derived from string.
Definition: word.H:59
Finite volume models.
Definition: fvModels.H:60
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
virtual tmp< volScalarField::Internal > dTilda(const volScalarField::Internal &chi, const volScalarField::Internal &fv1, const volTensorField::Internal &gradU) const
Length scale.
phi
Definition: correctPhi.H:3
dimensioned< scalar > magSqr(const dimensioned< Type > &)
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
tmp< volScalarField::Internal > r(const volScalarField::Internal &nur, const volScalarField::Internal &Omega, const volScalarField::Internal &dTilda) const
static autoPtr< dictionary > New(Istream &)
Construct top-level dictionary on freestore from Istream.
Definition: dictionaryIO.C:96
tmp< volScalarField > fv1(const volScalarField &chi) const
tmp< volScalarField > DnuTildaEff() const
Return the effective diffusivity for nuTilda.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
bool constrain(fvMatrix< Type > &eqn) const
Apply constraints to an equation.
Internal & ref()
Return a reference to the dimensioned internal field.
Foam::fvModels & fvModels
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
BasicMomentumTransportModel::alphaField alphaField
label patchi
U
Definition: pEqn.H:72
volScalarField & bound(volScalarField &, const dimensionedScalar &lowerBound)
Bound the given scalar field if it has gone unbounded.
Definition: bound.C:33
Finite volume constraints.
Definition: fvConstraints.H:57
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:46
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
SpalartAllmarasDES(const alphaField &alpha, const rhoField &rho, const volVectorField &U, const surfaceScalarField &alphaRhoPhi, const surfaceScalarField &phi, const transportModel &transport, const word &type=typeName)
Construct from components.
tmp< volScalarField::Internal > Omega(const volTensorField::Internal &gradU) const
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
BasicMomentumTransportModel::transportModel transportModel
virtual void cacheLESRegion(const volScalarField::Internal &dTilda) const
Cache the LES region indicator field.
void correctBoundaryConditions()
Correct boundary field.
dimensionedScalar pow6(const dimensionedScalar &ds)
tmp< volScalarField::Internal > Stilda(const volScalarField::Internal &chi, const volScalarField::Internal &fv1, const volScalarField::Internal &Omega, const volScalarField::Internal &dTilda) const
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.
dimensioned< scalar > mag(const dimensioned< Type > &)
A class for managing temporary objects.
Definition: PtrList.H:53
const dimensionedVector & g
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
tmp< volScalarField > chi() const
virtual void correct()
Correct nuTilda and related properties.
tmp< volScalarField::Internal > fw(const volScalarField::Internal &Omega, const volScalarField::Internal &dTilda) const
tmp< volScalarField::Internal > fv2(const volScalarField::Internal &chi, const volScalarField::Internal &fv1) const
scalar nut
Namespace for OpenFOAM.
virtual tmp< volScalarField > k() const
Return SGS kinetic energy.