DAkOmegaSSTLM.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
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  dimensionedScalar("omegaRes", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
138  zeroGradientFvPatchField<scalar>::typeName),
139  k_(const_cast<volScalarField&>(
140  mesh_.thisDb().lookupObject<volScalarField>("k"))),
141  kRes_(
142  IOobject(
143  "kRes",
144  mesh.time().timeName(),
145  mesh,
146  IOobject::NO_READ,
147  IOobject::NO_WRITE),
148  mesh,
149  dimensionedScalar("kRes", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
150  zeroGradientFvPatchField<scalar>::typeName),
151  ReThetat_(const_cast<volScalarField&>(
152  mesh_.thisDb().lookupObject<volScalarField>("ReThetat"))),
153  ReThetatRes_(
154  IOobject(
155  "ReThetatRes",
156  mesh_.time().timeName(),
157  mesh_,
158  IOobject::NO_READ,
159  IOobject::NO_WRITE),
160  mesh_,
161  dimensionedScalar("ReThetatRes", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
162  zeroGradientFvPatchScalarField::typeName),
163  gammaInt_(const_cast<volScalarField&>(
164  mesh_.thisDb().lookupObject<volScalarField>("gammaInt"))),
165  gammaIntRes_(
166  IOobject(
167  "gammaIntRes",
168  mesh_.time().timeName(),
169  mesh_,
170  IOobject::NO_READ,
171  IOobject::NO_WRITE),
172  mesh_,
173  dimensionedScalar("gammaIntRes", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
174  zeroGradientFvPatchScalarField::typeName),
175  gammaIntEff_(const_cast<volScalarField::Internal&>(
176  mesh_.thisDb().lookupObject<volScalarField::Internal>("gammaIntEff"))),
177  y_(mesh_.thisDb().lookupObject<volScalarField>("yWall")),
178  GPtr_(nullptr),
179  betaFIK_(
180  IOobject(
181  "betaFIK",
182  mesh.time().timeName(),
183  mesh,
184  IOobject::READ_IF_PRESENT,
185  IOobject::AUTO_WRITE),
186  mesh,
187  dimensionedScalar("betaFIK", dimensionSet(0, 0, 0, 0, 0, 0, 0), 1.0),
188  "zeroGradient"),
189  betaFIOmega_(
190  IOobject(
191  "betaFIOmega",
192  mesh.time().timeName(),
193  mesh,
194  IOobject::READ_IF_PRESENT,
195  IOobject::AUTO_WRITE),
196  mesh,
197  dimensionedScalar("betaFIOmega", dimensionSet(0, 0, 0, 0, 0, 0, 0), 1.0),
198  "zeroGradient")
199 {
200 
201  if (turbModelType_ == "incompressible")
202  {
203  omegaRes_.dimensions().reset(dimensionSet(0, 0, -2, 0, 0, 0, 0));
204  kRes_.dimensions().reset(dimensionSet(0, 2, -3, 0, 0, 0, 0));
205  ReThetatRes_.dimensions().reset(dimensionSet(0, 0, -1, 0, 0, 0, 0));
206  gammaIntRes_.dimensions().reset(dimensionSet(0, 0, -1, 0, 0, 0, 0));
207  }
208  else if (turbModelType_ == "compressible")
209  {
210  omegaRes_.dimensions().reset(dimensionSet(1, -3, -2, 0, 0, 0, 0));
211  kRes_.dimensions().reset(dimensionSet(1, -1, -3, 0, 0, 0, 0));
212  ReThetatRes_.dimensions().reset(dimensionSet(1, -3, -1, 0, 0, 0, 0));
213  gammaIntRes_.dimensions().reset(dimensionSet(1, -3, -1, 0, 0, 0, 0));
214  }
215 
216  // calculate the size of omegaWallFunction faces
217  label nWallFaces = 0;
218  forAll(omega_.boundaryField(), patchI)
219  {
220  if (omega_.boundaryField()[patchI].type() == "omegaWallFunction"
221  && omega_.boundaryField()[patchI].size() > 0)
222  {
223  forAll(omega_.boundaryField()[patchI], faceI)
224  {
225  nWallFaces++;
226  }
227  }
228  }
229 
230  // initialize omegaNearWall
231  omegaNearWall_.setSize(nWallFaces);
232 
233  // initialize the G field
234  tmp<volTensorField> tgradU = fvc::grad(U_);
235  volScalarField S2(2 * magSqr(symm(tgradU())));
236  volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
237  GPtr_.reset(new volScalarField::Internal("kOmegaSSTLM:G", nut_ * GbyNu0));
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  return tmp<fvScalarMatrix>(
323  new fvScalarMatrix(
324  k_,
325  dimVolume * this->rhoDimensions() * k_.dimensions() / dimTime));
326 }
327 
328 tmp<fvScalarMatrix> DAkOmegaSSTLM::omegaSource() const
329 {
330  return tmp<fvScalarMatrix>(
331  new fvScalarMatrix(
332  omega_,
333  dimVolume * this->rhoDimensions() * omega_.dimensions() / dimTime));
334 }
335 
336 tmp<fvScalarMatrix> DAkOmegaSSTLM::Qsas(
337  const volScalarField::Internal& S2,
338  const volScalarField::Internal& gamma,
339  const volScalarField::Internal& beta) const
340 {
341  return tmp<fvScalarMatrix>(
342  new fvScalarMatrix(
343  omega_,
344  dimVolume * this->rhoDimensions() * omega_.dimensions() / dimTime));
345 }
346 
347 // SSTLM functions
348 tmp<volScalarField> DAkOmegaSSTLM::F1(
349  const volScalarField& CDkOmega) const
350 {
351  const volScalarField Ry(y_ * sqrt(k_) / this->nu());
352  const volScalarField F3(exp(-pow(Ry / 120.0, 8)));
353 
354  return max(this->F1SST(CDkOmega), F3);
355 }
356 
357 tmp<volScalarField::Internal> DAkOmegaSSTLM::Pk(
358  const volScalarField::Internal& G) const
359 {
360  return gammaIntEff_ * this->PkSST(G);
361 }
362 
363 tmp<volScalarField::Internal> DAkOmegaSSTLM::epsilonByk(
364  const volScalarField& F1,
365  const volTensorField& gradU) const
366 {
367  return min(max(gammaIntEff_, scalar(0.1)), scalar(1))
368  * this->epsilonBykSST(F1, gradU);
369 }
370 
371 tmp<volScalarField::Internal> DAkOmegaSSTLM::Fthetat(
372  const volScalarField::Internal& Us,
373  const volScalarField::Internal& Omega,
374  const volScalarField::Internal& nu) const
375 {
376  const volScalarField::Internal& omega = omega_();
377  const volScalarField::Internal& y = y_();
378 
379  const volScalarField::Internal delta(375 * Omega * nu * ReThetat_() * y / sqr(Us));
380  const volScalarField::Internal ReOmega(sqr(y) * omega / nu);
381  const volScalarField::Internal Fwake(exp(-sqr(ReOmega / 1e5)));
382 
383  return tmp<volScalarField::Internal>(
384  new volScalarField::Internal(
385  IOobject::groupName("Fthetat", U_.group()),
386  min(
387  max(
388  Fwake * exp(-pow4((y / delta))),
389  (1 - sqr((gammaInt_() - 1.0 / ce2_) / (1 - 1.0 / ce2_)))),
390  scalar(1))));
391 }
392 
393 tmp<volScalarField::Internal> DAkOmegaSSTLM::ReThetac() const
394 {
395 
396  tmp<volScalarField::Internal> tReThetac(
397  new volScalarField::Internal(
398  IOobject(
399  IOobject::groupName("ReThetac", U_.group()),
400  mesh_.time().timeName(),
401  mesh_),
402  mesh_,
403  dimless));
404  volScalarField::Internal& ReThetac = tReThetac.ref();
405 
406  forAll(ReThetac, celli)
407  {
408  const scalar ReThetat = ReThetat_[celli];
409 
410  ReThetac[celli] =
411  ReThetat <= 1870
412  ? scalar(ReThetat
413  - 396.035e-2
414  + 120.656e-4 * ReThetat
415  - 868.230e-6 * sqr(ReThetat)
416  + 696.506e-9 * pow3(ReThetat)
417  - 174.105e-12 * pow4(ReThetat))
418  : scalar(ReThetat - 593.11 - 0.482 * (ReThetat - 1870));
419  }
420 
421  return tReThetac;
422 }
423 
424 tmp<volScalarField::Internal> DAkOmegaSSTLM::Flength(
425  const volScalarField::Internal& nu) const
426 {
427 
428  tmp<volScalarField::Internal> tFlength(
429  new volScalarField::Internal(
430  IOobject(
431  IOobject::groupName("Flength", U_.group()),
432  mesh_.time().timeName(),
433  mesh_),
434  mesh_,
435  dimless));
436  volScalarField::Internal& Flength = tFlength.ref();
437 
438  const volScalarField::Internal& omega = omega_();
439  const volScalarField::Internal& y = y_();
440 
441  forAll(ReThetat_, celli)
442  {
443  const scalar ReThetat = ReThetat_[celli];
444 
445  if (ReThetat < 400)
446  {
447  Flength[celli] =
448  398.189e-1
449  - 119.270e-4 * ReThetat
450  - 132.567e-6 * sqr(ReThetat);
451  }
452  else if (ReThetat < 596)
453  {
454  Flength[celli] =
455  263.404
456  - 123.939e-2 * ReThetat
457  + 194.548e-5 * sqr(ReThetat)
458  - 101.695e-8 * pow3(ReThetat);
459  }
460  else if (ReThetat < 1200)
461  {
462  Flength[celli] = 0.5 - 3e-4 * (ReThetat - 596);
463  }
464  else
465  {
466  Flength[celli] = 0.3188;
467  }
468 
469  const scalar Fsublayer =
470  exp(-sqr(sqr(y[celli]) * omega[celli] / (200 * nu[celli])));
471 
472  Flength[celli] = Flength[celli] * (1 - Fsublayer) + 40 * Fsublayer;
473  }
474 
475  return tFlength;
476 }
477 
478 tmp<volScalarField::Internal> DAkOmegaSSTLM::Fonset(
479  const volScalarField::Internal& Rev,
480  const volScalarField::Internal& ReThetac,
481  const volScalarField::Internal& RT) const
482 {
483  const volScalarField::Internal Fonset1(Rev / (2.193 * ReThetac));
484 
485  const volScalarField::Internal Fonset2(
486  min(max(Fonset1, pow4(Fonset1)), scalar(2)));
487 
488  const volScalarField::Internal Fonset3(max(1 - pow3(RT / 2.5), scalar(0)));
489 
490  return tmp<volScalarField::Internal>(
491  new volScalarField::Internal(
492  IOobject::groupName("Fonset", U_.group()),
493  max(Fonset2 - Fonset3, scalar(0))));
494 }
495 
496 tmp<volScalarField::Internal> DAkOmegaSSTLM::ReThetat0(
497  const volScalarField::Internal& Us,
498  const volScalarField::Internal& dUsds,
499  const volScalarField::Internal& nu) const
500 {
501 
502  tmp<volScalarField::Internal> tReThetat0(
503  new volScalarField::Internal(
504  IOobject(
505  IOobject::groupName("ReThetat0", U_.group()),
506  mesh_.time().timeName(),
507  mesh_),
508  mesh_,
509  dimless));
510  volScalarField::Internal& ReThetat0 = tReThetat0.ref();
511 
512  const volScalarField& k = k_;
513 
514  label maxIter = 0;
515 
516  forAll(ReThetat0, celli)
517  {
518  const scalar Tu(
519  max(100 * sqrt((2.0 / 3.0) * k[celli]) / Us[celli], scalar(0.027)));
520 
521  // Initialize lambda to zero.
522  // If lambda were cached between time-steps convergence would be faster
523  // starting from the previous time-step value.
524  scalar lambda = 0;
525 
526  scalar lambdaErr;
527  scalar thetat;
528  label iter = 0;
529 
530  do
531  {
532  // Previous iteration lambda for convergence test
533  const scalar lambda0 = lambda;
534 
535  if (Tu <= 1.3)
536  {
537  const scalar Flambda =
538  dUsds[celli] <= 0
539  ? scalar(1
540  - (-12.986 * lambda
541  - 123.66 * sqr(lambda)
542  - 405.689 * pow3(lambda))
543  * exp(-pow(Tu / 1.5, 1.5)))
544  : scalar(1
545  + 0.275 * (1 - exp(-35 * lambda))
546  * exp(-Tu / 0.5));
547 
548  thetat =
549  (1173.51 - 589.428 * Tu + 0.2196 / sqr(Tu))
550  * Flambda * nu[celli]
551  / Us[celli];
552  }
553  else
554  {
555  const scalar Flambda =
556  dUsds[celli] <= 0
557  ? scalar(1
558  - (-12.986 * lambda
559  - 123.66 * sqr(lambda)
560  - 405.689 * pow3(lambda))
561  * exp(-pow(Tu / 1.5, 1.5)))
562  : scalar(1
563  + 0.275 * (1 - exp(-35 * lambda))
564  * exp(-2 * Tu));
565 
566  thetat =
567  331.50 * pow((Tu - 0.5658), -0.671)
568  * Flambda * nu[celli] / Us[celli];
569  }
570 
571  lambda = sqr(thetat) / nu[celli] * dUsds[celli];
572  lambda = max(min(lambda, 0.1), -0.1);
573 
574  lambdaErr = mag(lambda - lambda0);
575 
576  maxIter = max(maxIter, ++iter);
577 
578  } while (lambdaErr > lambdaErr_);
579 
580  ReThetat0[celli] = max(thetat * Us[celli] / nu[celli], scalar(20));
581  }
582 
583  if (maxIter > maxLambdaIter_)
584  {
585  WarningInFunction
586  << "Number of lambda iterations exceeds maxLambdaIter("
587  << maxLambdaIter_ << ')' << endl;
588  }
589 
590  return tReThetat0;
591 }
592 
593 // Augmented functions
594 void DAkOmegaSSTLM::correctModelStates(wordList& modelStates) const
595 {
596  /*
597  Description:
598  Update the name in modelStates based on the selected physical model at runtime
599 
600  Example:
601  In DAStateInfo, if the modelStates reads:
602 
603  modelStates = {"nut"}
604 
605  then for the SA model, calling correctModelStates(modelStates) will give:
606 
607  modelStates={"nuTilda"}
608 
609  while calling correctModelStates(modelStates) for the SST model will give
610 
611  modelStates={"k","omega"}
612 
613  We don't udpate the names for the radiation model because users are
614  supposed to set modelStates={"G"}
615  */
616 
617  // For SST model, we need to replace nut with k, omega
618 
619  forAll(modelStates, idxI)
620  {
621  word stateName = modelStates[idxI];
622  if (stateName == "nut")
623  {
624  modelStates[idxI] = "omega";
625  modelStates.append("k");
626  modelStates.append("ReThetat");
627  modelStates.append("gammaInt");
628  }
629  }
630 }
631 
632 void DAkOmegaSSTLM::correctNut()
633 {
634  /*
635  Description:
636  Update nut based on other turbulence variables and update the BCs
637  Also update alphat if is present
638  */
639 
640  const volVectorField U = mesh_.thisDb().lookupObject<volVectorField>("U");
641  tmp<volTensorField> tgradU = fvc::grad(U);
642  volScalarField S2(2 * magSqr(symm(tgradU())));
643 
644  nut_ = a1_ * k_ / max(a1_ * omega_, b1_ * F23() * sqrt(S2));
645 
646  nut_.correctBoundaryConditions(); // nutkWallFunction: update wall face nut based on k
647 
648  // this is basically BasicTurbulenceModel::correctNut();
649  this->correctAlphat();
650 
651  return;
652 }
653 
654 void DAkOmegaSSTLM::correctBoundaryConditions()
655 {
656  /*
657  Description:
658  Update turbulence variable boundary values
659  */
660 
661  // correct the BCs for the perturbed fields
662  // kqWallFunction is a zero-gradient BC
663  k_.correctBoundaryConditions();
664 
665  ReThetat_.correctBoundaryConditions();
666  gammaInt_.correctBoundaryConditions();
667 }
668 
669 void DAkOmegaSSTLM::correctOmegaBoundaryConditions()
670 {
671  /*
672  Description:
673  this is a special treatment for omega BC because we cant directly call omega.
674  correctBoundaryConditions() because it will modify the internal omega and G that
675  are right next to walls. This will mess up adjoint Jacobians
676  To solve this issue,
677  1. we store the near wall omega before calling omega.correctBoundaryConditions()
678  2. call omega.correctBoundaryConditions()
679  3. Assign the stored near wall omega back
680  4. Apply a zeroGradient BC for omega at the wall patches
681  *********** NOTE *************
682  this treatment will obviously downgrade the accuracy of adjoint derivative since it is
683  not 100% consistent with what is used for the flow solver; however, our observation
684  shows that the impact is not very large.
685  */
686 
687  // save the perturbed omega at the wall
688  this->saveOmegaNearWall();
689  // correct omega boundary conditions, this includes updating wall face and near wall omega values,
690  // updating the inter-proc BCs
691  omega_.correctBoundaryConditions();
692  // reset the corrected omega near wall cell to its perturbed value
693  this->setOmegaNearWall();
694 }
695 
696 void DAkOmegaSSTLM::saveOmegaNearWall()
697 {
698  /*
699  Description:
700  Save the current near wall omega values to omegaNearWall_
701  */
702  label counterI = 0;
703  forAll(omega_.boundaryField(), patchI)
704  {
705  if (omega_.boundaryField()[patchI].type() == "omegaWallFunction"
706  and omega_.boundaryField()[patchI].size() > 0)
707  {
708  const UList<label>& faceCells = mesh_.boundaryMesh()[patchI].faceCells();
709  forAll(faceCells, faceI)
710  {
711  //Info<<"Near Wall cellI: "<<faceCells[faceI]<<endl;
712  omegaNearWall_[counterI] = omega_[faceCells[faceI]];
713  counterI++;
714  }
715  }
716  }
717  return;
718 }
719 
720 void DAkOmegaSSTLM::setOmegaNearWall()
721 {
722  /*
723  Description:
724  Set the current near wall omega values from omegaNearWall_
725  Here we also apply a zeroGradient BC to the wall faces
726  */
727  label counterI = 0;
728  forAll(omega_.boundaryField(), patchI)
729  {
730  if (omega_.boundaryField()[patchI].type() == "omegaWallFunction"
731  && omega_.boundaryField()[patchI].size() > 0)
732  {
733  const UList<label>& faceCells = mesh_.boundaryMesh()[patchI].faceCells();
734  forAll(faceCells, faceI)
735  {
736  omega_[faceCells[faceI]] = omegaNearWall_[counterI];
737  // zeroGradient BC
738  omega_.boundaryFieldRef()[patchI][faceI] = omega_[faceCells[faceI]];
739  counterI++;
740  }
741  }
742  }
743  return;
744 }
745 
746 void DAkOmegaSSTLM::updateIntermediateVariables()
747 {
748  /*
749  Description:
750  Update nut based on nuTilda. Note: we need to update nut and its BC since we
751  may have perturbed other turbulence vars that affect the nut values
752  */
753 
754  this->correctNut();
755 
756  // for SSTLM also update gammaIntEff_
757  // NOTE: this is not implemented yet!!!
758 }
759 
760 void DAkOmegaSSTLM::correctStateResidualModelCon(List<List<word>>& stateCon) const
761 {
762  /*
763  Description:
764  Update the original variable connectivity for the adjoint state
765  residuals in stateCon. Basically, we modify/add state variables based on the
766  original model variables defined in stateCon.
767 
768  Input:
769 
770  stateResCon: the connectivity levels for a state residual, defined in Foam::DAJacCon
771 
772  Example:
773  If stateCon reads:
774  stateCon=
775  {
776  {"U", "p", "nut"},
777  {"p"}
778  }
779 
780  For the SA turbulence model, calling this function for will get a new stateCon
781  stateCon=
782  {
783  {"U", "p", "nuTilda"},
784  {"p"}
785  }
786 
787  For the SST turbulence model, calling this function will give
788  stateCon=
789  {
790  {"U", "p", "k", "omega"},
791  {"p", "U"}
792  }
793  ***NOTE***: we add a extra level of U connectivity because nut is
794  related to grad(U), k, and omega in SST!
795  */
796 
797  label stateConSize = stateCon.size();
798  forAll(stateCon, idxI)
799  {
800  label addUCon = 0;
801  forAll(stateCon[idxI], idxJ)
802  {
803  word conStateName = stateCon[idxI][idxJ];
804  if (conStateName == "nut")
805  {
806  stateCon[idxI][idxJ] = "omega"; // replace nut with omega
807  stateCon[idxI].append("k"); // also add k for that level
808  stateCon[idxI].append("ReThetat");
809  stateCon[idxI].append("gammaInt");
810  addUCon = 1;
811  }
812  }
813  // add U for the current level and level+1 if it is not there yet
814  label isU;
815  if (addUCon == 1)
816  {
817  // first add U for the current level
818  isU = 0;
819  forAll(stateCon[idxI], idxJ)
820  {
821  word conStateName = stateCon[idxI][idxJ];
822  if (conStateName == "U")
823  {
824  isU = 1;
825  }
826  }
827  if (!isU)
828  {
829  stateCon[idxI].append("U");
830  }
831 
832  // now add U for level+1 if idxI is not the largest level
833  // if idxI is already the largest level, we have a problem
834  if (idxI != stateConSize - 1)
835  {
836  isU = 0;
837  forAll(stateCon[idxI + 1], idxJ)
838  {
839  word conStateName = stateCon[idxI + 1][idxJ];
840  if (conStateName == "U")
841  {
842  isU = 1;
843  }
844  }
845  if (!isU)
846  {
847  stateCon[idxI + 1].append("U");
848  }
849  }
850  else
851  {
852  FatalErrorIn(
853  "In DAStateInfo, nut shows in the largest connectivity level! "
854  "This is not supported!")
855  << exit(FatalError);
856  }
857  }
858  }
859 }
860 
861 void DAkOmegaSSTLM::addModelResidualCon(HashTable<List<List<word>>>& allCon) const
862 {
863  /*
864  Description:
865  Add the connectivity levels for all physical model residuals to allCon
866 
867  Input:
868  allCon: the connectivity levels for all state residual, defined in DAJacCon
869 
870  Example:
871  If stateCon reads:
872  allCon=
873  {
874  "URes":
875  {
876  {"U", "p", "nut"},
877  {"p"}
878  }
879  }
880 
881  For the SA turbulence model, calling this function for will get a new stateCon,
882  something like this:
883  allCon=
884  {
885  "URes":
886  {
887  {"U", "p", "nuTilda"},
888  {"p"}
889  },
890  "nuTildaRes":
891  {
892  {"U", "phi", "nuTilda"},
893  {"U"}
894  }
895  }
896 
897  */
898 
899  word pName;
900 
901  if (mesh_.thisDb().foundObject<volScalarField>("p"))
902  {
903  pName = "p";
904  }
905  else if (mesh_.thisDb().foundObject<volScalarField>("p_rgh"))
906  {
907  pName = "p_rgh";
908  }
909  else
910  {
911  FatalErrorIn(
912  "Neither p nor p_rgh was found in mesh.thisDb()!"
913  "addModelResidualCon failed to setup turbulence residuals!")
914  << exit(FatalError);
915  }
916 
917  // NOTE: for compressible flow, it depends on rho so we need to add T and p
918  if (turbModelType_ == "incompressible")
919  {
920  allCon.set(
921  "omegaRes",
922  {
923  {"U", "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
924  {"U", "omega", "k", "ReThetat", "gammaInt"}, // lv1
925  {"U", "omega", "k", "ReThetat", "gammaInt"} // lv2
926  });
927  allCon.set(
928  "kRes",
929  {
930  {"U", "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
931  {"U", "omega", "k", "ReThetat", "gammaInt"}, // lv1
932  {"U", "omega", "k", "ReThetat", "gammaInt"} // lv2
933  });
934 
935  allCon.set(
936  "ReThetatRes",
937  {
938  {"U", "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
939  {"U", "omega", "k", "ReThetat", "gammaInt"}, // lv1
940  {"U", "omega", "k", "ReThetat", "gammaInt"} // lv2
941  });
942 
943  allCon.set(
944  "gammaIntRes",
945  {
946  {"U", "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
947  {"U", "omega", "k", "ReThetat", "gammaInt"}, // lv1
948  {"U", "omega", "k", "ReThetat", "gammaInt"} // lv2
949  });
950  }
951  else if (turbModelType_ == "compressible")
952  {
953  allCon.set(
954  "omegaRes",
955  {
956  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
957  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"}, // lv1
958  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"} // lv2
959  });
960  allCon.set(
961  "kRes",
962  {
963  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
964  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"}, // lv1
965  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"} // lv2
966  });
967 
968  allCon.set(
969  "ReThetatRes",
970  {
971  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
972  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"}, // lv1
973  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"} // lv2
974  });
975 
976  allCon.set(
977  "gammaIntRes",
978  {
979  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt", "phi"}, // lv0
980  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"}, // lv1
981  {"U", "T", pName, "omega", "k", "ReThetat", "gammaInt"} // lv2
982  });
983  }
984 }
985 
986 void DAkOmegaSSTLM::correct(label printToScreen)
987 {
988  /*
989  Descroption:
990  Solve the residual equations and update the state. This function will be called
991  by the DASolver. It is needed because we want to control the output frequency
992  of the residual convergence every 100 steps. If using the correct from turbulence
993  it will output residual every step which will be too much of information.
994  */
995 
996  // We set the flag solveTurbState_ to 1 such that in the calcResiduals function
997  // we will solve and update nuTilda
998  solveTurbState_ = 1;
999  dictionary dummyOptions;
1000  dummyOptions.set("printToScreen", printToScreen);
1001  this->calcResiduals(dummyOptions);
1002  // after it, we reset solveTurbState_ = 0 such that calcResiduals will not
1003  // update nuTilda when calling from the adjoint class, i.e., solveAdjoint from DASolver.
1004  solveTurbState_ = 0;
1005 }
1006 
1007 void DAkOmegaSSTLM::calcResiduals(const dictionary& options)
1008 {
1009  /*
1010  Descroption:
1011  If solveTurbState_ == 1, this function solve and update k and omega, and
1012  is the same as calling turbulence.correct(). If solveTurbState_ == 0,
1013  this function compute residuals for turbulence variables, e.g., nuTildaRes_
1014 
1015  Input:
1016  options.isPC: 1 means computing residuals for preconditioner matrix.
1017  This essentially use the first order scheme for div(phi,nuTilda)
1018 
1019  p_, U_, phi_, etc: State variables in OpenFOAM
1020 
1021  Output:
1022  kRes_/omegaRes_: If solveTurbState_ == 0, update the residual field variable
1023 
1024  k_/omega_: If solveTurbState_ == 1, update them
1025  */
1026 
1027  // Copy and modify based on the "correct" function
1028 
1029  label printToScreen = options.lookupOrDefault("printToScreen", 0);
1030 
1031  word divKScheme = "div(phi,k)";
1032  word divOmegaScheme = "div(phi,omega)";
1033  word divReThetatScheme = "div(phi,ReThetat)";
1034  word divGammaIntScheme = "div(phi,gammaInt)";
1035 
1036  volScalarField rho = this->rho();
1037 
1038  label isPC = 0;
1039 
1040  if (!solveTurbState_)
1041  {
1042  isPC = options.getLabel("isPC");
1043 
1044  if (isPC)
1045  {
1046  divKScheme = "div(pc)";
1047  divOmegaScheme = "div(pc)";
1048  divReThetatScheme = "div(pc)";
1049  divGammaIntScheme = "div(pc)";
1050  }
1051  }
1052  // ************ SST part ********
1053  {
1054  // Note: for compressible flow, the "this->phi()" function divides phi by fvc:interpolate(rho),
1055  // while for the incompresssible "this->phi()" returns phi only
1056  // see src/TurbulenceModels/compressible/compressibleTurbulenceModel.C line 62 to 73
1057  volScalarField::Internal divU(fvc::div(fvc::absolute(phi_ / fvc::interpolate(rho), U_)));
1058 
1059  tmp<volTensorField> tgradU = fvc::grad(U_);
1060  volScalarField S2(2 * magSqr(symm(tgradU())));
1061  volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
1062  volScalarField::Internal& G = const_cast<volScalarField::Internal&>(GPtr_());
1063  G = nut_() * GbyNu0;
1064 
1065  if (solveTurbState_)
1066  {
1067  // Update omega and G at the wall
1068  omega_.boundaryFieldRef().updateCoeffs();
1069  }
1070  else
1071  {
1072  // NOTE instead of calling omega_.boundaryFieldRef().updateCoeffs();
1073  // here we call our self-defined boundary conditions
1074  this->correctOmegaBoundaryConditions();
1075  }
1076 
1077  volScalarField CDkOmega(
1078  (scalar(2) * alphaOmega2_) * (fvc::grad(k_) & fvc::grad(omega_)) / omega_);
1079 
1080  volScalarField F1(this->F1(CDkOmega));
1081  volScalarField F23(this->F23());
1082 
1083  {
1084 
1085  volScalarField::Internal gamma(this->gamma(F1));
1086  volScalarField::Internal beta(this->beta(F1));
1087 
1088  // Turbulent frequency equation
1089  tmp<fvScalarMatrix> omegaEqn(
1090  fvm::ddt(phase_, rho, omega_)
1091  + fvm::div(phaseRhoPhi_, omega_, divOmegaScheme)
1092  - fvm::laplacian(phase_ * rho * DomegaEff(F1), omega_)
1093  == phase_() * rho() * gamma * GbyNu(GbyNu0, F23(), S2()) * betaFIOmega_()
1094  - fvm::SuSp((2.0 / 3.0) * phase_() * rho() * gamma * divU, omega_)
1095  - fvm::Sp(phase_() * rho() * beta * omega_(), omega_)
1096  - fvm::SuSp(
1097  phase_() * rho() * (F1() - scalar(1)) * CDkOmega() / omega_(),
1098  omega_)
1099  + Qsas(S2(), gamma, beta)
1100  + omegaSource()
1101 
1102  );
1103 
1104  omegaEqn.ref().relax();
1105  omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
1106 
1107  if (solveTurbState_)
1108  {
1109  // get the solver performance info such as initial
1110  // and final residuals
1111  SolverPerformance<scalar> solverOmega = solve(omegaEqn);
1112  DAUtility::primalResidualControl(solverOmega, printToScreen, "omega", daGlobalVar_.primalMaxRes);
1113 
1114  DAUtility::boundVar(allOptions_, omega_, printToScreen);
1115  }
1116  else
1117  {
1118  // reset the corrected omega near wall cell to its perturbed value
1119  this->setOmegaNearWall();
1120 
1121  // calculate residuals
1122  omegaRes_ = omegaEqn() & omega_;
1123  // need to normalize residuals
1124  normalizeResiduals(omegaRes);
1125  }
1126  }
1127 
1128  // Turbulent kinetic energy equation
1129  tmp<fvScalarMatrix> kEqn(
1130  fvm::ddt(phase_, rho, k_)
1131  + fvm::div(phaseRhoPhi_, k_, divKScheme)
1132  - fvm::laplacian(phase_ * rho * DkEff(F1), k_)
1133  == phase_() * rho() * Pk(G) * betaFIK_()
1134  - fvm::SuSp((2.0 / 3.0) * phase_() * rho() * divU, k_)
1135  - fvm::Sp(phase_() * rho() * epsilonByk(F1, tgradU()), k_)
1136  + kSource());
1137 
1138  tgradU.clear();
1139 
1140  kEqn.ref().relax();
1141 
1142  if (solveTurbState_)
1143  {
1144 
1145  // get the solver performance info such as initial
1146  // and final residuals
1147  SolverPerformance<scalar> solverK = solve(kEqn);
1148  DAUtility::primalResidualControl(solverK, printToScreen, "k", daGlobalVar_.primalMaxRes);
1149 
1150  DAUtility::boundVar(allOptions_, k_, printToScreen);
1151 
1152  this->correctNut();
1153  }
1154  else
1155  {
1156  // calculate residuals
1157  kRes_ = kEqn() & k_;
1158  // need to normalize residuals
1159  normalizeResiduals(kRes);
1160  }
1161  }
1162 
1163  // ************ LM part ********
1164  {
1165  // we need to bound before computing residuals
1166  // this will avoid having NaN residuals
1167  DAUtility::boundVar(allOptions_, ReThetat_, printToScreen);
1168 
1169  // Local references
1170  const tmp<volScalarField> tnu = this->nu();
1171  const volScalarField::Internal& nu = tnu()();
1172  const volScalarField::Internal& y = y_();
1173 
1174  // Fields derived from the velocity gradient
1175  tmp<volTensorField> tgradULM = fvc::grad(U_);
1176  const volScalarField::Internal Omega(sqrt(2 * magSqr(skew(tgradULM()()))));
1177  const volScalarField::Internal S(sqrt(2 * magSqr(symm(tgradULM()()))));
1178  const volScalarField::Internal Us(max(mag(U_()), deltaU_));
1179  const volScalarField::Internal dUsds((U_() & (U_() & tgradULM()())) / sqr(Us));
1180  tgradULM.clear();
1181 
1182  const volScalarField::Internal Fthetat(this->Fthetat(Us, Omega, nu));
1183 
1184  {
1185  const volScalarField::Internal t(500 * nu / sqr(Us));
1186  const volScalarField::Internal Pthetat(
1187  phase_() * rho() * (cThetat_ / t) * (1 - Fthetat));
1188 
1189  // Transition onset momentum-thickness Reynolds number equation
1190  tmp<fvScalarMatrix> ReThetatEqn(
1191  fvm::ddt(phase_, rho, ReThetat_)
1192  + fvm::div(phaseRhoPhi_, ReThetat_, divReThetatScheme)
1193  - fvm::laplacian(phase_ * rho * DReThetatEff(), ReThetat_)
1194  == Pthetat * ReThetat0(Us, dUsds, nu) - fvm::Sp(Pthetat, ReThetat_));
1195 
1196  ReThetatEqn.ref().relax();
1197 
1198  if (solveTurbState_)
1199  {
1200 
1201  // get the solver performance info such as initial
1202  // and final residuals
1203  SolverPerformance<scalar> solverReThetat = solve(ReThetatEqn);
1204  DAUtility::primalResidualControl(solverReThetat, printToScreen, "ReThetat", daGlobalVar_.primalMaxRes);
1205 
1206  DAUtility::boundVar(allOptions_, ReThetat_, printToScreen);
1207  }
1208  else
1209  {
1210  // calculate residuals
1211  ReThetatRes_ = ReThetatEqn() & ReThetat_;
1212  // need to normalize residuals
1213  normalizeResiduals(ReThetatRes);
1214  }
1215  }
1216 
1217  // we need to bound before computing residuals
1218  // this will avoid having NaN residuals
1219  DAUtility::boundVar(allOptions_, gammaInt_, printToScreen);
1220 
1221  const volScalarField::Internal ReThetac(this->ReThetac());
1222  const volScalarField::Internal Rev(sqr(y) * S / nu);
1223  const volScalarField::Internal RT(k_() / (nu * omega_()));
1224 
1225  {
1226  const volScalarField::Internal Pgamma(
1227  phase_() * rho()
1228  * ca1_ * Flength(nu) * S * sqrt(gammaInt_() * Fonset(Rev, ReThetac, RT)));
1229 
1230  const volScalarField::Internal Fturb(exp(-pow4(0.25 * RT)));
1231 
1232  const volScalarField::Internal Egamma(
1233  phase_() * rho() * ca2_ * Omega * Fturb * gammaInt_());
1234 
1235  // Intermittency equation
1236  tmp<fvScalarMatrix> gammaIntEqn(
1237  fvm::ddt(phase_, rho, gammaInt_)
1238  + fvm::div(phaseRhoPhi_, gammaInt_, divGammaIntScheme)
1239  - fvm::laplacian(phase_ * rho * DgammaIntEff(), gammaInt_)
1240  == Pgamma - fvm::Sp(ce1_ * Pgamma, gammaInt_)
1241  + Egamma - fvm::Sp(ce2_ * Egamma, gammaInt_));
1242 
1243  gammaIntEqn.ref().relax();
1244 
1245  if (solveTurbState_)
1246  {
1247 
1248  // get the solver performance info such as initial
1249  // and final residuals
1250  SolverPerformance<scalar> solverGammaInt = solve(gammaIntEqn);
1251  DAUtility::primalResidualControl(solverGammaInt, printToScreen, "gammaInt", daGlobalVar_.primalMaxRes);
1252 
1253  DAUtility::boundVar(allOptions_, gammaInt_, printToScreen);
1254 
1255  const volScalarField::Internal Freattach(exp(-pow4(RT / 20.0)));
1256  const volScalarField::Internal gammaSep(
1257  min(2 * max(Rev / (3.235 * ReThetac) - 1, scalar(0)) * Freattach, scalar(2))
1258  * Fthetat);
1259 
1260  gammaIntEff_ = max(gammaInt_(), gammaSep);
1261  }
1262  else
1263  {
1264  // calculate residuals
1265  gammaIntRes_ = gammaIntEqn() & gammaInt_;
1266  // need to normalize residuals
1267  normalizeResiduals(gammaIntRes);
1268  }
1269  }
1270  }
1271 
1272  return;
1273 }
1274 
1275 void DAkOmegaSSTLM::getFvMatrixFields(
1276  const word varName,
1277  scalarField& diag,
1278  scalarField& upper,
1279  scalarField& lower)
1280 {
1281  /*
1282  Description:
1283  return the diag(), upper(), and lower() scalarFields from the turbulence model's fvMatrix
1284  this will be use to compute the preconditioner matrix
1285  */
1286 
1287  // NOTE: there is a bug in this func. The calcPCMatWithFvMatrix returns an error
1288  // saying try to insert a value that is out of range...
1289  Info << "Warning!!!!!! this child class is not implemented!" << endl;
1290 
1291  /*
1292  if (varName != "k" && varName != "omega" && varName != "ReThetat" && varName != "gammaInt")
1293  {
1294  FatalErrorIn(
1295  "varName not valid. It has to be k, omega, ReThetat, or gamma")
1296  << exit(FatalError);
1297  }
1298 
1299  volScalarField rho = this->rho();
1300 
1301  // Note: for compressible flow, the "this->phi()" function divides phi by fvc:interpolate(rho),
1302  // while for the incompresssible "this->phi()" returns phi only
1303  // see src/TurbulenceModels/compressible/compressibleTurbulenceModel.C line 62 to 73
1304  volScalarField::Internal divU(fvc::div(fvc::absolute(phi_ / fvc::interpolate(rho), U_)));
1305 
1306  // NOTE instead of calling omega_.boundaryFieldRef().updateCoeffs();
1307  // here we call our self-defined boundary conditions
1308  this->correctOmegaBoundaryConditions();
1309 
1310  volScalarField CDkOmega(
1311  (scalar(2) * alphaOmega2_) * (fvc::grad(k_) & fvc::grad(omega_)) / omega_);
1312 
1313  volScalarField F1(this->F1(CDkOmega));
1314  volScalarField F23(this->F23());
1315 
1316  if (varName == "omega" || varName == "k")
1317  {
1318  // Note: for compressible flow, the "this->phi()" function divides phi by fvc:interpolate(rho),
1319  // while for the incompresssible "this->phi()" returns phi only
1320  // see src/TurbulenceModels/compressible/compressibleTurbulenceModel.C line 62 to 73
1321  volScalarField::Internal divU(fvc::div(fvc::absolute(phi_ / fvc::interpolate(rho), U_)));
1322 
1323  tmp<volTensorField> tgradU = fvc::grad(U_);
1324  volScalarField S2(2 * magSqr(symm(tgradU())));
1325  volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
1326  volScalarField::Internal& G = const_cast<volScalarField::Internal&>(GPtr_());
1327  G = nut_() * GbyNu0;
1328 
1329  // NOTE instead of calling omega_.boundaryFieldRef().updateCoeffs();
1330  // here we call our self-defined boundary conditions
1331  this->correctOmegaBoundaryConditions();
1332 
1333  volScalarField CDkOmega(
1334  (scalar(2) * alphaOmega2_) * (fvc::grad(k_) & fvc::grad(omega_)) / omega_);
1335 
1336  volScalarField F1(this->F1(CDkOmega));
1337  volScalarField F23(this->F23());
1338 
1339  if (varName == "omega")
1340  {
1341 
1342  volScalarField::Internal gamma(this->gamma(F1));
1343  volScalarField::Internal beta(this->beta(F1));
1344 
1345  // Turbulent frequency equation
1346  fvScalarMatrix omegaEqn(
1347  fvm::ddt(phase_, rho, omega_)
1348  + fvm::div(phaseRhoPhi_, omega_, "div(pc)")
1349  - fvm::laplacian(phase_ * rho * DomegaEff(F1), omega_)
1350  == phase_() * rho() * gamma * GbyNu(GbyNu0, F23(), S2()) * betaFIOmega_()
1351  - fvm::SuSp((2.0 / 3.0) * phase_() * rho() * gamma * divU, omega_)
1352  - fvm::Sp(phase_() * rho() * beta * omega_(), omega_)
1353  - fvm::SuSp(
1354  phase_() * rho() * (F1() - scalar(1)) * CDkOmega() / omega_(),
1355  omega_)
1356  + Qsas(S2(), gamma, beta)
1357  + omegaSource()
1358 
1359  );
1360 
1361  omegaEqn.relax();
1362 
1363  // reset the corrected omega near wall cell to its perturbed value
1364  this->setOmegaNearWall();
1365 
1366  diag = omegaEqn.D();
1367  upper = omegaEqn.upper();
1368  lower = omegaEqn.lower();
1369  }
1370 
1371  if (varName == "k")
1372  {
1373  // Turbulent kinetic energy equation
1374  fvScalarMatrix kEqn(
1375  fvm::ddt(phase_, rho, k_)
1376  + fvm::div(phaseRhoPhi_, k_, "div(pc)")
1377  - fvm::laplacian(phase_ * rho * DkEff(F1), k_)
1378  == phase_() * rho() * Pk(G)
1379  - fvm::SuSp((2.0 / 3.0) * phase_() * rho() * divU, k_)
1380  - fvm::Sp(phase_() * rho() * epsilonByk(F1, tgradU()), k_)
1381  + kSource());
1382 
1383  tgradU.clear();
1384 
1385  kEqn.relax();
1386 
1387  diag = kEqn.D();
1388  upper = kEqn.upper();
1389  lower = kEqn.lower();
1390  }
1391  }
1392  else if (varName == "gammaInt" || varName == "ReThetat")
1393  {
1394  // we need to bound before computing residuals
1395  // this will avoid having NaN residuals
1396  DAUtility::boundVar(allOptions_, ReThetat_, 0);
1397 
1398  // Local references
1399  const tmp<volScalarField> tnu = this->nu();
1400  const volScalarField::Internal& nu = tnu()();
1401  const volScalarField::Internal& y = y_();
1402 
1403  // Fields derived from the velocity gradient
1404  tmp<volTensorField> tgradULM = fvc::grad(U_);
1405  const volScalarField::Internal Omega(sqrt(2 * magSqr(skew(tgradULM()()))));
1406  const volScalarField::Internal S(sqrt(2 * magSqr(symm(tgradULM()()))));
1407  const volScalarField::Internal Us(max(mag(U_()), deltaU_));
1408  const volScalarField::Internal dUsds((U_() & (U_() & tgradULM()())) / sqr(Us));
1409  tgradULM.clear();
1410 
1411  const volScalarField::Internal Fthetat(this->Fthetat(Us, Omega, nu));
1412 
1413  if (varName == "ReThetat")
1414  {
1415  const volScalarField::Internal t(500 * nu / sqr(Us));
1416  const volScalarField::Internal Pthetat(
1417  phase_() * rho() * (cThetat_ / t) * (1 - Fthetat));
1418 
1419  // Transition onset momentum-thickness Reynolds number equation
1420  fvScalarMatrix ReThetatEqn(
1421  fvm::ddt(phase_, rho, ReThetat_)
1422  + fvm::div(phaseRhoPhi_, ReThetat_, "div(pc)")
1423  - fvm::laplacian(phase_ * rho * DReThetatEff(), ReThetat_)
1424  == Pthetat * ReThetat0(Us, dUsds, nu) - fvm::Sp(Pthetat, ReThetat_));
1425 
1426  ReThetatEqn.relax();
1427 
1428  diag = ReThetatEqn.D();
1429  upper = ReThetatEqn.upper();
1430  lower = ReThetatEqn.lower();
1431  }
1432 
1433  if (varName == "gammaInt")
1434  {
1435  // we need to bound before computing residuals
1436  // this will avoid having NaN residuals
1437  DAUtility::boundVar(allOptions_, gammaInt_, 0);
1438 
1439  const volScalarField::Internal ReThetac(this->ReThetac());
1440  const volScalarField::Internal Rev(sqr(y) * S / nu);
1441  const volScalarField::Internal RT(k_() / (nu * omega_()));
1442 
1443  const volScalarField::Internal Pgamma(
1444  phase_() * rho()
1445  * ca1_ * Flength(nu) * S * sqrt(gammaInt_() * Fonset(Rev, ReThetac, RT)));
1446 
1447  const volScalarField::Internal Fturb(exp(-pow4(0.25 * RT)));
1448 
1449  const volScalarField::Internal Egamma(
1450  phase_() * rho() * ca2_ * Omega * Fturb * gammaInt_());
1451 
1452  // Intermittency equation
1453  fvScalarMatrix gammaIntEqn(
1454  fvm::ddt(phase_, rho, gammaInt_)
1455  + fvm::div(phaseRhoPhi_, gammaInt_, "div(pc)")
1456  - fvm::laplacian(phase_ * rho * DgammaIntEff(), gammaInt_)
1457  == Pgamma - fvm::Sp(ce1_ * Pgamma, gammaInt_)
1458  + Egamma - fvm::Sp(ce2_ * Egamma, gammaInt_));
1459 
1460  gammaIntEqn.relax();
1461 
1462  diag = gammaIntEqn.D();
1463  upper = gammaIntEqn.upper();
1464  lower = gammaIntEqn.lower();
1465  }
1466  }
1467 */
1468 }
1469 
1470 void DAkOmegaSSTLM::getTurbProdOverDestruct(volScalarField& PoD) const
1471 {
1472  /*
1473  Description:
1474  Return the value of the production over destruction term from the turbulence model
1475  */
1476  tmp<volTensorField> tgradU = fvc::grad(U_);
1477  volScalarField S2(2 * magSqr(symm(tgradU())));
1478  volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
1479  volScalarField::Internal& G = const_cast<volScalarField::Internal&>(GPtr_());
1480  G = nut_() * GbyNu0;
1481 
1482  volScalarField rho = this->rho();
1483 
1484  volScalarField CDkOmega(
1485  (scalar(2) * alphaOmega2_) * (fvc::grad(k_) & fvc::grad(omega_)) / omega_);
1486 
1487  volScalarField F1(this->F1(CDkOmega));
1488 
1489  volScalarField::Internal P = phase_() * rho() * Pk(G);
1490  volScalarField::Internal D = phase_() * rho() * epsilonByk(F1, tgradU()) * k_();
1491 
1492  forAll(P, cellI)
1493  {
1494  PoD[cellI] = P[cellI] / (D[cellI] + P[cellI] + 1e-16);
1495  }
1496 }
1497 
1498 void DAkOmegaSSTLM::getTurbConvOverProd(volScalarField& CoP) const
1499 {
1500  /*
1501  Description:
1502  Return the value of the convective over production term from the turbulence model
1503  */
1504 
1505  tmp<volTensorField> tgradU = fvc::grad(U_);
1506  volScalarField S2(2 * magSqr(symm(tgradU())));
1507  volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
1508  volScalarField::Internal& G = const_cast<volScalarField::Internal&>(GPtr_());
1509  G = nut_() * GbyNu0;
1510 
1511  volScalarField rho = this->rho();
1512 
1513  volScalarField CDkOmega(
1514  (scalar(2) * alphaOmega2_) * (fvc::grad(k_) & fvc::grad(omega_)) / omega_);
1515 
1516  volScalarField F1(this->F1(CDkOmega));
1517 
1518  volScalarField::Internal P = phase_() * rho() * Pk(G);
1519  volScalarField C = fvc::div(phaseRhoPhi_, k_);
1520 
1521  forAll(P, cellI)
1522  {
1523  CoP[cellI] = C[cellI] / (P[cellI] + C[cellI] + 1e-16);
1524  }
1525 }
1526 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1527 
1528 } // End namespace Foam
1529 
1530 // ************************************************************************* //
Foam::DAkOmegaSSTLM::omega_
volScalarField & omega_
Definition: DAkOmegaSSTLM.H:201
Foam::DAkOmegaSSTLM::ca2_
dimensionedScalar ca2_
Definition: DAkOmegaSSTLM.H:79
U
volVectorField & U
Definition: createRefsPimpleDyM.H:7
Foam::DAkOmegaSSTLM::omegaNearWall_
scalarList omegaNearWall_
Definition: DAkOmegaSSTLM.H:229
Foam::DAkOmegaSSTLM::ce2_
dimensionedScalar ce2_
Definition: DAkOmegaSSTLM.H:81
Foam::DAkOmegaSSTLM::epsilonByk
tmp< volScalarField::Internal > epsilonByk(const volScalarField &F1, const volTensorField &gradU) const
Definition: DAkOmegaSSTLM.C:363
Foam::DATurbulenceModel::phi_
surfaceScalarField & phi_
face flux
Definition: DATurbulenceModel.H:83
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:720
Foam::DAkOmegaSSTLM::y_
const volScalarField & y_
3D wall distance
Definition: DAkOmegaSSTLM.H:215
solve
nuTilda1Eqn solve(solverDictNuTilda)
Foam::DAkOmegaSSTLM::addModelResidualCon
virtual void addModelResidualCon(HashTable< List< List< word >>> &allCon) const
add the model residual connectivity to stateCon
Definition: DAkOmegaSSTLM.C:861
D
volVectorField & D
Definition: createRefsSolidDisplacement.H:8
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::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
Foam::DAkOmegaSSTLM::F2
tmp< volScalarField > F2() const
Definition: DAkOmegaSSTLM.C:262
Foam::DAkOmegaSSTLM::betaFIOmega_
volScalarField betaFIOmega_
Definition: DAkOmegaSSTLM.H:223
Foam::DAkOmegaSSTLM::DomegaEff
tmp< volScalarField > DomegaEff(const volScalarField &F1) const
Definition: DAkOmegaSSTLM.H:259
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:1007
Foam::DAkOmegaSSTLM::a1_
dimensionedScalar a1_
Definition: DAkOmegaSSTLM.H:69
Foam::DAkOmegaSSTLM::Pk
tmp< volScalarField::Internal > Pk(const volScalarField::Internal &G) const
Definition: DAkOmegaSSTLM.C:357
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:269
Foam::DATurbulenceModel::phaseRhoPhi_
surfaceScalarField & phaseRhoPhi_
phase*phi*density field
Definition: DATurbulenceModel.H:89
Foam::DATurbulenceModel::allOptions_
const dictionary & allOptions_
all DAFoam options
Definition: DATurbulenceModel.H:71
Foam::DATurbulenceModel::rho
tmp< volScalarField > rho() const
get the density field
Definition: DATurbulenceModel.C:294
Foam::DAkOmegaSSTLM::updateIntermediateVariables
virtual void updateIntermediateVariables()
update any intermediate variables that are dependent on state variables and are used in calcResiduals
Definition: DAkOmegaSSTLM.C:746
Foam::DAkOmegaSSTLM::correctNut
virtual void correctNut()
update nut based on other turbulence variables and update the BCs
Definition: DAkOmegaSSTLM.C:632
Foam::DAkOmegaSSTLM::correctBoundaryConditions
virtual void correctBoundaryConditions()
update turbulence variable boundary values
Definition: DAkOmegaSSTLM.C:654
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:393
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:475
DAkOmegaSSTLM.H
Foam::DAkOmegaSSTLM::c1_
dimensionedScalar c1_
Definition: DAkOmegaSSTLM.H:71
Foam::DAUtility::primalResidualControl
static void primalResidualControl(const SolverPerformance< scalar > &solverP, const label printToScreen, const word varName, scalar &primalMaxRes)
control when to print the residual and also compute the maxInitRes
Definition: DAUtility.C:734
Foam::DATurbulenceModel::daGlobalVar_
DAGlobalVar & daGlobalVar_
global variable
Definition: DATurbulenceModel.H:74
Foam::DATurbulenceModel::U_
volVectorField & U_
velocity
Definition: DATurbulenceModel.H:80
Foam::DATurbulenceModel::nut_
volScalarField & nut_
turbulence viscosity
Definition: DATurbulenceModel.H:77
Foam::DAkOmegaSSTLM::GPtr_
autoPtr< volScalarField::Internal > GPtr_
Definition: DAkOmegaSSTLM.H:219
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::DAGlobalVar::primalMaxRes
scalar primalMaxRes
the maximal residual for primal this var will be used in multiple classes
Definition: DAGlobalVar.H:62
Foam::DAkOmegaSSTLM::solveTurbState_
label solveTurbState_
whether to solve for turb states
Definition: DAkOmegaSSTLM.H:232
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
Foam::DAkOmegaSSTLM::gammaInt_
volScalarField & gammaInt_
Definition: DAkOmegaSSTLM.H:207
k
dimensionedScalar & k
Definition: createRefsHeatTransfer.H:6
Foam
Definition: checkGeometry.C:32
Foam::DAkOmegaSSTLM::Fthetat
tmp< volScalarField::Internal > Fthetat(const volScalarField::Internal &Us, const volScalarField::Internal &Omega, const volScalarField::Internal &nu) const
Definition: DAkOmegaSSTLM.C:371
Foam::DAkOmegaSSTLM::correct
virtual void correct(label printToScreen)
solve the residual equations and update the state
Definition: DAkOmegaSSTLM.C:986
Foam::DAkOmegaSSTLM::Qsas
tmp< fvScalarMatrix > Qsas(const volScalarField::Internal &S2, const volScalarField::Internal &gamma, const volScalarField::Internal &beta) const
Definition: DAkOmegaSSTLM.C:336
Foam::DAkOmegaSSTLM::Fonset
tmp< volScalarField::Internal > Fonset(const volScalarField::Internal &Rev, const volScalarField::Internal &ReThetac, const volScalarField::Internal &RT) const
Definition: DAkOmegaSSTLM.C:478
Foam::DATurbulenceModel::rhoDimensions
dimensionSet rhoDimensions() const
return the dimension of rho
Definition: DATurbulenceModel.C:321
Foam::DATurbulenceModel::nu
tmp< volScalarField > nu() const
get the nu field
Definition: DATurbulenceModel.C:267
Foam::DAkOmegaSSTLM::deltaU_
const dimensionedScalar deltaU_
Definition: DAkOmegaSSTLM.H:86
Foam::DAkOmegaSSTLM::alphaOmega2_
dimensionedScalar alphaOmega2_
Definition: DAkOmegaSSTLM.H:59
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
Foam::DAkOmegaSSTLM::getTurbProdOverDestruct
virtual void getTurbProdOverDestruct(volScalarField &PoD) const
return the value of the destruction term from the turbulence model
Definition: DAkOmegaSSTLM.C:1470
Foam::DATurbulenceModel::correctAlphat
void correctAlphat()
update alphat
Definition: DATurbulenceModel.C:203
Foam::DAkOmegaSSTLM::getTurbConvOverProd
virtual void getTurbConvOverProd(volScalarField &CoP) const
return the value of the convective over production term from the turbulence model
Definition: DAkOmegaSSTLM.C:1498
Foam::DAkOmegaSSTLM::saveOmegaNearWall
void saveOmegaNearWall()
save near wall omega values to omegaNearWall_
Definition: DAkOmegaSSTLM.C:696
Foam::DAkOmegaSSTLM::omegaSource
tmp< fvScalarMatrix > omegaSource() const
Definition: DAkOmegaSSTLM.C:328
Foam::DATurbulenceModel
Definition: DATurbulenceModel.H:52
Foam::DAkOmegaSSTLM::F1
tmp< volScalarField > F1(const volScalarField &CDkOmega) const
Definition: DAkOmegaSSTLM.C:348
Foam::DAkOmegaSSTLM::kRes_
volScalarField kRes_
Definition: DAkOmegaSSTLM.H:204
Foam::DAkOmegaSSTLM::DkEff
tmp< volScalarField > DkEff(const volScalarField &F1) const
Definition: DAkOmegaSSTLM.H:252
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::DAkOmegaSSTLM::DgammaIntEff
tmp< volScalarField > DgammaIntEff() const
Definition: DAkOmegaSSTLM.H:277
Foam::DATurbulenceModel::mesh_
const fvMesh & mesh_
fvMesh
Definition: DATurbulenceModel.H:65
lambda
volScalarField & lambda
Definition: createRefsSolidDisplacement.H:7
Foam::DAkOmegaSSTLM::ReThetat0
tmp< volScalarField::Internal > ReThetat0(const volScalarField::Internal &Us, const volScalarField::Internal &dUsds, const volScalarField::Internal &nu) const
Definition: DAkOmegaSSTLM.C:496
Foam::DAkOmegaSSTLM::getFvMatrixFields
virtual void getFvMatrixFields(const word varName, scalarField &diag, scalarField &upper, scalarField &lower)
return the diag(), upper(), and lower() scalarFields from the turbulence model's fvMatrix
Definition: DAkOmegaSSTLM.C:1275
Foam::DAkOmegaSSTLM::gammaIntRes_
volScalarField gammaIntRes_
Definition: DAkOmegaSSTLM.H:208
Foam::DAkOmegaSSTLM::betaFIK_
volScalarField betaFIK_
beta field for field inversion
Definition: DAkOmegaSSTLM.H:222
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:424
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:760
Foam::DAkOmegaSSTLM::correctOmegaBoundaryConditions
void correctOmegaBoundaryConditions()
specially treatment to correct epsilon BC
Definition: DAkOmegaSSTLM.C:669
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::turbModelType_
word turbModelType_
whether the turbulence model is incompressible or compressible
Definition: DATurbulenceModel.H:95
Foam::DATurbulenceModel::phase_
volScalarField phase_
phase field
Definition: DATurbulenceModel.H:86
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:594
Foam::DAkOmegaSSTLM::F1SST
tmp< volScalarField > F1SST(const volScalarField &CDkOmega) const
Definition: DAkOmegaSSTLM.C:243

Generated by   doxygen 1.8.17