kOmegaSSTLM.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) 2016-2018 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 "kOmegaSSTLM.H"
27 
28 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32 namespace RASModels
33 {
34 
35 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
36 
37 template<class BasicTurbulenceModel>
38 tmp<volScalarField> kOmegaSSTLM<BasicTurbulenceModel>::F1
39 (
40  const volScalarField& CDkOmega
41 ) const
42 {
43  const volScalarField Ry(this->y_*sqrt(this->k_)/this->nu());
44  const volScalarField F3(exp(-pow(Ry/120.0, 8)));
45 
46  return max(kOmegaSST<BasicTurbulenceModel>::F1(CDkOmega), F3);
47 }
48 
49 
50 template<class BasicTurbulenceModel>
52 (
54 ) const
55 {
56  return gammaIntEff_*kOmegaSST<BasicTurbulenceModel>::Pk(G);
57 }
58 
59 
60 template<class BasicTurbulenceModel>
62 (
63  const volScalarField::Internal& F1,
64  const volScalarField::Internal& F2
65 ) const
66 {
67  return
68  min(max(gammaIntEff_, scalar(0.1)), scalar(1))
70 }
71 
72 
73 template<class BasicTurbulenceModel>
75 (
76  const volScalarField::Internal& Us,
77  const volScalarField::Internal& Omega,
78  const volScalarField::Internal& nu
79 ) const
80 {
81  const volScalarField::Internal& omega = this->omega_();
82  const volScalarField::Internal& y = this->y_();
83 
84  const volScalarField::Internal delta(375*Omega*nu*ReThetat_()*y/sqr(Us));
85  const volScalarField::Internal ReOmega(sqr(y)*omega/nu);
86  const volScalarField::Internal Fwake(exp(-sqr(ReOmega/1e5)));
87 
89  (
91  (
92  IOobject::groupName("Fthetat", this->alphaRhoPhi_.group()),
93  min
94  (
95  max
96  (
97  Fwake*exp(-pow4((y/delta))),
98  (1 - sqr((gammaInt_() - 1.0/ce2_)/(1 - 1.0/ce2_)))
99  ),
100  scalar(1)
101  )
102  )
103  );
104 }
105 
106 
107 template<class BasicTurbulenceModel>
110 {
112  (
114  (
115  IOobject
116  (
117  IOobject::groupName("ReThetac", this->alphaRhoPhi_.group()),
118  this->runTime_.timeName(),
119  this->mesh_
120  ),
121  this->mesh_,
122  dimless
123  )
124  );
125  volScalarField::Internal& ReThetac = tReThetac.ref();
126 
127  forAll(ReThetac, celli)
128  {
129  const scalar ReThetat = ReThetat_[celli];
130 
131  ReThetac[celli] =
132  ReThetat <= 1870
133  ?
134  ReThetat
135  - 396.035e-2
136  + 120.656e-4*ReThetat
137  - 868.230e-6*sqr(ReThetat)
138  + 696.506e-9*pow3(ReThetat)
139  - 174.105e-12*pow4(ReThetat)
140  :
141  ReThetat - 593.11 - 0.482*(ReThetat - 1870);
142  }
143 
144  return tReThetac;
145 }
146 
147 
148 template<class BasicTurbulenceModel>
150 (
151  const volScalarField::Internal& nu
152 ) const
153 {
155  (
157  (
158  IOobject
159  (
160  IOobject::groupName("Flength", this->alphaRhoPhi_.group()),
161  this->runTime_.timeName(),
162  this->mesh_
163  ),
164  this->mesh_,
165  dimless
166  )
167  );
168  volScalarField::Internal& Flength = tFlength.ref();
169 
170  const volScalarField::Internal& omega = this->omega_();
171  const volScalarField::Internal& y = this->y_();
172 
173  forAll(ReThetat_, celli)
174  {
175  const scalar ReThetat = ReThetat_[celli];
176 
177  if (ReThetat < 400)
178  {
179  Flength[celli] =
180  398.189e-1
181  - 119.270e-4*ReThetat
182  - 132.567e-6*sqr(ReThetat);
183  }
184  else if (ReThetat < 596)
185  {
186  Flength[celli] =
187  263.404
188  - 123.939e-2*ReThetat
189  + 194.548e-5*sqr(ReThetat)
190  - 101.695e-8*pow3(ReThetat);
191  }
192  else if (ReThetat < 1200)
193  {
194  Flength[celli] = 0.5 - 3e-4*(ReThetat - 596);
195  }
196  else
197  {
198  Flength[celli] = 0.3188;
199  }
200 
201  const scalar Fsublayer =
202  exp(-sqr(sqr(y[celli])*omega[celli]/(200*nu[celli])));
203 
204  Flength[celli] = Flength[celli]*(1 - Fsublayer) + 40*Fsublayer;
205  }
206 
207  return tFlength;
208 }
209 
210 
211 template<class BasicTurbulenceModel>
213 (
214  const volScalarField::Internal& Us,
215  const volScalarField::Internal& dUsds,
216  const volScalarField::Internal& nu
217 ) const
218 {
220  (
222  (
223  IOobject
224  (
225  IOobject::groupName("ReThetat0", this->alphaRhoPhi_.group()),
226  this->runTime_.timeName(),
227  this->mesh_
228  ),
229  this->mesh_,
230  dimless
231  )
232  );
233  volScalarField::Internal& ReThetat0 = tReThetat0.ref();
234 
235  const volScalarField& k = this->k_;
236 
237  label maxIter = 0;
238 
239  forAll(ReThetat0, celli)
240  {
241  const scalar Tu
242  (
243  max(100*sqrt((2.0/3.0)*k[celli])/Us[celli], scalar(0.027))
244  );
245 
246  // Initialize lambda to zero.
247  // If lambda were cached between time-steps convergence would be faster
248  // starting from the previous time-step value.
249  scalar lambda = 0;
250 
251  scalar lambdaErr;
252  scalar thetat;
253  label iter = 0;
254 
255  do
256  {
257  // Previous iteration lambda for convergence test
258  const scalar lambda0 = lambda;
259 
260  if (Tu <= 1.3)
261  {
262  const scalar Flambda =
263  dUsds[celli] <= 0
264  ?
265  1
266  - (
267  - 12.986*lambda
268  - 123.66*sqr(lambda)
269  - 405.689*pow3(lambda)
270  )*exp(-pow(Tu/1.5, 1.5))
271  :
272  1
273  + 0.275*(1 - exp(-35*lambda))
274  *exp(-Tu/0.5);
275 
276  thetat =
277  (1173.51 - 589.428*Tu + 0.2196/sqr(Tu))
278  *Flambda*nu[celli]
279  /Us[celli];
280  }
281  else
282  {
283  const scalar Flambda =
284  dUsds[celli] <= 0
285  ?
286  1
287  - (
288  -12.986*lambda
289  -123.66*sqr(lambda)
290  -405.689*pow3(lambda)
291  )*exp(-pow(Tu/1.5, 1.5))
292  :
293  1
294  + 0.275*(1 - exp(-35*lambda))
295  *exp(-2*Tu);
296 
297  thetat =
298  331.50*pow((Tu - 0.5658), -0.671)
299  *Flambda*nu[celli]/Us[celli];
300  }
301 
302  lambda = sqr(thetat)/nu[celli]*dUsds[celli];
303  lambda = max(min(lambda, 0.1), -0.1);
304 
305  lambdaErr = mag(lambda - lambda0);
306 
307  maxIter = max(maxIter, ++iter);
308 
309  } while (lambdaErr > lambdaErr_);
310 
311  ReThetat0[celli] = max(thetat*Us[celli]/nu[celli], scalar(20));
312  }
313 
314  if (maxIter > maxLambdaIter_)
315  {
317  << "Number of lambda iterations exceeds maxLambdaIter("
318  << maxLambdaIter_ << ')'<< endl;
319  }
320 
321  return tReThetat0;
322 }
323 
324 
325 template<class BasicTurbulenceModel>
327 (
328  const volScalarField::Internal& Rev,
329  const volScalarField::Internal& ReThetac,
330  const volScalarField::Internal& RT
331 ) const
332 {
333  const volScalarField::Internal Fonset1(Rev/(2.193*ReThetac));
334 
335  const volScalarField::Internal Fonset2
336  (
337  min(max(Fonset1, pow4(Fonset1)), scalar(2))
338  );
339 
340  const volScalarField::Internal Fonset3(max(1 - pow3(RT/2.5), scalar(0)));
341 
343  (
345  (
346  IOobject::groupName("Fonset", this->alphaRhoPhi_.group()),
347  max(Fonset2 - Fonset3, scalar(0))
348  )
349  );
350 }
351 
352 
353 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
354 
355 template<class BasicTurbulenceModel>
357 (
358  const alphaField& alpha,
359  const rhoField& rho,
360  const volVectorField& U,
361  const surfaceScalarField& alphaRhoPhi,
362  const surfaceScalarField& phi,
363  const transportModel& transport,
364  const word& propertiesName,
365  const word& type
366 )
367 :
369  (
370  alpha,
371  rho,
372  U,
373  alphaRhoPhi,
374  phi,
375  transport,
376  propertiesName
377  ),
378 
379  ca1_
380  (
382  (
383  "ca1",
384  this->coeffDict_,
385  2
386  )
387  ),
388  ca2_
389  (
391  (
392  "ca2",
393  this->coeffDict_,
394  0.06
395  )
396  ),
397  ce1_
398  (
400  (
401  "ce1",
402  this->coeffDict_,
403  1
404  )
405  ),
406  ce2_
407  (
409  (
410  "ce2",
411  this->coeffDict_,
412  50
413  )
414  ),
415  cThetat_
416  (
418  (
419  "cThetat",
420  this->coeffDict_,
421  0.03
422  )
423  ),
424  sigmaThetat_
425  (
427  (
428  "sigmaThetat",
429  this->coeffDict_,
430  2
431  )
432  ),
433  lambdaErr_
434  (
435  this->coeffDict_.lookupOrDefault("lambdaErr", 1e-6)
436  ),
437  maxLambdaIter_
438  (
439  this->coeffDict_.lookupOrDefault("maxLambdaIter", 10)
440  ),
441  deltaU_("deltaU", dimVelocity, small),
442 
443  ReThetat_
444  (
445  IOobject
446  (
447  IOobject::groupName("ReThetat", alphaRhoPhi.group()),
448  this->runTime_.timeName(),
449  this->mesh_,
452  ),
453  this->mesh_
454  ),
455 
456  gammaInt_
457  (
458  IOobject
459  (
460  IOobject::groupName("gammaInt", alphaRhoPhi.group()),
461  this->runTime_.timeName(),
462  this->mesh_,
465  ),
466  this->mesh_
467  ),
468 
469  gammaIntEff_
470  (
471  IOobject
472  (
473  IOobject::groupName("gammaIntEff", alphaRhoPhi.group()),
474  this->runTime_.timeName(),
475  this->mesh_
476  ),
477  this->mesh_,
478  dimensionedScalar("0", dimless, 0)
479  )
480 {}
481 
482 
483 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
484 
485 template<class BasicTurbulenceModel>
487 {
489  {
490  ca1_.readIfPresent(this->coeffDict());
491  ca2_.readIfPresent(this->coeffDict());
492  ce1_.readIfPresent(this->coeffDict());
493  ce2_.readIfPresent(this->coeffDict());
494  sigmaThetat_.readIfPresent(this->coeffDict());
495  cThetat_.readIfPresent(this->coeffDict());
496  this->coeffDict().readIfPresent("lambdaErr", lambdaErr_);
497  this->coeffDict().readIfPresent("maxLambdaIter", maxLambdaIter_);
498 
499  return true;
500  }
501  else
502  {
503  return false;
504  }
505 }
506 
507 
508 template<class BasicTurbulenceModel>
510 {
511  // Local references
512  const alphaField& alpha = this->alpha_;
513  const rhoField& rho = this->rho_;
514  const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
515  const volVectorField& U = this->U_;
516  const volScalarField& k = this->k_;
517  const volScalarField& omega = this->omega_;
518  const tmp<volScalarField> tnu = this->nu();
519  const volScalarField::Internal& nu = tnu()();
520  const volScalarField::Internal& y = this->y_();
521  fv::options& fvOptions(fv::options::New(this->mesh_));
522 
523  // Fields derived from the velocity gradient
524  tmp<volTensorField> tgradU = fvc::grad(U);
525  const volScalarField::Internal Omega(sqrt(2*magSqr(skew(tgradU()()))));
526  const volScalarField::Internal S(sqrt(2*magSqr(symm(tgradU()()))));
527  const volScalarField::Internal Us(max(mag(U()), deltaU_));
528  const volScalarField::Internal dUsds((U() & (U() & tgradU()()))/sqr(Us));
529  tgradU.clear();
530 
531  const volScalarField::Internal Fthetat(this->Fthetat(Us, Omega, nu));
532 
533  {
534  const volScalarField::Internal t(500*nu/sqr(Us));
535  const volScalarField::Internal Pthetat
536  (
537  alpha()*rho()*(cThetat_/t)*(1 - Fthetat)
538  );
539 
540  // Transition onset momentum-thickness Reynolds number equation
541  tmp<fvScalarMatrix> ReThetatEqn
542  (
543  fvm::ddt(alpha, rho, ReThetat_)
544  + fvm::div(alphaRhoPhi, ReThetat_)
545  - fvm::laplacian(alpha*rho*DReThetatEff(), ReThetat_)
546  ==
547  Pthetat*ReThetat0(Us, dUsds, nu) - fvm::Sp(Pthetat, ReThetat_)
548  + fvOptions(alpha, rho, ReThetat_)
549  );
550 
551  ReThetatEqn.ref().relax();
552  fvOptions.constrain(ReThetatEqn.ref());
553  solve(ReThetatEqn);
554  fvOptions.correct(ReThetat_);
555  bound(ReThetat_, 0);
556  }
557 
558  const volScalarField::Internal ReThetac(this->ReThetac());
559  const volScalarField::Internal Rev(sqr(y)*S/nu);
560  const volScalarField::Internal RT(k()/(nu*omega()));
561 
562  {
563  const volScalarField::Internal Pgamma
564  (
565  alpha()*rho()
566  *ca1_*Flength(nu)*S*sqrt(gammaInt_()*Fonset(Rev, ReThetac, RT))
567  );
568 
569  const volScalarField::Internal Fturb(exp(-pow4(0.25*RT)));
570 
571  const volScalarField::Internal Egamma
572  (
573  alpha()*rho()*ca2_*Omega*Fturb*gammaInt_()
574  );
575 
576  // Intermittency equation
577  tmp<fvScalarMatrix> gammaIntEqn
578  (
579  fvm::ddt(alpha, rho, gammaInt_)
580  + fvm::div(alphaRhoPhi, gammaInt_)
581  - fvm::laplacian(alpha*rho*DgammaIntEff(), gammaInt_)
582  ==
583  Pgamma - fvm::Sp(ce1_*Pgamma, gammaInt_)
584  + Egamma - fvm::Sp(ce2_*Egamma, gammaInt_)
585  + fvOptions(alpha, rho, gammaInt_)
586  );
587 
588  gammaIntEqn.ref().relax();
589  fvOptions.constrain(gammaIntEqn.ref());
590  solve(gammaIntEqn);
591  fvOptions.correct(gammaInt_);
592  bound(gammaInt_, 0);
593  }
594 
595  const volScalarField::Internal Freattach(exp(-pow4(RT/20.0)));
596  const volScalarField::Internal gammaSep
597  (
598  min(2*max(Rev/(3.235*ReThetac) - 1, scalar(0))*Freattach, scalar(2))
599  *Fthetat
600  );
601 
602  gammaIntEff_ = max(gammaInt_(), gammaSep);
603 }
604 
605 
606 template<class BasicTurbulenceModel>
608 {
609  if (!this->turbulence_)
610  {
611  return;
612  }
613 
614  // Correct ReThetat and gammaInt
615  correctReThetatGammaInt();
616 
617  // Correct k and omega
619 }
620 
621 
622 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
623 
624 } // End namespace RASModels
625 } // End namespace Foam
626 
627 // ************************************************************************* //
scalar delta
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
BasicTurbulenceModel::transportModel transportModel
Definition: kOmegaSST.H:70
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
static dimensioned< Type > lookupOrAddToDict(const word &, dictionary &, const dimensionSet &dims=dimless, const Type &defaultValue=pTraits< Type >::zero)
Construct from dictionary, with default value.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:230
fv::options & fvOptions
BasicTurbulenceModel::rhoField rhoField
Definition: kOmegaSST.H:69
surfaceScalarField & phi
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
tmp< volScalarField::Internal > ReThetat0(const volScalarField::Internal &Us, const volScalarField::Internal &dUsds, const volScalarField::Internal &nu) const
Return the transition onset momentum-thickness Reynolds number.
Definition: kOmegaSSTLM.C:213
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedTensor skew(const dimensionedTensor &dt)
void correct(GeometricField< Type, fvPatchField, volMesh > &field)
Apply correction to field.
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
label k
Boltzmann constant.
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
Finite-volume options.
Definition: fvOptions.H:52
DimensionedField< scalar, volMesh > Internal
Type of the internal field from which this GeometricField is derived.
virtual tmp< volScalarField::Internal > epsilonByk(const volScalarField::Internal &F1, const volScalarField::Internal &F2) const
Modified form of the k-omega SST epsilon/k.
Definition: kOmegaSSTLM.C:62
tmp< volScalarField::Internal > ReThetac() const
Empirical correlation for critical Reynolds number where the.
Definition: kOmegaSSTLM.C:109
scalar y
dimensionedScalar exp(const dimensionedScalar &ds)
#define F3(B, C, D)
Definition: SHA1.C:175
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)
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
Definition: kOmegaSSTLM.C:607
Specialisation for RAS of the generic kOmegaSSTBase base class. For more information, see Description of kOmegaSSTBase.H.
Definition: kOmegaSST.H:57
void constrain(fvMatrix< Type > &eqn)
Apply constraints to equation.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
virtual bool read()
Re-read model coefficients if they have changed.
Definition: kOmegaSSTLM.C:486
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
tmp< volScalarField::Internal > Fonset(const volScalarField::Internal &Rev, const volScalarField::Internal &ReThetac, const volScalarField::Internal &RT) const
Transition onset location control function.
Definition: kOmegaSSTLM.C:327
tmp< volScalarField::Internal > Flength(const volScalarField::Internal &nu) const
Empirical correlation that controls the length of the.
Definition: kOmegaSSTLM.C:150
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
virtual tmp< volScalarField::Internal > Pk(const volScalarField::Internal &G) const
Return k production rate.
tensor Ry(const scalar &omega)
Rotational transformation tensor about the y-axis by omega radians.
Definition: transform.H:95
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
U
Definition: pEqn.H:72
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
#define WarningInFunction
Report a warning using Foam::Warning.
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
Langtry-Menter 4-equation transitional SST model based on the k-omega-SST RAS model.
Definition: kOmegaSSTLM.H:101
dimensionedScalar pow4(const dimensionedScalar &ds)
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
dimensioned< scalar > mag(const dimensioned< Type > &)
void correctReThetatGammaInt()
Solve the turbulence equations and correct the turbulence viscosity.
Definition: kOmegaSSTLM.C:509
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:98
virtual tmp< volScalarField::Internal > epsilonByk(const volScalarField::Internal &F1, const volScalarField::Internal &F2) const
Return epsilon/k which for standard RAS is betaStar*omega.
tmp< volScalarField::Internal > Fthetat(const volScalarField::Internal &Us, const volScalarField::Internal &Omega, const volScalarField::Internal &nu) const
Freestream blending-function.
Definition: kOmegaSSTLM.C:75
A class for managing temporary objects.
Definition: PtrList.H:53
static options & New(const fvMesh &mesh)
Construct fvOptions and register to datbase if not present.
Definition: fvOptions.C:101
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
BasicTurbulenceModel::alphaField alphaField
Definition: kOmegaSST.H:68
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
virtual tmp< volScalarField > F1(const volScalarField &CDkOmega) const
Modified form of the k-omega SST F1 function.
Definition: kOmegaSSTLM.C:39
volScalarField & nu
Namespace for OpenFOAM.
const dimensionSet dimVelocity
virtual tmp< volScalarField::Internal > Pk(const volScalarField::Internal &G) const
Modified form of the k-omega SST k production rate.
Definition: kOmegaSSTLM.C:52