DAkOmegaSSTLM.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v3
5 
6  This file is modified from OpenFOAM's source code
7  src/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTLM/kOmegaSSTLM.C
8 
9  OpenFOAM: The Open Source CFD Toolbox
10 
11  Copyright (C): 2011-2016 OpenFOAM Foundation
12 
13  OpenFOAM License:
14 
15  OpenFOAM is free software: you can redistribute it and/or modify it
16  under the terms of the GNU General Public License as published by
17  the Free Software Foundation, either version 3 of the License, or
18  (at your option) any later version.
19 
20  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
21  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23  for more details.
24 
25  You should have received a copy of the GNU General Public License
26  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
27 
28 \*---------------------------------------------------------------------------*/
29 
30 #include "DAkOmegaSSTLM.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 defineTypeNameAndDebug(DAkOmegaSSTLM, 0);
38 addToRunTimeSelectionTable(DATurbulenceModel, DAkOmegaSSTLM, dictionary);
39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40 
41 DAkOmegaSSTLM::DAkOmegaSSTLM(
42  const word modelType,
43  const fvMesh& mesh,
44  const DAOption& daOption)
45  : DATurbulenceModel(modelType, mesh, daOption),
46  // SST parameters
47  alphaK1_(dimensioned<scalar>::lookupOrAddToDict(
48  "alphaK1",
49  this->coeffDict_,
50  0.85)),
51  alphaK2_(dimensioned<scalar>::lookupOrAddToDict(
52  "alphaK2",
53  this->coeffDict_,
54  1.0)),
55  alphaOmega1_(dimensioned<scalar>::lookupOrAddToDict(
56  "alphaOmega1",
57  this->coeffDict_,
58  0.5)),
59  alphaOmega2_(dimensioned<scalar>::lookupOrAddToDict(
60  "alphaOmega2",
61  this->coeffDict_,
62  0.856)),
63  gamma1_(dimensioned<scalar>::lookupOrAddToDict(
64  "gamma1",
65  this->coeffDict_,
66  5.0 / 9.0)),
67  gamma2_(dimensioned<scalar>::lookupOrAddToDict(
68  "gamma2",
69  this->coeffDict_,
70  0.44)),
71  beta1_(dimensioned<scalar>::lookupOrAddToDict(
72  "beta1",
73  this->coeffDict_,
74  0.075)),
75  beta2_(dimensioned<scalar>::lookupOrAddToDict(
76  "beta2",
77  this->coeffDict_,
78  0.0828)),
79  betaStar_(dimensioned<scalar>::lookupOrAddToDict(
80  "betaStar",
81  this->coeffDict_,
82  0.09)),
83  a1_(dimensioned<scalar>::lookupOrAddToDict(
84  "a1",
85  this->coeffDict_,
86  0.31)),
87  b1_(dimensioned<scalar>::lookupOrAddToDict(
88  "b1",
89  this->coeffDict_,
90  1.0)),
91  c1_(dimensioned<scalar>::lookupOrAddToDict(
92  "c1",
93  this->coeffDict_,
94  10.0)),
95  F3_(Switch::lookupOrAddToDict(
96  "F3",
97  this->coeffDict_,
98  false)),
99  ca1_(dimensionedScalar::lookupOrAddToDict(
100  "ca1",
101  this->coeffDict_,
102  2)),
103  ca2_(dimensionedScalar::lookupOrAddToDict(
104  "ca2",
105  this->coeffDict_,
106  0.06)),
107  ce1_(dimensionedScalar::lookupOrAddToDict(
108  "ce1",
109  this->coeffDict_,
110  1)),
111  ce2_(dimensionedScalar::lookupOrAddToDict(
112  "ce2",
113  this->coeffDict_,
114  50)),
115  cThetat_(dimensionedScalar::lookupOrAddToDict(
116  "cThetat",
117  this->coeffDict_,
118  0.03)),
119  sigmaThetat_(dimensionedScalar::lookupOrAddToDict(
120  "sigmaThetat",
121  this->coeffDict_,
122  2)),
123  lambdaErr_(this->coeffDict_.lookupOrDefault("lambdaErr", 1e-6)),
124  maxLambdaIter_(this->coeffDict_.lookupOrDefault("maxLambdaIter", 10)),
125  deltaU_("deltaU", dimVelocity, SMALL),
126  // Augmented variables
127  omega_(const_cast<volScalarField&>(
128  mesh_.thisDb().lookupObject<volScalarField>("omega"))),
129  omegaRes_(
130  IOobject(
131  "omegaRes",
132  mesh.time().timeName(),
133  mesh,
134  IOobject::NO_READ,
135  IOobject::NO_WRITE),
136  mesh,
137 #ifdef CompressibleFlow
138  dimensionedScalar("omegaRes", dimensionSet(1, -3, -2, 0, 0, 0, 0), 0.0),
139 #endif
140 #ifdef IncompressibleFlow
141  dimensionedScalar("omegaRes", dimensionSet(0, 0, -2, 0, 0, 0, 0), 0.0),
142 #endif
143  zeroGradientFvPatchField<scalar>::typeName),
144  k_(const_cast<volScalarField&>(
145  mesh_.thisDb().lookupObject<volScalarField>("k"))),
146  kRes_(
147  IOobject(
148  "kRes",
149  mesh.time().timeName(),
150  mesh,
151  IOobject::NO_READ,
152  IOobject::NO_WRITE),
153  mesh,
154 #ifdef CompressibleFlow
155  dimensionedScalar("kRes", dimensionSet(1, -1, -3, 0, 0, 0, 0), 0.0),
156 #endif
157 #ifdef IncompressibleFlow
158  dimensionedScalar("kRes", dimensionSet(0, 2, -3, 0, 0, 0, 0), 0.0),
159 #endif
160  zeroGradientFvPatchField<scalar>::typeName),
161  ReThetat_(const_cast<volScalarField&>(
162  mesh_.thisDb().lookupObject<volScalarField>("ReThetat"))),
163  ReThetatRes_(
164  IOobject(
165  "ReThetatRes",
166  mesh_.time().timeName(),
167  mesh_,
168  IOobject::NO_READ,
169  IOobject::NO_WRITE),
170  mesh_,
171 #ifdef CompressibleFlow
172  dimensionedScalar("ReThetatRes", dimensionSet(1, -3, -1, 0, 0, 0, 0), 0.0),
173 #endif
174 #ifdef IncompressibleFlow
175  dimensionedScalar("ReThetatRes", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0.0),
176 #endif
177  zeroGradientFvPatchScalarField::typeName),
178  gammaInt_(const_cast<volScalarField&>(
179  mesh_.thisDb().lookupObject<volScalarField>("gammaInt"))),
180  gammaIntRes_(
181  IOobject(
182  "gammaIntRes",
183  mesh_.time().timeName(),
184  mesh_,
185  IOobject::NO_READ,
186  IOobject::NO_WRITE),
187  mesh_,
188 #ifdef CompressibleFlow
189  dimensionedScalar("gammaIntRes", dimensionSet(1, -3, -1, 0, 0, 0, 0), 0.0),
190 #endif
191 #ifdef IncompressibleFlow
192  dimensionedScalar("gammaIntRes", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0.0),
193 #endif
194  zeroGradientFvPatchScalarField::typeName),
195  gammaIntEff_(const_cast<volScalarField::Internal&>(
196  mesh_.thisDb().lookupObject<volScalarField::Internal>("gammaIntEff"))),
197  y_(mesh_.thisDb().lookupObject<volScalarField>("yWall"))
198 {
199 
200  // initialize printInterval_ we need to check whether it is a steady state
201  // or unsteady primal solver
202  IOdictionary fvSchemes(
203  IOobject(
204  "fvSchemes",
205  mesh.time().system(),
206  mesh,
207  IOobject::MUST_READ,
208  IOobject::NO_WRITE,
209  false));
210  word ddtScheme = word(fvSchemes.subDict("ddtSchemes").lookup("default"));
211  if (ddtScheme == "steadyState")
212  {
213  printInterval_ =
214  daOption.getAllOptions().lookupOrDefault<label>("printInterval", 100);
215  }
216  else
217  {
218  printInterval_ =
219  daOption.getAllOptions().lookupOrDefault<label>("printIntervalUnsteady", 500);
220  }
221 
222  // calculate the size of omegaWallFunction faces
223  label nWallFaces = 0;
224  forAll(omega_.boundaryField(), patchI)
225  {
226  if (omega_.boundaryField()[patchI].type() == "omegaWallFunction"
227  && omega_.boundaryField()[patchI].size() > 0)
228  {
229  forAll(omega_.boundaryField()[patchI], faceI)
230  {
231  nWallFaces++;
232  }
233  }
234  }
235 
236  // initialize omegaNearWall
237  omegaNearWall_.setSize(nWallFaces);
238 }
239 
240 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
241 
242 // SA member functions. these functions are copied from
243 tmp<volScalarField> DAkOmegaSSTLM::F1SST(
244  const volScalarField& CDkOmega) const
245 {
246 
247  tmp<volScalarField> CDkOmegaPlus = max(
248  CDkOmega,
249  dimensionedScalar("1.0e-10", dimless / sqr(dimTime), 1.0e-10));
250 
251  tmp<volScalarField> arg1 = min(
252  min(
253  max(
254  (scalar(1) / betaStar_) * sqrt(k_) / (omega_ * y_),
255  scalar(500) * (this->nu()) / (sqr(y_) * omega_)),
256  (scalar(4) * alphaOmega2_) * k_ / (CDkOmegaPlus * sqr(y_))),
257  scalar(10));
258 
259  return tanh(pow4(arg1));
260 }
261 
262 tmp<volScalarField> DAkOmegaSSTLM::F2() const
263 {
264 
265  tmp<volScalarField> arg2 = min(
266  max(
267  (scalar(2) / betaStar_) * sqrt(k_) / (omega_ * y_),
268  scalar(500) * (this->nu()) / (sqr(y_) * omega_)),
269  scalar(100));
270 
271  return tanh(sqr(arg2));
272 }
273 
274 tmp<volScalarField> DAkOmegaSSTLM::F3() const
275 {
276 
277  tmp<volScalarField> arg3 = min(
278  150 * (this->nu()) / (omega_ * sqr(y_)),
279  scalar(10));
280 
281  return 1 - tanh(pow4(arg3));
282 }
283 
284 tmp<volScalarField> DAkOmegaSSTLM::F23() const
285 {
286  tmp<volScalarField> f23(F2());
287 
288  if (F3_)
289  {
290  f23.ref() *= F3();
291  }
292 
293  return f23;
294 }
295 
296 tmp<volScalarField::Internal> DAkOmegaSSTLM::GbyNu(
297  const volScalarField::Internal& GbyNu0,
298  const volScalarField::Internal& F2,
299  const volScalarField::Internal& S2) const
300 {
301  return min(
302  GbyNu0,
303  (c1_ / a1_) * betaStar_ * omega_()
304  * max(a1_ * omega_(), b1_ * F2 * sqrt(S2)));
305 }
306 
307 tmp<volScalarField::Internal> DAkOmegaSSTLM::PkSST(
308  const volScalarField::Internal& G) const
309 {
310  return min(G, (c1_ * betaStar_) * k_() * omega_());
311 }
312 
313 tmp<volScalarField::Internal> DAkOmegaSSTLM::epsilonBykSST(
314  const volScalarField& F1,
315  const volTensorField& gradU) const
316 {
317  return betaStar_ * omega_();
318 }
319 
320 tmp<fvScalarMatrix> DAkOmegaSSTLM::kSource() const
321 {
322  const volScalarField& rho = rho_;
323  return tmp<fvScalarMatrix>(
324  new fvScalarMatrix(
325  k_,
326  dimVolume * rho.dimensions() * k_.dimensions() / dimTime));
327 }
328 
329 tmp<fvScalarMatrix> DAkOmegaSSTLM::omegaSource() const
330 {
331  const volScalarField& rho = rho_;
332  return tmp<fvScalarMatrix>(
333  new fvScalarMatrix(
334  omega_,
335  dimVolume * rho.dimensions() * omega_.dimensions() / dimTime));
336 }
337 
338 tmp<fvScalarMatrix> DAkOmegaSSTLM::Qsas(
339  const volScalarField::Internal& S2,
340  const volScalarField::Internal& gamma,
341  const volScalarField::Internal& beta) const
342 {
343  const volScalarField& rho = rho_;
344  return tmp<fvScalarMatrix>(
345  new fvScalarMatrix(
346  omega_,
347  dimVolume * rho.dimensions() * omega_.dimensions() / dimTime));
348 }
349 
350 // SSTLM functions
351 tmp<volScalarField> DAkOmegaSSTLM::F1(
352  const volScalarField& CDkOmega) const
353 {
354  const volScalarField Ry(y_ * sqrt(k_) / this->nu());
355  const volScalarField F3(exp(-pow(Ry / 120.0, 8)));
356 
357  return max(this->F1SST(CDkOmega), F3);
358 }
359 
360 tmp<volScalarField::Internal> DAkOmegaSSTLM::Pk(
361  const volScalarField::Internal& G) const
362 {
363  return gammaIntEff_ * this->PkSST(G);
364 }
365 
366 tmp<volScalarField::Internal> DAkOmegaSSTLM::epsilonByk(
367  const volScalarField& F1,
368  const volTensorField& gradU) const
369 {
370  return min(max(gammaIntEff_, scalar(0.1)), scalar(1))
371  * this->epsilonBykSST(F1, gradU);
372 }
373 
374 tmp<volScalarField::Internal> DAkOmegaSSTLM::Fthetat(
375  const volScalarField::Internal& Us,
376  const volScalarField::Internal& Omega,
377  const volScalarField::Internal& nu) const
378 {
379  const volScalarField::Internal& omega = omega_();
380  const volScalarField::Internal& y = y_();
381 
382  const volScalarField::Internal delta(375 * Omega * nu * ReThetat_() * y / sqr(Us));
383  const volScalarField::Internal ReOmega(sqr(y) * omega / nu);
384  const volScalarField::Internal Fwake(exp(-sqr(ReOmega / 1e5)));
385 
386  return tmp<volScalarField::Internal>(
387  new volScalarField::Internal(
388  IOobject::groupName("Fthetat", U_.group()),
389  min(
390  max(
391  Fwake * exp(-pow4((y / delta))),
392  (1 - sqr((gammaInt_() - 1.0 / ce2_) / (1 - 1.0 / ce2_)))),
393  scalar(1))));
394 }
395 
396 tmp<volScalarField::Internal> DAkOmegaSSTLM::ReThetac() const
397 {
398 
399  tmp<volScalarField::Internal> tReThetac(
400  new volScalarField::Internal(
401  IOobject(
402  IOobject::groupName("ReThetac", U_.group()),
403  mesh_.time().timeName(),
404  mesh_),
405  mesh_,
406  dimless));
407  volScalarField::Internal& ReThetac = tReThetac.ref();
408 
409  forAll(ReThetac, celli)
410  {
411  const scalar ReThetat = ReThetat_[celli];
412 
413  ReThetac[celli] =
414  ReThetat <= 1870
415  ? scalar(ReThetat
416  - 396.035e-2
417  + 120.656e-4 * ReThetat
418  - 868.230e-6 * sqr(ReThetat)
419  + 696.506e-9 * pow3(ReThetat)
420  - 174.105e-12 * pow4(ReThetat))
421  : scalar(ReThetat - 593.11 - 0.482 * (ReThetat - 1870));
422  }
423 
424  return tReThetac;
425 }
426 
427 tmp<volScalarField::Internal> DAkOmegaSSTLM::Flength(
428  const volScalarField::Internal& nu) const
429 {
430 
431  tmp<volScalarField::Internal> tFlength(
432  new volScalarField::Internal(
433  IOobject(
434  IOobject::groupName("Flength", U_.group()),
435  mesh_.time().timeName(),
436  mesh_),
437  mesh_,
438  dimless));
439  volScalarField::Internal& Flength = tFlength.ref();
440 
441  const volScalarField::Internal& omega = omega_();
442  const volScalarField::Internal& y = y_();
443 
444  forAll(ReThetat_, celli)
445  {
446  const scalar ReThetat = ReThetat_[celli];
447 
448  if (ReThetat < 400)
449  {
450  Flength[celli] =
451  398.189e-1
452  - 119.270e-4 * ReThetat
453  - 132.567e-6 * sqr(ReThetat);
454  }
455  else if (ReThetat < 596)
456  {
457  Flength[celli] =
458  263.404
459  - 123.939e-2 * ReThetat
460  + 194.548e-5 * sqr(ReThetat)
461  - 101.695e-8 * pow3(ReThetat);
462  }
463  else if (ReThetat < 1200)
464  {
465  Flength[celli] = 0.5 - 3e-4 * (ReThetat - 596);
466  }
467  else
468  {
469  Flength[celli] = 0.3188;
470  }
471 
472  const scalar Fsublayer =
473  exp(-sqr(sqr(y[celli]) * omega[celli] / (200 * nu[celli])));
474 
475  Flength[celli] = Flength[celli] * (1 - Fsublayer) + 40 * Fsublayer;
476  }
477 
478  return tFlength;
479 }
480 
481 tmp<volScalarField::Internal> DAkOmegaSSTLM::Fonset(
482  const volScalarField::Internal& Rev,
483  const volScalarField::Internal& ReThetac,
484  const volScalarField::Internal& RT) const
485 {
486  const volScalarField::Internal Fonset1(Rev / (2.193 * ReThetac));
487 
488  const volScalarField::Internal Fonset2(
489  min(max(Fonset1, pow4(Fonset1)), scalar(2)));
490 
491  const volScalarField::Internal Fonset3(max(1 - pow3(RT / 2.5), scalar(0)));
492 
493  return tmp<volScalarField::Internal>(
494  new volScalarField::Internal(
495  IOobject::groupName("Fonset", U_.group()),
496  max(Fonset2 - Fonset3, scalar(0))));
497 }
498 
499 tmp<volScalarField::Internal> DAkOmegaSSTLM::ReThetat0(
500  const volScalarField::Internal& Us,
501  const volScalarField::Internal& dUsds,
502  const volScalarField::Internal& nu) const
503 {
504 
505  tmp<volScalarField::Internal> tReThetat0(
506  new volScalarField::Internal(
507  IOobject(
508  IOobject::groupName("ReThetat0", U_.group()),
509  mesh_.time().timeName(),
510  mesh_),
511  mesh_,
512  dimless));
513  volScalarField::Internal& ReThetat0 = tReThetat0.ref();
514 
515  const volScalarField& k = k_;
516 
517  label maxIter = 0;
518 
519  forAll(ReThetat0, celli)
520  {
521  const scalar Tu(
522  max(100 * sqrt((2.0 / 3.0) * k[celli]) / Us[celli], scalar(0.027)));
523 
524  // Initialize lambda to zero.
525  // If lambda were cached between time-steps convergence would be faster
526  // starting from the previous time-step value.
527  scalar lambda = 0;
528 
529  scalar lambdaErr;
530  scalar thetat;
531  label iter = 0;
532 
533  do
534  {
535  // Previous iteration lambda for convergence test
536  const scalar lambda0 = lambda;
537 
538  if (Tu <= 1.3)
539  {
540  const scalar Flambda =
541  dUsds[celli] <= 0
542  ? scalar(1
543  - (-12.986 * lambda
544  - 123.66 * sqr(lambda)
545  - 405.689 * pow3(lambda))
546  * exp(-pow(Tu / 1.5, 1.5)))
547  : scalar(1
548  + 0.275 * (1 - exp(-35 * lambda))
549  * exp(-Tu / 0.5));
550 
551  thetat =
552  (1173.51 - 589.428 * Tu + 0.2196 / sqr(Tu))
553  * Flambda * nu[celli]
554  / Us[celli];
555  }
556  else
557  {
558  const scalar Flambda =
559  dUsds[celli] <= 0
560  ? scalar(1
561  - (-12.986 * lambda
562  - 123.66 * sqr(lambda)
563  - 405.689 * pow3(lambda))
564  * exp(-pow(Tu / 1.5, 1.5)))
565  : scalar(1
566  + 0.275 * (1 - exp(-35 * lambda))
567  * exp(-2 * Tu));
568 
569  thetat =
570  331.50 * pow((Tu - 0.5658), -0.671)
571  * Flambda * nu[celli] / Us[celli];
572  }
573 
574  lambda = sqr(thetat) / nu[celli] * dUsds[celli];
575  lambda = max(min(lambda, 0.1), -0.1);
576 
577  lambdaErr = mag(lambda - lambda0);
578 
579  maxIter = max(maxIter, ++iter);
580 
581  } while (lambdaErr > lambdaErr_);
582 
583  ReThetat0[celli] = max(thetat * Us[celli] / nu[celli], scalar(20));
584  }
585 
586  if (maxIter > maxLambdaIter_)
587  {
588  WarningInFunction
589  << "Number of lambda iterations exceeds maxLambdaIter("
590  << maxLambdaIter_ << ')' << endl;
591  }
592 
593  return tReThetat0;
594 }
595 
596 // Augmented functions
597 void DAkOmegaSSTLM::correctModelStates(wordList& modelStates) const
598 {
599  /*
600  Description:
601  Update the name in modelStates based on the selected physical model at runtime
602 
603  Example:
604  In DAStateInfo, if the modelStates reads:
605 
606  modelStates = {"nut"}
607 
608  then for the SA model, calling correctModelStates(modelStates) will give:
609 
610  modelStates={"nuTilda"}
611 
612  while calling correctModelStates(modelStates) for the SST model will give
613 
614  modelStates={"k","omega"}
615 
616  We don't udpate the names for the radiation model because users are
617  supposed to set modelStates={"G"}
618  */
619 
620  // For SST model, we need to replace nut with k, omega
621 
622  forAll(modelStates, idxI)
623  {
624  word stateName = modelStates[idxI];
625  if (stateName == "nut")
626  {
627  modelStates[idxI] = "omega";
628  modelStates.append("k");
629  modelStates.append("ReThetat");
630  modelStates.append("gammaInt");
631  }
632  }
633 }
634 
635 void DAkOmegaSSTLM::correctNut()
636 {
637  /*
638  Description:
639  Update nut based on other turbulence variables and update the BCs
640  Also update alphat if is present
641  */
642 
643  const volVectorField U = mesh_.thisDb().lookupObject<volVectorField>("U");
644  tmp<volTensorField> tgradU = fvc::grad(U);
645  volScalarField S2(2 * magSqr(symm(tgradU())));
646 
647  nut_ = a1_ * k_ / max(a1_ * omega_, b1_ * F23() * sqrt(S2));
648 
649  nut_.correctBoundaryConditions(); // nutkWallFunction: update wall face nut based on k
650 
651  // this is basically BasicTurbulenceModel::correctNut();
652  this->correctAlphat();
653 
654  return;
655 }
656 
657 void DAkOmegaSSTLM::correctBoundaryConditions()
658 {
659  /*
660  Description:
661  Update turbulence variable boundary values
662  */
663 
664  // correct the BCs for the perturbed fields
665  // kqWallFunction is a zero-gradient BC
666  k_.correctBoundaryConditions();
667 
668  ReThetat_.correctBoundaryConditions();
669  gammaInt_.correctBoundaryConditions();
670 }
671 
672 void DAkOmegaSSTLM::correctOmegaBoundaryConditions()
673 {
674  /*
675  Description:
676  this is a special treatment for omega BC because we cant directly call omega.
677  correctBoundaryConditions() because it will modify the internal omega and G that
678  are right next to walls. This will mess up adjoint Jacobians
679  To solve this issue,
680  1. we store the near wall omega before calling omega.correctBoundaryConditions()
681  2. call omega.correctBoundaryConditions()
682  3. Assign the stored near wall omega back
683  4. Apply a zeroGradient BC for omega at the wall patches
684  *********** NOTE *************
685  this treatment will obviously downgrade the accuracy of adjoint derivative since it is
686  not 100% consistent with what is used for the flow solver; however, our observation
687  shows that the impact is not very large.
688  */
689 
690  // save the perturbed omega at the wall
691  this->saveOmegaNearWall();
692  // correct omega boundary conditions, this includes updating wall face and near wall omega values,
693  // updating the inter-proc BCs
694  omega_.correctBoundaryConditions();
695  // reset the corrected omega near wall cell to its perturbed value
696  this->setOmegaNearWall();
697 }
698 
699 void DAkOmegaSSTLM::saveOmegaNearWall()
700 {
701  /*
702  Description:
703  Save the current near wall omega values to omegaNearWall_
704  */
705  label counterI = 0;
706  forAll(omega_.boundaryField(), patchI)
707  {
708  if (omega_.boundaryField()[patchI].type() == "omegaWallFunction"
709  and omega_.boundaryField()[patchI].size() > 0)
710  {
711  const UList<label>& faceCells = mesh_.boundaryMesh()[patchI].faceCells();
712  forAll(faceCells, faceI)
713  {
714  //Info<<"Near Wall cellI: "<<faceCells[faceI]<<endl;
715  omegaNearWall_[counterI] = omega_[faceCells[faceI]];
716  counterI++;
717  }
718  }
719  }
720  return;
721 }
722 
723 void DAkOmegaSSTLM::setOmegaNearWall()
724 {
725  /*
726  Description:
727  Set the current near wall omega values from omegaNearWall_
728  Here we also apply a zeroGradient BC to the wall faces
729  */
730  label counterI = 0;
731  forAll(omega_.boundaryField(), patchI)
732  {
733  if (omega_.boundaryField()[patchI].type() == "omegaWallFunction"
734  && omega_.boundaryField()[patchI].size() > 0)
735  {
736  const UList<label>& faceCells = mesh_.boundaryMesh()[patchI].faceCells();
737  forAll(faceCells, faceI)
738  {
739  omega_[faceCells[faceI]] = omegaNearWall_[counterI];
740  // zeroGradient BC
741  omega_.boundaryFieldRef()[patchI][faceI] = omega_[faceCells[faceI]];
742  counterI++;
743  }
744  }
745  }
746  return;
747 }
748 
749 void DAkOmegaSSTLM::updateIntermediateVariables()
750 {
751  /*
752  Description:
753  Update nut based on nuTilda. Note: we need to update nut and its BC since we
754  may have perturbed other turbulence vars that affect the nut values
755  */
756 
757  this->correctNut();
758 
759  // for SSTLM also update gammaIntEff_
760  // NOTE: this is not implemented yet!!!
761 }
762 
763 void DAkOmegaSSTLM::correctStateResidualModelCon(List<List<word>>& stateCon) const
764 {
765  /*
766  Description:
767  Update the original variable connectivity for the adjoint state
768  residuals in stateCon. Basically, we modify/add state variables based on the
769  original model variables defined in stateCon.
770 
771  Input:
772 
773  stateResCon: the connectivity levels for a state residual, defined in Foam::DAJacCon
774 
775  Example:
776  If stateCon reads:
777  stateCon=
778  {
779  {"U", "p", "nut"},
780  {"p"}
781  }
782 
783  For the SA turbulence model, calling this function for will get a new stateCon
784  stateCon=
785  {
786  {"U", "p", "nuTilda"},
787  {"p"}
788  }
789 
790  For the SST turbulence model, calling this function will give
791  stateCon=
792  {
793  {"U", "p", "k", "omega"},
794  {"p", "U"}
795  }
796  ***NOTE***: we add a extra level of U connectivity because nut is
797  related to grad(U), k, and omega in SST!
798  */
799 
800  label stateConSize = stateCon.size();
801  forAll(stateCon, idxI)
802  {
803  label addUCon = 0;
804  forAll(stateCon[idxI], idxJ)
805  {
806  word conStateName = stateCon[idxI][idxJ];
807  if (conStateName == "nut")
808  {
809  stateCon[idxI][idxJ] = "omega"; // replace nut with omega
810  stateCon[idxI].append("k"); // also add k for that level
811  stateCon[idxI].append("ReThetat");
812  stateCon[idxI].append("gammaInt");
813  addUCon = 1;
814  }
815  }
816  // add U for the current level and level+1 if it is not there yet
817  label isU;
818  if (addUCon == 1)
819  {
820  // first add U for the current level
821  isU = 0;
822  forAll(stateCon[idxI], idxJ)
823  {
824  word conStateName = stateCon[idxI][idxJ];
825  if (conStateName == "U")
826  {
827  isU = 1;
828  }
829  }
830  if (!isU)
831  {
832  stateCon[idxI].append("U");
833  }
834 
835  // now add U for level+1 if idxI is not the largest level
836  // if idxI is already the largest level, we have a problem
837  if (idxI != stateConSize - 1)
838  {
839  isU = 0;
840  forAll(stateCon[idxI + 1], idxJ)
841  {
842  word conStateName = stateCon[idxI + 1][idxJ];
843  if (conStateName == "U")
844  {
845  isU = 1;
846  }
847  }
848  if (!isU)
849  {
850  stateCon[idxI + 1].append("U");
851  }
852  }
853  else
854  {
855  FatalErrorIn(
856  "In DAStateInfo, nut shows in the largest connectivity level! "
857  "This is not supported!")
858  << exit(FatalError);
859  }
860  }
861  }
862 }
863 
864 void DAkOmegaSSTLM::addModelResidualCon(HashTable<List<List<word>>>& allCon) const
865 {
866  /*
867  Description:
868  Add the connectivity levels for all physical model residuals to allCon
869 
870  Input:
871  allCon: the connectivity levels for all state residual, defined in DAJacCon
872 
873  Example:
874  If stateCon reads:
875  allCon=
876  {
877  "URes":
878  {
879  {"U", "p", "nut"},
880  {"p"}
881  }
882  }
883 
884  For the SA turbulence model, calling this function for will get a new stateCon,
885  something like this:
886  allCon=
887  {
888  "URes":
889  {
890  {"U", "p", "nuTilda"},
891  {"p"}
892  },
893  "nuTildaRes":
894  {
895  {"U", "phi", "nuTilda"},
896  {"U"}
897  }
898  }
899 
900  */
901 
902  word pName;
903 
904  if (mesh_.thisDb().foundObject<volScalarField>("p"))
905  {
906  pName = "p";
907  }
908  else if (mesh_.thisDb().foundObject<volScalarField>("p_rgh"))
909  {
910  pName = "p_rgh";
911  }
912  else
913  {
914  FatalErrorIn(
915  "Neither p nor p_rgh was found in mesh.thisDb()!"
916  "addModelResidualCon failed to setup turbulence residuals!")
917  << exit(FatalError);
918  }
919 
920  // NOTE: for compressible flow, it depends on rho so we need to add T and p
921 #ifdef IncompressibleFlow
922  allCon.set(
923  "omegaRes",
924  {
925  {"U", "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
926  {"U", "omega", "k", "ReThetat", "gammaInt"}, // lv1
927  {"U", "omega", "k", "ReThetat", "gammaInt"} // lv2
928  });
929  allCon.set(
930  "kRes",
931  {
932  {"U", "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
933  {"U", "omega", "k", "ReThetat", "gammaInt"}, // lv1
934  {"U", "omega", "k", "ReThetat", "gammaInt"} // lv2
935  });
936 
937  allCon.set(
938  "ReThetatRes",
939  {
940  {"U", "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
941  {"U", "omega", "k", "ReThetat", "gammaInt"}, // lv1
942  {"U", "omega", "k", "ReThetat", "gammaInt"} // lv2
943  });
944 
945  allCon.set(
946  "gammaIntRes",
947  {
948  {"U", "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
949  {"U", "omega", "k", "ReThetat", "gammaInt"}, // lv1
950  {"U", "omega", "k", "ReThetat", "gammaInt"} // lv2
951  });
952 #endif
953 
954 #ifdef CompressibleFlow
955  allCon.set(
956  "omegaRes",
957  {
958  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
959  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"}, // lv1
960  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"} // lv2
961  });
962  allCon.set(
963  "kRes",
964  {
965  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
966  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"}, // lv1
967  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"} // lv2
968  });
969 
970  allCon.set(
971  "ReThetatRes",
972  {
973  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
974  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"}, // lv1
975  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"} // lv2
976  });
977 
978  allCon.set(
979  "gammaIntRes",
980  {
981  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
982  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"}, // lv1
983  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"} // lv2
984  });
985 #endif
986 }
987 
988 void DAkOmegaSSTLM::correct()
989 {
990  /*
991  Descroption:
992  Solve the residual equations and update the state. This function will be called
993  by the DASolver. It is needed because we want to control the output frequency
994  of the residual convergence every 100 steps. If using the correct from turbulence
995  it will output residual every step which will be too much of information.
996  */
997 
998  // We set the flag solveTurbState_ to 1 such that in the calcResiduals function
999  // we will solve and update nuTilda
1000  solveTurbState_ = 1;
1001  dictionary dummyOptions;
1002  this->calcResiduals(dummyOptions);
1003  // after it, we reset solveTurbState_ = 0 such that calcResiduals will not
1004  // update nuTilda when calling from the adjoint class, i.e., solveAdjoint from DASolver.
1005  solveTurbState_ = 0;
1006 }
1007 
1008 void DAkOmegaSSTLM::calcResiduals(const dictionary& options)
1009 {
1010  /*
1011  Descroption:
1012  If solveTurbState_ == 1, this function solve and update k and omega, and
1013  is the same as calling turbulence.correct(). If solveTurbState_ == 0,
1014  this function compute residuals for turbulence variables, e.g., nuTildaRes_
1015 
1016  Input:
1017  options.isPC: 1 means computing residuals for preconditioner matrix.
1018  This essentially use the first order scheme for div(phi,nuTilda)
1019 
1020  p_, U_, phi_, etc: State variables in OpenFOAM
1021 
1022  Output:
1023  kRes_/omegaRes_: If solveTurbState_ == 0, update the residual field variable
1024 
1025  k_/omega_: If solveTurbState_ == 1, update them
1026  */
1027 
1028  // Copy and modify based on the "correct" function
1029 
1030  label printToScreen = this->isPrintTime(mesh_.time(), printInterval_);
1031 
1032  word divKScheme = "div(phi,k)";
1033  word divOmegaScheme = "div(phi,omega)";
1034  word divReThetatScheme = "div(phi,ReThetat)";
1035  word divGammaIntScheme = "div(phi,gammaInt)";
1036 
1037  label isPC = 0;
1038 
1039  if (!solveTurbState_)
1040  {
1041  isPC = options.getLabel("isPC");
1042 
1043  if (isPC)
1044  {
1045  divKScheme = "div(pc)";
1046  divOmegaScheme = "div(pc)";
1047  divReThetatScheme = "div(pc)";
1048  divGammaIntScheme = "div(pc)";
1049  }
1050  }
1051  // ************ SST part ********
1052  {
1053  // Note: for compressible flow, the "this->phi()" function divides phi by fvc:interpolate(rho),
1054  // while for the incompresssible "this->phi()" returns phi only
1055  // see src/TurbulenceModels/compressible/compressibleTurbulenceModel.C line 62 to 73
1056  volScalarField::Internal divU(fvc::div(fvc::absolute(phi_ / fvc::interpolate(rho_), U_)));
1057 
1058  tmp<volTensorField> tgradU = fvc::grad(U_);
1059  volScalarField S2(2 * magSqr(symm(tgradU())));
1060  volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
1061  volScalarField::Internal G("kOmegaSSTLM:G", nut_ * GbyNu0);
1062 
1063  if (solveTurbState_)
1064  {
1065  // Update omega and G at the wall
1066  omega_.boundaryFieldRef().updateCoeffs();
1067  }
1068  else
1069  {
1070  // NOTE instead of calling omega_.boundaryFieldRef().updateCoeffs();
1071  // here we call our self-defined boundary conditions
1072  this->correctOmegaBoundaryConditions();
1073  }
1074 
1075  volScalarField CDkOmega(
1076  (scalar(2) * alphaOmega2_) * (fvc::grad(k_) & fvc::grad(omega_)) / omega_);
1077 
1078  volScalarField F1(this->F1(CDkOmega));
1079  volScalarField F23(this->F23());
1080 
1081  {
1082 
1083  volScalarField::Internal gamma(this->gamma(F1));
1084  volScalarField::Internal beta(this->beta(F1));
1085 
1086  // Turbulent frequency equation
1087  tmp<fvScalarMatrix> omegaEqn(
1088  fvm::ddt(phase_, rho_, omega_)
1089  + fvm::div(phaseRhoPhi_, omega_, divOmegaScheme)
1090  - fvm::laplacian(phase_ * rho_ * DomegaEff(F1), omega_)
1091  == phase_() * rho_() * gamma * GbyNu(GbyNu0, F23(), S2())
1092  - fvm::SuSp((2.0 / 3.0) * phase_() * rho_() * gamma * divU, omega_)
1093  - fvm::Sp(phase_() * rho_() * beta * omega_(), omega_)
1094  - fvm::SuSp(
1095  phase_() * rho_() * (F1() - scalar(1)) * CDkOmega() / omega_(),
1096  omega_)
1097  + Qsas(S2(), gamma, beta)
1098  + omegaSource()
1099 
1100  );
1101 
1102  omegaEqn.ref().relax();
1103  omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
1104 
1105  if (solveTurbState_)
1106  {
1107 
1108  // get the solver performance info such as initial
1109  // and final residuals
1110  SolverPerformance<scalar> solverOmega = solve(omegaEqn);
1111  if (printToScreen)
1112  {
1113  Info << "omega Initial residual: " << solverOmega.initialResidual() << endl
1114  << " Final residual: " << solverOmega.finalResidual() << endl;
1115  }
1116 
1117  DAUtility::boundVar(allOptions_, omega_, printToScreen);
1118  }
1119  else
1120  {
1121  // reset the corrected omega near wall cell to its perturbed value
1122  this->setOmegaNearWall();
1123 
1124  // calculate residuals
1125  omegaRes_ = omegaEqn() & omega_;
1126  // need to normalize residuals
1127  normalizeResiduals(omegaRes);
1128  }
1129  }
1130 
1131  // Turbulent kinetic energy equation
1132  tmp<fvScalarMatrix> kEqn(
1133  fvm::ddt(phase_, rho_, k_)
1134  + fvm::div(phaseRhoPhi_, k_, divKScheme)
1135  - fvm::laplacian(phase_ * rho_ * DkEff(F1), k_)
1136  == phase_() * rho_() * Pk(G)
1137  - fvm::SuSp((2.0 / 3.0) * phase_() * rho_() * divU, k_)
1138  - fvm::Sp(phase_() * rho_() * epsilonByk(F1, tgradU()), k_)
1139  + kSource());
1140 
1141  tgradU.clear();
1142 
1143  kEqn.ref().relax();
1144 
1145  if (solveTurbState_)
1146  {
1147 
1148  // get the solver performance info such as initial
1149  // and final residuals
1150  SolverPerformance<scalar> solverK = solve(kEqn);
1151  if (printToScreen)
1152  {
1153  Info << "k Initial residual: " << solverK.initialResidual() << endl
1154  << " Final residual: " << solverK.finalResidual() << endl;
1155  }
1156 
1157  DAUtility::boundVar(allOptions_, k_, printToScreen);
1158 
1159  this->correctNut();
1160  }
1161  else
1162  {
1163  // calculate residuals
1164  kRes_ = kEqn() & k_;
1165  // need to normalize residuals
1166  normalizeResiduals(kRes);
1167  }
1168  }
1169 
1170  // ************ LM part ********
1171  {
1172  // we need to bound before computing residuals
1173  // this will avoid having NaN residuals
1174  DAUtility::boundVar(allOptions_, ReThetat_, printToScreen);
1175 
1176  // Local references
1177  const tmp<volScalarField> tnu = this->nu();
1178  const volScalarField::Internal& nu = tnu()();
1179  const volScalarField::Internal& y = y_();
1180 
1181  // Fields derived from the velocity gradient
1182  tmp<volTensorField> tgradULM = fvc::grad(U_);
1183  const volScalarField::Internal Omega(sqrt(2 * magSqr(skew(tgradULM()()))));
1184  const volScalarField::Internal S(sqrt(2 * magSqr(symm(tgradULM()()))));
1185  const volScalarField::Internal Us(max(mag(U_()), deltaU_));
1186  const volScalarField::Internal dUsds((U_() & (U_() & tgradULM()())) / sqr(Us));
1187  tgradULM.clear();
1188 
1189  const volScalarField::Internal Fthetat(this->Fthetat(Us, Omega, nu));
1190 
1191  {
1192  const volScalarField::Internal t(500 * nu / sqr(Us));
1193  const volScalarField::Internal Pthetat(
1194  phase_() * rho_() * (cThetat_ / t) * (1 - Fthetat));
1195 
1196  // Transition onset momentum-thickness Reynolds number equation
1197  tmp<fvScalarMatrix> ReThetatEqn(
1198  fvm::ddt(phase_, rho_, ReThetat_)
1199  + fvm::div(phaseRhoPhi_, ReThetat_, divReThetatScheme)
1200  - fvm::laplacian(phase_ * rho_ * DReThetatEff(), ReThetat_)
1201  == Pthetat * ReThetat0(Us, dUsds, nu) - fvm::Sp(Pthetat, ReThetat_));
1202 
1203  ReThetatEqn.ref().relax();
1204 
1205  if (solveTurbState_)
1206  {
1207 
1208  // get the solver performance info such as initial
1209  // and final residuals
1210  SolverPerformance<scalar> solverReThetat = solve(ReThetatEqn);
1211  if (printToScreen)
1212  {
1213  Info << "ReThetat Initial residual: " << solverReThetat.initialResidual() << endl
1214  << " Final residual: " << solverReThetat.finalResidual() << endl;
1215  }
1216 
1217  DAUtility::boundVar(allOptions_, ReThetat_, printToScreen);
1218  }
1219  else
1220  {
1221  // calculate residuals
1222  ReThetatRes_ = ReThetatEqn() & ReThetat_;
1223  // need to normalize residuals
1224  normalizeResiduals(ReThetatRes);
1225  }
1226  }
1227 
1228  // we need to bound before computing residuals
1229  // this will avoid having NaN residuals
1230  DAUtility::boundVar(allOptions_, gammaInt_, printToScreen);
1231 
1232  const volScalarField::Internal ReThetac(this->ReThetac());
1233  const volScalarField::Internal Rev(sqr(y) * S / nu);
1234  const volScalarField::Internal RT(k_() / (nu * omega_()));
1235 
1236  {
1237  const volScalarField::Internal Pgamma(
1238  phase_() * rho_()
1239  * ca1_ * Flength(nu) * S * sqrt(gammaInt_() * Fonset(Rev, ReThetac, RT)));
1240 
1241  const volScalarField::Internal Fturb(exp(-pow4(0.25 * RT)));
1242 
1243  const volScalarField::Internal Egamma(
1244  phase_() * rho_() * ca2_ * Omega * Fturb * gammaInt_());
1245 
1246  // Intermittency equation
1247  tmp<fvScalarMatrix> gammaIntEqn(
1248  fvm::ddt(phase_, rho_, gammaInt_)
1249  + fvm::div(phaseRhoPhi_, gammaInt_, divGammaIntScheme)
1250  - fvm::laplacian(phase_ * rho_ * DgammaIntEff(), gammaInt_)
1251  == Pgamma - fvm::Sp(ce1_ * Pgamma, gammaInt_)
1252  + Egamma - fvm::Sp(ce2_ * Egamma, gammaInt_));
1253 
1254  gammaIntEqn.ref().relax();
1255 
1256  if (solveTurbState_)
1257  {
1258 
1259  // get the solver performance info such as initial
1260  // and final residuals
1261  SolverPerformance<scalar> solverGammaInt = solve(gammaIntEqn);
1262  if (printToScreen)
1263  {
1264  Info << "gammaInt Initial residual: " << solverGammaInt.initialResidual() << endl
1265  << " Final residual: " << solverGammaInt.finalResidual() << endl;
1266  }
1267 
1268  DAUtility::boundVar(allOptions_, gammaInt_, printToScreen);
1269 
1270  const volScalarField::Internal Freattach(exp(-pow4(RT / 20.0)));
1271  const volScalarField::Internal gammaSep(
1272  min(2 * max(Rev / (3.235 * ReThetac) - 1, scalar(0)) * Freattach, scalar(2))
1273  * Fthetat);
1274 
1275  gammaIntEff_ = max(gammaInt_(), gammaSep);
1276  }
1277  else
1278  {
1279  // calculate residuals
1280  gammaIntRes_ = gammaIntEqn() & gammaInt_;
1281  // need to normalize residuals
1282  normalizeResiduals(gammaIntRes);
1283  }
1284  }
1285  }
1286 
1287  return;
1288 }
1289 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1290 
1291 } // End namespace Foam
1292 
1293 // ************************************************************************* //
Foam::DAkOmegaSSTLM::omega_
volScalarField & omega_
Definition: DAkOmegaSSTLM.H:201
Foam::DAkOmegaSSTLM::ca2_
dimensionedScalar ca2_
Definition: DAkOmegaSSTLM.H:79
U
U
Definition: pEqnPimpleDyM.H:60
Foam::DAkOmegaSSTLM::omegaNearWall_
scalarList omegaNearWall_
Definition: DAkOmegaSSTLM.H:221
Foam::DAkOmegaSSTLM::ce2_
dimensionedScalar ce2_
Definition: DAkOmegaSSTLM.H:81
forAll
forAll(pseudoP.boundaryField(), patchI)
Definition: solvePseudoPEqn.H:10
Foam::DAkOmegaSSTLM::epsilonByk
tmp< volScalarField::Internal > epsilonByk(const volScalarField &F1, const volTensorField &gradU) const
Definition: DAkOmegaSSTLM.C:366
Foam::DATurbulenceModel::phi_
surfaceScalarField & phi_
face flux
Definition: DATurbulenceModel.H:89
Foam::DAkOmegaSSTLM::gammaIntEff_
volScalarField::Internal & gammaIntEff_
Effective intermittency.
Definition: DAkOmegaSSTLM.H:212
Foam::DAkOmegaSSTLM::setOmegaNearWall
void setOmegaNearWall()
set omegaNearWall_ to near wall omega values
Definition: DAkOmegaSSTLM.C:723
Foam::DAkOmegaSSTLM::y_
const volScalarField & y_
3D wall distance
Definition: DAkOmegaSSTLM.H:215
Foam::DAkOmegaSSTLM::addModelResidualCon
virtual void addModelResidualCon(HashTable< List< List< word >>> &allCon) const
add the model residual connectivity to stateCon
Definition: DAkOmegaSSTLM.C:864
Foam::DAkOmegaSSTLM::ca1_
dimensionedScalar ca1_
Definition: DAkOmegaSSTLM.H:78
Foam::DAkOmegaSSTLM::beta
tmp< volScalarField::Internal > beta(const volScalarField::Internal &F1) const
Definition: DAkOmegaSSTLM.H:122
Foam::DAkOmegaSSTLM::ce1_
dimensionedScalar ce1_
Definition: DAkOmegaSSTLM.H:80
Foam::DAOption
Definition: DAOption.H:29
Foam::DAkOmegaSSTLM::F2
tmp< volScalarField > F2() const
Definition: DAkOmegaSSTLM.C:262
daOption
DAOption daOption(mesh, pyOptions_)
Foam::DAkOmegaSSTLM::DomegaEff
tmp< volScalarField > DomegaEff(const volScalarField &F1) const
Definition: DAkOmegaSSTLM.H:254
Foam::DAkOmegaSSTLM::F23
tmp< volScalarField > F23() const
Definition: DAkOmegaSSTLM.C:284
Foam::DAkOmegaSSTLM::calcResiduals
virtual void calcResiduals(const dictionary &options)
compute the turbulence residuals
Definition: DAkOmegaSSTLM.C:1008
Foam::DAkOmegaSSTLM::a1_
dimensionedScalar a1_
Definition: DAkOmegaSSTLM.H:69
Foam::DAkOmegaSSTLM::Pk
tmp< volScalarField::Internal > Pk(const volScalarField::Internal &G) const
Definition: DAkOmegaSSTLM.C:360
normalizeResiduals
#define normalizeResiduals(resName)
Definition: DAMacroFunctions.H:28
Foam::DAkOmegaSSTLM::omegaRes_
volScalarField omegaRes_
Definition: DAkOmegaSSTLM.H:202
Foam::DAkOmegaSSTLM::F3_
Switch F3_
Definition: DAkOmegaSSTLM.H:73
Foam::DAkOmegaSSTLM::DReThetatEff
tmp< volScalarField > DReThetatEff() const
Definition: DAkOmegaSSTLM.H:264
Foam::DATurbulenceModel::phaseRhoPhi_
surfaceScalarField & phaseRhoPhi_
phase*phi*density field
Definition: DATurbulenceModel.H:95
Foam::DATurbulenceModel::allOptions_
const dictionary & allOptions_
all DAFoam options
Definition: DATurbulenceModel.H:80
Foam::DAkOmegaSSTLM::updateIntermediateVariables
virtual void updateIntermediateVariables()
update any intermediate variables that are dependent on state variables and are used in calcResiduals
Definition: DAkOmegaSSTLM.C:749
Foam::DAkOmegaSSTLM::correctNut
virtual void correctNut()
update nut based on other turbulence variables and update the BCs
Definition: DAkOmegaSSTLM.C:635
Foam::DAkOmegaSSTLM::correctBoundaryConditions
virtual void correctBoundaryConditions()
update turbulence variable boundary values
Definition: DAkOmegaSSTLM.C:657
Foam::DAkOmegaSSTLM::ReThetat_
volScalarField & ReThetat_
Definition: DAkOmegaSSTLM.H:205
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam::DAkOmegaSSTLM::ReThetac
tmp< volScalarField::Internal > ReThetac() const
Definition: DAkOmegaSSTLM.C:396
Foam::DAUtility::boundVar
static void boundVar(const dictionary &allOptions, volScalarField &var, const label printToScreen)
bound a volScalar variable based on parametes defined in DAOption::allOptions_
Definition: DAUtility.C:419
DAkOmegaSSTLM.H
Foam::DAkOmegaSSTLM::c1_
dimensionedScalar c1_
Definition: DAkOmegaSSTLM.H:71
Foam::DATurbulenceModel::U_
volVectorField & U_
velocity
Definition: DATurbulenceModel.H:86
Foam::DATurbulenceModel::nut_
volScalarField & nut_
turbulence viscosity
Definition: DATurbulenceModel.H:83
Foam::DAkOmegaSSTLM::gamma
tmp< volScalarField::Internal > gamma(const volScalarField::Internal &F1) const
Definition: DAkOmegaSSTLM.H:128
Foam::DAkOmegaSSTLM::maxLambdaIter_
label maxLambdaIter_
Definition: DAkOmegaSSTLM.H:85
Foam::DAkOmegaSSTLM::solveTurbState_
label solveTurbState_
whether to solve for turb states
Definition: DAkOmegaSSTLM.H:224
Foam::DAkOmegaSSTLM::gammaInt_
volScalarField & gammaInt_
Definition: DAkOmegaSSTLM.H:207
k
dimensionedScalar & k
Definition: createRefsHeatTransfer.H:6
Foam::DAkOmegaSSTLM::printInterval_
label printInterval_
time step interval to print residual
Definition: DAkOmegaSSTLM.H:227
Foam
Definition: multiFreqScalarFvPatchField.C:144
Foam::DAkOmegaSSTLM::correct
virtual void correct()
solve the residual equations and update the state
Definition: DAkOmegaSSTLM.C:988
Foam::DAkOmegaSSTLM::Fthetat
tmp< volScalarField::Internal > Fthetat(const volScalarField::Internal &Us, const volScalarField::Internal &Omega, const volScalarField::Internal &nu) const
Definition: DAkOmegaSSTLM.C:374
Foam::DAkOmegaSSTLM::Qsas
tmp< fvScalarMatrix > Qsas(const volScalarField::Internal &S2, const volScalarField::Internal &gamma, const volScalarField::Internal &beta) const
Definition: DAkOmegaSSTLM.C:338
Foam::DAkOmegaSSTLM::Fonset
tmp< volScalarField::Internal > Fonset(const volScalarField::Internal &Rev, const volScalarField::Internal &ReThetac, const volScalarField::Internal &RT) const
Definition: DAkOmegaSSTLM.C:481
Foam::DATurbulenceModel::nu
tmp< volScalarField > nu() const
get the nu field
Definition: DATurbulenceModel.C:262
Foam::DAkOmegaSSTLM::deltaU_
const dimensionedScalar deltaU_
Definition: DAkOmegaSSTLM.H:86
Foam::DAkOmegaSSTLM::alphaOmega2_
dimensionedScalar alphaOmega2_
Definition: DAkOmegaSSTLM.H:59
Foam::DATurbulenceModel::correctAlphat
void correctAlphat()
update alphat
Definition: DATurbulenceModel.C:206
Foam::DAkOmegaSSTLM::saveOmegaNearWall
void saveOmegaNearWall()
save near wall omega values to omegaNearWall_
Definition: DAkOmegaSSTLM.C:699
Foam::DAkOmegaSSTLM::omegaSource
tmp< fvScalarMatrix > omegaSource() const
Definition: DAkOmegaSSTLM.C:329
solve
pseudoPEqn solve(solverDictP_)
Foam::DATurbulenceModel
Definition: DATurbulenceModel.H:61
Foam::DAkOmegaSSTLM::F1
tmp< volScalarField > F1(const volScalarField &CDkOmega) const
Definition: DAkOmegaSSTLM.C:351
Foam::DAkOmegaSSTLM::kRes_
volScalarField kRes_
Definition: DAkOmegaSSTLM.H:204
Foam::DAkOmegaSSTLM::DkEff
tmp< volScalarField > DkEff(const volScalarField &F1) const
Definition: DAkOmegaSSTLM.H:247
Foam::DAkOmegaSSTLM::ReThetatRes_
volScalarField ReThetatRes_
Definition: DAkOmegaSSTLM.H:206
Foam::DAkOmegaSSTLM::DAkOmegaSSTLM
DAkOmegaSSTLM(const word modelType, const fvMesh &mesh, const DAOption &daOption)
Definition: DAkOmegaSSTLM.C:41
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFvSource, 0)
Foam::DAkOmegaSSTLM::DgammaIntEff
tmp< volScalarField > DgammaIntEff() const
Definition: DAkOmegaSSTLM.H:272
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFvSource, DAFvSourceActuatorDisk, dictionary)
Foam::DATurbulenceModel::mesh_
const fvMesh & mesh_
fvMesh
Definition: DATurbulenceModel.H:74
lambda
volScalarField & lambda
Definition: createRefsSolidDisplacement.H:7
Foam::DATurbulenceModel::isPrintTime
label isPrintTime(const Time &runTime, const label printInterval) const
Definition: DATurbulenceModel.C:464
Foam::DAkOmegaSSTLM::ReThetat0
tmp< volScalarField::Internal > ReThetat0(const volScalarField::Internal &Us, const volScalarField::Internal &dUsds, const volScalarField::Internal &nu) const
Definition: DAkOmegaSSTLM.C:499
Foam::DAkOmegaSSTLM::gammaIntRes_
volScalarField gammaIntRes_
Definition: DAkOmegaSSTLM.H:208
Foam::DAkOmegaSSTLM::epsilonBykSST
tmp< volScalarField::Internal > epsilonBykSST(const volScalarField &F1, const volTensorField &gradU) const
Definition: DAkOmegaSSTLM.C:313
Foam::DAkOmegaSSTLM::Flength
tmp< volScalarField::Internal > Flength(const volScalarField::Internal &nu) const
Definition: DAkOmegaSSTLM.C:427
Foam::DAkOmegaSSTLM::lambdaErr_
scalar lambdaErr_
Definition: DAkOmegaSSTLM.H:84
Foam::DAkOmegaSSTLM::b1_
dimensionedScalar b1_
Definition: DAkOmegaSSTLM.H:70
Foam::DAkOmegaSSTLM::betaStar_
dimensionedScalar betaStar_
Definition: DAkOmegaSSTLM.H:67
Foam::DAkOmegaSSTLM::kSource
tmp< fvScalarMatrix > kSource() const
Definition: DAkOmegaSSTLM.C:320
Foam::DAkOmegaSSTLM::PkSST
tmp< volScalarField::Internal > PkSST(const volScalarField::Internal &G) const
Definition: DAkOmegaSSTLM.C:307
Foam::DAkOmegaSSTLM::correctStateResidualModelCon
virtual void correctStateResidualModelCon(List< List< word >> &stateCon) const
update the original variable connectivity for the adjoint state residuals in stateCon
Definition: DAkOmegaSSTLM.C:763
Foam::DAkOmegaSSTLM::correctOmegaBoundaryConditions
void correctOmegaBoundaryConditions()
specially treatment to correct epsilon BC
Definition: DAkOmegaSSTLM.C:672
Foam::DAkOmegaSSTLM::k_
volScalarField & k_
Definition: DAkOmegaSSTLM.H:203
Foam::DAkOmegaSSTLM::F3
tmp< volScalarField > F3() const
Definition: DAkOmegaSSTLM.C:274
Foam::DAkOmegaSSTLM::cThetat_
dimensionedScalar cThetat_
Definition: DAkOmegaSSTLM.H:82
Foam::DATurbulenceModel::phase_
volScalarField phase_
phase field
Definition: DATurbulenceModel.H:92
Foam::DAkOmegaSSTLM::GbyNu
tmp< volScalarField::Internal > GbyNu(const volScalarField::Internal &GbyNu0, const volScalarField::Internal &F2, const volScalarField::Internal &S2) const
Definition: DAkOmegaSSTLM.C:296
Foam::DAkOmegaSSTLM::correctModelStates
virtual void correctModelStates(wordList &modelStates) const
update the turbulence state for DAStateInfo::regStates_
Definition: DAkOmegaSSTLM.C:597
Foam::DAkOmegaSSTLM::F1SST
tmp< volScalarField > F1SST(const volScalarField &CDkOmega) const
Definition: DAkOmegaSSTLM.C:243
rho
volScalarField & rho
Definition: createRefsRhoSimpleC.H:8

Generated by   doxygen 1.8.17