FreeFOAM The Cross-Platform CFD Toolkit
NonlinearKEShih.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 "NonlinearKEShih.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 namespace incompressible
35 {
36 namespace RASModels
37 {
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
41 defineTypeNameAndDebug(NonlinearKEShih, 0);
42 addToRunTimeSelectionTable(RASModel, NonlinearKEShih, dictionary);
43 
44 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
45 
47 (
48  const volVectorField& U,
49  const surfaceScalarField& phi,
50  transportModel& lamTransportModel
51 )
52 :
53  RASModel(typeName, U, phi, lamTransportModel),
54 
55  C1_
56  (
58  (
59  "C1",
60  coeffDict_,
61  1.44
62  )
63  ),
64  C2_
65  (
67  (
68  "C2",
69  coeffDict_,
70  1.92
71  )
72  ),
73  sigmak_
74  (
76  (
77  "sigmak",
78  coeffDict_,
79  1.0
80  )
81  ),
82  sigmaEps_
83  (
85  (
86  "sigmaEps",
87  coeffDict_,
88  1.3
89  )
90  ),
91  A1_
92  (
94  (
95  "A1",
96  coeffDict_,
97  1.25
98  )
99  ),
100  A2_
101  (
103  (
104  "A2",
105  coeffDict_,
106  1000.0
107  )
108  ),
109  Ctau1_
110  (
112  (
113  "Ctau1",
114  coeffDict_,
115  -4.0
116  )
117  ),
118  Ctau2_
119  (
121  (
122  "Ctau2",
123  coeffDict_,
124  13.0
125  )
126  ),
127  Ctau3_
128  (
130  (
131  "Ctau3",
132  coeffDict_,
133  -2.0
134  )
135  ),
136  alphaKsi_
137  (
139  (
140  "alphaKsi",
141  coeffDict_,
142  0.9
143  )
144  ),
145 
146  kappa_
147  (
149  (
150  "kappa_",
151  coeffDict_,
152  0.41
153  )
154  ),
155  E_
156  (
158  (
159  "E",
160  coeffDict_,
161  9.8
162  )
163  ),
164 
165  k_
166  (
167  IOobject
168  (
169  "k",
170  runTime_.timeName(),
171  mesh_,
174  ),
175  mesh_
176  ),
177 
178  epsilon_
179  (
180  IOobject
181  (
182  "epsilon",
183  runTime_.timeName(),
184  mesh_,
187  ),
188  mesh_
189  ),
190 
191  gradU_(fvc::grad(U)),
192  eta_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))),
193  ksi_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))),
194  Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
195  fEta_(A2_ + pow(eta_, 3.0)),
196 
197  nut_("nut", Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_)),
198 
199  nonlinearStress_
200  (
201  "nonlinearStress",
202  symm
203  (
204  pow(k_, 3.0)/sqr(epsilon_)
205  *(
206  Ctau1_/fEta_
207  *(
208  (gradU_ & gradU_)
209  + (gradU_ & gradU_)().T()
210  )
211  + Ctau2_/fEta_*(gradU_ & gradU_.T())
212  + Ctau3_/fEta_*(gradU_.T() & gradU_)
213  )
214  )
215  )
216 {
218 
219  printCoeffs();
220 }
221 
222 
223 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
224 
226 {
227  volTensorField gradU_ = fvc::grad(U_);
228 
230  (
232  (
233  IOobject
234  (
235  "R",
236  runTime_.timeName(),
237  mesh_,
240  ),
241  ((2.0/3.0)*I)*k_ - nut_*twoSymm(gradU_) + nonlinearStress_,
242  k_.boundaryField().types()
243  )
244  );
245 }
246 
247 
249 {
251  (
253  (
254  IOobject
255  (
256  "devRhoReff",
257  runTime_.timeName(),
258  mesh_,
261  ),
262  -nuEff()*dev(twoSymm(fvc::grad(U_))) + nonlinearStress_
263  )
264  );
265 }
266 
267 
269 {
270  return
271  (
272  fvc::div(nonlinearStress_)
273  - fvm::laplacian(nuEff(), U)
274  - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
275  );
276 }
277 
278 
280 {
281  if (RASModel::read())
282  {
283  C1_.readIfPresent(coeffDict());
284  C2_.readIfPresent(coeffDict());
285  sigmak_.readIfPresent(coeffDict());
286  sigmaEps_.readIfPresent(coeffDict());
287  A1_.readIfPresent(coeffDict());
288  A2_.readIfPresent(coeffDict());
289  Ctau1_.readIfPresent(coeffDict());
290  Ctau2_.readIfPresent(coeffDict());
291  Ctau3_.readIfPresent(coeffDict());
292  alphaKsi_.readIfPresent(coeffDict());
293 
294  kappa_.readIfPresent(coeffDict());
295  E_.readIfPresent(coeffDict());
296 
297  return true;
298  }
299  else
300  {
301  return false;
302  }
303 }
304 
305 
307 {
309 
310  if (!turbulence_)
311  {
312  return;
313  }
314 
315  gradU_ = fvc::grad(U_);
316 
317  // generation term
318  volScalarField S2 = symm(gradU_) && gradU_;
319 
321  (
322  "RASModel::G",
323  Cmu_*sqr(k_)/epsilon_*S2
324  - (nonlinearStress_ && gradU_)
325  );
326 
328 
329  // Dissipation equation
330  tmp<fvScalarMatrix> epsEqn
331  (
332  fvm::ddt(epsilon_)
333  + fvm::div(phi_, epsilon_)
334  - fvm::laplacian(DepsilonEff(), epsilon_)
335  ==
336  C1_*G*epsilon_/k_
337  - fvm::Sp(C2_*epsilon_/k_, epsilon_)
338  );
339 
340  epsEqn().relax();
341 
343 
344  solve(epsEqn);
345  bound(epsilon_, epsilon0_);
346 
347 
348  // Turbulent kinetic energy equation
349 
351  (
352  fvm::ddt(k_)
353  + fvm::div(phi_, k_)
354  - fvm::laplacian(DkEff(), k_)
355  ==
356  G
357  - fvm::Sp(epsilon_/k_, k_)
358  );
359 
360  kEqn().relax();
361  solve(kEqn);
362  bound(k_, k0_);
363 
364 
365  // Re-calculate viscosity
366 
367  eta_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())));
368  ksi_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())));
369  Cmu_ = 2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_));
370  fEta_ = A2_ + pow(eta_, 3.0);
371 
372  nut_ = Cmu_*sqr(k_)/epsilon_;
373 
375 
376  nonlinearStress_ = symm
377  (
378  pow(k_, 3.0)/sqr(epsilon_)
379  *(
380  Ctau1_/fEta_
381  *(
382  (gradU_ & gradU_)
383  + (gradU_ & gradU_)().T()
384  )
385  + Ctau2_/fEta_*(gradU_ & gradU_.T())
386  + Ctau3_/fEta_*(gradU_.T() & gradU_)
387  )
388  );
389 }
390 
391 
392 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
393 
394 } // End namespace RASModels
395 } // End namespace incompressible
396 } // End namespace Foam
397 
398 // ************************ vim: set sw=4 sts=4 et: ************************ //