FreeFOAM The Cross-Platform CFD Toolkit
realizableKE.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 "realizableKE.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(realizableKE, 0);
43 addToRunTimeSelectionTable(RASModel, realizableKE, dictionary);
44 
45 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
46 
47 tmp<volScalarField> realizableKE::rCmu
48 (
49  const volTensorField& gradU,
50  const volScalarField& S2,
51  const volScalarField& magS
52 )
53 {
54  tmp<volSymmTensorField> tS = dev(symm(gradU));
55  const volSymmTensorField& S = tS();
56 
57  volScalarField W =
58  (2*sqrt(2.0))*((S&S)&&S)
59  /(
60  magS*S2
61  + dimensionedScalar("small", dimensionSet(0, 0, -3, 0, 0), SMALL)
62  );
63 
64  tS.clear();
65 
66  volScalarField phis =
67  (1.0/3.0)*acos(min(max(sqrt(6.0)*W, -scalar(1)), scalar(1)));
68  volScalarField As = sqrt(6.0)*cos(phis);
69  volScalarField Us = sqrt(S2/2.0 + magSqr(skew(gradU)));
70 
71  return 1.0/(A0_ + As*Us*k_/(epsilon_ + epsilonSmall_));
72 }
73 
74 
75 tmp<volScalarField> realizableKE::rCmu
76 (
77  const volTensorField& gradU
78 )
79 {
80  volScalarField S2 = 2*magSqr(dev(symm(gradU)));
81  volScalarField magS = sqrt(S2);
82  return rCmu(gradU, S2, magS);
83 }
84 
85 
86 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
87 
88 realizableKE::realizableKE
89 (
90  const volScalarField& rho,
91  const volVectorField& U,
92  const surfaceScalarField& phi,
93  const basicThermo& thermophysicalModel
94 )
95 :
96  RASModel(typeName, rho, U, phi, thermophysicalModel),
97 
98  Cmu_
99  (
101  (
102  "Cmu",
103  coeffDict_,
104  0.09
105  )
106  ),
107  A0_
108  (
110  (
111  "A0",
112  coeffDict_,
113  4.0
114  )
115  ),
116  C2_
117  (
119  (
120  "C2",
121  coeffDict_,
122  1.9
123  )
124  ),
125  sigmak_
126  (
128  (
129  "sigmak",
130  coeffDict_,
131  1.0
132  )
133  ),
134  sigmaEps_
135  (
137  (
138  "sigmaEps",
139  coeffDict_,
140  1.2
141  )
142  ),
143  Prt_
144  (
146  (
147  "Prt",
148  coeffDict_,
149  1.0
150  )
151  ),
152 
153  k_
154  (
155  IOobject
156  (
157  "k",
158  runTime_.timeName(),
159  mesh_,
162  ),
163  autoCreateK("k", mesh_)
164  ),
165  epsilon_
166  (
167  IOobject
168  (
169  "epsilon",
170  runTime_.timeName(),
171  mesh_,
174  ),
175  autoCreateEpsilon("epsilon", mesh_)
176  ),
177  mut_
178  (
179  IOobject
180  (
181  "mut",
182  runTime_.timeName(),
183  mesh_,
186  ),
187  autoCreateMut("mut", mesh_)
188  ),
189  alphat_
190  (
191  IOobject
192  (
193  "alphat",
194  runTime_.timeName(),
195  mesh_,
198  ),
199  autoCreateAlphat("alphat", mesh_)
200  )
201 {
202  bound(k_, k0_);
203  bound(epsilon_, epsilon0_);
204 
205  mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
206  mut_.correctBoundaryConditions();
207 
208  alphat_ = mut_/Prt_;
209  alphat_.correctBoundaryConditions();
210 
211  printCoeffs();
212 }
213 
214 
215 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
216 
217 tmp<volSymmTensorField> realizableKE::R() const
218 {
220  (
222  (
223  IOobject
224  (
225  "R",
226  runTime_.timeName(),
227  mesh_,
230  ),
231  ((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
232  k_.boundaryField().types()
233  )
234  );
235 }
236 
237 
238 tmp<volSymmTensorField> realizableKE::devRhoReff() const
239 {
241  (
243  (
244  IOobject
245  (
246  "devRhoReff",
247  runTime_.timeName(),
248  mesh_,
251  ),
252  -muEff()*dev(twoSymm(fvc::grad(U_)))
253  )
254  );
255 }
256 
257 
258 tmp<fvVectorMatrix> realizableKE::divDevRhoReff(volVectorField& U) const
259 {
260  return
261  (
262  - fvm::laplacian(muEff(), U) - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
263  );
264 }
265 
266 
267 bool realizableKE::read()
268 {
269  if (RASModel::read())
270  {
271  Cmu_.readIfPresent(coeffDict());
272  A0_.readIfPresent(coeffDict());
273  C2_.readIfPresent(coeffDict());
274  sigmak_.readIfPresent(coeffDict());
275  sigmaEps_.readIfPresent(coeffDict());
276  Prt_.readIfPresent(coeffDict());
277 
278  return true;
279  }
280  else
281  {
282  return false;
283  }
284 }
285 
286 
288 {
289  if (!turbulence_)
290  {
291  // Re-calculate viscosity
292  mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/epsilon_;
293  mut_.correctBoundaryConditions();
294 
295  // Re-calculate thermal diffusivity
296  alphat_ = mut_/Prt_;
297  alphat_.correctBoundaryConditions();
298 
299  return;
300  }
301 
303 
304  volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));
305 
306  if (mesh_.moving())
307  {
308  divU += fvc::div(mesh_.phi());
309  }
310 
311  volTensorField gradU = fvc::grad(U_);
312  volScalarField S2 = 2*magSqr(dev(symm(gradU)));
313  volScalarField magS = sqrt(S2);
314 
315  volScalarField eta = magS*k_/epsilon_;
316  volScalarField C1 = max(eta/(scalar(5) + eta), scalar(0.43));
317 
318  volScalarField G("RASModel::G", mut_*(gradU && dev(twoSymm(gradU))));
319 
320  // Update espsilon and G at the wall
321  epsilon_.boundaryField().updateCoeffs();
322 
323  // Dissipation equation
324  tmp<fvScalarMatrix> epsEqn
325  (
326  fvm::ddt(rho_, epsilon_)
327  + fvm::div(phi_, epsilon_)
328  - fvm::laplacian(DepsilonEff(), epsilon_)
329  ==
330  C1*rho_*magS*epsilon_
331  - fvm::Sp
332  (
333  C2_*rho_*epsilon_/(k_ + sqrt((mu()/rho_)*epsilon_)),
334  epsilon_
335  )
336  );
337 
338  epsEqn().relax();
339 
340  epsEqn().boundaryManipulate(epsilon_.boundaryField());
341 
342  solve(epsEqn);
343  bound(epsilon_, epsilon0_);
344 
345 
346  // Turbulent kinetic energy equation
347 
349  (
350  fvm::ddt(rho_, k_)
351  + fvm::div(phi_, k_)
352  - fvm::laplacian(DkEff(), k_)
353  ==
354  G - fvm::SuSp(2.0/3.0*rho_*divU, k_)
355  - fvm::Sp(rho_*(epsilon_)/k_, k_)
356  );
357 
358  kEqn().relax();
359  solve(kEqn);
360  bound(k_, k0_);
361 
362  // Re-calculate viscosity
363  mut_ = rCmu(gradU, S2, magS)*rho_*sqr(k_)/epsilon_;
364  mut_.correctBoundaryConditions();
365 
366  // Re-calculate thermal diffusivity
367  alphat_ = mut_/Prt_;
368  alphat_.correctBoundaryConditions();
369 }
370 
371 
372 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
373 
374 } // End namespace RASModels
375 } // End namespace compressible
376 } // End namespace Foam
377 
378 // ************************ vim: set sw=4 sts=4 et: ************************ //