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