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 incompressible
37 {
38 namespace RASModels
39 {
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
43 defineTypeNameAndDebug(LaunderGibsonRSTM, 0);
44 addToRunTimeSelectionTable(RASModel, LaunderGibsonRSTM, dictionary);
45 
46 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
47 
49 (
50  const volVectorField& U,
51  const surfaceScalarField& phi,
52  transportModel& lamTransportModel
53 )
54 :
55  RASModel(typeName, U, phi, lamTransportModel),
56 
57  Cmu_
58  (
60  (
61  "Cmu",
62  coeffDict_,
63  0.09
64  )
65  ),
66  kappa_
67  (
69  (
70  "kappa",
71  coeffDict_,
72  0.41
73  )
74  ),
75  Clg1_
76  (
78  (
79  "Clg1",
80  coeffDict_,
81  1.8
82  )
83  ),
84  Clg2_
85  (
87  (
88  "Clg2",
89  coeffDict_,
90  0.6
91  )
92  ),
93  C1_
94  (
96  (
97  "C1",
98  coeffDict_,
99  1.44
100  )
101  ),
102  C2_
103  (
105  (
106  "C2",
107  coeffDict_,
108  1.92
109  )
110  ),
111  Cs_
112  (
114  (
115  "Cs",
116  coeffDict_,
117  0.25
118  )
119  ),
120  Ceps_
121  (
123  (
124  "Ceps",
125  coeffDict_,
126  0.15
127  )
128  ),
129  sigmaR_
130  (
132  (
133  "sigmaR",
134  coeffDict_,
135  0.81967
136  )
137  ),
138  sigmaEps_
139  (
141  (
142  "sigmaEps",
143  coeffDict_,
144  1.3
145  )
146  ),
147  C1Ref_
148  (
150  (
151  "C1Ref",
152  coeffDict_,
153  0.5
154  )
155  ),
156  C2Ref_
157  (
159  (
160  "C2Ref",
161  coeffDict_,
162  0.3
163  )
164  ),
165  couplingFactor_
166  (
168  (
169  "couplingFactor",
170  coeffDict_,
171  0.0
172  )
173  ),
174 
175  yr_(mesh_),
176 
177  R_
178  (
179  IOobject
180  (
181  "R",
182  runTime_.timeName(),
183  mesh_,
186  ),
187  autoCreateR("R", mesh_)
188  ),
189  k_
190  (
191  IOobject
192  (
193  "k",
194  runTime_.timeName(),
195  mesh_,
198  ),
199  autoCreateK("k", mesh_)
200  ),
201  epsilon_
202  (
203  IOobject
204  (
205  "epsilon",
206  runTime_.timeName(),
207  mesh_,
210  ),
211  autoCreateEpsilon("epsilon", mesh_)
212  ),
213  nut_
214  (
215  IOobject
216  (
217  "nut",
218  runTime_.timeName(),
219  mesh_,
222  ),
223  autoCreateNut("nut", mesh_)
224  )
225 {
226  nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
227  nut_.correctBoundaryConditions();
228 
229  if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
230  {
232  (
233  "LaunderGibsonRSTM::LaunderGibsonRSTM"
234  "(const volVectorField& U, const surfaceScalarField& phi,"
235  "transportModel& lamTransportModel)"
236  ) << "couplingFactor = " << couplingFactor_
237  << " is not in range 0 - 1" << nl
238  << exit(FatalError);
239  }
240 
241  printCoeffs();
242 }
243 
244 
245 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
246 
248 {
250  (
252  (
253  IOobject
254  (
255  "devRhoReff",
256  runTime_.timeName(),
257  mesh_,
260  ),
261  R_ - nu()*dev(twoSymm(fvc::grad(U_)))
262  )
263  );
264 }
265 
266 
268 {
269  if (couplingFactor_.value() > 0.0)
270  {
271  return
272  (
273  fvc::div(R_ + couplingFactor_*nut_*fvc::grad(U), "div(R)")
274  + fvc::laplacian((1.0-couplingFactor_)*nut_, U, "laplacian(nuEff,U)")
275  - fvm::laplacian(nuEff(), U)
276  );
277  }
278  else
279  {
280  return
281  (
282  fvc::div(R_)
283  + fvc::laplacian(nut_, U, "laplacian(nuEff,U)")
284  - fvm::laplacian(nuEff(), U)
285  );
286  }
287 }
288 
289 
291 {
292  if (RASModel::read())
293  {
294  Cmu_.readIfPresent(coeffDict());
295  Clg1_.readIfPresent(coeffDict());
296  Clg2_.readIfPresent(coeffDict());
297  C1_.readIfPresent(coeffDict());
298  C2_.readIfPresent(coeffDict());
299  Cs_.readIfPresent(coeffDict());
300  Ceps_.readIfPresent(coeffDict());
301  sigmaR_.readIfPresent(coeffDict());
302  sigmaEps_.readIfPresent(coeffDict());
303  C1Ref_.readIfPresent(coeffDict());
304  C2Ref_.readIfPresent(coeffDict());
305 
306  couplingFactor_.readIfPresent(coeffDict());
307 
308  if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
309  {
310  FatalErrorIn("LaunderGibsonRSTM::read()")
311  << "couplingFactor = " << couplingFactor_
312  << " is not in range 0 - 1"
313  << exit(FatalError);
314  }
315 
316  return true;
317  }
318  else
319  {
320  return false;
321  }
322 }
323 
324 
326 {
328 
329  if (!turbulence_)
330  {
331  return;
332  }
333 
334  if (mesh_.changing())
335  {
336  yr_.correct();
337  }
338 
340  volScalarField G("RASModel::G", 0.5*mag(tr(P)));
341 
342  // Update espsilon and G at the wall
343  epsilon_.boundaryField().updateCoeffs();
344 
345  // Dissipation equation
346  tmp<fvScalarMatrix> epsEqn
347  (
348  fvm::ddt(epsilon_)
349  + fvm::div(phi_, epsilon_)
350  - fvm::Sp(fvc::div(phi_), epsilon_)
351  //- fvm::laplacian(Ceps*(k_/epsilon_)*R_, epsilon_)
352  - fvm::laplacian(DepsilonEff(), epsilon_)
353  ==
354  C1_*G*epsilon_/k_
355  - fvm::Sp(C2_*epsilon_/k_, epsilon_)
356  );
357 
358  epsEqn().relax();
359 
360  epsEqn().boundaryManipulate(epsilon_.boundaryField());
361 
362  solve(epsEqn);
363  bound(epsilon_, epsilon0_);
364 
365 
366  // Reynolds stress equation
367 
368  const fvPatchList& patches = mesh_.boundary();
369 
370  forAll(patches, patchi)
371  {
372  const fvPatch& curPatch = patches[patchi];
373 
374  if (isA<wallFvPatch>(curPatch))
375  {
376  forAll(curPatch, facei)
377  {
378  label faceCelli = curPatch.faceCells()[facei];
379  P[faceCelli] *=
380  min(G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL), 1.0);
381  }
382  }
383  }
384 
385  volSymmTensorField reflect = C1Ref_*epsilon_/k_*R_ - C2Ref_*Clg2_*dev(P);
386 
388  (
389  fvm::ddt(R_)
390  + fvm::div(phi_, R_)
391  - fvm::Sp(fvc::div(phi_), R_)
392  //- fvm::laplacian(Cs*(k_/epsilon_)*R_, R_)
393  - fvm::laplacian(DREff(), R_)
394  + fvm::Sp(Clg1_*epsilon_/k_, R_)
395  ==
396  P
397  + (2.0/3.0*(Clg1_ - 1)*I)*epsilon_
398  - Clg2_*dev(P)
399 
400  // wall reflection terms
401  + symm
402  (
403  I*((yr_.n() & reflect) & yr_.n())
404  - 1.5*(yr_.n()*(reflect & yr_.n())
405  + (yr_.n() & reflect)*yr_.n())
406  )*pow(Cmu_, 0.75)*pow(k_, 1.5)/(kappa_*yr_*epsilon_)
407  );
408 
409  REqn().relax();
410  solve(REqn);
411 
412  R_.max
413  (
415  (
416  "zero",
417  R_.dimensions(),
418  symmTensor
419  (
420  k0_.value(), -GREAT, -GREAT,
421  k0_.value(), -GREAT,
422  k0_.value()
423  )
424  )
425  );
426 
427  k_ == 0.5*tr(R_);
428  bound(k_, k0_);
429 
430 
431  // Re-calculate turbulent viscosity
432  nut_ = Cmu_*sqr(k_)/epsilon_;
434 
435 
436  // Correct wall shear stresses
437 
438  forAll(patches, patchi)
439  {
440  const fvPatch& curPatch = patches[patchi];
441 
442  if (isA<wallFvPatch>(curPatch))
443  {
445 
446  const scalarField& nutw = nut_.boundaryField()[patchi];
447 
448  vectorField snGradU = U_.boundaryField()[patchi].snGrad();
449 
450  const vectorField& faceAreas
451  = mesh_.Sf().boundaryField()[patchi];
452 
453  const scalarField& magFaceAreas
455 
456  forAll(curPatch, facei)
457  {
458  // Calculate near-wall velocity gradient
459  tensor gradUw
460  = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
461 
462  // Calculate near-wall shear-stress tensor
463  tensor tauw = -nutw[facei]*2*symm(gradUw);
464 
465  // Reset the shear components of the stress tensor
466  Rw[facei].xy() = tauw.xy();
467  Rw[facei].xz() = tauw.xz();
468  Rw[facei].yz() = tauw.yz();
469  }
470  }
471  }
472 }
473 
474 
475 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
476 
477 } // End namespace RASModels
478 } // End namespace incompressible
479 } // End namespace Foam
480 
481 // ************************ vim: set sw=4 sts=4 et: ************************ //