SpalartAllmaras.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 "SpalartAllmaras.H"
27 #include "fvModels.H"
28 #include "fvConstraints.H"
29 #include "bound.H"
30 #include "wallDist.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 namespace RASModels
37 {
38 
39 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
40 
41 template<class BasicMomentumTransportModel>
43 {
44  return volScalarField::New(modelName("chi"), nuTilda_/this->nu());
45 }
46 
47 
48 template<class BasicMomentumTransportModel>
50 (
51  const volScalarField& chi
52 ) const
53 {
54  const volScalarField chi3(modelName("chi3"), pow3(chi));
55  return volScalarField::New(modelName("fv1"), chi3/(chi3 + pow3(Cv1_)));
56 }
57 
58 
59 template<class BasicMomentumTransportModel>
61 (
62  const volScalarField::Internal& chi,
63  const volScalarField::Internal& fv1
64 ) const
65 {
67  (
68  modelName("fv2"),
69  1.0 - chi/(1.0 + chi*fv1)
70  );
71 }
72 
73 
74 template<class BasicMomentumTransportModel>
77 (
78  const volScalarField::Internal& chi,
79  const volScalarField::Internal& fv1
80 ) const
81 {
82  const volScalarField::Internal Omega
83  (
84  modelName("Omega"),
85  ::sqrt(2.0)*mag(skew(fvc::grad(this->U_)().v()))
86  );
87 
89  (
90  modelName("Stilda"),
91  (
92  max
93  (
94  Omega
95  + fv2(chi, fv1)*nuTilda_/sqr(kappa_*y_),
96  Cs_*Omega
97  )
98  )
99  );
100 }
101 
102 
103 template<class BasicMomentumTransportModel>
105 (
106  const volScalarField::Internal& Stilda
107 ) const
108 {
110  (
111  modelName("r"),
112  min
113  (
114  nuTilda_()
115  /(
116  max
117  (
118  Stilda,
119  dimensionedScalar(Stilda.dimensions(), small)
120  )
121  *sqr(kappa_*y_)
122  ),
123  scalar(10.0)
124  )
125  );
126 
127  const volScalarField::Internal g(modelName("g"), r + Cw2_*(pow6(r) - r));
128 
130  (
131  modelName("fw"),
132  g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0)
133  );
134 }
135 
136 
137 template<class BasicMomentumTransportModel>
139 (
140  const volScalarField& fv1
141 )
142 {
143  this->nut_ = nuTilda_*fv1;
144  this->nut_.correctBoundaryConditions();
145  fvConstraints::New(this->mesh_).constrain(this->nut_);
146 }
147 
148 
149 template<class BasicMomentumTransportModel>
151 {
152  correctNut(fv1(this->chi()));
153 }
154 
155 
156 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
157 
158 template<class BasicMomentumTransportModel>
160 (
161  const alphaField& alpha,
162  const rhoField& rho,
163  const volVectorField& U,
164  const surfaceScalarField& alphaRhoPhi,
165  const surfaceScalarField& phi,
166  const transportModel& transport,
167  const word& type
168 )
169 :
171  (
172  type,
173  alpha,
174  rho,
175  U,
176  alphaRhoPhi,
177  phi,
178  transport
179  ),
180 
181  sigmaNut_
182  (
184  (
185  "sigmaNut",
186  this->coeffDict_,
187  0.66666
188  )
189  ),
190  kappa_
191  (
193  (
194  "kappa",
195  this->coeffDict_,
196  0.41
197  )
198  ),
199 
200  Cb1_
201  (
203  (
204  "Cb1",
205  this->coeffDict_,
206  0.1355
207  )
208  ),
209  Cb2_
210  (
212  (
213  "Cb2",
214  this->coeffDict_,
215  0.622
216  )
217  ),
218  Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
219  Cw2_
220  (
222  (
223  "Cw2",
224  this->coeffDict_,
225  0.3
226  )
227  ),
228  Cw3_
229  (
231  (
232  "Cw3",
233  this->coeffDict_,
234  2.0
235  )
236  ),
237  Cv1_
238  (
240  (
241  "Cv1",
242  this->coeffDict_,
243  7.1
244  )
245  ),
246  Cs_
247  (
249  (
250  "Cs",
251  this->coeffDict_,
252  0.3
253  )
254  ),
255 
256  nuTilda_
257  (
258  IOobject
259  (
260  "nuTilda",
261  this->runTime_.timeName(),
262  this->mesh_,
265  ),
266  this->mesh_
267  ),
268 
269  y_(wallDist::New(this->mesh_).y())
270 {
271  if (type == typeName)
272  {
273  this->printCoeffs(type);
274  }
275 }
276 
277 
278 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
279 
280 template<class BasicMomentumTransportModel>
282 {
284  {
285  sigmaNut_.readIfPresent(this->coeffDict());
286  kappa_.readIfPresent(this->coeffDict());
287 
288  Cb1_.readIfPresent(this->coeffDict());
289  Cb2_.readIfPresent(this->coeffDict());
290  Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
291  Cw2_.readIfPresent(this->coeffDict());
292  Cw3_.readIfPresent(this->coeffDict());
293  Cv1_.readIfPresent(this->coeffDict());
294  Cs_.readIfPresent(this->coeffDict());
295 
296  return true;
297  }
298  else
299  {
300  return false;
301  }
302 }
303 
304 
305 template<class BasicMomentumTransportModel>
308 {
309  return volScalarField::New
310  (
311  "DnuTildaEff",
312  (nuTilda_ + this->nu())/sigmaNut_
313  );
314 }
315 
316 
317 template<class BasicMomentumTransportModel>
319 {
320  return volScalarField::New
321  (
322  "k",
323  this->mesh_,
324  dimensionedScalar(dimensionSet(0, 2, -2, 0, 0), 0)
325  );
326 }
327 
328 
329 template<class BasicMomentumTransportModel>
332 {
334  << "Turbulence kinetic energy dissipation rate not defined for "
335  << "Spalart-Allmaras model. Returning zero field"
336  << endl;
337 
338  return volScalarField::New
339  (
340  "epsilon",
341  this->mesh_,
342  dimensionedScalar(dimensionSet(0, 2, -3, 0, 0), 0)
343  );
344 }
345 
346 
347 template<class BasicMomentumTransportModel>
349 {
350  if (!this->turbulence_)
351  {
352  return;
353  }
354 
355  // Local references
356  const alphaField& alpha = this->alpha_;
357  const rhoField& rho = this->rho_;
358  const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
359  const Foam::fvModels& fvModels(Foam::fvModels::New(this->mesh_));
361  (
362  Foam::fvConstraints::New(this->mesh_)
363  );
364 
366 
367  const volScalarField chi(this->chi());
368  const volScalarField fv1(this->fv1(chi));
369 
370  const volScalarField::Internal Stilda(this->Stilda(chi, fv1));
371 
372  tmp<fvScalarMatrix> nuTildaEqn
373  (
374  fvm::ddt(alpha, rho, nuTilda_)
375  + fvm::div(alphaRhoPhi, nuTilda_)
376  - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
377  - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
378  ==
379  Cb1_*alpha()*rho()*Stilda*nuTilda_()
380  - fvm::Sp(Cw1_*alpha()*rho()*fw(Stilda)*nuTilda_()/sqr(y_), nuTilda_)
381  + fvModels.source(alpha, rho, nuTilda_)
382  );
383 
384  nuTildaEqn.ref().relax();
385  fvConstraints.constrain(nuTildaEqn.ref());
386  solve(nuTildaEqn);
387  fvConstraints.constrain(nuTilda_);
388  bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), 0));
389  nuTilda_.correctBoundaryConditions();
390 
391  correctNut(fv1);
392 }
393 
394 
395 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
396 
397 } // End namespace RASModels
398 } // End namespace Foam
399 
400 // ************************************************************************* //
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
static tmp< DimensionedField< Type, GeoMesh > > New(const word &name, const Mesh &mesh, const dimensionSet &)
Return a temporary field constructed from name, mesh.
virtual tmp< volScalarField > epsilon() const
Return the turbulence kinetic energy dissipation rate.
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)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
tmp< volScalarField::Internal > fw(const volScalarField::Internal &Stilda) const
bool readIfPresent(const dictionary &)
Update the value of dimensioned<Type> if found in the dictionary.
Generic dimensioned Type class.
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
Eddy viscosity turbulence model base class.
Definition: eddyViscosity.H:49
tmp< volScalarField > chi() const
DimensionedField< scalar, volMesh > Internal
Type of the internal field from which this GeometricField is derived.
BasicMomentumTransportModel::alphaField alphaField
Definition: RASModel.H:88
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
Templated abstract base class for RAS turbulence models.
Definition: RASModel.H:49
tmp< volScalarField::Internal > Stilda(const volScalarField::Internal &chi, const volScalarField::Internal &fv1) const
scalar y
SpalartAllmaras(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.
Dimension set for the base types.
Definition: dimensionSet.H:120
bool read(const char *, int32_t &)
Definition: int32IO.C:85
BasicMomentumTransportModel::transportModel transportModel
Definition: RASModel.H:90
Foam::fvConstraints & fvConstraints
A class for handling words, derived from string.
Definition: word.H:59
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-9/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
Definition: pEqn.H:68
Finite volume models.
Definition: fvModels.H:60
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
phi
Definition: correctPhi.H:3
tmp< volScalarField > DnuTildaEff() const
Return the effective diffusivity for nuTilda.
virtual bool read()
Read RASProperties dictionary.
tmp< volScalarField::Internal > fv2(const volScalarField::Internal &chi, const volScalarField::Internal &fv1) const
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
Bound the given scalar field if it has gone unbounded.
static autoPtr< dictionary > New(Istream &)
Construct top-level dictionary on freestore from Istream.
Definition: dictionaryIO.C:96
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
bool constrain(fvMatrix< Type > &eqn) const
Apply constraints to an equation.
Foam::fvModels & fvModels
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
BasicMomentumTransportModel::rhoField rhoField
Definition: RASModel.H:89
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
#define WarningInFunction
Report a warning using Foam::Warning.
Finite volume constraints.
Definition: fvConstraints.H:57
tmp< volScalarField > fv1(const volScalarField &chi) const
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
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
virtual tmp< volScalarField > k() const
Return the turbulence kinetic energy.
void correctBoundaryConditions()
Correct boundary field.
dimensionedScalar pow6(const dimensionedScalar &ds)
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
Namespace for OpenFOAM.