FreeFOAM The Cross-Platform CFD Toolkit
LienLeschzinerLowRe.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 "LienLeschzinerLowRe.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(LienLeschzinerLowRe, 0);
44 addToRunTimeSelectionTable(RASModel, LienLeschzinerLowRe, 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  C1_
58  (
60  (
61  "C1",
62  coeffDict_,
63  1.44
64  )
65  ),
66  C2_
67  (
69  (
70  "C2",
71  coeffDict_,
72  1.92
73  )
74  ),
75  sigmak_
76  (
78  (
79  "sigmak",
80  coeffDict_,
81  1.0
82  )
83  ),
84  sigmaEps_
85  (
87  (
88  "sigmaEps",
89  coeffDict_,
90  1.3
91  )
92  ),
93  Cmu_
94  (
96  (
97  "Cmu",
98  coeffDict_,
99  0.09
100  )
101  ),
102  kappa_
103  (
105  (
106  "kappa",
107  coeffDict_,
108  0.41
109  )
110  ),
111  Am_
112  (
114  (
115  "Am",
116  coeffDict_,
117  0.016
118  )
119  ),
120  Aepsilon_
121  (
123  (
124  "Aepsilon",
125  coeffDict_,
126  0.263
127  )
128  ),
129  Amu_
130  (
132  (
133  "Amu",
134  coeffDict_,
135  0.00222
136  )
137  ),
138 
139  k_
140  (
141  IOobject
142  (
143  "k",
144  runTime_.timeName(),
145  mesh_,
148  ),
149  mesh_
150  ),
151 
152  epsilon_
153  (
154  IOobject
155  (
156  "epsilon",
157  runTime_.timeName(),
158  mesh_,
161  ),
162  mesh_
163  ),
164 
165  y_(mesh_),
166 
167  yStar_(sqrt(k_)*y_/nu() + SMALL),
168 
169  nut_
170  (
171  IOobject
172  (
173  "nut",
174  runTime_.timeName(),
175  mesh_,
178  ),
179  autoCreateLowReNut("nut", mesh_)
180  )
181 {
182  nut_ = Cmu_*(scalar(1) - exp(-Am_*yStar_))
183  /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL)*sqr(k_)
184  /(epsilon_ + epsilonSmall_);
185  nut_.correctBoundaryConditions();
186 
187  printCoeffs();
188 }
189 
190 
191 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
192 
194 {
195  volTensorField gradU = fvc::grad(U_);
196 
198  (
200  (
201  IOobject
202  (
203  "R",
204  runTime_.timeName(),
205  mesh_,
208  ),
209  ((2.0/3.0)*I)*k_ - nut_*twoSymm(gradU),
210  k_.boundaryField().types()
211  )
212  );
213 }
214 
215 
217 {
219  (
221  (
222  IOobject
223  (
224  "devRhoReff",
225  runTime_.timeName(),
226  mesh_,
229  ),
231  )
232  );
233 }
234 
235 
237 {
238  return
239  (
240  - fvm::laplacian(nuEff(), U)
241  //- (fvc::grad(U) & fvc::grad(nuEff()))
242  - fvc::div(nuEff()*fvc::grad(U)().T())
243  );
244 }
245 
246 
248 {
249  if (RASModel::read())
250  {
251  C1_.readIfPresent(coeffDict());
252  C2_.readIfPresent(coeffDict());
253  sigmak_.readIfPresent(coeffDict());
254  sigmaEps_.readIfPresent(coeffDict());
255  Cmu_.readIfPresent(coeffDict());
256  kappa_.readIfPresent(coeffDict());
257  Am_.readIfPresent(coeffDict());
258  Aepsilon_.readIfPresent(coeffDict());
259  Amu_.readIfPresent(coeffDict());
260 
261  return true;
262  }
263  else
264  {
265  return false;
266  }
267 }
268 
269 
271 {
273 
274  if (!turbulence_)
275  {
276  return;
277  }
278 
279  if (mesh_.changing())
280  {
281  y_.correct();
282  }
283 
284  scalar Cmu75 = pow(Cmu_.value(), 0.75);
285 
286  volTensorField gradU = fvc::grad(U_);
287 
288  // generation term
289  volScalarField S2 = symm(gradU) && gradU;
290 
291  yStar_ = sqrt(k_)*y_/nu() + SMALL;
292  volScalarField Rt = sqr(k_)/(nu()*epsilon_);
293 
294  volScalarField fMu =
295  (scalar(1) - exp(-Am_*yStar_))
296  /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL);
297 
298  volScalarField f2 = scalar(1) - 0.3*exp(-sqr(Rt));
299 
300  volScalarField G("RASModel::G", Cmu_*fMu*sqr(k_)/epsilon_*S2);
301 
302 
303  // Dissipation equation
304  tmp<fvScalarMatrix> epsEqn
305  (
306  fvm::ddt(epsilon_)
307  + fvm::div(phi_, epsilon_)
308  - fvm::laplacian(DepsilonEff(), epsilon_)
309  ==
310  C1_*G*epsilon_/k_
311  // E-term
312  + C2_*f2*Cmu75*sqrt(k_)
313  /(kappa_*y_*(scalar(1) - exp(-Aepsilon_*yStar_)))
314  *exp(-Amu_*sqr(yStar_))*epsilon_
315  - fvm::Sp(C2_*f2*epsilon_/k_, epsilon_)
316  );
317 
318  epsEqn().relax();
319 
322 
323  solve(epsEqn);
324  bound(epsilon_, epsilon0_);
325 
326 
327  // Turbulent kinetic energy equation
328 
330  (
331  fvm::ddt(k_)
332  + fvm::div(phi_, k_)
333  - fvm::laplacian(DkEff(), k_)
334  ==
335  G
336  - fvm::Sp(epsilon_/k_, k_)
337  );
338 
339  kEqn().relax();
340  solve(kEqn);
341  bound(k_, k0_);
342 
343 
344  // Re-calculate viscosity
345  nut_ = Cmu_*fMu*sqr(k_)/epsilon_;
346 }
347 
348 
349 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
350 
351 } // End namespace RASModels
352 } // End namespace incompressible
353 } // End namespace Foam
354 
355 // ************************ vim: set sw=4 sts=4 et: ************************ //