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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2019-2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "kkLOmega.H"
30#include "bound.H"
31#include "wallDist.H"
33
34// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35
36namespace Foam
37{
38namespace incompressible
39{
40namespace RASModels
41{
42
43// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
44
47
48// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
49
50tmp<volScalarField> kkLOmega::fv(const volScalarField& Ret) const
51{
52 return(1.0 - exp(-sqrt(Ret)/Av_));
53}
54
55
56tmp<volScalarField> kkLOmega::fINT() const
57{
58 return
59 (
60 min
61 (
62 kt_/(Cint_*(kl_ + kt_ + kMin_)),
63 dimensionedScalar("1.0", dimless, 1.0)
64 )
65 );
66}
67
68
69tmp<volScalarField> kkLOmega::fSS(const volScalarField& Omega) const
70{
71 return(exp(-sqr(Css_*nu()*Omega/(kt_ + kMin_))));
72}
73
74
75tmp<volScalarField> kkLOmega::Cmu(const volScalarField& S) const
76{
77 return(1.0/(A0_ + As_*(S/(omega_ + omegaMin_))));
78}
79
80
81tmp<volScalarField> kkLOmega::BetaTS(const volScalarField& ReOmega) const
82{
83 return(scalar(1) - exp(-sqr(max(ReOmega - CtsCrit_, scalar(0)))/Ats_));
84}
85
86
87tmp<volScalarField> kkLOmega::fTaul
88(
89 const volScalarField& lambdaEff,
90 const volScalarField& ktL,
91 const volScalarField& Omega
92) const
93{
94 return
95 (
96 scalar(1)
97 - exp
98 (
99 -CtauL_*ktL
100 /
101 (
102 sqr
103 (
104 lambdaEff*Omega
106 (
107 "ROOTVSMALL",
109 ROOTVSMALL
110 )
111 )
112 )
113 )
114 );
115}
116
117
118tmp<volScalarField> kkLOmega::alphaT
119(
120 const volScalarField& lambdaEff,
121 const volScalarField& fv,
122 const volScalarField& ktS
123) const
124{
125 return(fv*CmuStd_*sqrt(ktS)*lambdaEff);
126}
127
128
129tmp<volScalarField> kkLOmega::fOmega
130(
131 const volScalarField& lambdaEff,
132 const volScalarField& lambdaT
133) const
134{
135 return
136 (
137 scalar(1)
138 - exp
139 (
140 -0.41
141 *pow4
142 (
143 lambdaEff
144 / (
145 lambdaT
147 (
148 "ROTVSMALL",
149 lambdaT.dimensions(),
150 ROOTVSMALL
151 )
152 )
153 )
154 )
155 );
156}
157
158
159tmp<volScalarField> kkLOmega::phiBP(const volScalarField& Omega) const
160{
161 return
162 (
163 min
164 (
165 max
166 (
167 kt_/nu()
168 / (
169 Omega
171 (
172 "ROTVSMALL",
173 Omega.dimensions(),
174 ROOTVSMALL
175 )
176 )
177 - CbpCrit_,
178 scalar(0)
179 ),
180 scalar(50)
181 )
182 );
183}
184
185
186tmp<volScalarField> kkLOmega::phiNAT
187(
188 const volScalarField& ReOmega,
189 const volScalarField& fNatCrit
190) const
191{
192 return
193 (
194 max
195 (
196 ReOmega
197 - CnatCrit_
198 / (
199 fNatCrit + dimensionedScalar("ROTVSMALL", dimless, ROOTVSMALL)
200 ),
201 scalar(0)
202 )
203 );
204}
205
206
207tmp<volScalarField> kkLOmega::D(const volScalarField& k) const
208{
209 return nu()*magSqr(fvc::grad(sqrt(k)));
210}
211
212
213// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
214
216{
217 // Currently this function is not implemented due to the complexity of
218 // evaluating nut. Better calculate nut at the end of correct()
220}
221
222
223// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
224
226(
228 const geometricOneField& rho,
229 const volVectorField& U,
230 const surfaceScalarField& alphaRhoPhi,
231 const surfaceScalarField& phi,
232 const transportModel& transport,
233 const word& propertiesName,
234 const word& type
235)
236:
237 eddyViscosity<incompressible::RASModel>
238 (
239 type,
240 alpha,
241 rho,
242 U,
243 alphaRhoPhi,
244 phi,
245 transport,
246 propertiesName
247 ),
248
249 A0_
250 (
251 dimensioned<scalar>::getOrAddToDict
252 (
253 "A0",
254 coeffDict_,
255 4.04
256 )
257 ),
258 As_
259 (
260 dimensioned<scalar>::getOrAddToDict
261 (
262 "As",
263 coeffDict_,
264 2.12
265 )
266 ),
267 Av_
268 (
269 dimensioned<scalar>::getOrAddToDict
270 (
271 "Av",
272 coeffDict_,
273 6.75
274 )
275 ),
276 Abp_
277 (
278 dimensioned<scalar>::getOrAddToDict
279 (
280 "Abp",
281 coeffDict_,
282 0.6
283 )
284 ),
285 Anat_
286 (
287 dimensioned<scalar>::getOrAddToDict
288 (
289 "Anat",
290 coeffDict_,
291 200
292 )
293 ),
294 Ats_
295 (
296 dimensioned<scalar>::getOrAddToDict
297 (
298 "Ats",
299 coeffDict_,
300 200
301 )
302 ),
303 CbpCrit_
304 (
305 dimensioned<scalar>::getOrAddToDict
306 (
307 "CbpCrit",
308 coeffDict_,
309 1.2
310 )
311 ),
312 Cnc_
313 (
314 dimensioned<scalar>::getOrAddToDict
315 (
316 "Cnc",
317 coeffDict_,
318 0.1
319 )
320 ),
321 CnatCrit_
322 (
323 dimensioned<scalar>::getOrAddToDict
324 (
325 "CnatCrit",
326 coeffDict_,
327 1250
328 )
329 ),
330 Cint_
331 (
332 dimensioned<scalar>::getOrAddToDict
333 (
334 "Cint",
335 coeffDict_,
336 0.75
337 )
338 ),
339 CtsCrit_
340 (
341 dimensioned<scalar>::getOrAddToDict
342 (
343 "CtsCrit",
344 coeffDict_,
345 1000
346 )
347 ),
348 CrNat_
349 (
350 dimensioned<scalar>::getOrAddToDict
351 (
352 "CrNat",
353 coeffDict_,
354 0.02
355 )
356 ),
357 C11_
358 (
359 dimensioned<scalar>::getOrAddToDict
360 (
361 "C11",
362 coeffDict_,
363 3.4e-6
364 )
365 ),
366 C12_
367 (
368 dimensioned<scalar>::getOrAddToDict
369 (
370 "C12",
371 coeffDict_,
372 1.0e-10
373 )
374 ),
375 CR_
376 (
377 dimensioned<scalar>::getOrAddToDict
378 (
379 "CR",
380 coeffDict_,
381 0.12
382 )
383 ),
384 CalphaTheta_
385 (
386 dimensioned<scalar>::getOrAddToDict
387 (
388 "CalphaTheta",
389 coeffDict_,
390 0.035
391 )
392 ),
393 Css_
394 (
395 dimensioned<scalar>::getOrAddToDict
396 (
397 "Css",
398 coeffDict_,
399 1.5
400 )
401 ),
402 CtauL_
403 (
404 dimensioned<scalar>::getOrAddToDict
405 (
406 "CtauL",
407 coeffDict_,
408 4360
409 )
410 ),
411 Cw1_
412 (
413 dimensioned<scalar>::getOrAddToDict
414 (
415 "Cw1",
416 coeffDict_,
417 0.44
418 )
419 ),
420 Cw2_
421 (
422 dimensioned<scalar>::getOrAddToDict
423 (
424 "Cw2",
425 coeffDict_,
426 0.92
427 )
428 ),
429 Cw3_
430 (
431 dimensioned<scalar>::getOrAddToDict
432 (
433 "Cw3",
434 coeffDict_,
435 0.3
436 )
437 ),
438 CwR_
439 (
440 dimensioned<scalar>::getOrAddToDict
441 (
442 "CwR",
443 coeffDict_,
444 1.5
445 )
446 ),
447 Clambda_
448 (
449 dimensioned<scalar>::getOrAddToDict
450 (
451 "Clambda",
452 coeffDict_,
453 2.495
454 )
455 ),
456 CmuStd_
457 (
458 dimensioned<scalar>::getOrAddToDict
459 (
460 "CmuStd",
461 coeffDict_,
462 0.09
463 )
464 ),
465 Prtheta_
466 (
467 dimensioned<scalar>::getOrAddToDict
468 (
469 "Prtheta",
470 coeffDict_,
471 0.85
472 )
473 ),
474 Sigmak_
475 (
476 dimensioned<scalar>::getOrAddToDict
477 (
478 "Sigmak",
479 coeffDict_,
480 1
481 )
482 ),
483 Sigmaw_
484 (
485 dimensioned<scalar>::getOrAddToDict
486 (
487 "Sigmaw",
488 coeffDict_,
489 1.17
490 )
491 ),
492 kt_
493 (
495 (
496 IOobject::groupName("kt", alphaRhoPhi.group()),
497 runTime_.timeName(),
498 mesh_,
499 IOobject::MUST_READ,
500 IOobject::AUTO_WRITE
501 ),
502 mesh_
503 ),
504 kl_
505 (
507 (
508 IOobject::groupName("kl", alphaRhoPhi.group()),
509 runTime_.timeName(),
510 mesh_,
511 IOobject::MUST_READ,
512 IOobject::AUTO_WRITE
513 ),
514 mesh_
515 ),
516 omega_
517 (
519 (
520 IOobject::groupName("omega", alphaRhoPhi.group()),
521 runTime_.timeName(),
522 mesh_,
523 IOobject::MUST_READ,
524 IOobject::AUTO_WRITE
525 ),
526 mesh_
527 ),
528 epsilon_
529 (
531 (
532 "epsilon",
533 runTime_.timeName(),
534 mesh_
535 ),
536 kt_*omega_ + D(kl_) + D(kt_)
537 ),
538 y_(wallDist::New(mesh_).y())
539{
540 bound(kt_, kMin_);
541 bound(kl_, kMin_);
542 bound(omega_, omegaMin_);
543 bound(epsilon_, epsilonMin_);
544
545 if (type == typeName)
546 {
547 // Evaluating nut_ is complex so start from the field read from file
549
550 printCoeffs(type);
551 }
552}
553
554
555// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
556
558{
560 {
561 A0_.readIfPresent(coeffDict());
562 As_.readIfPresent(coeffDict());
563 Av_.readIfPresent(coeffDict());
564 Abp_.readIfPresent(coeffDict());
565 Anat_.readIfPresent(coeffDict());
566 Abp_.readIfPresent(coeffDict());
567 Ats_.readIfPresent(coeffDict());
568 CbpCrit_.readIfPresent(coeffDict());
569 Cnc_.readIfPresent(coeffDict());
570 CnatCrit_.readIfPresent(coeffDict());
571 Cint_.readIfPresent(coeffDict());
572 CtsCrit_.readIfPresent(coeffDict());
573 CrNat_.readIfPresent(coeffDict());
574 C11_.readIfPresent(coeffDict());
575 C12_.readIfPresent(coeffDict());
576 CR_.readIfPresent(coeffDict());
577 CalphaTheta_.readIfPresent(coeffDict());
578 Css_.readIfPresent(coeffDict());
579 CtauL_.readIfPresent(coeffDict());
580 Cw1_.readIfPresent(coeffDict());
581 Cw2_.readIfPresent(coeffDict());
582 Cw3_.readIfPresent(coeffDict());
583 CwR_.readIfPresent(coeffDict());
584 Clambda_.readIfPresent(coeffDict());
585 CmuStd_.readIfPresent(coeffDict());
586 Prtheta_.readIfPresent(coeffDict());
587 Sigmak_.readIfPresent(coeffDict());
588 Sigmaw_.readIfPresent(coeffDict());
589
590 return true;
591 }
592
593 return false;
594}
595
596
598{}
599
600
602{
604
605 if (!turbulence_)
606 {
607 return;
608 }
609
610 const volScalarField lambdaT(sqrt(kt_)/(omega_ + omegaMin_));
611
612 const volScalarField lambdaEff(min(Clambda_*y_, lambdaT));
613
614 const volScalarField fw
615 (
616 pow
617 (
618 lambdaEff
619 /(lambdaT + dimensionedScalar("SMALL", dimLength, ROOTVSMALL)),
620 2.0/3.0
621 )
622 );
623
624 tmp<volTensorField> tgradU(fvc::grad(U_));
625 const volTensorField& gradU = tgradU();
626
627 const volScalarField Omega(sqrt(2.0)*mag(skew(gradU)));
628
629 const volScalarField S2(2.0*magSqr(dev(symm(gradU))));
630
631 const volScalarField ktS(fSS(Omega)*fw*kt_);
632
633 const volScalarField nuts
634 (
635 fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_))
636 *fINT()
637 *Cmu(sqrt(S2))*sqrt(ktS)*lambdaEff
638 );
639 volScalarField Pkt(this->GName(), nuts*S2);
640
641 const volScalarField ktL(kt_ - ktS);
642 const volScalarField ReOmega(sqr(y_)*Omega/nu());
643 const volScalarField nutl
644 (
645 min
646 (
647 C11_*fTaul(lambdaEff, ktL, Omega)*Omega*sqr(lambdaEff)
648 *sqrt(ktL)*lambdaEff/nu()
649 + C12_*BetaTS(ReOmega)*ReOmega*sqr(y_)*Omega
650 ,
651 0.5*(kl_ + ktL)/(sqrt(S2) + omegaMin_)
652 )
653 );
654
655 const volScalarField Pkl(nutl*S2);
656
657 const volScalarField alphaTEff
658 (
659 alphaT(lambdaEff, fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_)), ktS)
660 );
661
662 // By pass source term divided by kl_
663
664 const dimensionedScalar fwMin("SMALL", dimless, ROOTVSMALL);
665
666 const volScalarField Rbp
667 (
668 CR_*(1.0 - exp(-phiBP(Omega)()/Abp_))*omega_
669 /(fw + fwMin)
670 );
671
672 const volScalarField fNatCrit(1.0 - exp(-Cnc_*sqrt(kl_)*y_/nu()));
673
674 // Natural source term divided by kl_
675 const volScalarField Rnat
676 (
677 CrNat_*(1.0 - exp(-phiNAT(ReOmega, fNatCrit)/Anat_))*Omega
678 );
679
680
682
683 // Turbulence specific dissipation rate equation
684 tmp<fvScalarMatrix> omegaEqn
685 (
687 + fvm::div(phi_, omega_)
688 - fvm::laplacian(DomegaEff(alphaTEff), omega_)
689 ==
690 Cw1_*Pkt*omega_/(kt_ + kMin_)
691 - fvm::SuSp
692 (
693 (1.0 - CwR_/(fw + fwMin))*kl_*(Rbp + Rnat)/(kt_ + kMin_)
694 , omega_
695 )
697 + (
698 Cw3_*fOmega(lambdaEff, lambdaT)*alphaTEff*sqr(fw)*sqrt(kt_)
699 )()()/pow3(y_())
700 );
701
702 omegaEqn.ref().relax();
703 omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
704
705 solve(omegaEqn);
706 bound(omega_, omegaMin_);
707
708
709 const volScalarField Dl(D(kl_));
710
711 // Laminar kinetic energy equation
713 (
715 + fvm::div(phi_, kl_)
716 - fvm::laplacian(nu(), kl_)
717 ==
718 Pkl
719 - fvm::Sp(Rbp + Rnat + Dl/(kl_ + kMin_), kl_)
720 );
721
722 klEqn.ref().relax();
723 klEqn.ref().boundaryManipulate(kl_.boundaryFieldRef());
724
725 solve(klEqn);
726 bound(kl_, kMin_);
727
728
729 const volScalarField Dt(D(kt_));
730
731 // Turbulent kinetic energy equation
733 (
735 + fvm::div(phi_, kt_)
736 - fvm::laplacian(DkEff(alphaTEff), kt_)
737 ==
738 Pkt
739 + (Rbp + Rnat)*kl_
740 - fvm::Sp(omega_ + Dt/(kt_+ kMin_), kt_)
741 );
742
743 ktEqn.ref().relax();
744 ktEqn.ref().boundaryManipulate(kt_.boundaryFieldRef());
745
746 solve(ktEqn);
747 bound(kt_, kMin_);
748
749
750 // Update total fluctuation kinetic energy dissipation rate
751 epsilon_ = kt_*omega_ + Dl + Dt;
752 bound(epsilon_, epsilonMin_);
753
754
755 // Re-calculate turbulent viscosity
756 nut_ = nuts + nutl;
758}
759
760
761// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
762
763} // End namespace RASModels
764} // End namespace incompressible
765} // End namespace Foam
766
767// ************************************************************************* //
scalar y
label k
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
Bound the given scalar field if it has gone unbounded.
surfaceScalarField & phi
void updateCoeffs()
Update the boundary condition coefficients.
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
void correctBoundaryConditions()
Correct boundary field.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Generic dimensioned Type class.
bool readIfPresent(const dictionary &dict)
Eddy viscosity turbulence model base class.
Definition: eddyViscosity.H:58
virtual void correct()=0
Solve the turbulence equations and correct the turbulence viscosity.
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
Low Reynolds-number k-kl-omega turbulence model for incompressible flows.
Definition: kkLOmega.H:113
tmp< volScalarField > DkEff(const volScalarField &alphaT) const
Return the effective diffusivity for k.
Definition: kkLOmega.H:242
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
Definition: kkLOmega.C:601
virtual void validate()
Validate the turbulence fields after construction.
Definition: kkLOmega.C:597
const volScalarField & y_
Wall distance.
Definition: kkLOmega.H:202
tmp< volScalarField > DomegaEff(const volScalarField &alphaT) const
Return the effective diffusivity for omega.
Definition: kkLOmega.H:251
virtual tmp< volScalarField > k() const
Return the total fluctuation kinetic energy.
Definition: kkLOmega.H:278
virtual bool read()
Re-read model coefficients if they have changed.
Definition: kkLOmega.C:557
BasicTurbulenceModel::transportModel transportModel
A class for managing temporary objects.
Definition: tmp.H:65
T & ref() const
Definition: tmpI.H:227
Interface to run-time selectable methods to calculate the distance-to-wall and normal-to-wall fields.
Definition: wallDist.H:78
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
U
Definition: pEqn.H:72
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:517
word timeName
Definition: getTimeIndex.H:3
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:54
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:48
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:48
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:48
zeroField SuSp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
zeroField Sp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
autoPtr< BasicCompressibleTurbulenceModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &alphaRhoPhi, const surfaceScalarField &phi, const typename BasicCompressibleTurbulenceModel::transportModel &transport, const word &propertiesName)
RASModel< turbulenceModel > RASModel
Namespace for OpenFOAM.
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
dimensionedScalar exp(const dimensionedScalar &ds)
const dimensionSet dimless
Dimensionless.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
dimensionedScalar pow3(const dimensionedScalar &ds)
volScalarField & bound(volScalarField &, const dimensionedScalar &lowerBound)
Bound the given scalar field if it has gone unbounded.
Definition: bound.C:35
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:598
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
dimensionedScalar pow4(const dimensionedScalar &ds)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
dimensionedTensor skew(const dimensionedTensor &dt)
labelList fv(nPoints)
volScalarField & nu
volScalarField & alpha
const dimensionedScalar & D
CEqn solve()