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-2024 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  typedName("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  typedName("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  typedName("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  typedName("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  typedName("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  typedName("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(typedName("g"), r + Cw2_*(pow6(r) - r));
157 
159  (
160  typedName("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  typedName("dTilda"),
178  min(CDES_*this->delta()(), this->y())
179  );
180 }
181 
182 
183 template<class BasicMomentumTransportModel>
185 (
186  const volScalarField::Internal& dTilda
187 ) const
188 {
189  if (this->mesh_.cacheTemporaryObject(typedName("LESRegion")))
190  {
192  (
193  typedName("LESRegion"),
194  neg(dTilda - this->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 viscosity& viscosity,
230  const word& type
231 )
232 :
233  LESeddyViscosity<BasicMomentumTransportModel>
234  (
235  type,
236  alpha,
237  rho,
238  U,
239  alphaRhoPhi,
240  phi,
241  viscosity
242  ),
243 
244  sigmaNut_
245  (
246  dimensioned<scalar>::lookupOrAddToDict
247  (
248  "sigmaNut",
249  this->coeffDict_,
250  0.66666
251  )
252  ),
253  kappa_
254  (
255  dimensioned<scalar>::lookupOrAddToDict
256  (
257  "kappa",
258  this->coeffDict_,
259  0.41
260  )
261  ),
262  Cb1_
263  (
264  dimensioned<scalar>::lookupOrAddToDict
265  (
266  "Cb1",
267  this->coeffDict_,
268  0.1355
269  )
270  ),
271  Cb2_
272  (
273  dimensioned<scalar>::lookupOrAddToDict
274  (
275  "Cb2",
276  this->coeffDict_,
277  0.622
278  )
279  ),
280  Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
281  Cw2_
282  (
283  dimensioned<scalar>::lookupOrAddToDict
284  (
285  "Cw2",
286  this->coeffDict_,
287  0.3
288  )
289  ),
290  Cw3_
291  (
292  dimensioned<scalar>::lookupOrAddToDict
293  (
294  "Cw3",
295  this->coeffDict_,
296  2.0
297  )
298  ),
299  Cv1_
300  (
301  dimensioned<scalar>::lookupOrAddToDict
302  (
303  "Cv1",
304  this->coeffDict_,
305  7.1
306  )
307  ),
308  Cs_
309  (
310  dimensioned<scalar>::lookupOrAddToDict
311  (
312  "Cs",
313  this->coeffDict_,
314  0.3
315  )
316  ),
317  CDES_
318  (
319  dimensioned<scalar>::lookupOrAddToDict
320  (
321  "CDES",
322  this->coeffDict_,
323  0.65
324  )
325  ),
326  ck_
327  (
328  dimensioned<scalar>::lookupOrAddToDict
329  (
330  "ck",
331  this->coeffDict_,
332  0.07
333  )
334  ),
335 
336  nuTilda_
337  (
338  IOobject
339  (
340  "nuTilda",
341  this->runTime_.name(),
342  this->mesh_,
343  IOobject::MUST_READ,
344  IOobject::AUTO_WRITE
345  ),
346  this->mesh_
347  )
348 {
349  if (type == typeName)
350  {
351  this->printCoeffs(type);
352  }
353 }
354 
355 
356 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
357 
358 template<class BasicMomentumTransportModel>
360 {
362  {
363  sigmaNut_.readIfPresent(this->coeffDict());
364  kappa_.readIfPresent(*this);
365 
366  Cb1_.readIfPresent(this->coeffDict());
367  Cb2_.readIfPresent(this->coeffDict());
368  Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
369  Cw2_.readIfPresent(this->coeffDict());
370  Cw3_.readIfPresent(this->coeffDict());
371  Cv1_.readIfPresent(this->coeffDict());
372  Cs_.readIfPresent(this->coeffDict());
373 
374  CDES_.readIfPresent(this->coeffDict());
375  ck_.readIfPresent(this->coeffDict());
376 
377  return true;
378  }
379  else
380  {
381  return false;
382  }
383 }
384 
385 
386 template<class BasicMomentumTransportModel>
388 DnuTildaEff() const
389 {
390  return volScalarField::New
391  (
392  "DnuTildaEff",
393  (nuTilda_ + this->nu())/sigmaNut_
394  );
395 }
396 
397 
398 template<class BasicMomentumTransportModel>
400 {
401  const volScalarField chi(this->chi());
402  const volScalarField fv1(this->fv1(chi));
403 
404  volScalarField dTildaExtrapolated
405  (
406  IOobject
407  (
408  "dTildaExtrapolated",
409  this->mesh_.time().name(),
410  this->mesh_
411  ),
412  this->mesh_,
413  dimLength,
414  zeroGradientFvPatchScalarField::typeName
415  );
416  dTildaExtrapolated.internalFieldRef() =
417  dTilda(chi, fv1, fvc::grad(this->U_));
418  dTildaExtrapolated.correctBoundaryConditions();
419 
421  (
423  (
424  typedName("k"),
425  sqr(this->nut()/ck_/dTildaExtrapolated)
426  )
427  );
428 
429  const fvPatchList& patches = this->mesh_.boundary();
430  volScalarField::Boundary& kBf = tk.ref().boundaryFieldRef();
431 
433  {
434  if (isA<wallFvPatch>(patches[patchi]))
435  {
436  kBf[patchi] = 0;
437  }
438  }
439 
440  return tk;
441 }
442 
443 
444 template<class BasicMomentumTransportModel>
446 {
447  if (!this->turbulence_)
448  {
449  return;
450  }
451 
452  // Local references
453  const alphaField& alpha = this->alpha_;
454  const rhoField& rho = this->rho_;
455  const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
456  const volVectorField& U = this->U_;
457  const Foam::fvModels& fvModels(Foam::fvModels::New(this->mesh_));
459  (
460  Foam::fvConstraints::New(this->mesh_)
461  );
462 
464 
465  const volScalarField chi(this->chi());
466  const volScalarField fv1(this->fv1(chi));
467 
468  tmp<volTensorField> tgradU = fvc::grad(U);
469  const volScalarField::Internal Omega(this->Omega(tgradU()));
470  const volScalarField::Internal dTilda(this->dTilda(chi, fv1, tgradU()));
471  const volScalarField::Internal Stilda
472  (
473  this->Stilda(chi, fv1, Omega, dTilda)
474  );
475  tgradU.clear();
476 
477  tmp<fvScalarMatrix> nuTildaEqn
478  (
479  fvm::ddt(alpha, rho, nuTilda_)
480  + fvm::div(alphaRhoPhi, nuTilda_)
481  - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
482  - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
483  ==
484  Cb1_*alpha()*rho()*Stilda*nuTilda_()
485  - fvm::Sp
486  (
487  Cw1_*alpha()*rho()*fw(Stilda, dTilda)*nuTilda_()/sqr(dTilda),
488  nuTilda_
489  )
490  + fvModels.source(alpha, rho, nuTilda_)
491  );
492 
493  nuTildaEqn.ref().relax();
494  fvConstraints.constrain(nuTildaEqn.ref());
495  solve(nuTildaEqn);
496  fvConstraints.constrain(nuTilda_);
497  bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), 0));
498  nuTilda_.correctBoundaryConditions();
499 
500  correctNut();
501 
502  // Optionally cache the LESRegion field
503  cacheLESRegion(dTilda);
504 }
505 
506 
507 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
508 
509 } // End namespace LESModels
510 } // End namespace Foam
511 
512 // ************************************************************************* //
scalar y
scalar delta
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
const dimensionSet & dimensions() const
Return dimensions.
static tmp< DimensionedField< Type, GeoMesh > > New(const word &name, const Mesh &mesh, const dimensionSet &, const Field< Type > &)
Return a temporary field constructed from name, mesh,.
Generic GeometricBoundaryField class.
Generic GeometricField class.
Internal & internalFieldRef()
Return a reference to the dimensioned internal field.
void correctBoundaryConditions()
Correct boundary field.
static tmp< GeometricField< Type, PatchField, GeoMesh > > New(const word &name, const Internal &, const PtrList< PatchField< Type >> &, const HashPtrTable< Source > &=HashPtrTable< Source >())
Return a temporary field constructed from name,.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
Eddy viscosity LES SGS model base class.
BasicMomentumTransportModel::alphaField alphaField
BasicMomentumTransportModel::rhoField rhoField
tmp< volScalarField > chi() const
tmp< volScalarField::Internal > Stilda(const volScalarField::Internal &chi, const volScalarField::Internal &fv1, const volScalarField::Internal &Omega, const volScalarField::Internal &dTilda) const
tmp< volScalarField::Internal > Omega(const volTensorField::Internal &gradU) const
virtual tmp< volScalarField > k() const
Return SGS kinetic energy.
virtual void correct()
Correct nuTilda and related properties.
tmp< volScalarField > DnuTildaEff() const
Return the effective diffusivity for nuTilda.
tmp< volScalarField > fv1(const volScalarField &chi) const
tmp< volScalarField::Internal > fv2(const volScalarField::Internal &chi, const volScalarField::Internal &fv1) const
virtual tmp< volScalarField::Internal > dTilda(const volScalarField::Internal &chi, const volScalarField::Internal &fv1, const volTensorField::Internal &gradU) const
Length scale.
SpalartAllmarasDES(const alphaField &alpha, const rhoField &rho, const volVectorField &U, const surfaceScalarField &alphaRhoPhi, const surfaceScalarField &phi, const viscosity &viscosity, const word &type=typeName)
Construct from components.
virtual void cacheLESRegion(const volScalarField::Internal &dTilda) const
Cache the LES region indicator field.
tmp< volScalarField::Internal > fw(const volScalarField::Internal &Omega, const volScalarField::Internal &dTilda) const
tmp< volScalarField::Internal > r(const volScalarField::Internal &nur, const volScalarField::Internal &Omega, const volScalarField::Internal &dTilda) const
virtual bool read()
Read model coefficients if they have changed.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: PtrList.H:75
static autoPtr< dictionary > New(Istream &)
Construct top-level dictionary on freestore from Istream.
Definition: dictionaryIO.C:100
Generic dimensioned Type class.
bool readIfPresent(const dictionary &)
Update the value of dimensioned<Type> if found in the dictionary.
virtual void correct()=0
Solve the turbulence equations and correct the turbulence viscosity.
Finite volume constraints.
Definition: fvConstraints.H:67
bool constrain(fvMatrix< Type > &eqn) const
Apply constraints to an equation.
Finite volume models.
Definition: fvModels.H:65
tmp< fvMatrix< Type > > source(const VolField< Type > &field) const
Return source for an equation.
A class for managing temporary objects.
Definition: tmp.H:55
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:237
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
Abstract base class for all fluid physical properties.
Definition: viscosity.H:50
A class for handling words, derived from string.
Definition: word.H:62
Foam::fvConstraints & fvConstraints(Foam::fvConstraints::New(mesh))
Foam::fvModels & fvModels(Foam::fvModels::New(mesh))
const scalar nut
label patchi
const fvPatchList & patches
U
Definition: pEqn.H:72
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
tmp< VolField< typename outerProduct< vector, Type >::type > > grad(const SurfaceField< Type > &ssf)
Definition: fvcGrad.C:46
tmp< fvMatrix< Type > > laplacian(const VolField< Type > &vf, const word &name)
Definition: fvmLaplacian.C:47
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const VolField< Type > &vf, const word &name)
Definition: fvmDiv.C:48
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const VolField< Type > &)
tmp< fvMatrix< Type > > ddt(const VolField< Type > &vf)
Definition: fvmDdt.C:46
Namespace for OpenFOAM.
dimensionedScalar pow6(const dimensionedScalar &ds)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar pow3(const dimensionedScalar &ds)
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
const dimensionSet dimLength
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
word typedName(Name name)
Return the name of the object within the given type.
Definition: typeInfo.H:176
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
bool bound(volScalarField &, const dimensionedScalar &min)
Bound the given scalar field where it is below the specified min value.
Definition: bound.C:31
dimensionedScalar neg(const dimensionedScalar &ds)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
dimensionedTensor skew(const dimensionedTensor &dt)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.