DAkOmegaSST.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/kOmegaSST/kOmegaSST.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 "DAkOmegaSST.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 defineTypeNameAndDebug(DAkOmegaSST, 0);
38 addToRunTimeSelectionTable(DATurbulenceModel, DAkOmegaSST, dictionary);
39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40 
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  // Augmented variables
100  omega_(const_cast<volScalarField&>(
101  mesh_.thisDb().lookupObject<volScalarField>("omega"))),
102  omegaRes_(
103  IOobject(
104  "omegaRes",
105  mesh.time().timeName(),
106  mesh,
107  IOobject::NO_READ,
108  IOobject::NO_WRITE),
109  mesh,
110  dimensionedScalar("omegaRes", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
111  zeroGradientFvPatchField<scalar>::typeName),
112  k_(
113  const_cast<volScalarField&>(
114  mesh_.thisDb().lookupObject<volScalarField>("k"))),
115  kRes_(
116  IOobject(
117  "kRes",
118  mesh.time().timeName(),
119  mesh,
120  IOobject::NO_READ,
121  IOobject::NO_WRITE),
122  mesh,
123  dimensionedScalar("kRes", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
124  zeroGradientFvPatchField<scalar>::typeName),
125  y_(mesh_.thisDb().lookupObject<volScalarField>("yWall")),
126  GPtr_(nullptr),
127  betaFIK_(
128  IOobject(
129  "betaFIK",
130  mesh.time().timeName(),
131  mesh,
132  IOobject::READ_IF_PRESENT,
133  IOobject::AUTO_WRITE),
134  mesh,
135  dimensionedScalar("betaFIK", dimensionSet(0, 0, 0, 0, 0, 0, 0), 1.0),
136  "zeroGradient"),
137  betaFIOmega_(
138  IOobject(
139  "betaFIOmega",
140  mesh.time().timeName(),
141  mesh,
142  IOobject::READ_IF_PRESENT,
143  IOobject::AUTO_WRITE),
144  mesh,
145  dimensionedScalar("betaFIOmega", dimensionSet(0, 0, 0, 0, 0, 0, 0), 1.0),
146  "zeroGradient")
147 {
148 
149  if (turbModelType_ == "incompressible")
150  {
151  omegaRes_.dimensions().reset(dimensionSet(0, 0, -2, 0, 0, 0, 0));
152  kRes_.dimensions().reset(dimensionSet(0, 2, -3, 0, 0, 0, 0));
153  }
154  else if (turbModelType_ == "compressible")
155  {
156  omegaRes_.dimensions().reset(dimensionSet(1, -3, -2, 0, 0, 0, 0));
157  kRes_.dimensions().reset(dimensionSet(1, -1, -3, 0, 0, 0, 0));
158  }
159 
160  // calculate the size of omegaWallFunction faces
161  label nWallFaces = 0;
162  forAll(omega_.boundaryField(), patchI)
163  {
164  if (omega_.boundaryField()[patchI].type() == "omegaWallFunction"
165  && omega_.boundaryField()[patchI].size() > 0)
166  {
167  forAll(omega_.boundaryField()[patchI], faceI)
168  {
169  nWallFaces++;
170  }
171  }
172  }
173 
174  // initialize omegaNearWall
175  omegaNearWall_.setSize(nWallFaces);
176 
177  // initialize the G field
178  tmp<volTensorField> tgradU = fvc::grad(U_);
179  volScalarField S2(2 * magSqr(symm(tgradU())));
180  volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
181  GPtr_.reset(new volScalarField::Internal("kOmegaSST:G", nut_ * GbyNu0));
182 }
183 
184 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
185 
186 // SA member functions. these functions are copied from
187 tmp<volScalarField> DAkOmegaSST::F1(
188  const volScalarField& CDkOmega) const
189 {
190 
191  tmp<volScalarField> CDkOmegaPlus = max(
192  CDkOmega,
193  dimensionedScalar("1.0e-10", dimless / sqr(dimTime), 1.0e-10));
194 
195  tmp<volScalarField> arg1 = min(
196  min(
197  max(
198  (scalar(1) / betaStar_) * sqrt(k_) / (omega_ * y_),
199  scalar(500) * (this->nu()) / (sqr(y_) * omega_)),
200  (scalar(4) * alphaOmega2_) * k_ / (CDkOmegaPlus * sqr(y_))),
201  scalar(10));
202 
203  return tanh(pow4(arg1));
204 }
205 
206 tmp<volScalarField> DAkOmegaSST::F2() const
207 {
208 
209  tmp<volScalarField> arg2 = min(
210  max(
211  (scalar(2) / betaStar_) * sqrt(k_) / (omega_ * y_),
212  scalar(500) * (this->nu()) / (sqr(y_) * omega_)),
213  scalar(100));
214 
215  return tanh(sqr(arg2));
216 }
217 
218 tmp<volScalarField> DAkOmegaSST::F3() const
219 {
220 
221  tmp<volScalarField> arg3 = min(
222  150 * (this->nu()) / (omega_ * sqr(y_)),
223  scalar(10));
224 
225  return 1 - tanh(pow4(arg3));
226 }
227 
228 tmp<volScalarField> DAkOmegaSST::F23() const
229 {
230  tmp<volScalarField> f23(F2());
231 
232  if (F3_)
233  {
234  f23.ref() *= F3();
235  }
236 
237  return f23;
238 }
239 
240 tmp<volScalarField::Internal> DAkOmegaSST::GbyNu(
241  const volScalarField::Internal& GbyNu0,
242  const volScalarField::Internal& F2,
243  const volScalarField::Internal& S2) const
244 {
245  return min(
246  GbyNu0,
247  (c1_ / a1_) * betaStar_ * omega_()
248  * max(a1_ * omega_(), b1_ * F2 * sqrt(S2)));
249 }
250 
251 tmp<volScalarField::Internal> DAkOmegaSST::Pk(
252  const volScalarField::Internal& G) const
253 {
254  return min(G, (c1_ * betaStar_) * k_() * omega_());
255 }
256 
257 tmp<volScalarField::Internal> DAkOmegaSST::epsilonByk(
258  const volScalarField& F1,
259  const volTensorField& gradU) const
260 {
261  return betaStar_ * omega_();
262 }
263 
264 tmp<fvScalarMatrix> DAkOmegaSST::kSource() const
265 {
266  return tmp<fvScalarMatrix>(
267  new fvScalarMatrix(
268  k_,
269  dimVolume * this->rhoDimensions() * k_.dimensions() / dimTime));
270 }
271 
272 tmp<fvScalarMatrix> DAkOmegaSST::omegaSource() const
273 {
274  return tmp<fvScalarMatrix>(
275  new fvScalarMatrix(
276  omega_,
277  dimVolume * this->rhoDimensions() * omega_.dimensions() / dimTime));
278 }
279 
280 tmp<fvScalarMatrix> DAkOmegaSST::Qsas(
281  const volScalarField::Internal& S2,
282  const volScalarField::Internal& gamma,
283  const volScalarField::Internal& beta) const
284 {
285  return tmp<fvScalarMatrix>(
286  new fvScalarMatrix(
287  omega_,
288  dimVolume * this->rhoDimensions() * omega_.dimensions() / dimTime));
289 }
290 
291 // Augmented functions
292 void DAkOmegaSST::correctModelStates(wordList& modelStates) const
293 {
294  /*
295  Description:
296  Update the name in modelStates based on the selected physical model at runtime
297 
298  Example:
299  In DAStateInfo, if the modelStates reads:
300 
301  modelStates = {"nut"}
302 
303  then for the SA model, calling correctModelStates(modelStates) will give:
304 
305  modelStates={"nuTilda"}
306 
307  while calling correctModelStates(modelStates) for the SST model will give
308 
309  modelStates={"k","omega"}
310 
311  We don't udpate the names for the radiation model because users are
312  supposed to set modelStates={"G"}
313  */
314 
315  // For SST model, we need to replace nut with k, omega
316 
317  forAll(modelStates, idxI)
318  {
319  word stateName = modelStates[idxI];
320  if (stateName == "nut")
321  {
322  modelStates[idxI] = "omega";
323  modelStates.append("k");
324  }
325  }
326 }
327 
329 {
330  /*
331  Description:
332  Update nut based on other turbulence variables and update the BCs
333  Also update alphat if is present
334  */
335 
336  const volVectorField U = mesh_.thisDb().lookupObject<volVectorField>("U");
337  tmp<volTensorField> tgradU = fvc::grad(U);
338  volScalarField S2(2 * magSqr(symm(tgradU())));
339 
340  nut_ = a1_ * k_ / max(a1_ * omega_, b1_ * F23() * sqrt(S2));
341 
342  nut_.correctBoundaryConditions(); // nutkWallFunction: update wall face nut based on k
343 
344  // this is basically BasicTurbulenceModel::correctNut();
345  this->correctAlphat();
346 
347  return;
348 }
349 
351 {
352  /*
353  Description:
354  Update turbulence variable boundary values
355  */
356 
357  // correct the BCs for the perturbed fields
358  // kqWallFunction is a zero-gradient BC
359  k_.correctBoundaryConditions();
360 }
361 
363 {
364  /*
365  Description:
366  this is a special treatment for omega BC because we cant directly call omega.
367  correctBoundaryConditions() because it will modify the internal omega and G that
368  are right next to walls. This will mess up adjoint Jacobians
369  To solve this issue,
370  1. we store the near wall omega before calling omega.correctBoundaryConditions()
371  2. call omega.correctBoundaryConditions()
372  3. Assign the stored near wall omega back
373  4. Apply a zeroGradient BC for omega at the wall patches
374  *********** NOTE *************
375  this treatment will obviously downgrade the accuracy of adjoint derivative since it is
376  not 100% consistent with what is used for the flow solver; however, our observation
377  shows that the impact is not very large.
378  */
379 
380  // save the perturbed omega at the wall
381  this->saveOmegaNearWall();
382  // correct omega boundary conditions, this includes updating wall face and near wall omega values,
383  // updating the inter-proc BCs
384  omega_.correctBoundaryConditions();
385  // reset the corrected omega near wall cell to its perturbed value
386  this->setOmegaNearWall();
387 }
388 
390 {
391  /*
392  Description:
393  Save the current near wall omega values to omegaNearWall_
394  */
395  label counterI = 0;
396  forAll(omega_.boundaryField(), patchI)
397  {
398  if (omega_.boundaryField()[patchI].type() == "omegaWallFunction"
399  and omega_.boundaryField()[patchI].size() > 0)
400  {
401  const UList<label>& faceCells = mesh_.boundaryMesh()[patchI].faceCells();
402  forAll(faceCells, faceI)
403  {
404  //Info<<"Near Wall cellI: "<<faceCells[faceI]<<endl;
405  omegaNearWall_[counterI] = omega_[faceCells[faceI]];
406  counterI++;
407  }
408  }
409  }
410  return;
411 }
412 
414 {
415  /*
416  Description:
417  Set the current near wall omega values from omegaNearWall_
418  Here we also apply a zeroGradient BC to the wall faces
419  */
420  label counterI = 0;
421  forAll(omega_.boundaryField(), patchI)
422  {
423  if (omega_.boundaryField()[patchI].type() == "omegaWallFunction"
424  && omega_.boundaryField()[patchI].size() > 0)
425  {
426  const UList<label>& faceCells = mesh_.boundaryMesh()[patchI].faceCells();
427  forAll(faceCells, faceI)
428  {
429  omega_[faceCells[faceI]] = omegaNearWall_[counterI];
430  // zeroGradient BC
431  omega_.boundaryFieldRef()[patchI][faceI] = omega_[faceCells[faceI]];
432  counterI++;
433  }
434  }
435  }
436  return;
437 }
438 
440 {
441  /*
442  Description:
443  Update nut based on nuTilda. Note: we need to update nut and its BC since we
444  may have perturbed other turbulence vars that affect the nut values
445  */
446 
447  this->correctNut();
448 }
449 
450 void DAkOmegaSST::correctStateResidualModelCon(List<List<word>>& stateCon) const
451 {
452  /*
453  Description:
454  Update the original variable connectivity for the adjoint state
455  residuals in stateCon. Basically, we modify/add state variables based on the
456  original model variables defined in stateCon.
457 
458  Input:
459 
460  stateResCon: the connectivity levels for a state residual, defined in Foam::DAJacCon
461 
462  Example:
463  If stateCon reads:
464  stateCon=
465  {
466  {"U", "p", "nut"},
467  {"p"}
468  }
469 
470  For the SA turbulence model, calling this function for will get a new stateCon
471  stateCon=
472  {
473  {"U", "p", "nuTilda"},
474  {"p"}
475  }
476 
477  For the SST turbulence model, calling this function will give
478  stateCon=
479  {
480  {"U", "p", "k", "omega"},
481  {"p", "U"}
482  }
483  ***NOTE***: we add a extra level of U connectivity because nut is
484  related to grad(U), k, and omega in SST!
485  */
486 
487  label stateConSize = stateCon.size();
488  forAll(stateCon, idxI)
489  {
490  label addUCon = 0;
491  forAll(stateCon[idxI], idxJ)
492  {
493  word conStateName = stateCon[idxI][idxJ];
494  if (conStateName == "nut")
495  {
496  stateCon[idxI][idxJ] = "omega"; // replace nut with omega
497  stateCon[idxI].append("k"); // also add k for that level
498  addUCon = 1;
499  }
500  }
501  // add U for the current level and level+1 if it is not there yet
502  label isU;
503  if (addUCon == 1)
504  {
505  // first add U for the current level
506  isU = 0;
507  forAll(stateCon[idxI], idxJ)
508  {
509  word conStateName = stateCon[idxI][idxJ];
510  if (conStateName == "U")
511  {
512  isU = 1;
513  }
514  }
515  if (!isU)
516  {
517  stateCon[idxI].append("U");
518  }
519 
520  // now add U for level+1 if idxI is not the largest level
521  // if idxI is already the largest level, we have a problem
522  if (idxI != stateConSize - 1)
523  {
524  isU = 0;
525  forAll(stateCon[idxI + 1], idxJ)
526  {
527  word conStateName = stateCon[idxI + 1][idxJ];
528  if (conStateName == "U")
529  {
530  isU = 1;
531  }
532  }
533  if (!isU)
534  {
535  stateCon[idxI + 1].append("U");
536  }
537  }
538  else
539  {
540  FatalErrorIn(
541  "In DAStateInfo, nut shows in the largest connectivity level! "
542  "This is not supported!")
543  << exit(FatalError);
544  }
545  }
546  }
547 }
548 
549 void DAkOmegaSST::addModelResidualCon(HashTable<List<List<word>>>& allCon) const
550 {
551  /*
552  Description:
553  Add the connectivity levels for all physical model residuals to allCon
554 
555  Input:
556  allCon: the connectivity levels for all state residual, defined in DAJacCon
557 
558  Example:
559  If stateCon reads:
560  allCon=
561  {
562  "URes":
563  {
564  {"U", "p", "nut"},
565  {"p"}
566  }
567  }
568 
569  For the SA turbulence model, calling this function for will get a new stateCon,
570  something like this:
571  allCon=
572  {
573  "URes":
574  {
575  {"U", "p", "nuTilda"},
576  {"p"}
577  },
578  "nuTildaRes":
579  {
580  {"U", "phi", "nuTilda"},
581  {"U"}
582  }
583  }
584 
585  */
586 
587  word pName;
588 
589  if (mesh_.thisDb().foundObject<volScalarField>("p"))
590  {
591  pName = "p";
592  }
593  else if (mesh_.thisDb().foundObject<volScalarField>("p_rgh"))
594  {
595  pName = "p_rgh";
596  }
597  else
598  {
599  FatalErrorIn(
600  "Neither p nor p_rgh was found in mesh.thisDb()!"
601  "addModelResidualCon failed to setup turbulence residuals!")
602  << exit(FatalError);
603  }
604 
605  // NOTE: for compressible flow, it depends on rho so we need to add T and p
606  if (turbModelType_ == "incompressible")
607  {
608  allCon.set(
609  "omegaRes",
610  {
611  {"U", "omega", "k", "phi"}, // lv0
612  {"U", "omega", "k"}, // lv1
613  {"U", "omega", "k"} // lv2
614  });
615  allCon.set(
616  "kRes",
617  {
618  {"U", "omega", "k", "phi"}, // lv0
619  {"U", "omega", "k"}, // lv1
620  {"U", "omega", "k"} // lv2
621  });
622  }
623  else if (turbModelType_ == "compressible")
624  {
625  allCon.set(
626  "omegaRes",
627  {
628  {"U", "T", pName, "omega", "k", "phi"}, // lv0
629  {"U", "T", pName, "omega", "k"}, // lv1
630  {"U", "T", pName, "omega", "k"} // lv2
631  });
632  allCon.set(
633  "kRes",
634  {
635  {"U", "T", pName, "omega", "k", "phi"}, // lv0
636  {"U", "T", pName, "omega", "k"}, // lv1
637  {"U", "T", pName, "omega", "k"} // lv2
638  });
639  }
640 }
641 
642 void DAkOmegaSST::correct(label printToScreen)
643 {
644  /*
645  Descroption:
646  Solve the residual equations and update the state. This function will be called
647  by the DASolver. It is needed because we want to control the output frequency
648  of the residual convergence every 100 steps. If using the correct from turbulence
649  it will output residual every step which will be too much of information.
650  */
651 
652  // We set the flag solveTurbState_ to 1 such that in the calcResiduals function
653  // we will solve and update nuTilda
654  solveTurbState_ = 1;
655  dictionary dummyOptions;
656  dummyOptions.set("printToScreen", printToScreen);
657  this->calcResiduals(dummyOptions);
658  // after it, we reset solveTurbState_ = 0 such that calcResiduals will not
659  // update nuTilda when calling from the adjoint class, i.e., solveAdjoint from DASolver.
660  solveTurbState_ = 0;
661 }
662 
663 void DAkOmegaSST::calcResiduals(const dictionary& options)
664 {
665  /*
666  Descroption:
667  If solveTurbState_ == 1, this function solve and update k and omega, and
668  is the same as calling turbulence.correct(). If solveTurbState_ == 0,
669  this function compute residuals for turbulence variables, e.g., nuTildaRes_
670 
671  Input:
672  options.isPC: 1 means computing residuals for preconditioner matrix.
673  This essentially use the first order scheme for div(phi,nuTilda)
674 
675  p_, U_, phi_, etc: State variables in OpenFOAM
676 
677  Output:
678  kRes_/omegaRes_: If solveTurbState_ == 0, update the residual field variable
679 
680  k_/omega_: If solveTurbState_ == 1, update them
681  */
682 
683  // Copy and modify based on the "correct" function
684 
685  word divKScheme = "div(phi,k)";
686  word divOmegaScheme = "div(phi,omega)";
687 
688  label isPC = 0;
689  label printToScreen = options.lookupOrDefault("printToScreen", 0);
690 
691  if (!solveTurbState_)
692  {
693  isPC = options.getLabel("isPC");
694 
695  if (isPC)
696  {
697  divKScheme = "div(pc)";
698  divOmegaScheme = "div(pc)";
699  }
700  }
701 
702  volScalarField rho = this->rho();
703 
704  // Note: for compressible flow, the "this->phi()" function divides phi by fvc:interpolate(rho),
705  // while for the incompresssible "this->phi()" returns phi only
706  // see src/TurbulenceModels/compressible/compressibleTurbulenceModel.C line 62 to 73
707  volScalarField::Internal divU(fvc::div(fvc::absolute(phi_ / fvc::interpolate(rho), U_)));
708 
709  tmp<volTensorField> tgradU = fvc::grad(U_);
710  volScalarField S2(2 * magSqr(symm(tgradU())));
711  volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
712  volScalarField::Internal& G = const_cast<volScalarField::Internal&>(GPtr_());
713  G = nut_() * GbyNu0;
714 
715  if (solveTurbState_)
716  {
717  // Update omega and G at the wall
718  omega_.boundaryFieldRef().updateCoeffs();
719  }
720  else
721  {
722  // NOTE instead of calling omega_.boundaryFieldRef().updateCoeffs();
723  // here we call our self-defined boundary conditions
725  }
726 
727  volScalarField CDkOmega(
728  (scalar(2) * alphaOmega2_) * (fvc::grad(k_) & fvc::grad(omega_)) / omega_);
729 
730  volScalarField F1(this->F1(CDkOmega));
731  volScalarField F23(this->F23());
732 
733  {
734 
735  volScalarField::Internal gamma(this->gamma(F1));
736  volScalarField::Internal beta(this->beta(F1));
737 
738  // Turbulent frequency equation
739  tmp<fvScalarMatrix> omegaEqn(
740  fvm::ddt(phase_, rho, omega_)
741  + fvm::div(phaseRhoPhi_, omega_, divOmegaScheme)
742  - fvm::laplacian(phase_ * rho * DomegaEff(F1), omega_)
743  == phase_() * rho() * gamma * GbyNu(GbyNu0, F23(), S2()) * betaFIOmega_()
744  - fvm::SuSp((2.0 / 3.0) * phase_() * rho() * gamma * divU, omega_)
745  - fvm::Sp(phase_() * rho() * beta * omega_(), omega_)
746  - fvm::SuSp(
747  phase_() * rho() * (F1() - scalar(1)) * CDkOmega() / omega_(),
748  omega_)
749  + Qsas(S2(), gamma, beta)
750  + omegaSource()
751 
752  );
753 
754  omegaEqn.ref().relax();
755  omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
756 
757  if (solveTurbState_)
758  {
759 
760  // get the solver performance info such as initial
761  // and final residuals
762  SolverPerformance<scalar> solverOmega = solve(omegaEqn);
763  DAUtility::primalResidualControl(solverOmega, printToScreen, "omega", daGlobalVar_.primalMaxRes);
764 
765  DAUtility::boundVar(allOptions_, omega_, printToScreen);
766  }
767  else
768  {
769  // reset the corrected omega near wall cell to its perturbed value
770  this->setOmegaNearWall();
771 
772  // calculate residuals
773  omegaRes_ = omegaEqn() & omega_;
774  // need to normalize residuals
775  normalizeResiduals(omegaRes);
776  }
777  }
778 
779  // Turbulent kinetic energy equation
780  tmp<fvScalarMatrix> kEqn(
781  fvm::ddt(phase_, rho, k_)
782  + fvm::div(phaseRhoPhi_, k_, divKScheme)
783  - fvm::laplacian(phase_ * rho * DkEff(F1), k_)
784  == phase_() * rho() * Pk(G) * betaFIK_()
785  - fvm::SuSp((2.0 / 3.0) * phase_() * rho() * divU, k_)
786  - fvm::Sp(phase_() * rho() * epsilonByk(F1, tgradU()), k_)
787  + kSource());
788 
789  tgradU.clear();
790 
791  kEqn.ref().relax();
792 
793  if (solveTurbState_)
794  {
795 
796  // get the solver performance info such as initial
797  // and final residuals
798  SolverPerformance<scalar> solverK = solve(kEqn);
799  DAUtility::primalResidualControl(solverK, printToScreen, "k", daGlobalVar_.primalMaxRes);
800 
801  DAUtility::boundVar(allOptions_, k_, printToScreen);
802 
803  this->correctNut();
804  }
805  else
806  {
807  // calculate residuals
808  kRes_ = kEqn() & k_;
809  // need to normalize residuals
810  normalizeResiduals(kRes);
811  }
812 
813  return;
814 }
815 
817  const word varName,
818  scalarField& diag,
819  scalarField& upper,
820  scalarField& lower)
821 {
822  /*
823  Description:
824  return the diag(), upper(), and lower() scalarFields from the turbulence model's fvMatrix
825  this will be use to compute the preconditioner matrix
826  */
827 
828  if (varName != "k" && varName != "omega")
829  {
830  FatalErrorIn(
831  "varName not valid. It has to be k or omega")
832  << exit(FatalError);
833  }
834 
835  volScalarField rho = this->rho();
836 
837  // Note: for compressible flow, the "this->phi()" function divides phi by fvc:interpolate(rho),
838  // while for the incompresssible "this->phi()" returns phi only
839  // see src/TurbulenceModels/compressible/compressibleTurbulenceModel.C line 62 to 73
840  volScalarField::Internal divU(fvc::div(fvc::absolute(phi_ / fvc::interpolate(rho), U_)));
841 
842  tmp<volTensorField> tgradU = fvc::grad(U_);
843  volScalarField S2(2 * magSqr(symm(tgradU())));
844  volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
845  volScalarField::Internal& G = const_cast<volScalarField::Internal&>(GPtr_());
846  G = nut_() * GbyNu0;
847 
848  // NOTE instead of calling omega_.boundaryFieldRef().updateCoeffs();
849  // here we call our self-defined boundary conditions
851 
852  volScalarField CDkOmega(
853  (scalar(2) * alphaOmega2_) * (fvc::grad(k_) & fvc::grad(omega_)) / omega_);
854 
855  volScalarField F1(this->F1(CDkOmega));
856  volScalarField F23(this->F23());
857 
858  if (varName == "omega")
859  {
860  volScalarField::Internal gamma(this->gamma(F1));
861  volScalarField::Internal beta(this->beta(F1));
862 
863  // Turbulent frequency equation
864  fvScalarMatrix omegaEqn(
865  fvm::ddt(phase_, rho, omega_)
866  + fvm::div(phaseRhoPhi_, omega_, "div(pc)")
867  - fvm::laplacian(phase_ * rho * DomegaEff(F1), omega_)
868  == phase_() * rho() * gamma * GbyNu(GbyNu0, F23(), S2()) * betaFIOmega_()
869  - fvm::SuSp((2.0 / 3.0) * phase_() * rho() * gamma * divU, omega_)
870  - fvm::Sp(phase_() * rho() * beta * omega_(), omega_)
871  - fvm::SuSp(
872  phase_() * rho() * (F1() - scalar(1)) * CDkOmega() / omega_(),
873  omega_)
874  + Qsas(S2(), gamma, beta)
875  + omegaSource()
876 
877  );
878 
879  omegaEqn.relax();
880 
881  // reset the corrected omega near wall cell to its perturbed value
882  this->setOmegaNearWall();
883 
884  diag = omegaEqn.D();
885  upper = omegaEqn.upper();
886  lower = omegaEqn.lower();
887  }
888  else if (varName == "k")
889  {
890  // Turbulent kinetic energy equation
891  fvScalarMatrix kEqn(
892  fvm::ddt(phase_, rho, k_)
893  + fvm::div(phaseRhoPhi_, k_, "div(pc)")
894  - fvm::laplacian(phase_ * rho * DkEff(F1), k_)
895  == phase_() * rho() * Pk(G) * betaFIK_()
896  - fvm::SuSp((2.0 / 3.0) * phase_() * rho() * divU, k_)
897  - fvm::Sp(phase_() * rho() * epsilonByk(F1, tgradU()), k_)
898  + kSource());
899 
900  kEqn.relax();
901 
902  diag = kEqn.D();
903  upper = kEqn.upper();
904  lower = kEqn.lower();
905  }
906 }
907 
908 void DAkOmegaSST::getTurbProdOverDestruct(volScalarField& PoD) const
909 {
910  /*
911  Description:
912  Return the value of the production over destruction term from the turbulence model
913  */
914  tmp<volTensorField> tgradU = fvc::grad(U_);
915  volScalarField S2(2 * magSqr(symm(tgradU())));
916  volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
917  volScalarField::Internal& G = const_cast<volScalarField::Internal&>(GPtr_());
918  G = nut_() * GbyNu0;
919 
920  volScalarField rho = this->rho();
921 
922  volScalarField CDkOmega(
923  (scalar(2) * alphaOmega2_) * (fvc::grad(k_) & fvc::grad(omega_)) / omega_);
924 
925  volScalarField F1(this->F1(CDkOmega));
926 
927  volScalarField::Internal P = phase_() * rho() * Pk(G);
928  volScalarField::Internal D = phase_() * rho() * epsilonByk(F1, tgradU()) * k_();
929 
930  forAll(P, cellI)
931  {
932  PoD[cellI] = P[cellI] / (D[cellI] + P[cellI] + 1e-16);
933  }
934 }
935 
936 void DAkOmegaSST::getTurbConvOverProd(volScalarField& CoP) const
937 {
938  /*
939  Description:
940  Return the value of the convective over production term from the turbulence model
941  */
942 
943  tmp<volTensorField> tgradU = fvc::grad(U_);
944  volScalarField S2(2 * magSqr(symm(tgradU())));
945  volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
946  volScalarField::Internal& G = const_cast<volScalarField::Internal&>(GPtr_());
947  G = nut_() * GbyNu0;
948 
949  volScalarField rho = this->rho();
950 
951  volScalarField CDkOmega(
952  (scalar(2) * alphaOmega2_) * (fvc::grad(k_) & fvc::grad(omega_)) / omega_);
953 
954  volScalarField F1(this->F1(CDkOmega));
955 
956  volScalarField::Internal P = phase_() * rho() * Pk(G);
957  volScalarField C = fvc::div(phaseRhoPhi_, k_);
958 
959  forAll(P, cellI)
960  {
961  CoP[cellI] = C[cellI] / (P[cellI] + C[cellI] + 1e-16);
962  }
963 }
964 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
965 
966 } // End namespace Foam
967 
968 // ************************************************************************* //
Foam::DAkOmegaSST::c1_
dimensionedScalar c1_
Definition: DAkOmegaSST.H:71
U
volVectorField & U
Definition: createRefsPimpleDyM.H:7
Foam::DAkOmegaSST::epsilonByk
tmp< volScalarField::Internal > epsilonByk(const volScalarField &F1, const volTensorField &gradU) const
Definition: DAkOmegaSST.C:257
Foam::DAkOmegaSST::omegaNearWall_
scalarList omegaNearWall_
Definition: DAkOmegaSST.H:185
Foam::DAkOmegaSST::correctModelStates
virtual void correctModelStates(wordList &modelStates) const
update the turbulence state for DAStateInfo::regStates_
Definition: DAkOmegaSST.C:292
Foam::DATurbulenceModel::phi_
surfaceScalarField & phi_
face flux
Definition: DATurbulenceModel.H:83
solve
nuTilda1Eqn solve(solverDictNuTilda)
D
volVectorField & D
Definition: createRefsSolidDisplacement.H:8
Foam::DAkOmegaSST::betaStar_
dimensionedScalar betaStar_
Definition: DAkOmegaSST.H:67
Foam::DAkOmegaSST::kSource
tmp< fvScalarMatrix > kSource() const
Definition: DAkOmegaSST.C:264
Foam::DAOption
Definition: DAOption.H:29
Foam::DAkOmegaSST::F1
tmp< volScalarField > F1(const volScalarField &CDkOmega) const
Definition: DAkOmegaSST.C:187
Foam::DAkOmegaSST::Qsas
tmp< fvScalarMatrix > Qsas(const volScalarField::Internal &S2, const volScalarField::Internal &gamma, const volScalarField::Internal &beta) const
Definition: DAkOmegaSST.C:280
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
Foam::DAkOmegaSST::k_
volScalarField & k_
Definition: DAkOmegaSST.H:166
Foam::DAkOmegaSST::Pk
tmp< volScalarField::Internal > Pk(const volScalarField::Internal &G) const
Definition: DAkOmegaSST.C:251
Foam::DAkOmegaSST::DAkOmegaSST
DAkOmegaSST(const word modelType, const fvMesh &mesh, const DAOption &daOption)
Definition: DAkOmegaSST.C:41
normalizeResiduals
#define normalizeResiduals(resName)
Definition: DAMacroFunctions.H:28
Foam::DATurbulenceModel::phaseRhoPhi_
surfaceScalarField & phaseRhoPhi_
phase*phi*density field
Definition: DATurbulenceModel.H:89
Foam::DAkOmegaSST::F3_
Switch F3_
Definition: DAkOmegaSST.H:73
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::DAkOmegaSST::omegaSource
tmp< fvScalarMatrix > omegaSource() const
Definition: DAkOmegaSST.C:272
Foam::DAkOmegaSST::solveTurbState_
label solveTurbState_
whether to solve for turb states
Definition: DAkOmegaSST.H:188
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam::DAkOmegaSST::omega_
volScalarField & omega_
Definition: DAkOmegaSST.H:164
Foam::DAkOmegaSST::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: DAkOmegaSST.C:816
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
Foam::DAkOmegaSST::saveOmegaNearWall
void saveOmegaNearWall()
save near wall omega values to omegaNearWall_
Definition: DAkOmegaSST.C:389
Foam::DAkOmegaSST::correctOmegaBoundaryConditions
void correctOmegaBoundaryConditions()
specially treatment to correct epsilon BC
Definition: DAkOmegaSST.C:362
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::DAkOmegaSST::F23
tmp< volScalarField > F23() const
Definition: DAkOmegaSST.C:228
Foam::DAkOmegaSST::DomegaEff
tmp< volScalarField > DomegaEff(const volScalarField &F1) const
Definition: DAkOmegaSST.H:129
Foam::DATurbulenceModel::daGlobalVar_
DAGlobalVar & daGlobalVar_
global variable
Definition: DATurbulenceModel.H:74
Foam::DAkOmegaSST::betaFIOmega_
volScalarField betaFIOmega_
Definition: DAkOmegaSST.H:179
Foam::DAkOmegaSST::GbyNu
tmp< volScalarField::Internal > GbyNu(const volScalarField::Internal &GbyNu0, const volScalarField::Internal &F2, const volScalarField::Internal &S2) const
Definition: DAkOmegaSST.C:240
Foam::DATurbulenceModel::U_
volVectorField & U_
velocity
Definition: DATurbulenceModel.H:80
Foam::DATurbulenceModel::nut_
volScalarField & nut_
turbulence viscosity
Definition: DATurbulenceModel.H:77
Foam::DAGlobalVar::primalMaxRes
scalar primalMaxRes
the maximal residual for primal this var will be used in multiple classes
Definition: DAGlobalVar.H:62
Foam::DAkOmegaSST::a1_
dimensionedScalar a1_
Definition: DAkOmegaSST.H:69
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
Foam::DAkOmegaSST::getTurbProdOverDestruct
virtual void getTurbProdOverDestruct(volScalarField &PoD) const
return the value of the destruction term from the turbulence model
Definition: DAkOmegaSST.C:908
Foam::DAkOmegaSST::y_
const volScalarField & y_
3D wall distance
Definition: DAkOmegaSST.H:171
Foam::DAkOmegaSST::betaFIK_
volScalarField betaFIK_
beta field for field inversion
Definition: DAkOmegaSST.H:178
Foam
Definition: checkGeometry.C:32
Foam::DAkOmegaSST::gamma
tmp< volScalarField::Internal > gamma(const volScalarField::Internal &F1) const
Definition: DAkOmegaSST.H:115
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
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
Foam::DAkOmegaSST::F2
tmp< volScalarField > F2() const
Definition: DAkOmegaSST.C:206
Foam::DAkOmegaSST::b1_
dimensionedScalar b1_
Definition: DAkOmegaSST.H:70
Foam::DATurbulenceModel::correctAlphat
void correctAlphat()
update alphat
Definition: DATurbulenceModel.C:203
Foam::DATurbulenceModel
Definition: DATurbulenceModel.H:52
Foam::DAkOmegaSST::updateIntermediateVariables
virtual void updateIntermediateVariables()
update any intermediate variables that are dependent on state variables and are used in calcResiduals
Definition: DAkOmegaSST.C:439
Foam::DAkOmegaSST::GPtr_
autoPtr< volScalarField::Internal > GPtr_
Definition: DAkOmegaSST.H:175
Foam::DATurbulenceModel::mesh_
const fvMesh & mesh_
fvMesh
Definition: DATurbulenceModel.H:65
Foam::DAkOmegaSST::DkEff
tmp< volScalarField > DkEff(const volScalarField &F1) const
Definition: DAkOmegaSST.H:122
Foam::DAkOmegaSST::correct
virtual void correct(label printToScreen)
solve the residual equations and update the state
Definition: DAkOmegaSST.C:642
Foam::DAkOmegaSST::correctStateResidualModelCon
virtual void correctStateResidualModelCon(List< List< word >> &stateCon) const
update the original variable connectivity for the adjoint state residuals in stateCon
Definition: DAkOmegaSST.C:450
Foam::DAkOmegaSST::calcResiduals
virtual void calcResiduals(const dictionary &options)
compute the turbulence residuals
Definition: DAkOmegaSST.C:663
Foam::DAkOmegaSST::setOmegaNearWall
void setOmegaNearWall()
set omegaNearWall_ to near wall omega values
Definition: DAkOmegaSST.C:413
Foam::DAkOmegaSST::getTurbConvOverProd
virtual void getTurbConvOverProd(volScalarField &CoP) const
return the value of the convective over production term from the turbulence model
Definition: DAkOmegaSST.C:936
DAkOmegaSST.H
Foam::DAkOmegaSST::kRes_
volScalarField kRes_
Definition: DAkOmegaSST.H:167
Foam::DAkOmegaSST::alphaOmega2_
dimensionedScalar alphaOmega2_
Definition: DAkOmegaSST.H:59
Foam::DAkOmegaSST::correctBoundaryConditions
virtual void correctBoundaryConditions()
update turbulence variable boundary values
Definition: DAkOmegaSST.C:350
Foam::DAkOmegaSST::F3
tmp< volScalarField > F3() const
Definition: DAkOmegaSST.C:218
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::DAkOmegaSST::addModelResidualCon
virtual void addModelResidualCon(HashTable< List< List< word >>> &allCon) const
add the model residual connectivity to stateCon
Definition: DAkOmegaSST.C:549
Foam::DAkOmegaSST::beta
tmp< volScalarField::Internal > beta(const volScalarField::Internal &F1) const
Definition: DAkOmegaSST.H:109
Foam::DAkOmegaSST::correctNut
virtual void correctNut()
update nut based on other turbulence variables and update the BCs
Definition: DAkOmegaSST.C:328
Foam::DAkOmegaSST::omegaRes_
volScalarField omegaRes_
Definition: DAkOmegaSST.H:165