FreeFOAM The Cross-Platform CFD Toolkit
LRR.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 "LRR.H"
29 
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 namespace compressible
37 {
38 namespace RASModels
39 {
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
44 addToRunTimeSelectionTable(RASModel, LRR, dictionary);
45 
46 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
47 
48 LRR::LRR
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  Clrr1_
68  (
70  (
71  "Clrr1",
72  coeffDict_,
73  1.8
74  )
75  ),
76  Clrr2_
77  (
79  (
80  "Clrr2",
81  coeffDict_,
82  0.6
83  )
84  ),
85  C1_
86  (
88  (
89  "C1",
90  coeffDict_,
91  1.44
92  )
93  ),
94  C2_
95  (
97  (
98  "C2",
99  coeffDict_,
100  1.92
101  )
102  ),
103  Cs_
104  (
106  (
107  "Cs",
108  coeffDict_,
109  0.25
110  )
111  ),
112  Ceps_
113  (
115  (
116  "Ceps",
117  coeffDict_,
118  0.15
119  )
120  ),
121  couplingFactor_
122  (
124  (
125  "couplingFactor",
126  coeffDict_,
127  0.0
128  )
129  ),
130  sigmaR_
131  (
133  (
134  "sigmaR",
135  coeffDict_,
136  0.81967
137  )
138  ),
139  sigmaEps_
140  (
142  (
143  "sigmaEps",
144  coeffDict_,
145  1.3
146  )
147  ),
148  Prt_
149  (
151  (
152  "Prt",
153  coeffDict_,
154  1.0
155  )
156  ),
157 
158  R_
159  (
160  IOobject
161  (
162  "R",
163  runTime_.timeName(),
164  mesh_,
167  ),
168  autoCreateR("R", mesh_)
169  ),
170  k_
171  (
172  IOobject
173  (
174  "k",
175  runTime_.timeName(),
176  mesh_,
179  ),
180  autoCreateK("k", mesh_)
181  ),
182  epsilon_
183  (
184  IOobject
185  (
186  "epsilon",
187  runTime_.timeName(),
188  mesh_,
191  ),
192  autoCreateEpsilon("epsilon", mesh_)
193  ),
194  mut_
195  (
196  IOobject
197  (
198  "mut",
199  runTime_.timeName(),
200  mesh_,
203  ),
204  autoCreateMut("mut", mesh_)
205  ),
206  alphat_
207  (
208  IOobject
209  (
210  "alphat",
211  runTime_.timeName(),
212  mesh_,
215  ),
216  autoCreateAlphat("alphat", mesh_)
217  )
218 {
219  if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
220  {
222  (
223  "LRR::LRR"
224  "( const volScalarField&, const volVectorField&"
225  ", const surfaceScalarField&, incompressibleTransportModel&)"
226  ) << "couplingFactor = " << couplingFactor_
227  << " is not in range 0 - 1" << nl
228  << exit(FatalError);
229  }
230 
231  mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
232  mut_.correctBoundaryConditions();
233 
234  alphat_ = mut_/Prt_;
235  alphat_.correctBoundaryConditions();
236 
237  printCoeffs();
238 }
239 
240 
241 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
242 
243 tmp<volSymmTensorField> LRR::devRhoReff() const
244 {
246  (
248  (
249  IOobject
250  (
251  "devRhoReff",
252  runTime_.timeName(),
253  mesh_,
256  ),
257  rho_*R_ - mu()*dev(twoSymm(fvc::grad(U_)))
258  )
259  );
260 }
261 
262 
263 tmp<fvVectorMatrix> LRR::divDevRhoReff(volVectorField& U) const
264 {
265  if (couplingFactor_.value() > 0.0)
266  {
267  return
268  (
269  fvc::div(rho_*R_ + couplingFactor_*mut_*fvc::grad(U))
270  + fvc::laplacian((1.0 - couplingFactor_)*mut_, U)
271  - fvm::laplacian(muEff(), U)
272  - fvc::div(mu()*dev2(fvc::grad(U)().T()))
273  );
274  }
275  else
276  {
277  return
278  (
279  fvc::div(rho_*R_)
280  + fvc::laplacian(mut_, U)
281  - fvm::laplacian(muEff(), U)
282  - fvc::div(mu()*dev2(fvc::grad(U)().T()))
283  );
284  }
285 }
286 
287 
288 bool LRR::read()
289 {
290  if (RASModel::read())
291  {
292  Cmu_.readIfPresent(coeffDict());
293  Clrr1_.readIfPresent(coeffDict());
294  Clrr2_.readIfPresent(coeffDict());
295  C1_.readIfPresent(coeffDict());
296  C2_.readIfPresent(coeffDict());
297  Cs_.readIfPresent(coeffDict());
298  Ceps_.readIfPresent(coeffDict());
299  sigmaR_.readIfPresent(coeffDict());
300  sigmaEps_.readIfPresent(coeffDict());
301  Prt_.readIfPresent(coeffDict());
302  couplingFactor_.readIfPresent(coeffDict());
303 
304  if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
305  {
306  FatalErrorIn("LRR::read()")
307  << "couplingFactor = " << couplingFactor_
308  << " is not in range 0 - 1" << nl
309  << exit(FatalError);
310  }
311 
312  return true;
313  }
314  else
315  {
316  return false;
317  }
318 }
319 
320 
322 {
323  if (!turbulence_)
324  {
325  // Re-calculate viscosity
326  mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
327  mut_.correctBoundaryConditions();
328 
329  // Re-calculate thermal diffusivity
330  alphat_ = mut_/Prt_;
331  alphat_.correctBoundaryConditions();
332 
333  return;
334  }
335 
337 
338  volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
339  volScalarField G("RASModel::G", 0.5*mag(tr(P)));
340 
341  // Update espsilon and G at the wall
342  epsilon_.boundaryField().updateCoeffs();
343 
344  // Dissipation equation
345  tmp<fvScalarMatrix> epsEqn
346  (
347  fvm::ddt(rho_, epsilon_)
348  + fvm::div(phi_, epsilon_)
349  //- fvm::laplacian(Ceps*rho_*(k_/epsilon_)*R_, epsilon_)
350  - fvm::laplacian(DepsilonEff(), epsilon_)
351  ==
352  C1_*rho_*G*epsilon_/k_
353  - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
354  );
355 
356  epsEqn().relax();
357 
358  epsEqn().boundaryManipulate(epsilon_.boundaryField());
359 
360  solve(epsEqn);
361  bound(epsilon_, epsilon0_);
362 
363 
364  // Reynolds stress equation
365 
366  const fvPatchList& patches = mesh_.boundary();
367 
368  forAll(patches, patchi)
369  {
370  const fvPatch& curPatch = patches[patchi];
371 
372  if (isA<wallFvPatch>(curPatch))
373  {
374  forAll(curPatch, facei)
375  {
376  label faceCelli = curPatch.faceCells()[facei];
377  P[faceCelli]
378  *= min(G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL), 100.0);
379  }
380  }
381  }
382 
383 
385  (
386  fvm::ddt(rho_, R_)
387  + fvm::div(phi_, R_)
388  //- fvm::laplacian(Cs*rho_*(k_/epsilon_)*R_, R_)
389  - fvm::laplacian(DREff(), R_)
390  + fvm::Sp(Clrr1_*rho_*epsilon_/k_, R_)
391  ==
392  rho_*P
393  - (2.0/3.0*(1 - Clrr1_)*I)*rho_*epsilon_
394  - Clrr2_*rho_*dev(P)
395  );
396 
397  REqn().relax();
398  solve(REqn);
399 
400  R_.max
401  (
403  (
404  "zero",
405  R_.dimensions(),
406  symmTensor
407  (
408  k0_.value(), -GREAT, -GREAT,
409  k0_.value(), -GREAT,
410  k0_.value()
411  )
412  )
413  );
414 
415  k_ = 0.5*tr(R_);
416  bound(k_, k0_);
417 
418 
419  // Re-calculate viscosity
420  mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
421  mut_.correctBoundaryConditions();
422 
423  // Re-calculate thermal diffusivity
424  alphat_ = mut_/Prt_;
425  alphat_.correctBoundaryConditions();
426 
427 
428  // Correct wall shear stresses
429 
430  forAll(patches, patchi)
431  {
432  const fvPatch& curPatch = patches[patchi];
433 
434  if (isA<wallFvPatch>(curPatch))
435  {
436  symmTensorField& Rw = R_.boundaryField()[patchi];
437 
438  const scalarField& rhow = rho_.boundaryField()[patchi];
439  const scalarField& mutw = mut_.boundaryField()[patchi];
440 
441  vectorField snGradU = U_.boundaryField()[patchi].snGrad();
442 
443  const vectorField& faceAreas
444  = mesh_.Sf().boundaryField()[patchi];
445 
446  const scalarField& magFaceAreas
447  = mesh_.magSf().boundaryField()[patchi];
448 
449  forAll(curPatch, facei)
450  {
451  // Calculate near-wall velocity gradient
452  tensor gradUw
453  = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
454 
455  // Calculate near-wall shear-stress tensor
456  tensor tauw = -(mutw[facei]/rhow[facei])*2*dev(symm(gradUw));
457 
458  // Reset the shear components of the stress tensor
459  Rw[facei].xy() = tauw.xy();
460  Rw[facei].xz() = tauw.xz();
461  Rw[facei].yz() = tauw.yz();
462  }
463  }
464  }
465 }
466 
467 
468 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
469 
470 } // End namespace RASModels
471 } // End namespace compressible
472 } // End namespace Foam
473 
474 // ************************ vim: set sw=4 sts=4 et: ************************ //