kkLOmega.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 "kkLOmega.H"
27 #include "bound.H"
28 #include "wallDist.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 namespace incompressible
36 {
37 namespace RASModels
38 {
39 
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 
44 
45 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
46 
47 tmp<volScalarField> kkLOmega::fv(const volScalarField& Ret) const
48 {
49  return(1 - exp(-sqrt(Ret)/Av_));
50 }
51 
52 
53 tmp<volScalarField> kkLOmega::fINT() const
54 {
55  return
56  (
57  min
58  (
59  kt_/(Cint_*(kl_ + kt_ + kMin_)),
61  )
62  );
63 }
64 
65 
66 tmp<volScalarField> kkLOmega::fSS(const volScalarField& Omega) const
67 {
68  return(exp(-sqr(Css_*nu()*Omega/(kt_ + kMin_))));
69 }
70 
71 
72 tmp<volScalarField> kkLOmega::Cmu(const volScalarField& S) const
73 {
74  return(1/(A0_ + As_*(S/(omega_ + omegaMin_))));
75 }
76 
77 
78 tmp<volScalarField> kkLOmega::BetaTS(const volScalarField& ReOmega) const
79 {
80  return(scalar(1) - exp(-sqr(max(ReOmega - CtsCrit_, scalar(0)))/Ats_));
81 }
82 
83 
84 tmp<volScalarField> kkLOmega::fTaul
85 (
86  const volScalarField& lambdaEff,
87  const volScalarField& ktL,
88  const volScalarField& Omega
89 ) const
90 {
91  return
92  (
93  scalar(1)
94  - exp
95  (
96  -CtauL_*ktL
97  /(
98  sqr(lambdaEff*Omega)
100  (
101  "rootVSmall",
103  rootVSmall
104  )
105  )
106  )
107  );
108 }
109 
110 
111 tmp<volScalarField> kkLOmega::alphaT
112 (
113  const volScalarField& lambdaEff,
114  const volScalarField& fv,
115  const volScalarField& ktS
116 ) const
117 {
118  return(fv*CmuStd_*sqrt(ktS)*lambdaEff);
119 }
120 
121 
122 tmp<volScalarField> kkLOmega::fOmega
123 (
124  const volScalarField& lambdaEff,
125  const volScalarField& lambdaT
126 ) const
127 {
128  return
129  (
130  scalar(1)
131  - exp
132  (
133  -0.41
134  *pow4
135  (
136  lambdaEff
137  / (
138  lambdaT
140  (
141  "ROTvSmall",
142  lambdaT.dimensions(),
143  rootVSmall
144  )
145  )
146  )
147  )
148  );
149 }
150 
151 
152 tmp<volScalarField> kkLOmega::phiBP(const volScalarField& Omega) const
153 {
154  return
155  (
156  min
157  (
158  max
159  (
160  kt_/nu()
161  / (
162  Omega
164  (
165  "ROTvSmall",
166  Omega.dimensions(),
167  rootVSmall
168  )
169  )
170  - CbpCrit_,
171  scalar(0)
172  ),
173  scalar(50)
174  )
175  );
176 }
177 
178 
179 tmp<volScalarField> kkLOmega::phiNAT
180 (
181  const volScalarField& ReOmega,
182  const volScalarField& fNatCrit
183 ) const
184 {
185  return
186  (
187  max
188  (
189  ReOmega
190  - CnatCrit_
191  / (
192  fNatCrit + dimensionedScalar(dimless, rootVSmall)
193  ),
194  scalar(0)
195  )
196  );
197 }
198 
199 
200 tmp<volScalarField> kkLOmega::D(const volScalarField& k) const
201 {
202  return nu()*magSqr(fvc::grad(sqrt(k)));
203 }
204 
205 
206 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
207 
209 {
210  // Currently this function is not implemented due to the complexity of
211  // evaluating nut. Better calculate nut at the end of correct()
213 }
214 
215 
216 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
217 
219 (
220  const geometricOneField& alpha,
221  const geometricOneField& rho,
222  const volVectorField& U,
223  const surfaceScalarField& alphaRhoPhi,
224  const surfaceScalarField& phi,
225  const transportModel& transport,
226  const word& type
227 )
228 :
230  (
231  type,
232  alpha,
233  rho,
234  U,
235  alphaRhoPhi,
236  phi,
237  transport
238  ),
239 
240  A0_
241  (
243  (
244  "A0",
245  coeffDict_,
246  4.04
247  )
248  ),
249  As_
250  (
252  (
253  "As",
254  coeffDict_,
255  2.12
256  )
257  ),
258  Av_
259  (
261  (
262  "Av",
263  coeffDict_,
264  6.75
265  )
266  ),
267  Abp_
268  (
270  (
271  "Abp",
272  coeffDict_,
273  0.6
274  )
275  ),
276  Anat_
277  (
279  (
280  "Anat",
281  coeffDict_,
282  200
283  )
284  ),
285  Ats_
286  (
288  (
289  "Ats",
290  coeffDict_,
291  200
292  )
293  ),
294  CbpCrit_
295  (
297  (
298  "CbpCrit",
299  coeffDict_,
300  1.2
301  )
302  ),
303  Cnc_
304  (
306  (
307  "Cnc",
308  coeffDict_,
309  0.1
310  )
311  ),
312  CnatCrit_
313  (
315  (
316  "CnatCrit",
317  coeffDict_,
318  1250
319  )
320  ),
321  Cint_
322  (
324  (
325  "Cint",
326  coeffDict_,
327  0.75
328  )
329  ),
330  CtsCrit_
331  (
333  (
334  "CtsCrit",
335  coeffDict_,
336  1000
337  )
338  ),
339  CrNat_
340  (
342  (
343  "CrNat",
344  coeffDict_,
345  0.02
346  )
347  ),
348  C11_
349  (
351  (
352  "C11",
353  coeffDict_,
354  3.4e-6
355  )
356  ),
357  C12_
358  (
360  (
361  "C12",
362  coeffDict_,
363  1.0e-10
364  )
365  ),
366  CR_
367  (
369  (
370  "CR",
371  coeffDict_,
372  0.12
373  )
374  ),
376  (
378  (
379  "CalphaTheta",
380  coeffDict_,
381  0.035
382  )
383  ),
384  Css_
385  (
387  (
388  "Css",
389  coeffDict_,
390  1.5
391  )
392  ),
393  CtauL_
394  (
396  (
397  "CtauL",
398  coeffDict_,
399  4360
400  )
401  ),
402  Cw1_
403  (
405  (
406  "Cw1",
407  coeffDict_,
408  0.44
409  )
410  ),
411  Cw2_
412  (
414  (
415  "Cw2",
416  coeffDict_,
417  0.92
418  )
419  ),
420  Cw3_
421  (
423  (
424  "Cw3",
425  coeffDict_,
426  0.3
427  )
428  ),
429  CwR_
430  (
432  (
433  "CwR",
434  coeffDict_,
435  1.5
436  )
437  ),
438  Clambda_
439  (
441  (
442  "Clambda",
443  coeffDict_,
444  2.495
445  )
446  ),
447  CmuStd_
448  (
450  (
451  "CmuStd",
452  coeffDict_,
453  0.09
454  )
455  ),
456  Prtheta_
457  (
459  (
460  "Prtheta",
461  coeffDict_,
462  0.85
463  )
464  ),
465  Sigmak_
466  (
468  (
469  "Sigmak",
470  coeffDict_,
471  1
472  )
473  ),
474  Sigmaw_
475  (
477  (
478  "Sigmaw",
479  coeffDict_,
480  1.17
481  )
482  ),
483  kt_
484  (
485  IOobject
486  (
487  IOobject::groupName("kt", alphaRhoPhi.group()),
488  runTime_.timeName(),
489  mesh_,
492  ),
493  mesh_
494  ),
495  kl_
496  (
497  IOobject
498  (
499  IOobject::groupName("kl", alphaRhoPhi.group()),
500  runTime_.timeName(),
501  mesh_,
504  ),
505  mesh_
506  ),
507  omega_
508  (
509  IOobject
510  (
511  IOobject::groupName("omega", alphaRhoPhi.group()),
512  runTime_.timeName(),
513  mesh_,
516  ),
517  mesh_
518  ),
519  epsilon_
520  (
521  IOobject
522  (
523  "epsilon",
524  runTime_.timeName(),
525  mesh_
526  ),
527  kt_*omega_ + D(kl_) + D(kt_)
528  ),
529  y_(wallDist::New(mesh_).y())
530 {
531  bound(kt_, kMin_);
532  bound(kl_, kMin_);
533  bound(omega_, omegaMin_);
534  bound(epsilon_, epsilonMin_);
535 
536  if (type == typeName)
537  {
538  // Evaluating nut_ is complex so start from the field read from file
540 
541  printCoeffs(type);
542  }
543 }
544 
545 
546 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
547 
549 {
551  {
552  A0_.readIfPresent(coeffDict());
553  As_.readIfPresent(coeffDict());
554  Av_.readIfPresent(coeffDict());
555  Abp_.readIfPresent(coeffDict());
556  Anat_.readIfPresent(coeffDict());
557  Abp_.readIfPresent(coeffDict());
558  Ats_.readIfPresent(coeffDict());
559  CbpCrit_.readIfPresent(coeffDict());
560  Cnc_.readIfPresent(coeffDict());
561  CnatCrit_.readIfPresent(coeffDict());
562  Cint_.readIfPresent(coeffDict());
563  CtsCrit_.readIfPresent(coeffDict());
564  CrNat_.readIfPresent(coeffDict());
565  C11_.readIfPresent(coeffDict());
566  C12_.readIfPresent(coeffDict());
567  CR_.readIfPresent(coeffDict());
568  CalphaTheta_.readIfPresent(coeffDict());
569  Css_.readIfPresent(coeffDict());
570  CtauL_.readIfPresent(coeffDict());
571  Cw1_.readIfPresent(coeffDict());
572  Cw2_.readIfPresent(coeffDict());
573  Cw3_.readIfPresent(coeffDict());
574  CwR_.readIfPresent(coeffDict());
575  Clambda_.readIfPresent(coeffDict());
576  CmuStd_.readIfPresent(coeffDict());
577  Prtheta_.readIfPresent(coeffDict());
578  Sigmak_.readIfPresent(coeffDict());
579  Sigmaw_.readIfPresent(coeffDict());
580 
581  return true;
582  }
583  else
584  {
585  return false;
586  }
587 }
588 
589 
591 {}
592 
593 
595 {
597 
598  if (!turbulence_)
599  {
600  return;
601  }
602 
603  const volScalarField lambdaT(sqrt(kt_)/(omega_ + omegaMin_));
604 
605  const volScalarField lambdaEff(min(Clambda_*y_, lambdaT));
606 
607  const volScalarField fw
608  (
609  pow
610  (
611  lambdaEff
612  /(lambdaT + dimensionedScalar(dimLength, rootVSmall)),
613  2.0/3.0
614  )
615  );
616 
617  tmp<volTensorField> tgradU(fvc::grad(U_));
618  const volTensorField& gradU = tgradU();
619 
620  const volScalarField Omega(sqrt(2.0)*mag(skew(gradU)));
621 
622  const volScalarField S2(2*magSqr(dev(symm(gradU))));
623 
624  const volScalarField ktS(fSS(Omega)*fw*kt_);
625 
626  const volScalarField nuts
627  (
628  fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_))
629  *fINT()
630  *Cmu(sqrt(S2))*sqrt(ktS)*lambdaEff
631  );
632  const volScalarField Pkt(this->GName(), nuts*S2);
633 
634  const volScalarField ktL(kt_ - ktS);
635  const volScalarField ReOmega(sqr(y_)*Omega/nu());
636 
637  const volScalarField nutl
638  (
639  min
640  (
641  C11_*fTaul(lambdaEff, ktL, Omega)*Omega*sqr(lambdaEff)
642  *sqrt(ktL)*lambdaEff/nu()
643  + C12_*BetaTS(ReOmega)*sqr(sqr(lambdaEff/Clambda_)*Omega)/nu()
644  , 0.5*(kl_ + ktL)/(sqrt(S2) + omegaMin_)
645  )
646  );
647 
648  const volScalarField Pkl(nutl*S2);
649 
650  const volScalarField alphaTEff
651  (
652  alphaT(lambdaEff, fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_)), ktS)
653  );
654 
655  // By pass source term divided by kl_
656 
657  const dimensionedScalar fwMin("small", dimless, rootVSmall);
658 
659  const volScalarField Rbp
660  (
661  CR_*(1 - exp(-phiBP(Omega)()/Abp_))*omega_
662  /(fw + fwMin)
663  );
664 
665  const volScalarField fNatCrit(1 - exp(-Cnc_*sqrt(kl_)*y_/nu()));
666 
667  // Natural source term divided by kl_
668  const volScalarField Rnat
669  (
670  CrNat_*(1 - exp(-phiNAT(ReOmega, fNatCrit)/Anat_))*Omega
671  );
672 
673 
674  omega_.boundaryFieldRef().updateCoeffs();
675 
676  // Turbulence specific dissipation rate equation
677  tmp<fvScalarMatrix> omegaEqn
678  (
680  + fvm::div(phi_, omega_)
681  - fvm::laplacian(DomegaEff(alphaTEff), omega_)
682  ==
683  Cw1_*Pkt*omega_/(kt_ + kMin_)
684  - fvm::SuSp
685  (
686  (1 - CwR_/(fw + fwMin))*kl_*(Rbp + Rnat)/(kt_ + kMin_)
687  , omega_
688  )
689  - fvm::Sp(Cw2_*sqr(fw)*omega_, omega_)
690  + (
691  Cw3_*fOmega(lambdaEff, lambdaT)*alphaTEff*sqr(fw)*sqrt(kt_)
692  )()()/pow3(y_())
693  );
694 
695  omegaEqn.ref().relax();
696  omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
697 
698  solve(omegaEqn);
699  bound(omega_, omegaMin_);
700 
701 
702  const volScalarField Dl(D(kl_));
703 
704  // Laminar kinetic energy equation
705  tmp<fvScalarMatrix> klEqn
706  (
707  fvm::ddt(kl_)
708  + fvm::div(phi_, kl_)
709  - fvm::laplacian(nu(), kl_)
710  ==
711  Pkl
712  - fvm::Sp(Rbp + Rnat + Dl/(kl_ + kMin_), kl_)
713  );
714 
715  klEqn.ref().relax();
716  klEqn.ref().boundaryManipulate(kl_.boundaryFieldRef());
717 
718  solve(klEqn);
719  bound(kl_, kMin_);
720 
721 
722  const volScalarField Dt(D(kt_));
723 
724  // Turbulent kinetic energy equation
725  tmp<fvScalarMatrix> ktEqn
726  (
727  fvm::ddt(kt_)
728  + fvm::div(phi_, kt_)
729  - fvm::laplacian(DkEff(alphaTEff), kt_)
730  ==
731  Pkt
732  + (Rbp + Rnat)*kl_
733  - fvm::Sp(omega_ + Dt/(kt_+ kMin_), kt_)
734  );
735 
736  ktEqn.ref().relax();
737  ktEqn.ref().boundaryManipulate(kt_.boundaryFieldRef());
738 
739  solve(ktEqn);
740  bound(kt_, kMin_);
741 
742 
743  // Update total fluctuation kinetic energy dissipation rate
744  epsilon_ = kt_*omega_ + Dl + Dt;
745  bound(epsilon_, epsilonMin_);
746 
747 
748  // Re-calculate turbulent viscosity
749  nut_ = nuts + nutl;
751 }
752 
753 
754 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
755 
756 } // End namespace RASModels
757 } // End namespace incompressible
758 } // End namespace Foam
759 
760 // ************************************************************************* //
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
tmp< fvMatrix< Type > > SuSp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
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)
kkLOmega(const geometricOneField &alpha, const geometricOneField &rho, const volVectorField &U, const surfaceScalarField &alphaRhoPhi, const surfaceScalarField &phi, const transportModel &transport, const word &type=typeName)
Construct from components.
Definition: kkLOmega.C:219
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
const volScalarField & y() const
Return reference to cached distance-to-wall field.
Definition: wallDist.H:139
dimensionedScalar sqrt(const dimensionedScalar &ds)
bool readIfPresent(const dictionary &)
Update the value of dimensioned<Type> if found in the dictionary.
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
Definition: kkLOmega.C:594
const dimensionSet dimless
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
Macros for easy insertion into run-time selection tables.
const dimensionSet dimLength
virtual void correct()=0
Solve the turbulence equations and correct the turbulence viscosity.
const dimensionSet dimTime
const dimensionSet & dimensions() const
Return dimensions.
Low Reynolds-number k-kl-omega turbulence model for incompressible flows.
Definition: kkLOmega.H:113
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
dimensionedScalar exp(const dimensionedScalar &ds)
tmp< volScalarField > DkEff(const volScalarField &alphaT) const
Return the effective diffusivity for k.
Definition: kkLOmega.H:245
incompressible::RASModel ::transportModel transportModel
Definition: eddyViscosity.H:72
A class for handling words, derived from string.
Definition: word.H:59
static word groupName(Name name, const word &group)
virtual bool read()
Read RASProperties dictionary.
Definition: kkLOmega.C:548
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
phi
Definition: correctPhi.H:3
RASModel< momentumTransportModel > RASModel
virtual void validate()
Validate the turbulence fields after construction.
Definition: kkLOmega.C:590
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.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
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
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
addToRunTimeSelectionTable(RASModel, kkLOmega, dictionary)
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
dimensionedScalar pow4(const dimensionedScalar &ds)
void correctBoundaryConditions()
Correct boundary field.
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.
dimensioned< scalar > mag(const dimensioned< Type > &)
const volScalarField & y_
Wall distance.
Definition: kkLOmega.H:205
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
A class for managing temporary objects.
Definition: PtrList.H:53
tmp< volScalarField > DomegaEff(const volScalarField &alphaT) const
Return the effective diffusivity for omega.
Definition: kkLOmega.H:255
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:370
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Namespace for OpenFOAM.