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 incompressible
37 {
38 namespace RASModels
39 {
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
44 addToRunTimeSelectionTable(RASModel, LRR, 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  Clrr1_
67  (
69  (
70  "Clrr1",
71  coeffDict_,
72  1.8
73  )
74  ),
75  Clrr2_
76  (
78  (
79  "Clrr2",
80  coeffDict_,
81  0.6
82  )
83  ),
84  C1_
85  (
87  (
88  "C1",
89  coeffDict_,
90  1.44
91  )
92  ),
93  C2_
94  (
96  (
97  "C2",
98  coeffDict_,
99  1.92
100  )
101  ),
102  Cs_
103  (
105  (
106  "Cs",
107  coeffDict_,
108  0.25
109  )
110  ),
111  Ceps_
112  (
114  (
115  "Ceps",
116  coeffDict_,
117  0.15
118  )
119  ),
120  sigmaEps_
121  (
123  (
124  "sigmaEps",
125  coeffDict_,
126  1.3
127  )
128  ),
129  couplingFactor_
130  (
132  (
133  "couplingFactor",
134  coeffDict_,
135  0.0
136  )
137  ),
138 
139  R_
140  (
141  IOobject
142  (
143  "R",
144  runTime_.timeName(),
145  mesh_,
148  ),
149  autoCreateR("R", mesh_)
150  ),
151  k_
152  (
153  IOobject
154  (
155  "k",
156  runTime_.timeName(),
157  mesh_,
160  ),
161  autoCreateK("k", mesh_)
162  ),
163  epsilon_
164  (
165  IOobject
166  (
167  "epsilon",
168  runTime_.timeName(),
169  mesh_,
172  ),
173  autoCreateEpsilon("epsilon", mesh_)
174  ),
175  nut_
176  (
177  IOobject
178  (
179  "nut",
180  runTime_.timeName(),
181  mesh_,
184  ),
185  autoCreateNut("nut", mesh_)
186  )
187 {
188  nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
189  nut_.correctBoundaryConditions();
190 
191  if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
192  {
194  (
195  "LRR::LRR"
196  "(const volVectorField& U, const surfaceScalarField& phi,"
197  "transportModel& lamTransportModel)"
198  ) << "couplingFactor = " << couplingFactor_
199  << " is not in range 0 - 1" << nl
200  << exit(FatalError);
201  }
202 
203  printCoeffs();
204 }
205 
206 
207 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
208 
210 {
212  (
214  (
215  IOobject
216  (
217  "devRhoReff",
218  runTime_.timeName(),
219  mesh_,
222  ),
223  R_ - nu()*dev(twoSymm(fvc::grad(U_)))
224  )
225  );
226 }
227 
228 
230 {
231  if (couplingFactor_.value() > 0.0)
232  {
233  return
234  (
235  fvc::div(R_ + couplingFactor_*nut_*fvc::grad(U), "div(R)")
237  (
238  (1.0 - couplingFactor_)*nut_,
239  U,
240  "laplacian(nuEff,U)"
241  )
242  - fvm::laplacian(nuEff(), U)
243  );
244  }
245  else
246  {
247  return
248  (
249  fvc::div(R_)
250  + fvc::laplacian(nut_, U, "laplacian(nuEff,U)")
251  - fvm::laplacian(nuEff(), U)
252  );
253  }
254 }
255 
256 
257 bool LRR::read()
258 {
259  if (RASModel::read())
260  {
261  Cmu_.readIfPresent(coeffDict());
262  Clrr1_.readIfPresent(coeffDict());
263  Clrr2_.readIfPresent(coeffDict());
264  C1_.readIfPresent(coeffDict());
265  C2_.readIfPresent(coeffDict());
266  Cs_.readIfPresent(coeffDict());
267  Ceps_.readIfPresent(coeffDict());
268  sigmaEps_.readIfPresent(coeffDict());
269 
270  couplingFactor_.readIfPresent(coeffDict());
271 
272  if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
273  {
274  FatalErrorIn("LRR::read()")
275  << "couplingFactor = " << couplingFactor_
276  << " is not in range 0 - 1"
277  << exit(FatalError);
278  }
279 
280  return true;
281  }
282  else
283  {
284  return false;
285  }
286 }
287 
288 
290 {
292 
293  if (!turbulence_)
294  {
295  return;
296  }
297 
299  volScalarField G("RASModel::G", 0.5*mag(tr(P)));
300 
301  // Update espsilon and G at the wall
302  epsilon_.boundaryField().updateCoeffs();
303 
304  // Dissipation equation
305  tmp<fvScalarMatrix> epsEqn
306  (
307  fvm::ddt(epsilon_)
308  + fvm::div(phi_, epsilon_)
309  - fvm::Sp(fvc::div(phi_), epsilon_)
310  //- fvm::laplacian(Ceps*(K/epsilon_)*R, epsilon_)
311  - fvm::laplacian(DepsilonEff(), epsilon_)
312  ==
313  C1_*G*epsilon_/k_
314  - fvm::Sp(C2_*epsilon_/k_, epsilon_)
315  );
316 
317  epsEqn().relax();
318 
319  epsEqn().boundaryManipulate(epsilon_.boundaryField());
320 
321  solve(epsEqn);
322  bound(epsilon_, epsilon0_);
323 
324 
325  // Reynolds stress equation
326 
327  const fvPatchList& patches = mesh_.boundary();
328 
329  forAll(patches, patchi)
330  {
331  const fvPatch& curPatch = patches[patchi];
332 
333  if (isA<wallFvPatch>(curPatch))
334  {
335  forAll(curPatch, facei)
336  {
337  label faceCelli = curPatch.faceCells()[facei];
338  P[faceCelli]
339  *= min(G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL), 1.0);
340  }
341  }
342  }
343 
344 
346  (
347  fvm::ddt(R_)
348  + fvm::div(phi_, R_)
349  - fvm::Sp(fvc::div(phi_), R_)
350  //- fvm::laplacian(Cs*(k_/epsilon_)*R_, R_)
351  - fvm::laplacian(DREff(), R_)
352  + fvm::Sp(Clrr1_*epsilon_/k_, R_)
353  ==
354  P
355  - (2.0/3.0*(1 - Clrr1_)*I)*epsilon_
356  - Clrr2_*dev(P)
357  );
358 
359  REqn().relax();
360  solve(REqn);
361 
362  R_.max
363  (
365  (
366  "zero",
367  R_.dimensions(),
368  symmTensor
369  (
370  k0_.value(), -GREAT, -GREAT,
371  k0_.value(), -GREAT,
372  k0_.value()
373  )
374  )
375  );
376 
377  k_ = 0.5*tr(R_);
378  bound(k_, k0_);
379 
380 
381  // Re-calculate viscosity
382  nut_ = Cmu_*sqr(k_)/epsilon_;
384 
385 
386  // Correct wall shear stresses
387 
388  forAll(patches, patchi)
389  {
390  const fvPatch& curPatch = patches[patchi];
391 
392  if (isA<wallFvPatch>(curPatch))
393  {
395 
396  const scalarField& nutw = nut_.boundaryField()[patchi];
397 
398  vectorField snGradU = U_.boundaryField()[patchi].snGrad();
399 
400  const vectorField& faceAreas
401  = mesh_.Sf().boundaryField()[patchi];
402 
403  const scalarField& magFaceAreas
405 
406  forAll(curPatch, facei)
407  {
408  // Calculate near-wall velocity gradient
409  tensor gradUw
410  = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
411 
412  // Calculate near-wall shear-stress tensor
413  tensor tauw = -nutw[facei]*2*symm(gradUw);
414 
415  // Reset the shear components of the stress tensor
416  Rw[facei].xy() = tauw.xy();
417  Rw[facei].xz() = tauw.xz();
418  Rw[facei].yz() = tauw.yz();
419  }
420  }
421  }
422 }
423 
424 
425 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
426 
427 } // End namespace RASModels
428 } // End namespace incompressible
429 } // End namespace Foam
430 
431 // ************************ vim: set sw=4 sts=4 et: ************************ //