kOmegaSST.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) 2011-2015 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 "kOmegaSST.H"
27 #include "bound.H"
28 #include "wallDist.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 namespace RASModels
35 {
36 
37 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
38 
39 template<class BasicTurbulenceModel>
41 (
42  const volScalarField& CDkOmega
43 ) const
44 {
45  tmp<volScalarField> CDkOmegaPlus = max
46  (
47  CDkOmega,
48  dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
49  );
50 
51  tmp<volScalarField> arg1 = min
52  (
53  min
54  (
55  max
56  (
57  (scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
58  scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_)
59  ),
60  (4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
61  ),
62  scalar(10)
63  );
64 
65  return tanh(pow4(arg1));
66 }
67 
68 template<class BasicTurbulenceModel>
69 tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F2() const
70 {
71  tmp<volScalarField> arg2 = min
72  (
73  max
74  (
75  (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
76  scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_)
77  ),
78  scalar(100)
79  );
80 
81  return tanh(sqr(arg2));
82 }
83 
84 template<class BasicTurbulenceModel>
85 tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F3() const
86 {
87  tmp<volScalarField> arg3 = min
88  (
89  150*(this->mu()/this->rho_)/(omega_*sqr(y_)),
90  scalar(10)
91  );
92 
93  return 1 - tanh(pow4(arg3));
94 }
95 
96 template<class BasicTurbulenceModel>
97 tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F23() const
98 {
99  tmp<volScalarField> f23(F2());
100 
101  if (F3_)
102  {
103  f23() *= F3();
104  }
105 
106  return f23;
107 }
108 
109 
110 template<class BasicTurbulenceModel>
112 {
113  this->nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
114  this->nut_.correctBoundaryConditions();
115 
116  BasicTurbulenceModel::correctNut();
117 }
118 
119 
120 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
121 
122 template<class BasicTurbulenceModel>
124 {
125  correctNut(2*magSqr(symm(fvc::grad(this->U_))));
126 }
127 
128 
129 template<class BasicTurbulenceModel>
131 {
132  return tmp<fvScalarMatrix>
133  (
134  new fvScalarMatrix
135  (
136  k_,
137  dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
138  )
139  );
140 }
141 
142 
143 template<class BasicTurbulenceModel>
145 {
146  return tmp<fvScalarMatrix>
147  (
148  new fvScalarMatrix
149  (
150  omega_,
151  dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
152  )
153  );
154 }
155 
156 
157 template<class BasicTurbulenceModel>
159 (
160  const volScalarField& S2,
161  const volScalarField& gamma,
162  const volScalarField& beta
163 ) const
164 {
165  return tmp<fvScalarMatrix>
166  (
167  new fvScalarMatrix
168  (
169  omega_,
170  dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
171  )
172  );
173 }
174 
175 
176 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
177 
178 template<class BasicTurbulenceModel>
180 (
181  const alphaField& alpha,
182  const rhoField& rho,
183  const volVectorField& U,
184  const surfaceScalarField& alphaRhoPhi,
185  const surfaceScalarField& phi,
186  const transportModel& transport,
187  const word& propertiesName,
188  const word& type
189 )
190 :
192  (
193  type,
194  alpha,
195  rho,
196  U,
197  alphaRhoPhi,
198  phi,
199  transport,
200  propertiesName
201  ),
202 
203  alphaK1_
204  (
206  (
207  "alphaK1",
208  this->coeffDict_,
209  0.85
210  )
211  ),
212  alphaK2_
213  (
215  (
216  "alphaK2",
217  this->coeffDict_,
218  1.0
219  )
220  ),
221  alphaOmega1_
222  (
224  (
225  "alphaOmega1",
226  this->coeffDict_,
227  0.5
228  )
229  ),
230  alphaOmega2_
231  (
233  (
234  "alphaOmega2",
235  this->coeffDict_,
236  0.856
237  )
238  ),
239  gamma1_
240  (
242  (
243  "gamma1",
244  this->coeffDict_,
245  5.0/9.0
246  )
247  ),
248  gamma2_
249  (
251  (
252  "gamma2",
253  this->coeffDict_,
254  0.44
255  )
256  ),
257  beta1_
258  (
260  (
261  "beta1",
262  this->coeffDict_,
263  0.075
264  )
265  ),
266  beta2_
267  (
269  (
270  "beta2",
271  this->coeffDict_,
272  0.0828
273  )
274  ),
275  betaStar_
276  (
278  (
279  "betaStar",
280  this->coeffDict_,
281  0.09
282  )
283  ),
284  a1_
285  (
287  (
288  "a1",
289  this->coeffDict_,
290  0.31
291  )
292  ),
293  b1_
294  (
296  (
297  "b1",
298  this->coeffDict_,
299  1.0
300  )
301  ),
302  c1_
303  (
305  (
306  "c1",
307  this->coeffDict_,
308  10.0
309  )
310  ),
311  F3_
312  (
314  (
315  "F3",
316  this->coeffDict_,
317  false
318  )
319  ),
320 
321  y_(wallDist::New(this->mesh_).y()),
322 
323  k_
324  (
325  IOobject
326  (
327  IOobject::groupName("k", U.group()),
328  this->runTime_.timeName(),
329  this->mesh_,
332  ),
333  this->mesh_
334  ),
335  omega_
336  (
337  IOobject
338  (
339  IOobject::groupName("omega", U.group()),
340  this->runTime_.timeName(),
341  this->mesh_,
344  ),
345  this->mesh_
346  )
347 {
348  bound(k_, this->kMin_);
349  bound(omega_, this->omegaMin_);
350 
351  if (type == typeName)
352  {
353  correctNut();
354  this->printCoeffs(type);
355  }
356 }
357 
358 
359 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
360 
361 template<class BasicTurbulenceModel>
363 {
365  {
366  alphaK1_.readIfPresent(this->coeffDict());
367  alphaK2_.readIfPresent(this->coeffDict());
368  alphaOmega1_.readIfPresent(this->coeffDict());
369  alphaOmega2_.readIfPresent(this->coeffDict());
370  gamma1_.readIfPresent(this->coeffDict());
371  gamma2_.readIfPresent(this->coeffDict());
372  beta1_.readIfPresent(this->coeffDict());
373  beta2_.readIfPresent(this->coeffDict());
374  betaStar_.readIfPresent(this->coeffDict());
375  a1_.readIfPresent(this->coeffDict());
376  b1_.readIfPresent(this->coeffDict());
377  c1_.readIfPresent(this->coeffDict());
378  F3_.readIfPresent("F3", this->coeffDict());
379 
380  return true;
381  }
382  else
383  {
384  return false;
385  }
386 }
387 
388 
389 template<class BasicTurbulenceModel>
391 {
392  if (!this->turbulence_)
393  {
394  return;
395  }
396 
397  // Local references
398  const alphaField& alpha = this->alpha_;
399  const rhoField& rho = this->rho_;
400  const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
401  const volVectorField& U = this->U_;
402  volScalarField& nut = this->nut_;
403 
405 
407 
408  tmp<volTensorField> tgradU = fvc::grad(U);
409  volScalarField S2(2*magSqr(symm(tgradU())));
410  volScalarField GbyNu((tgradU() && dev(twoSymm(tgradU()))));
411  volScalarField G(this->GName(), nut*GbyNu);
412  tgradU.clear();
413 
414  // Update omega and G at the wall
415  omega_.boundaryField().updateCoeffs();
416 
417  volScalarField CDkOmega
418  (
419  (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
420  );
421 
422  volScalarField F1(this->F1(CDkOmega));
423 
424  {
425  volScalarField gamma(this->gamma(F1));
426  volScalarField beta(this->beta(F1));
427 
428  // Turbulent frequency equation
429  tmp<fvScalarMatrix> omegaEqn
430  (
431  fvm::ddt(alpha, rho, omega_)
432  + fvm::div(alphaRhoPhi, omega_)
433  - fvm::laplacian(alpha*rho*DomegaEff(F1), omega_)
434  ==
435  alpha*rho*gamma
436  *min
437  (
438  GbyNu,
439  (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2))
440  )
441  - fvm::SuSp((2.0/3.0)*alpha*rho*gamma*divU, omega_)
442  - fvm::Sp(alpha*rho*beta*omega_, omega_)
443  - fvm::SuSp
444  (
445  alpha*rho*(F1 - scalar(1))*CDkOmega/omega_,
446  omega_
447  )
448  + Qsas(S2, gamma, beta)
449  + omegaSource()
450  );
451 
452  omegaEqn().relax();
453 
454  omegaEqn().boundaryManipulate(omega_.boundaryField());
455 
456  solve(omegaEqn);
457  bound(omega_, this->omegaMin_);
458  }
459 
460  // Turbulent kinetic energy equation
462  (
463  fvm::ddt(alpha, rho, k_)
464  + fvm::div(alphaRhoPhi, k_)
465  - fvm::laplacian(alpha*rho*DkEff(F1), k_)
466  ==
467  min(alpha*rho*G, (c1_*betaStar_)*alpha*rho*k_*omega_)
468  - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
469  - fvm::Sp(alpha*rho*betaStar_*omega_, k_)
470  + kSource()
471  );
472 
473  kEqn().relax();
474  solve(kEqn);
475  bound(k_, this->kMin_);
476 
477  correctNut(S2);
478 }
479 
480 
481 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
482 
483 } // End namespace RASModels
484 } // End namespace Foam
485 
486 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
virtual tmp< fvScalarMatrix > omegaSource() const
Definition: kOmegaSST.C:144
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
Definition: kOmegaSST.C:390
const dimensionedScalar mu
Atomic mass unit.
Definition: createFields.H:13
#define F2(B, C, D)
Definition: SHA1.C:176
volScalarField divU(fvc::div(fvc::absolute(phi, U)))
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
virtual void correctNut()
Definition: kOmegaSST.C:123
Implementation of the k-omega-SST turbulence model for incompressible and compressible flows...
Definition: kOmegaSST.H:120
Bound the given scalar field if it has gone unbounded.
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:174
dimensioned< scalar > magSqr(const dimensioned< Type > &)
static const wallDist & New(const fvMesh &mesh)
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
BasicTurbulenceModel::rhoField rhoField
Definition: kOmegaSST.H:225
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(psi *p+alphal *rhol0+((alphav *psiv+alphal *psil)-psi)*pSat, rhoMin);#1"/home/dm2/henry/OpenFOAM/OpenFOAM-3.0.x/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:74
volScalarField & bound(volScalarField &, const dimensionedScalar &lowerBound)
Bound the given scalar field if it has gone unbounded.
Definition: bound.C:33
static Switch lookupOrAddToDict(const word &, dictionary &, const Switch &defaultValue=false)
Construct from dictionary, supplying default value so that if the.
Definition: Switch.C:107
A class for handling words, derived from string.
Definition: word.H:59
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:46
BasicTurbulenceModel::transportModel transportModel
Definition: kOmegaSST.H:226
Templated abstract base class for RAS turbulence models.
Definition: RASModel.H:49
tmp< fvMatrix< Type > > Sp(const DimensionedField< scalar, volMesh > &, const GeometricField< Type, fvPatchField, volMesh > &)
Namespace for OpenFOAM.
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
#define F3(B, C, D)
Definition: SHA1.C:177
static word groupName(Name name, const word &group)
virtual tmp< fvScalarMatrix > Qsas(const volScalarField &S2, const volScalarField &gamma, const volScalarField &beta) const
Definition: kOmegaSST.C:159
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
Generic dimensioned Type class.
const double e
Elementary charge.
Definition: doubleFloat.H:78
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
Definition: fvcMeshPhi.C:187
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
const volScalarField & y() const
Return reference to cached distance-to-wall field.
Definition: wallDist.H:142
tmp< fvMatrix< Type > > SuSp(const DimensionedField< scalar, volMesh > &, const GeometricField< Type, fvPatchField, volMesh > &)
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
scalar nut
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Eddy viscosity turbulence model base class.
Definition: eddyViscosity.H:52
dimensionedScalar tanh(const dimensionedScalar &ds)
#define F1(B, C, D)
Definition: SHA1.C:175
fileName::Type type(const fileName &)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:589
surfaceScalarField & phi
virtual tmp< fvScalarMatrix > kSource() const
Definition: kOmegaSST.C:130
solverPerformance solve(fvMatrix< Type > &, const dictionary &)
Solve returning the solution statistics given convergence tolerance.
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
dimensionedScalar pow4(const dimensionedScalar &ds)
bool read(const char *, int32_t &)
Definition: int32IO.C:87
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
BasicTurbulenceModel::alphaField alphaField
Definition: kOmegaSST.H:224
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
const dimensionedScalar G
Newtonian constant of gravitation.
U
Definition: pEqn.H:82
virtual bool read()
Re-read model coefficients if they have changed.
Definition: kOmegaSST.C:362