kOmegaSSTBase.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 "kOmegaSSTBase.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 
37 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
38 
39 template<class MomentumTransportModel, class BasicMomentumTransportModel>
40 tmp<volScalarField>
42 (
43  const volScalarField& CDkOmega
44 ) const
45 {
46  tmp<volScalarField> CDkOmegaPlus = max
47  (
48  CDkOmega,
50  );
51 
52  tmp<volScalarField> arg1 = min
53  (
54  min
55  (
56  max
57  (
58  (scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
59  scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_)
60  ),
61  (4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
62  ),
63  scalar(10)
64  );
65 
66  return tanh(pow4(arg1));
67 }
68 
69 template<class MomentumTransportModel, class BasicMomentumTransportModel>
70 tmp<volScalarField>
72 F2() const
73 {
74  tmp<volScalarField> arg2 = min
75  (
76  max
77  (
78  (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
79  scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_)
80  ),
81  scalar(100)
82  );
83 
84  return tanh(sqr(arg2));
85 }
86 
87 template<class MomentumTransportModel, class BasicMomentumTransportModel>
88 tmp<volScalarField>
90 F3() const
91 {
92  tmp<volScalarField> arg3 = min
93  (
94  150*(this->mu()/this->rho_)/(omega_*sqr(y_)),
95  scalar(10)
96  );
97 
98  return 1 - tanh(pow4(arg3));
99 }
100 
101 template<class MomentumTransportModel, class BasicMomentumTransportModel>
102 tmp<volScalarField>
103 kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::kOmegaSST::
104 F23() const
105 {
106  tmp<volScalarField> f23(F2());
107 
108  if (F3_)
109  {
110  f23.ref() *= F3();
111  }
112 
113  return f23;
114 }
115 
116 
117 template<class MomentumTransportModel, class BasicMomentumTransportModel>
119 (
120  const volScalarField& S2,
121  const volScalarField& F2
122 )
123 {
124  this->nut_ = a1_*k_/max(a1_*omega_, b1_*F2*sqrt(S2));
125  this->nut_.correctBoundaryConditions();
126  fvConstraints::New(this->mesh_).constrain(this->nut_);
127 }
128 
129 
130 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
131 
132 template<class MomentumTransportModel, class BasicMomentumTransportModel>
135 {
136  correctNut(2*magSqr(symm(fvc::grad(this->U_))), F23());
137 }
138 
139 
140 template<class MomentumTransportModel, class BasicMomentumTransportModel>
143 (
144  const volScalarField::Internal& G
145 ) const
146 {
147  return min(G, (c1_*betaStar_)*this->k_()*this->omega_());
148 }
149 
150 
151 template<class MomentumTransportModel, class BasicMomentumTransportModel>
154 (
155  const volScalarField::Internal& F1,
156  const volScalarField::Internal& F2
157 ) const
158 {
159  return betaStar_*omega_();
160 }
161 
162 
163 template<class MomentumTransportModel, class BasicMomentumTransportModel>
166 {
167  return tmp<fvScalarMatrix>
168  (
169  new fvScalarMatrix
170  (
171  k_,
172  dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
173  )
174  );
175 }
176 
177 
178 template<class MomentumTransportModel, class BasicMomentumTransportModel>
181 omegaSource() const
182 {
183  return tmp<fvScalarMatrix>
184  (
185  new fvScalarMatrix
186  (
187  omega_,
188  dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
189  )
190  );
191 }
192 
193 
194 template<class MomentumTransportModel, class BasicMomentumTransportModel>
197 (
198  const volScalarField::Internal& S2,
199  const volScalarField::Internal& gamma,
200  const volScalarField::Internal& beta
201 ) const
202 {
203  return tmp<fvScalarMatrix>
204  (
205  new fvScalarMatrix
206  (
207  omega_,
208  dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
209  )
210  );
211 }
212 
213 
214 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
215 
216 template<class MomentumTransportModel, class BasicMomentumTransportModel>
218 (
219  const word& type,
220  const alphaField& alpha,
221  const rhoField& rho,
222  const volVectorField& U,
223  const surfaceScalarField& alphaRhoPhi,
224  const surfaceScalarField& phi,
225  const transportModel& transport
226 )
227 :
229  (
230  type,
231  alpha,
232  rho,
233  U,
234  alphaRhoPhi,
235  phi,
236  transport
237  ),
238 
239  alphaK1_
240  (
242  (
243  "alphaK1",
244  this->coeffDict_,
245  0.85
246  )
247  ),
248  alphaK2_
249  (
251  (
252  "alphaK2",
253  this->coeffDict_,
254  1.0
255  )
256  ),
257  alphaOmega1_
258  (
260  (
261  "alphaOmega1",
262  this->coeffDict_,
263  0.5
264  )
265  ),
266  alphaOmega2_
267  (
269  (
270  "alphaOmega2",
271  this->coeffDict_,
272  0.856
273  )
274  ),
275  gamma1_
276  (
278  (
279  "gamma1",
280  this->coeffDict_,
281  5.0/9.0
282  )
283  ),
284  gamma2_
285  (
287  (
288  "gamma2",
289  this->coeffDict_,
290  0.44
291  )
292  ),
293  beta1_
294  (
296  (
297  "beta1",
298  this->coeffDict_,
299  0.075
300  )
301  ),
302  beta2_
303  (
305  (
306  "beta2",
307  this->coeffDict_,
308  0.0828
309  )
310  ),
311  betaStar_
312  (
314  (
315  "betaStar",
316  this->coeffDict_,
317  0.09
318  )
319  ),
320  a1_
321  (
323  (
324  "a1",
325  this->coeffDict_,
326  0.31
327  )
328  ),
329  b1_
330  (
332  (
333  "b1",
334  this->coeffDict_,
335  1.0
336  )
337  ),
338  c1_
339  (
341  (
342  "c1",
343  this->coeffDict_,
344  10.0
345  )
346  ),
347  F3_
348  (
350  (
351  "F3",
352  this->coeffDict_,
353  false
354  )
355  ),
356 
357  y_(wallDist::New(this->mesh_).y()),
358 
359  k_
360  (
361  IOobject
362  (
363  IOobject::groupName("k", alphaRhoPhi.group()),
364  this->runTime_.timeName(),
365  this->mesh_,
368  ),
369  this->mesh_
370  ),
371  omega_
372  (
373  IOobject
374  (
375  IOobject::groupName("omega", alphaRhoPhi.group()),
376  this->runTime_.timeName(),
377  this->mesh_,
380  ),
381  this->mesh_
382  )
383 {
384  bound(k_, this->kMin_);
385  bound(omega_, this->omegaMin_);
386 }
387 
388 
389 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
390 
391 template<class MomentumTransportModel, class BasicMomentumTransportModel>
393 {
395  {
396  alphaK1_.readIfPresent(this->coeffDict());
397  alphaK2_.readIfPresent(this->coeffDict());
398  alphaOmega1_.readIfPresent(this->coeffDict());
399  alphaOmega2_.readIfPresent(this->coeffDict());
400  gamma1_.readIfPresent(this->coeffDict());
401  gamma2_.readIfPresent(this->coeffDict());
402  beta1_.readIfPresent(this->coeffDict());
403  beta2_.readIfPresent(this->coeffDict());
404  betaStar_.readIfPresent(this->coeffDict());
405  a1_.readIfPresent(this->coeffDict());
406  b1_.readIfPresent(this->coeffDict());
407  c1_.readIfPresent(this->coeffDict());
408  F3_.readIfPresent("F3", this->coeffDict());
409 
410  return true;
411  }
412  else
413  {
414  return false;
415  }
416 }
417 
418 
419 template<class MomentumTransportModel, class BasicMomentumTransportModel>
421 {
422  if (!this->turbulence_)
423  {
424  return;
425  }
426 
427  // Local references
428  const alphaField& alpha = this->alpha_;
429  const rhoField& rho = this->rho_;
430  const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
431  const volVectorField& U = this->U_;
432  volScalarField& nut = this->nut_;
433  const Foam::fvModels& fvModels(Foam::fvModels::New(this->mesh_));
435  (
436  Foam::fvConstraints::New(this->mesh_)
437  );
438 
440 
442  (
443  fvc::div(fvc::absolute(this->phi(), U))()()
444  );
445 
446  tmp<volTensorField> tgradU = fvc::grad(U);
447  volScalarField S2(2*magSqr(symm(tgradU())));
448  volScalarField::Internal GbyNu(dev(twoSymm(tgradU()())) && tgradU()());
449  volScalarField::Internal G(this->GName(), nut()*GbyNu);
450  tgradU.clear();
451 
452  // Update omega and G at the wall
453  omega_.boundaryFieldRef().updateCoeffs();
454 
455  volScalarField CDkOmega
456  (
457  (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
458  );
459 
460  volScalarField F1(this->F1(CDkOmega));
461  volScalarField F23(this->F23());
462 
463  {
464  volScalarField::Internal gamma(this->gamma(F1));
465  volScalarField::Internal beta(this->beta(F1));
466 
467  // Turbulent frequency equation
468  tmp<fvScalarMatrix> omegaEqn
469  (
470  fvm::ddt(alpha, rho, omega_)
471  + fvm::div(alphaRhoPhi, omega_)
472  - fvm::laplacian(alpha*rho*DomegaEff(F1), omega_)
473  ==
474  alpha()*rho()*gamma
475  *min
476  (
477  GbyNu,
478  (c1_/a1_)*betaStar_*omega_()
479  *max(a1_*omega_(), b1_*F23()*sqrt(S2()))
480  )
481  - fvm::SuSp((2.0/3.0)*alpha()*rho()*gamma*divU, omega_)
482  - fvm::Sp(alpha()*rho()*beta*omega_(), omega_)
483  - fvm::SuSp
484  (
485  alpha()*rho()*(F1() - scalar(1))*CDkOmega()/omega_(),
486  omega_
487  )
488  + Qsas(S2(), gamma, beta)
489  + omegaSource()
490  + fvModels.source(alpha, rho, omega_)
491  );
492 
493  omegaEqn.ref().relax();
494  fvConstraints.constrain(omegaEqn.ref());
495  omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
496  solve(omegaEqn);
497  fvConstraints.constrain(omega_);
498  bound(omega_, this->omegaMin_);
499  }
500 
501  // Turbulent kinetic energy equation
503  (
504  fvm::ddt(alpha, rho, k_)
505  + fvm::div(alphaRhoPhi, k_)
506  - fvm::laplacian(alpha*rho*DkEff(F1), k_)
507  ==
508  alpha()*rho()*Pk(G)
509  - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
510  - fvm::Sp(alpha()*rho()*epsilonByk(F1, F23), k_)
511  + kSource()
512  + fvModels.source(alpha, rho, k_)
513  );
514 
515  kEqn.ref().relax();
516  fvConstraints.constrain(kEqn.ref());
517  solve(kEqn);
518  fvConstraints.constrain(k_);
519  bound(k_, this->kMin_);
520 
521  correctNut(S2, F23);
522 }
523 
524 
525 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
526 
527 } // End namespace Foam
528 
529 // ************************************************************************* //
static word group(const word &name)
Return group (extension part of name)
Definition: IOobject.C:134
dimensionedScalar tanh(const dimensionedScalar &ds)
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
tmp< fvMatrix< Type > > SuSp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:237
kOmegaSST(const word &type, const alphaField &alpha, const rhoField &rho, const volVectorField &U, const surfaceScalarField &alphaRhoPhi, const surfaceScalarField &phi, const transportModel &transport)
Construct from components.
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define F1(B, C, D)
Definition: SHA1.C:173
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
dimensionedSymmTensor sqr(const dimensionedVector &dv)
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
dimensionedScalar sqrt(const dimensionedScalar &ds)
#define F2(B, C, D)
Definition: SHA1.C:174
virtual tmp< fvScalarMatrix > omegaSource() const
const dimensionSet dimless
Generic dimensioned Type class.
const dimensionedScalar G
Newtonian constant of gravitation.
void read(Istream &, label &, const dictionary &)
In-place read with dictionary lookup.
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
BasicMomentumTransportModel::alphaField alphaField
Definition: RASModel.H:88
static Switch lookupOrAddToDict(const word &, dictionary &, const Switch &defaultValue=false)
Construct from dictionary, supplying default value so that if the.
Definition: Switch.C:111
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
const dimensionSet dimTime
scalar y
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
BasicMomentumTransportModel::transportModel transportModel
Definition: RASModel.H:90
Foam::fvConstraints & fvConstraints
#define F3(B, C, D)
Definition: SHA1.C:175
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
static word groupName(Name name, const word &group)
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
dimensioned< scalar > magSqr(const dimensioned< Type > &)
virtual tmp< fvScalarMatrix > Qsas(const volScalarField::Internal &S2, const volScalarField::Internal &gamma, const volScalarField::Internal &beta) const
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.
const dimensionedScalar mu
Atomic mass unit.
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 > &)
virtual tmp< volScalarField::Internal > Pk(const volScalarField::Internal &G) const
Return k production rate.
tmp< volScalarField > divU
Definition: alphaSuSp.H:9
bool constrain(fvMatrix< Type > &eqn) const
Apply constraints to an equation.
Foam::fvModels & fvModels
BasicMomentumTransportModel::rhoField rhoField
Definition: RASModel.H:89
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
Definition: fvcMeshPhi.C:188
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
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
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
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
dimensionedScalar pow4(const dimensionedScalar &ds)
virtual tmp< fvScalarMatrix > kSource() const
const dimensionSet dimVolume
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
virtual tmp< volScalarField::Internal > epsilonByk(const volScalarField::Internal &F1, const volScalarField::Internal &F2) const
Return epsilon/k which for standard RAS is betaStar*omega.
A class for managing temporary objects.
Definition: PtrList.H:53
Templated abstract base class for turbulence models.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
virtual bool read()
Re-read model coefficients if they have changed.
scalar nut
Namespace for OpenFOAM.
virtual void correctNut()