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-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 "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()
216  notImplemented("kkLOmega::correctNut()");
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 {
599 
600  if (!turbulence_)
601  {
602  return;
603  }
604 
605  const volScalarField lambdaT(sqrt(kt_)/(omega_ + omegaMin_));
606 
607  const volScalarField lambdaEff(min(Clambda_*y_, lambdaT));
608 
609  const volScalarField fw
610  (
611  pow
612  (
613  lambdaEff
614  /(lambdaT + dimensionedScalar("SMALL", dimLength, ROOTVSMALL)),
615  2.0/3.0
616  )
617  );
618 
619  tmp<volTensorField> tgradU(fvc::grad(U_));
620  const volTensorField& gradU = tgradU();
621 
622  const volScalarField Omega(sqrt(2.0)*mag(skew(gradU)));
623 
624  const volScalarField S2(2.0*magSqr(dev(symm(gradU))));
625 
626  const volScalarField ktS(fSS(Omega)*fw*kt_);
627 
628  const volScalarField nuts
629  (
630  fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_))
631  *fINT()
632  *Cmu(sqrt(S2))*sqrt(ktS)*lambdaEff
633  );
634  const volScalarField Pkt(nuts*S2);
635 
636  const volScalarField ktL(kt_ - ktS);
637  const volScalarField ReOmega(sqr(y_)*Omega/nu());
638  const volScalarField nutl
639  (
640  min
641  (
642  C11_*fTaul(lambdaEff, ktL, Omega)*Omega*sqr(lambdaEff)
643  *sqrt(ktL)*lambdaEff/nu()
644  + C12_*BetaTS(ReOmega)*ReOmega*sqr(y_)*Omega
645  ,
646  0.5*(kl_ + ktL)/(sqrt(S2) + omegaMin_)
647  )
648  );
649 
650  const volScalarField Pkl(nutl*S2);
651 
652  const volScalarField alphaTEff
653  (
654  alphaT(lambdaEff, fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_)), ktS)
655  );
656 
657  // By pass source term divided by kl_
658 
659  const dimensionedScalar fwMin("SMALL", dimless, ROOTVSMALL);
660 
661  const volScalarField Rbp
662  (
663  CR_*(1.0 - exp(-phiBP(Omega)()/Abp_))*omega_
664  /(fw + fwMin)
665  );
666 
667  const volScalarField fNatCrit(1.0 - exp(-Cnc_*sqrt(kl_)*y_/nu()));
668 
669  // Natural source term divided by kl_
670  const volScalarField Rnat
671  (
672  CrNat_*(1.0 - exp(-phiNAT(ReOmega, fNatCrit)/Anat_))*Omega
673  );
674 
675 
676  omega_.boundaryField().updateCoeffs();
677 
678  // Turbulence specific dissipation rate equation
679  tmp<fvScalarMatrix> omegaEqn
680  (
682  + fvm::div(phi_, omega_)
683  - fvm::laplacian(DomegaEff(alphaTEff), omega_)
684  ==
685  Cw1_*Pkt*omega_/(kt_ + kMin_)
686  - fvm::SuSp
687  (
688  (1.0 - CwR_/(fw + fwMin))*kl_*(Rbp + Rnat)/(kt_ + kMin_)
689  , omega_
690  )
691  - fvm::Sp(Cw2_*sqr(fw)*omega_, omega_)
692  + (
693  Cw3_*fOmega(lambdaEff, lambdaT)*alphaTEff*sqr(fw)*sqrt(kt_)
694  )().dimensionedInternalField()/pow3(y_.dimensionedInternalField())
695  );
696 
697  omegaEqn().relax();
698  omegaEqn().boundaryManipulate(omega_.boundaryField());
699 
700  solve(omegaEqn);
701  bound(omega_, omegaMin_);
702 
703 
704  const volScalarField Dl(D(kl_));
705 
706  // Laminar kinetic energy equation
707  tmp<fvScalarMatrix> klEqn
708  (
709  fvm::ddt(kl_)
710  + fvm::div(phi_, kl_)
711  - fvm::laplacian(nu(), kl_)
712  ==
713  Pkl
714  - fvm::Sp(Rbp + Rnat + Dl/(kl_ + kMin_), kl_)
715  );
716 
717  klEqn().relax();
718  klEqn().boundaryManipulate(kl_.boundaryField());
719 
720  solve(klEqn);
721  bound(kl_, kMin_);
722 
723 
724  const volScalarField Dt(D(kt_));
725 
726  // Turbulent kinetic energy equation
727  tmp<fvScalarMatrix> ktEqn
728  (
729  fvm::ddt(kt_)
730  + fvm::div(phi_, kt_)
731  - fvm::laplacian(DkEff(alphaTEff), kt_)
732  ==
733  Pkt
734  + (Rbp + Rnat)*kl_
735  - fvm::Sp(omega_ + Dt/(kt_+ kMin_), kt_)
736  );
737 
738  ktEqn().relax();
739  ktEqn().boundaryManipulate(kt_.boundaryField());
740 
741  solve(ktEqn);
742  bound(kt_, kMin_);
743 
744 
745  // Update total fluctuation kinetic energy dissipation rate
746  epsilon_ = kt_*omega_ + Dl + Dt;
747  bound(epsilon_, epsilonMin_);
748 
749 
750  // Re-calculate turbulent viscosity
751  nut_ = nuts + nutl;
753 }
754 
755 
756 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
757 
758 } // End namespace RASModels
759 } // End namespace incompressible
760 } // End namespace Foam
761 
762 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
dimensionedTensor skew(const dimensionedTensor &dt)
dimensioned< scalar > mag(const dimensioned< Type > &)
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
Bound the given scalar field if it has gone unbounded.
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
volScalarField & bound(volScalarField &, const dimensionedScalar &lowerBound)
Bound the given scalar field if it has gone unbounded.
Definition: bound.C:33
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
virtual void correct()=0
Solve the turbulence equations and correct the turbulence viscosity.
dimensionedScalar exp(const dimensionedScalar &ds)
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
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
tmp< fvMatrix< Type > > Sp(const DimensionedField< scalar, volMesh > &, const GeometricField< Type, fvPatchField, volMesh > &)
Namespace for OpenFOAM.
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
static word groupName(Name name, const word &group)
const double e
Elementary charge.
Definition: doubleFloat.H:78
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
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 > &)
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: [].
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
bool readIfPresent(const dictionary &)
Update the value of dimensioned<Type> if found in the dictionary.
Macros for easy insertion into run-time selection tables.
const dimensionSet & dimensions() const
Return dimensions.
addToRunTimeSelectionTable(RASModel, kkLOmega, dictionary)
fileName::Type type(const fileName &)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:589
surfaceScalarField & phi
rDeltaT dimensionedInternalField()
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
Definition: kkLOmega.C:596
solverPerformance solve(fvMatrix< Type > &, const dictionary &)
Solve returning the solution statistics given convergence tolerance.
Low Reynolds-number k-kl-omega turbulence model for incompressible flows.
Definition: kkLOmega.H:110
void correctBoundaryConditions()
Correct boundary field.
incompressible::RASModel::transportModel transportModel
Definition: eddyViscosity.H:75
dimensionedScalar pow4(const dimensionedScalar &ds)
volScalarField & nu
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H: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)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
#define notImplemented(functionName)
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:356
virtual bool read()
Read RASProperties dictionary.
Definition: kkLOmega.C:554
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
tmp< volScalarField > DkEff(const volScalarField &alphaT) const
Return the effective diffusivity for k.
Definition: kkLOmega.H:243
RASModel< turbulenceModel > RASModel
U
Definition: pEqn.H:82
const volScalarField & y_
Wall distance.
Definition: kkLOmega.H:202