FreeFOAM The Cross-Platform CFD Toolkit
RNGkEpsilon.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 "RNGkEpsilon.H"
28 
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 namespace compressible
36 {
37 namespace RASModels
38 {
39 
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 
42 defineTypeNameAndDebug(RNGkEpsilon, 0);
43 addToRunTimeSelectionTable(RASModel, RNGkEpsilon, dictionary);
44 
45 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
46 
47 RNGkEpsilon::RNGkEpsilon
48 (
49  const volScalarField& rho,
50  const volVectorField& U,
51  const surfaceScalarField& phi,
52  const basicThermo& thermophysicalModel
53 )
54 :
55  RASModel(typeName, rho, U, phi, thermophysicalModel),
56 
57  Cmu_
58  (
60  (
61  "Cmu",
62  coeffDict_,
63  0.0845
64  )
65  ),
66  C1_
67  (
69  (
70  "C1",
71  coeffDict_,
72  1.42
73  )
74  ),
75  C2_
76  (
78  (
79  "C2",
80  coeffDict_,
81  1.68
82  )
83  ),
84  C3_
85  (
87  (
88  "C3",
89  coeffDict_,
90  -0.33
91  )
92  ),
93  sigmak_
94  (
96  (
97  "sigmak",
98  coeffDict_,
99  0.71942
100  )
101  ),
102  sigmaEps_
103  (
105  (
106  "sigmaEps",
107  coeffDict_,
108  0.71942
109  )
110  ),
111  Prt_
112  (
114  (
115  "Prt",
116  coeffDict_,
117  1.0
118  )
119  ),
120  eta0_
121  (
123  (
124  "eta0",
125  coeffDict_,
126  4.38
127  )
128  ),
129  beta_
130  (
132  (
133  "beta",
134  coeffDict_,
135  0.012
136  )
137  ),
138 
139  k_
140  (
141  IOobject
142  (
143  "k",
144  runTime_.timeName(),
145  mesh_,
148  ),
149  autoCreateK("k", mesh_)
150  ),
151  epsilon_
152  (
153  IOobject
154  (
155  "epsilon",
156  runTime_.timeName(),
157  mesh_,
160  ),
161  autoCreateEpsilon("epsilon", mesh_)
162  ),
163  mut_
164  (
165  IOobject
166  (
167  "mut",
168  runTime_.timeName(),
169  mesh_,
172  ),
173  autoCreateMut("mut", mesh_)
174  ),
175  alphat_
176  (
177  IOobject
178  (
179  "alphat",
180  runTime_.timeName(),
181  mesh_,
184  ),
185  autoCreateAlphat("alphat", mesh_)
186  )
187 {
188  mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
189  mut_.correctBoundaryConditions();
190 
191  alphat_ = mut_/Prt_;
192  alphat_.correctBoundaryConditions();
193 
194  printCoeffs();
195 }
196 
197 
198 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
199 
200 tmp<volSymmTensorField> RNGkEpsilon::R() const
201 {
203  (
205  (
206  IOobject
207  (
208  "R",
209  runTime_.timeName(),
210  mesh_,
213  ),
214  ((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
215  k_.boundaryField().types()
216  )
217  );
218 }
219 
220 
221 tmp<volSymmTensorField> RNGkEpsilon::devRhoReff() const
222 {
224  (
226  (
227  IOobject
228  (
229  "devRhoReff",
230  runTime_.timeName(),
231  mesh_,
234  ),
235  -muEff()*dev(twoSymm(fvc::grad(U_)))
236  )
237  );
238 }
239 
240 
241 tmp<fvVectorMatrix> RNGkEpsilon::divDevRhoReff(volVectorField& U) const
242 {
243  return
244  (
245  - fvm::laplacian(muEff(), U) - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
246  );
247 }
248 
249 
250 bool RNGkEpsilon::read()
251 {
252  if (RASModel::read())
253  {
254  Cmu_.readIfPresent(coeffDict());
255  C1_.readIfPresent(coeffDict());
256  C2_.readIfPresent(coeffDict());
257  C3_.readIfPresent(coeffDict());
258  Prt_.readIfPresent(coeffDict());
259  sigmaEps_.readIfPresent(coeffDict());
260  Prt_.readIfPresent(coeffDict());
261  eta0_.readIfPresent(coeffDict());
262  beta_.readIfPresent(coeffDict());
263 
264  return true;
265  }
266  else
267  {
268  return false;
269  }
270 }
271 
272 
274 {
275  if (!turbulence_)
276  {
277  // Re-calculate viscosity
278  mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
279  mut_.correctBoundaryConditions();
280 
281  // Re-calculate thermal diffusivity
282  alphat_ = mut_/Prt_;
283  alphat_.correctBoundaryConditions();
284 
285  return;
286  }
287 
289 
290  volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));
291 
292  if (mesh_.moving())
293  {
294  divU += fvc::div(mesh_.phi());
295  }
296 
297  tmp<volTensorField> tgradU = fvc::grad(U_);
298  volScalarField S2 = (tgradU() && dev(twoSymm(tgradU())));
299  tgradU.clear();
300 
301  volScalarField G("RASModel::G", mut_*S2);
302 
303  volScalarField eta = sqrt(mag(S2))*k_/epsilon_;
304  volScalarField eta3 = eta*sqr(eta);
305 
306  volScalarField R =
307  ((eta*(-eta/eta0_ + scalar(1)))/(beta_*eta3 + scalar(1)));
308 
309  // Update espsilon and G at the wall
310  epsilon_.boundaryField().updateCoeffs();
311 
312  // Dissipation equation
313  tmp<fvScalarMatrix> epsEqn
314  (
315  fvm::ddt(rho_, epsilon_)
316  + fvm::div(phi_, epsilon_)
317  - fvm::laplacian(DepsilonEff(), epsilon_)
318  ==
319  (C1_ - R)*G*epsilon_/k_
320  - fvm::SuSp(((2.0/3.0)*C1_ + C3_)*rho_*divU, epsilon_)
321  - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
322  );
323 
324  epsEqn().relax();
325 
326  epsEqn().boundaryManipulate(epsilon_.boundaryField());
327 
328  solve(epsEqn);
329  bound(epsilon_, epsilon0_);
330 
331 
332  // Turbulent kinetic energy equation
333 
335  (
336  fvm::ddt(rho_, k_)
337  + fvm::div(phi_, k_)
338  - fvm::laplacian(DkEff(), k_)
339  ==
340  G - fvm::SuSp(2.0/3.0*rho_*divU, k_)
341  - fvm::Sp(rho_*(epsilon_)/k_, k_)
342  );
343 
344  kEqn().relax();
345  solve(kEqn);
346  bound(k_, k0_);
347 
348 
349  // Re-calculate viscosity
350  mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
351  mut_.correctBoundaryConditions();
352 
353  // Re-calculate thermal diffusivity
354  alphat_ = mut_/Prt_;
355  alphat_.correctBoundaryConditions();
356 }
357 
358 
359 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
360 
361 } // End namespace RASModels
362 } // End namespace compressible
363 } // End namespace Foam
364 
365 // ************************ vim: set sw=4 sts=4 et: ************************ //