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 incompressible
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 
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 
89 (
90  const volVectorField& U,
91  const surfaceScalarField& phi,
92  transportModel& lamTransportModel
93 )
94 :
95  RASModel(typeName, U, phi, lamTransportModel),
96 
97  Cmu_
98  (
100  (
101  "Cmu",
102  coeffDict_,
103  0.09
104  )
105  ),
106  A0_
107  (
109  (
110  "A0",
111  coeffDict_,
112  4.0
113  )
114  ),
115  C2_
116  (
118  (
119  "C2",
120  coeffDict_,
121  1.9
122  )
123  ),
124  sigmak_
125  (
127  (
128  "sigmak",
129  coeffDict_,
130  1.0
131  )
132  ),
133  sigmaEps_
134  (
136  (
137  "sigmaEps",
138  coeffDict_,
139  1.2
140  )
141  ),
142 
143  k_
144  (
145  IOobject
146  (
147  "k",
148  runTime_.timeName(),
149  mesh_,
152  ),
153  autoCreateK("k", mesh_)
154  ),
155  epsilon_
156  (
157  IOobject
158  (
159  "epsilon",
160  runTime_.timeName(),
161  mesh_,
164  ),
165  autoCreateEpsilon("epsilon", mesh_)
166  ),
167  nut_
168  (
169  IOobject
170  (
171  "nut",
172  runTime_.timeName(),
173  mesh_,
176  ),
177  autoCreateNut("nut", mesh_)
178  )
179 {
180  bound(k_, k0_);
181  bound(epsilon_, epsilon0_);
182 
183  nut_ = rCmu(fvc::grad(U_))*sqr(k_)/(epsilon_ + epsilonSmall_);
184  nut_.correctBoundaryConditions();
185 
186  printCoeffs();
187 }
188 
189 
190 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
191 
193 {
195  (
197  (
198  IOobject
199  (
200  "R",
201  runTime_.timeName(),
202  mesh_,
205  ),
206  ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
207  k_.boundaryField().types()
208  )
209  );
210 }
211 
212 
214 {
216  (
218  (
219  IOobject
220  (
221  "devRhoReff",
222  runTime_.timeName(),
223  mesh_,
226  ),
228  )
229  );
230 }
231 
232 
234 {
235  return
236  (
237  - fvm::laplacian(nuEff(), U)
238  - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
239  );
240 }
241 
242 
244 {
245  if (RASModel::read())
246  {
247  Cmu_.readIfPresent(coeffDict());
248  A0_.readIfPresent(coeffDict());
249  C2_.readIfPresent(coeffDict());
250  sigmak_.readIfPresent(coeffDict());
251  sigmaEps_.readIfPresent(coeffDict());
252 
253  return true;
254  }
255  else
256  {
257  return false;
258  }
259 }
260 
261 
263 {
265 
266  if (!turbulence_)
267  {
268  return;
269  }
270 
271  volTensorField gradU = fvc::grad(U_);
272  volScalarField S2 = 2*magSqr(dev(symm(gradU)));
273  volScalarField magS = sqrt(S2);
274 
275  volScalarField eta = magS*k_/epsilon_;
276  volScalarField C1 = max(eta/(scalar(5) + eta), scalar(0.43));
277 
278  volScalarField G("RASModel::G", nut_*S2);
279 
280  // Update espsilon and G at the wall
281  epsilon_.boundaryField().updateCoeffs();
282 
283 
284  // Dissipation equation
285  tmp<fvScalarMatrix> epsEqn
286  (
287  fvm::ddt(epsilon_)
288  + fvm::div(phi_, epsilon_)
289  - fvm::Sp(fvc::div(phi_), epsilon_)
290  - fvm::laplacian(DepsilonEff(), epsilon_)
291  ==
292  C1*magS*epsilon_
293  - fvm::Sp
294  (
295  C2_*epsilon_/(k_ + sqrt(nu()*epsilon_)),
296  epsilon_
297  )
298  );
299 
300  epsEqn().relax();
301 
302  epsEqn().boundaryManipulate(epsilon_.boundaryField());
303 
304  solve(epsEqn);
305  bound(epsilon_, epsilon0_);
306 
307 
308  // Turbulent kinetic energy equation
310  (
311  fvm::ddt(k_)
312  + fvm::div(phi_, k_)
313  - fvm::Sp(fvc::div(phi_), k_)
314  - fvm::laplacian(DkEff(), k_)
315  ==
316  G - fvm::Sp(epsilon_/k_, k_)
317  );
318 
319  kEqn().relax();
320  solve(kEqn);
321  bound(k_, k0_);
322 
323 
324  // Re-calculate viscosity
325  nut_ = rCmu(gradU, S2, magS)*sqr(k_)/epsilon_;
327 }
328 
329 
330 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
331 
332 } // End namespace RASModels
333 } // End namespace incompressible
334 } // End namespace Foam
335 
336 // ************************ vim: set sw=4 sts=4 et: ************************ //