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