FreeFOAM The Cross-Platform CFD Toolkit
LienCubicKELowRe.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 "LienCubicKELowRe.H"
29 
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 namespace incompressible
37 {
38 namespace RASModels
39 {
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
43 defineTypeNameAndDebug(LienCubicKELowRe, 0);
44 addToRunTimeSelectionTable(RASModel, LienCubicKELowRe, dictionary);
45 
46 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
47 
49 (
50  const volVectorField& U,
51  const surfaceScalarField& phi,
52  transportModel& lamTransportModel
53 )
54 :
55  RASModel(typeName, U, phi, lamTransportModel),
56 
57  C1_
58  (
60  (
61  "C1",
62  coeffDict_,
63  1.44
64  )
65  ),
66  C2_
67  (
69  (
70  "C2",
71  coeffDict_,
72  1.92
73  )
74  ),
75  sigmak_
76  (
78  (
79  "sigmak",
80  coeffDict_,
81  1.0
82  )
83  ),
84  sigmaEps_
85  (
87  (
88  "sigmaEps",
89  coeffDict_,
90  1.3
91  )
92  ),
93  A1_
94  (
96  (
97  "A1",
98  coeffDict_,
99  1.25
100  )
101  ),
102  A2_
103  (
105  (
106  "A2",
107  coeffDict_,
108  1000.0
109  )
110  ),
111  Ctau1_
112  (
114  (
115  "Ctau1",
116  coeffDict_,
117  -4.0
118  )
119  ),
120  Ctau2_
121  (
123  (
124  "Ctau2",
125  coeffDict_,
126  13.0
127  )
128  ),
129  Ctau3_
130  (
132  (
133  "Ctau3",
134  coeffDict_,
135  -2.0
136  )
137  ),
138  alphaKsi_
139  (
141  (
142  "alphaKsi",
143  coeffDict_,
144  0.9
145  )
146  ),
147  CmuWall_
148  (
150  (
151  "Cmu",
152  coeffDict_,
153  0.09
154  )
155  ),
156  kappa_
157  (
159  (
160  "kappa",
161  coeffDict_,
162  0.41
163  )
164  ),
165  Am_
166  (
168  (
169  "Am",
170  coeffDict_,
171  0.016
172  )
173  ),
174  Aepsilon_
175  (
177  (
178  "Aepsilon",
179  coeffDict_,
180  0.263
181  )
182  ),
183  Amu_
184  (
186  (
187  "Amu",
188  coeffDict_,
189  0.00222
190  )
191  ),
192 
193  k_
194  (
195  IOobject
196  (
197  "k",
198  runTime_.timeName(),
199  mesh_,
202  ),
203  mesh_
204  ),
205 
206  epsilon_
207  (
208  IOobject
209  (
210  "epsilon",
211  runTime_.timeName(),
212  mesh_,
215  ),
216  mesh_
217  ),
218 
219  y_(mesh_),
220 
221  gradU_(fvc::grad(U)),
222  eta_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))),
223  ksi_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))),
224  Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
225  fEta_(A2_ + pow(eta_, 3.0)),
226 
227  C5viscosity_
228  (
229  -2.0*pow(Cmu_, 3.0)*pow(k_, 4.0)/pow(epsilon_, 3.0)
230  *(magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T()))
231  ),
232 
233  yStar_(sqrt(k_)*y_/nu() + SMALL),
234 
235  nut_
236  (
237  IOobject
238  (
239  "nut",
240  runTime_.timeName(),
241  mesh_,
244  ),
245  autoCreateLowReNut("nut", mesh_)
246  ),
247 
248  nonlinearStress_
249  (
250  "nonlinearStress",
251  symm
252  (
253  // quadratic terms
254  pow(k_, 3.0)/sqr(epsilon_)
255  *(
256  Ctau1_/fEta_
257  *(
258  (gradU_ & gradU_)
259  + (gradU_ & gradU_)().T()
260  )
261  + Ctau2_/fEta_*(gradU_ & gradU_.T())
262  + Ctau3_/fEta_*(gradU_.T() & gradU_)
263  )
264  // cubic term C4
265  - 20.0*pow(k_, 4.0)/pow(epsilon_, 3.0)
266  *pow(Cmu_, 3.0)
267  *(
268  ((gradU_ & gradU_) & gradU_.T())
269  + ((gradU_ & gradU_.T()) & gradU_.T())
270  - ((gradU_.T() & gradU_) & gradU_)
271  - ((gradU_.T() & gradU_.T()) & gradU_)
272  )
273  // cubic term C5, explicit part
274  + min
275  (
276  C5viscosity_,
277  dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
278  )*gradU_
279  )
280  )
281 {
282  nut_ = Cmu_
283  *(
284  scalar(1) - exp(-Am_*yStar_))
285  /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL
286  )
287  *sqr(k_)/(epsilon_ + epsilonSmall_)
288  // cubic term C5, implicit part
289  + max
290  (
291  C5viscosity_,
292  dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
293  );
294 
295  nut_.correctBoundaryConditions();
296 
297  printCoeffs();
298 }
299 
300 
301 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
302 
304 {
306  (
308  (
309  IOobject
310  (
311  "R",
312  runTime_.timeName(),
313  mesh_,
316  ),
317  ((2.0/3.0)*I)*k_ - nut_*twoSymm(gradU_) + nonlinearStress_,
318  k_.boundaryField().types()
319  )
320  );
321 }
322 
323 
325 {
327  (
329  (
330  IOobject
331  (
332  "devRhoReff",
333  runTime_.timeName(),
334  mesh_,
337  ),
338  -nuEff()*dev(twoSymm(fvc::grad(U_))) + nonlinearStress_
339  )
340  );
341 }
342 
343 
345 {
346  return
347  (
348  fvc::div(nonlinearStress_)
349  - fvm::laplacian(nuEff(), U)
350  - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
351  );
352 }
353 
354 
356 {
357  if (RASModel::read())
358  {
359  C1_.readIfPresent(coeffDict());
360  C2_.readIfPresent(coeffDict());
361  sigmak_.readIfPresent(coeffDict());
362  sigmaEps_.readIfPresent(coeffDict());
363  A1_.readIfPresent(coeffDict());
364  A2_.readIfPresent(coeffDict());
365  Ctau1_.readIfPresent(coeffDict());
366  Ctau2_.readIfPresent(coeffDict());
367  Ctau3_.readIfPresent(coeffDict());
368  alphaKsi_.readIfPresent(coeffDict());
369  CmuWall_.readIfPresent(coeffDict());
370  kappa_.readIfPresent(coeffDict());
371  Am_.readIfPresent(coeffDict());
372  Aepsilon_.readIfPresent(coeffDict());
373  Amu_.readIfPresent(coeffDict());
374 
375  return true;
376  }
377  else
378  {
379  return false;
380  }
381 }
382 
383 
385 {
387 
388  if (!turbulence_)
389  {
390  return;
391  }
392 
393  if (mesh_.changing())
394  {
395  y_.correct();
396  }
397 
398  gradU_ = fvc::grad(U_);
399 
400  // generation term
401  volScalarField S2 = symm(gradU_) && gradU_;
402 
403  yStar_ = sqrt(k_)*y_/nu() + SMALL;
404  volScalarField Rt = sqr(k_)/(nu()*epsilon_);
405 
406  volScalarField fMu =
407  (scalar(1) - exp(-Am_*yStar_))
408  /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL);
409 
410  volScalarField f2 = scalar(1) - 0.3*exp(-sqr(Rt));
411 
412  volScalarField G =
413  Cmu_*fMu*sqr(k_)/epsilon_*S2 - (nonlinearStress_ && gradU_);
414 
415  // Dissipation equation
416  tmp<fvScalarMatrix> epsEqn
417  (
418  fvm::ddt(epsilon_)
419  + fvm::div(phi_, epsilon_)
420  - fvm::laplacian(DepsilonEff(), epsilon_)
421  ==
422  C1_*G*epsilon_/k_
423  // E-term
424  + C2_*f2*pow(Cmu_, 0.75)*pow(k_, scalar(0.5))
425  /(kappa_*y_*(scalar(1) - exp(-Aepsilon_*yStar_)))
426  *exp(-Amu_*sqr(yStar_))*epsilon_
427  - fvm::Sp(C2_*f2*epsilon_/k_, epsilon_)
428  );
429 
430  epsEqn().relax();
431 
434 
435  solve(epsEqn);
436  bound(epsilon_, epsilon0_);
437 
438 
439  // Turbulent kinetic energy equation
440 
442  (
443  fvm::ddt(k_)
444  + fvm::div(phi_, k_)
445  - fvm::laplacian(DkEff(), k_)
446  ==
447  G
448  - fvm::Sp(epsilon_/k_, k_)
449  );
450 
451  kEqn().relax();
452  solve(kEqn);
453  bound(k_, k0_);
454 
455 
456  // Re-calculate viscosity
457 
458  eta_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())));
459  ksi_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())));
460  Cmu_ = 2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_));
461  fEta_ = A2_ + pow(eta_, 3.0);
462 
463  C5viscosity_ =
464  - 2.0*pow(Cmu_, 3.0)*pow(k_, 4.0)/pow(epsilon_, 3.0)
465  *(magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T()));
466 
467  nut_ =
468  Cmu_*fMu*sqr(k_)/epsilon_
469  // C5 term, implicit
470  + max
471  (
472  C5viscosity_,
473  dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
474  );
475 
476  nonlinearStress_ = symm
477  (
478  // quadratic terms
479  pow(k_, 3.0)/sqr(epsilon_)
480  *(
481  Ctau1_/fEta_
482  *(
483  (gradU_ & gradU_)
484  + (gradU_ & gradU_)().T()
485  )
486  + Ctau2_/fEta_*(gradU_ & gradU_.T())
487  + Ctau3_/fEta_*(gradU_.T() & gradU_)
488  )
489  // cubic term C4
490  - 20.0*pow(k_, 4.0)/pow(epsilon_, 3.0)
491  *pow(Cmu_, 3.0)
492  *(
493  ((gradU_ & gradU_) & gradU_.T())
494  + ((gradU_ & gradU_.T()) & gradU_.T())
495  - ((gradU_.T() & gradU_) & gradU_)
496  - ((gradU_.T() & gradU_.T()) & gradU_)
497  )
498  // cubic term C5, explicit part
499  + min
500  (
501  C5viscosity_,
502  dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
503  )*gradU_
504  );
505 }
506 
507 
508 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
509 
510 } // End namespace RASModels
511 } // End namespace incompressible
512 } // End namespace Foam
513 
514 // ************************ vim: set sw=4 sts=4 et: ************************ //