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-2020 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 "fvOptions.H"
28 #include "bound.H"
29 #include "wallDist.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 namespace RASModels
36 {
37 
38 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
39 
40 template<class BasicMomentumTransportModel>
42 {
43  return volScalarField::New(modelName("chi"), nuTilda_/this->nu());
44 }
45 
46 
47 template<class BasicMomentumTransportModel>
49 (
50  const volScalarField& chi
51 ) const
52 {
53  const volScalarField chi3(modelName("chi3"), pow3(chi));
54  return volScalarField::New(modelName("fv1"), chi3/(chi3 + pow3(Cv1_)));
55 }
56 
57 
58 template<class BasicMomentumTransportModel>
60 (
61  const volScalarField::Internal& chi,
62  const volScalarField::Internal& fv1
63 ) const
64 {
66  (
67  modelName("fv2"),
68  1.0 - chi/(1.0 + chi*fv1)
69  );
70 }
71 
72 
73 template<class BasicMomentumTransportModel>
76 (
77  const volScalarField::Internal& chi,
78  const volScalarField::Internal& fv1
79 ) const
80 {
81  const volScalarField::Internal Omega
82  (
83  modelName("Omega"),
84  ::sqrt(2.0)*mag(skew(fvc::grad(this->U_)().v()))
85  );
86 
88  (
89  modelName("Stilda"),
90  (
91  max
92  (
93  Omega
94  + fv2(chi, fv1)*nuTilda_/sqr(kappa_*y_),
95  Cs_*Omega
96  )
97  )
98  );
99 }
100 
101 
102 template<class BasicMomentumTransportModel>
104 (
105  const volScalarField::Internal& Stilda
106 ) const
107 {
109  (
110  modelName("r"),
111  min
112  (
113  nuTilda_()
114  /(
115  max
116  (
117  Stilda,
118  dimensionedScalar(Stilda.dimensions(), small)
119  )
120  *sqr(kappa_*y_)
121  ),
122  scalar(10.0)
123  )
124  );
125 
126  const volScalarField::Internal g(modelName("g"), r + Cw2_*(pow6(r) - r));
127 
129  (
130  modelName("fw"),
131  g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0)
132  );
133 }
134 
135 
136 template<class BasicMomentumTransportModel>
138 (
139  const volScalarField& fv1
140 )
141 {
142  this->nut_ = nuTilda_*fv1;
143  this->nut_.correctBoundaryConditions();
144  fv::options::New(this->mesh_).correct(this->nut_);
145 }
146 
147 
148 template<class BasicMomentumTransportModel>
150 {
151  correctNut(fv1(this->chi()));
152 }
153 
154 
155 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
156 
157 template<class BasicMomentumTransportModel>
159 (
160  const alphaField& alpha,
161  const rhoField& rho,
162  const volVectorField& U,
163  const surfaceScalarField& alphaRhoPhi,
164  const surfaceScalarField& phi,
165  const transportModel& transport,
166  const word& type
167 )
168 :
170  (
171  type,
172  alpha,
173  rho,
174  U,
175  alphaRhoPhi,
176  phi,
177  transport
178  ),
179 
180  sigmaNut_
181  (
183  (
184  "sigmaNut",
185  this->coeffDict_,
186  0.66666
187  )
188  ),
189  kappa_
190  (
192  (
193  "kappa",
194  this->coeffDict_,
195  0.41
196  )
197  ),
198 
199  Cb1_
200  (
202  (
203  "Cb1",
204  this->coeffDict_,
205  0.1355
206  )
207  ),
208  Cb2_
209  (
211  (
212  "Cb2",
213  this->coeffDict_,
214  0.622
215  )
216  ),
217  Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
218  Cw2_
219  (
221  (
222  "Cw2",
223  this->coeffDict_,
224  0.3
225  )
226  ),
227  Cw3_
228  (
230  (
231  "Cw3",
232  this->coeffDict_,
233  2.0
234  )
235  ),
236  Cv1_
237  (
239  (
240  "Cv1",
241  this->coeffDict_,
242  7.1
243  )
244  ),
245  Cs_
246  (
248  (
249  "Cs",
250  this->coeffDict_,
251  0.3
252  )
253  ),
254 
255  nuTilda_
256  (
257  IOobject
258  (
259  "nuTilda",
260  this->runTime_.timeName(),
261  this->mesh_,
264  ),
265  this->mesh_
266  ),
267 
268  y_(wallDist::New(this->mesh_).y())
269 {
270  if (type == typeName)
271  {
272  this->printCoeffs(type);
273  }
274 }
275 
276 
277 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
278 
279 template<class BasicMomentumTransportModel>
281 {
283  {
284  sigmaNut_.readIfPresent(this->coeffDict());
285  kappa_.readIfPresent(this->coeffDict());
286 
287  Cb1_.readIfPresent(this->coeffDict());
288  Cb2_.readIfPresent(this->coeffDict());
289  Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
290  Cw2_.readIfPresent(this->coeffDict());
291  Cw3_.readIfPresent(this->coeffDict());
292  Cv1_.readIfPresent(this->coeffDict());
293  Cs_.readIfPresent(this->coeffDict());
294 
295  return true;
296  }
297  else
298  {
299  return false;
300  }
301 }
302 
303 
304 template<class BasicMomentumTransportModel>
307 {
308  return volScalarField::New
309  (
310  "DnuTildaEff",
311  (nuTilda_ + this->nu())/sigmaNut_
312  );
313 }
314 
315 
316 template<class BasicMomentumTransportModel>
318 {
319  return volScalarField::New
320  (
321  "k",
322  this->mesh_,
323  dimensionedScalar(dimensionSet(0, 2, -2, 0, 0), 0)
324  );
325 }
326 
327 
328 template<class BasicMomentumTransportModel>
331 {
333  << "Turbulence kinetic energy dissipation rate not defined for "
334  << "Spalart-Allmaras model. Returning zero field"
335  << endl;
336 
337  return volScalarField::New
338  (
339  "epsilon",
340  this->mesh_,
341  dimensionedScalar(dimensionSet(0, 2, -3, 0, 0), 0)
342  );
343 }
344 
345 
346 template<class BasicMomentumTransportModel>
348 {
349  if (!this->turbulence_)
350  {
351  return;
352  }
353 
354  // Local references
355  const alphaField& alpha = this->alpha_;
356  const rhoField& rho = this->rho_;
357  const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
358  fv::options& fvOptions(fv::options::New(this->mesh_));
359 
361 
362  const volScalarField chi(this->chi());
363  const volScalarField fv1(this->fv1(chi));
364 
365  const volScalarField::Internal Stilda(this->Stilda(chi, fv1));
366 
367  tmp<fvScalarMatrix> nuTildaEqn
368  (
369  fvm::ddt(alpha, rho, nuTilda_)
370  + fvm::div(alphaRhoPhi, nuTilda_)
371  - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
372  - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
373  ==
374  Cb1_*alpha()*rho()*Stilda*nuTilda_()
375  - fvm::Sp(Cw1_*alpha()*rho()*fw(Stilda)*nuTilda_()/sqr(y_), nuTilda_)
376  + fvOptions(alpha, rho, nuTilda_)
377  );
378 
379  nuTildaEqn.ref().relax();
380  fvOptions.constrain(nuTildaEqn.ref());
381  solve(nuTildaEqn);
382  fvOptions.correct(nuTilda_);
383  bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), 0));
384  nuTilda_.correctBoundaryConditions();
385 
386  correctNut(fv1);
387 }
388 
389 
390 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
391 
392 } // End namespace RASModels
393 } // End namespace Foam
394 
395 // ************************************************************************* //
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
fv::options & fvOptions
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)
void correct(GeometricField< Type, fvPatchField, volMesh > &field)
Apply correction to field.
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
dimensionedSymmTensor sqr(const dimensionedVector &dv)
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
Finite-volume options.
Definition: fvOptions.H:52
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
phi
Definition: pEqn.H:104
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
static const wallDist & New(const fvMesh &mesh)
Definition: MeshObject.C:44
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
A class for handling words, derived from string.
Definition: word.H:59
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
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
void constrain(fvMatrix< Type > &eqn)
Apply constraints to equation.
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.
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-8/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
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
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.
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 > &)
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
A class for managing temporary objects.
Definition: PtrList.H:53
static options & New(const fvMesh &mesh)
Construct fvOptions and register to datbase if not present.
Definition: fvOptions.C:101
const dimensionedVector & g
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Namespace for OpenFOAM.