FreeFOAM The Cross-Platform CFD Toolkit
kEpsilon.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 "kEpsilon.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(kEpsilon, 0);
43 addToRunTimeSelectionTable(RASModel, kEpsilon, dictionary);
44 
45 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
46 
47 kEpsilon::kEpsilon
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.09
64  )
65  ),
66  C1_
67  (
69  (
70  "C1",
71  coeffDict_,
72  1.44
73  )
74  ),
75  C2_
76  (
78  (
79  "C2",
80  coeffDict_,
81  1.92
82  )
83  ),
84  C3_
85  (
87  (
88  "C3",
89  coeffDict_,
90  -0.33
91  )
92  ),
93  sigmak_
94  (
96  (
97  "sigmak",
98  coeffDict_,
99  1.0
100  )
101  ),
102  sigmaEps_
103  (
105  (
106  "sigmaEps",
107  coeffDict_,
108  1.3
109  )
110  ),
111  Prt_
112  (
114  (
115  "Prt",
116  coeffDict_,
117  1.0
118  )
119  ),
120 
121  k_
122  (
123  IOobject
124  (
125  "k",
126  runTime_.timeName(),
127  mesh_,
130  ),
131  autoCreateK("k", mesh_)
132  ),
133  epsilon_
134  (
135  IOobject
136  (
137  "epsilon",
138  runTime_.timeName(),
139  mesh_,
142  ),
143  autoCreateEpsilon("epsilon", mesh_)
144  ),
145  mut_
146  (
147  IOobject
148  (
149  "mut",
150  runTime_.timeName(),
151  mesh_,
154  ),
155  autoCreateMut("mut", mesh_)
156  ),
157  alphat_
158  (
159  IOobject
160  (
161  "alphat",
162  runTime_.timeName(),
163  mesh_,
166  ),
167  autoCreateAlphat("alphat", mesh_)
168  )
169 {
170  mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
171  mut_.correctBoundaryConditions();
172 
173  alphat_ = mut_/Prt_;
174  alphat_.correctBoundaryConditions();
175 
176  printCoeffs();
177 }
178 
179 
180 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
181 
182 tmp<volSymmTensorField> kEpsilon::R() const
183 {
185  (
187  (
188  IOobject
189  (
190  "R",
191  runTime_.timeName(),
192  mesh_,
195  ),
196  ((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
197  k_.boundaryField().types()
198  )
199  );
200 }
201 
202 
203 tmp<volSymmTensorField> kEpsilon::devRhoReff() const
204 {
206  (
208  (
209  IOobject
210  (
211  "devRhoReff",
212  runTime_.timeName(),
213  mesh_,
216  ),
217  -muEff()*dev(twoSymm(fvc::grad(U_)))
218  )
219  );
220 }
221 
222 
223 tmp<fvVectorMatrix> kEpsilon::divDevRhoReff(volVectorField& U) const
224 {
225  return
226  (
227  - fvm::laplacian(muEff(), U)
228  - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
229  );
230 }
231 
232 
233 bool kEpsilon::read()
234 {
235  if (RASModel::read())
236  {
237  Cmu_.readIfPresent(coeffDict());
238  C1_.readIfPresent(coeffDict());
239  C2_.readIfPresent(coeffDict());
240  C3_.readIfPresent(coeffDict());
241  sigmak_.readIfPresent(coeffDict());
242  sigmaEps_.readIfPresent(coeffDict());
243  Prt_.readIfPresent(coeffDict());
244 
245  return true;
246  }
247  else
248  {
249  return false;
250  }
251 }
252 
253 
255 {
256  if (!turbulence_)
257  {
258  // Re-calculate viscosity
259  mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
260  mut_.correctBoundaryConditions();
261 
262  // Re-calculate thermal diffusivity
263  alphat_ = mut_/Prt_;
264  alphat_.correctBoundaryConditions();
265 
266  return;
267  }
268 
270 
271  volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));
272 
273  if (mesh_.moving())
274  {
275  divU += fvc::div(mesh_.phi());
276  }
277 
278  tmp<volTensorField> tgradU = fvc::grad(U_);
279  volScalarField G("RASModel::G", mut_*(tgradU() && dev(twoSymm(tgradU()))));
280  tgradU.clear();
281 
282  // Update espsilon and G at the wall
283  epsilon_.boundaryField().updateCoeffs();
284 
285  // Dissipation equation
286  tmp<fvScalarMatrix> epsEqn
287  (
288  fvm::ddt(rho_, epsilon_)
289  + fvm::div(phi_, epsilon_)
290  - fvm::laplacian(DepsilonEff(), epsilon_)
291  ==
292  C1_*G*epsilon_/k_
293  - fvm::SuSp(((2.0/3.0)*C1_ + C3_)*rho_*divU, epsilon_)
294  - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
295  );
296 
297  epsEqn().relax();
298 
299  epsEqn().boundaryManipulate(epsilon_.boundaryField());
300 
301  solve(epsEqn);
302  bound(epsilon_, epsilon0_);
303 
304 
305  // Turbulent kinetic energy equation
306 
308  (
309  fvm::ddt(rho_, k_)
310  + fvm::div(phi_, k_)
311  - fvm::laplacian(DkEff(), k_)
312  ==
313  G
314  - fvm::SuSp((2.0/3.0)*rho_*divU, k_)
315  - fvm::Sp(rho_*epsilon_/k_, k_)
316  );
317 
318  kEqn().relax();
319  solve(kEqn);
320  bound(k_, k0_);
321 
322 
323  // Re-calculate viscosity
324  mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
325  mut_.correctBoundaryConditions();
326 
327  // Re-calculate thermal diffusivity
328  alphat_ = mut_/Prt_;
329  alphat_.correctBoundaryConditions();
330 }
331 
332 
333 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
334 
335 } // End namespace RASModels
336 } // End namespace compressible
337 } // End namespace Foam
338 
339 // ************************ vim: set sw=4 sts=4 et: ************************ //