FreeFOAM The Cross-Platform CFD Toolkit
LaunderGibsonRSTM.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "LaunderGibsonRSTM.H"
29 
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 namespace compressible
37 {
38 namespace RASModels
39 {
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
43 defineTypeNameAndDebug(LaunderGibsonRSTM, 0);
44 addToRunTimeSelectionTable(RASModel, LaunderGibsonRSTM, dictionary);
45 
46 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
47 
48 LaunderGibsonRSTM::LaunderGibsonRSTM
49 (
50  const volScalarField& rho,
51  const volVectorField& U,
52  const surfaceScalarField& phi,
53  const basicThermo& thermophysicalModel
54 )
55 :
56  RASModel(typeName, rho, U, phi, thermophysicalModel),
57 
58  Cmu_
59  (
61  (
62  "Cmu",
63  coeffDict_,
64  0.09
65  )
66  ),
67  kappa_
68  (
70  (
71  "kappa",
72  coeffDict_,
73  0.41
74  )
75  ),
76  Clg1_
77  (
79  (
80  "Clg1",
81  coeffDict_,
82  1.8
83  )
84  ),
85  Clg2_
86  (
88  (
89  "Clg2",
90  coeffDict_,
91  0.6
92  )
93  ),
94  C1_
95  (
97  (
98  "C1",
99  coeffDict_,
100  1.44
101  )
102  ),
103  C2_
104  (
106  (
107  "C2",
108  coeffDict_,
109  1.92
110  )
111  ),
112  Cs_
113  (
115  (
116  "Cs",
117  coeffDict_,
118  0.25
119  )
120  ),
121  Ceps_
122  (
124  (
125  "Ceps",
126  coeffDict_,
127  0.15
128  )
129  ),
130  C1Ref_
131  (
133  (
134  "C1Ref",
135  coeffDict_,
136  0.5
137  )
138  ),
139  C2Ref_
140  (
142  (
143  "C2Ref",
144  coeffDict_,
145  0.3
146  )
147  ),
148  couplingFactor_
149  (
151  (
152  "couplingFactor",
153  coeffDict_,
154  0.0
155  )
156  ),
157  sigmaR_
158  (
160  (
161  "sigmaR",
162  coeffDict_,
163  0.81967
164  )
165  ),
166  sigmaEps_
167  (
169  (
170  "sigmaEps",
171  coeffDict_,
172  1.3
173  )
174  ),
175  Prt_
176  (
178  (
179  "Prt",
180  coeffDict_,
181  1.0
182  )
183  ),
184 
185  y_(mesh_),
186 
187  R_
188  (
189  IOobject
190  (
191  "R",
192  runTime_.timeName(),
193  mesh_,
196  ),
197  autoCreateR("R", mesh_)
198  ),
199  k_
200  (
201  IOobject
202  (
203  "k",
204  runTime_.timeName(),
205  mesh_,
208  ),
209  autoCreateK("k", mesh_)
210  ),
211  epsilon_
212  (
213  IOobject
214  (
215  "epsilon",
216  runTime_.timeName(),
217  mesh_,
220  ),
221  autoCreateEpsilon("epsilon", mesh_)
222  ),
223  mut_
224  (
225  IOobject
226  (
227  "mut",
228  runTime_.timeName(),
229  mesh_,
232  ),
233  autoCreateMut("mut", mesh_)
234  ),
235  alphat_
236  (
237  IOobject
238  (
239  "alphat",
240  runTime_.timeName(),
241  mesh_,
244  ),
245  autoCreateAlphat("alphat", mesh_)
246  )
247 {
248  if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
249  {
251  (
252  "LaunderGibsonRSTM::LaunderGibsonRSTM"
253  "(const volScalarField&, const volVectorField&"
254  ", const surfaceScalarField&, basicThermo&)"
255  ) << "couplingFactor = " << couplingFactor_
256  << " is not in range 0 - 1" << nl
257  << exit(FatalError);
258  }
259 
260  mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
261  mut_.correctBoundaryConditions();
262 
263  alphat_ = mut_/Prt_;
264  alphat_.correctBoundaryConditions();
265 
266  printCoeffs();
267 }
268 
269 
270 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
271 
272 tmp<volSymmTensorField> LaunderGibsonRSTM::devRhoReff() const
273 {
275  (
277  (
278  IOobject
279  (
280  "devRhoReff",
281  runTime_.timeName(),
282  mesh_,
285  ),
286  rho_*R_ - mu()*dev(twoSymm(fvc::grad(U_)))
287  )
288  );
289 }
290 
291 
292 tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevRhoReff(volVectorField& U) const
293 {
294  if (couplingFactor_.value() > 0.0)
295  {
296  return
297  (
298  fvc::div(rho_*R_ + couplingFactor_*mut_*fvc::grad(U))
299  + fvc::laplacian((1.0 - couplingFactor_)*mut_, U)
300  - fvm::laplacian(muEff(), U)
301  - fvc::div(mu()*dev2(fvc::grad(U)().T()))
302  );
303  }
304  else
305  {
306  return
307  (
308  fvc::div(rho_*R_)
309  + fvc::laplacian(mut_, U)
310  - fvm::laplacian(muEff(), U)
311  - fvc::div(mu()*dev2(fvc::grad(U)().T()))
312  );
313  }
314 }
315 
316 
317 bool LaunderGibsonRSTM::read()
318 {
319  if (RASModel::read())
320  {
321  Cmu_.readIfPresent(coeffDict());
322  kappa_.readIfPresent(coeffDict());
323  Clg1_.readIfPresent(coeffDict());
324  Clg2_.readIfPresent(coeffDict());
325  C1_.readIfPresent(coeffDict());
326  C2_.readIfPresent(coeffDict());
327  Cs_.readIfPresent(coeffDict());
328  Ceps_.readIfPresent(coeffDict());
329  C1Ref_.readIfPresent(coeffDict());
330  C2Ref_.readIfPresent(coeffDict());
331  sigmaR_.readIfPresent(coeffDict());
332  sigmaEps_.readIfPresent(coeffDict());
333  Prt_.readIfPresent(coeffDict());
334 
335  couplingFactor_.readIfPresent(coeffDict());
336 
337  if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
338  {
339  FatalErrorIn("LaunderGibsonRSTM::read()")
340  << "couplingFactor = " << couplingFactor_
341  << " is not in range 0 - 1" << nl
342  << exit(FatalError);
343  }
344 
345  return true;
346  }
347  else
348  {
349  return false;
350  }
351 }
352 
353 
355 {
356  if (!turbulence_)
357  {
358  // Re-calculate viscosity
359  mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
360  mut_.correctBoundaryConditions();
361 
362  // Re-calculate thermal diffusivity
363  alphat_ = mut_/Prt_;
364  alphat_.correctBoundaryConditions();
365 
366  return;
367  }
368 
370 
371  if (mesh_.changing())
372  {
373  y_.correct();
374  }
375 
376  volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
377  volScalarField G("RASModel::G", 0.5*mag(tr(P)));
378 
379  // Update espsilon and G at the wall
380  epsilon_.boundaryField().updateCoeffs();
381 
382  // Dissipation equation
383  tmp<fvScalarMatrix> epsEqn
384  (
385  fvm::ddt(rho_, epsilon_)
386  + fvm::div(phi_, epsilon_)
387  //- fvm::laplacian(Ceps*rho_*(k_/epsilon_)*R_, epsilon_)
388  - fvm::laplacian(DepsilonEff(), epsilon_)
389  ==
390  C1_*rho_*G*epsilon_/k_
391  - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
392  );
393 
394  epsEqn().relax();
395 
396  epsEqn().boundaryManipulate(epsilon_.boundaryField());
397 
398  solve(epsEqn);
399  bound(epsilon_, epsilon0_);
400 
401 
402  // Reynolds stress equation
403 
404  const fvPatchList& patches = mesh_.boundary();
405 
406  forAll(patches, patchi)
407  {
408  const fvPatch& curPatch = patches[patchi];
409 
410  if (isA<wallFvPatch>(curPatch))
411  {
412  forAll(curPatch, facei)
413  {
414  label faceCelli = curPatch.faceCells()[facei];
415  P[faceCelli] *=
416  min(G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL), 100.0);
417  }
418  }
419  }
420 
421  volSymmTensorField reflect = C1Ref_*epsilon_/k_*R_ - C2Ref_*Clg2_*dev(P);
422 
424  (
425  fvm::ddt(rho_, R_)
426  + fvm::div(phi_, R_)
427  //- fvm::laplacian(Cs*rho_*(k_/epsilon_)*R_, R_)
428  - fvm::laplacian(DREff(), R_)
429  + fvm::Sp(Clg1_*rho_*epsilon_/k_, R_)
430  ==
431  rho_*P
432  + (2.0/3.0*(Clg1_ - 1)*I)*rho_*epsilon_
433  - Clg2_*rho_*dev(P)
434 
435  // wall reflection terms
436  + symm
437  (
438  I*((y_.n() & reflect) & y_.n())
439  - 1.5*(y_.n()*(reflect & y_.n())
440  + (y_.n() & reflect)*y_.n())
441  )*pow(Cmu_, 0.75)*rho_*pow(k_, 1.5)/(kappa_*y_*epsilon_)
442  );
443 
444  REqn().relax();
445  solve(REqn);
446 
447  R_.max
448  (
450  (
451  "zero",
452  R_.dimensions(),
453  symmTensor
454  (
455  k0_.value(), -GREAT, -GREAT,
456  k0_.value(), -GREAT,
457  k0_.value()
458  )
459  )
460  );
461 
462  k_ == 0.5*tr(R_);
463  bound(k_, k0_);
464 
465 
466  // Re-calculate turbulent viscosity
467  mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
468  mut_.correctBoundaryConditions();
469 
470  // Re-calculate thermal diffusivity
471  alphat_ = mut_/Prt_;
472  alphat_.correctBoundaryConditions();
473 
474  // Correct wall shear stresses
475 
476  forAll(patches, patchi)
477  {
478  const fvPatch& curPatch = patches[patchi];
479 
480  if (isA<wallFvPatch>(curPatch))
481  {
482  symmTensorField& Rw = R_.boundaryField()[patchi];
483 
484  const scalarField& mutw = mut_.boundaryField()[patchi];
485  const scalarField& rhow = rho_.boundaryField()[patchi];
486 
487  vectorField snGradU = U_.boundaryField()[patchi].snGrad();
488 
489  const vectorField& faceAreas
490  = mesh_.Sf().boundaryField()[patchi];
491 
492  const scalarField& magFaceAreas
493  = mesh_.magSf().boundaryField()[patchi];
494 
495  forAll(curPatch, facei)
496  {
497  // Calculate near-wall velocity gradient
498  tensor gradUw
499  = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
500 
501  // Calculate near-wall shear-stress tensor
502  tensor tauw = -(mutw[facei]/rhow[facei])*2*dev(symm(gradUw));
503 
504  // Reset the shear components of the stress tensor
505  Rw[facei].xy() = tauw.xy();
506  Rw[facei].xz() = tauw.xz();
507  Rw[facei].yz() = tauw.yz();
508  }
509  }
510  }
511 }
512 
513 
514 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
515 
516 } // End namespace RASModels
517 } // End namespace compressible
518 } // End namespace Foam
519 
520 // ************************ vim: set sw=4 sts=4 et: ************************ //