FreeFOAM The Cross-Platform CFD Toolkit
LienCubicKE.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 "LienCubicKE.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(LienCubicKE, 0);
43 addToRunTimeSelectionTable(RASModel, LienCubicKE, dictionary);
44 
45 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
46 
48 (
49  const volVectorField& U,
50  const surfaceScalarField& phi,
51  transportModel& lamTransportModel
52 )
53 :
54  RASModel(typeName, U, phi, lamTransportModel),
55 
56  C1_
57  (
59  (
60  "C1",
61  coeffDict_,
62  1.44
63  )
64  ),
65  C2_
66  (
68  (
69  "C2",
70  coeffDict_,
71  1.92
72  )
73  ),
74  sigmak_
75  (
77  (
78  "sigmak",
79  coeffDict_,
80  1.0
81  )
82  ),
83  sigmaEps_
84  (
86  (
87  "sigmaEps",
88  coeffDict_,
89  1.3
90  )
91  ),
92  A1_
93  (
95  (
96  "A1",
97  coeffDict_,
98  1.25
99  )
100  ),
101  A2_
102  (
104  (
105  "A2",
106  coeffDict_,
107  1000.0
108  )
109  ),
110  Ctau1_
111  (
113  (
114  "Ctau1",
115  coeffDict_,
116  -4.0
117  )
118  ),
119  Ctau2_
120  (
122  (
123  "Ctau2",
124  coeffDict_,
125  13.0
126  )
127  ),
128  Ctau3_
129  (
131  (
132  "Ctau3",
133  coeffDict_,
134  -2.0
135  )
136  ),
137  alphaKsi_
138  (
140  (
141  "alphaKsi",
142  coeffDict_,
143  0.9
144  )
145  ),
146 
147  k_
148  (
149  IOobject
150  (
151  "k",
152  runTime_.timeName(),
153  mesh_,
156  ),
157  autoCreateK("k", mesh_)
158  ),
159  epsilon_
160  (
161  IOobject
162  (
163  "epsilon",
164  runTime_.timeName(),
165  mesh_,
168  ),
169  autoCreateEpsilon("epsilon", mesh_)
170  ),
171 
172  gradU_(fvc::grad(U)),
173  eta_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))),
174  ksi_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))),
175  Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
176  fEta_(A2_ + pow(eta_, 3.0)),
177 
178  C5viscosity_
179  (
180  - 2.0*pow3(Cmu_)*pow4(k_)/pow3(epsilon_)
181  *(
182  magSqr(gradU_ + gradU_.T())
183  - magSqr(gradU_ - gradU_.T())
184  )
185  ),
186 
187  nut_
188  (
189  IOobject
190  (
191  "nut",
192  runTime_.timeName(),
193  mesh_,
196  ),
197  autoCreateNut("nut", mesh_)
198  ),
199 
200  nonlinearStress_
201  (
202  "nonlinearStress",
203  // quadratic terms
204  symm
205  (
206  pow(k_, 3.0)/sqr(epsilon_)
207  *(
208  Ctau1_/fEta_
209  *(
210  (gradU_ & gradU_)
211  + (gradU_ & gradU_)().T()
212  )
213  + Ctau2_/fEta_*(gradU_ & gradU_.T())
214  + Ctau3_/fEta_*(gradU_.T() & gradU_)
215  )
216  // cubic term C4
217  - 20.0*pow(k_, 4.0)/pow(epsilon_, 3.0)
218  *pow(Cmu_, 3.0)
219  *(
220  ((gradU_ & gradU_) & gradU_.T())
221  + ((gradU_ & gradU_.T()) & gradU_.T())
222  - ((gradU_.T() & gradU_) & gradU_)
223  - ((gradU_.T() & gradU_.T()) & gradU_)
224  )
225  )
226  )
227 {
228  nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_) + C5viscosity_;
229  nut_.correctBoundaryConditions();
230 
231  printCoeffs();
232 }
233 
234 
235 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
236 
238 {
240  (
242  (
243  IOobject
244  (
245  "R",
246  runTime_.timeName(),
247  mesh_,
250  ),
251  ((2.0/3.0)*I)*k_ - nut_*twoSymm(gradU_) + nonlinearStress_,
252  k_.boundaryField().types()
253  )
254  );
255 }
256 
257 
259 {
261  (
263  (
264  IOobject
265  (
266  "devRhoReff",
267  runTime_.timeName(),
268  mesh_,
271  ),
272  -nuEff()*dev(twoSymm(fvc::grad(U_))) + nonlinearStress_
273  )
274  );
275 }
276 
277 
279 {
280  return
281  (
282  fvc::div(nonlinearStress_)
283  - fvm::laplacian(nuEff(), U)
284  - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
285  );
286 }
287 
288 
290 {
291  if (RASModel::read())
292  {
293  C1_.readIfPresent(coeffDict());
294  C2_.readIfPresent(coeffDict());
295  sigmak_.readIfPresent(coeffDict());
296  sigmaEps_.readIfPresent(coeffDict());
297  A1_.readIfPresent(coeffDict());
298  A2_.readIfPresent(coeffDict());
299  Ctau1_.readIfPresent(coeffDict());
300  Ctau2_.readIfPresent(coeffDict());
301  Ctau3_.readIfPresent(coeffDict());
302  alphaKsi_.readIfPresent(coeffDict());
303 
304  return true;
305  }
306  else
307  {
308  return false;
309  }
310 }
311 
312 
314 {
316 
317  if (!turbulence_)
318  {
319  return;
320  }
321 
322  gradU_ = fvc::grad(U_);
323 
324  // generation term
325  volScalarField S2 = symm(gradU_) && gradU_;
326 
328  (
329  "RASModel::G",
330  Cmu_*sqr(k_)/epsilon_*S2 - (nonlinearStress_ && gradU_)
331  );
332 
333  // Update espsilon and G at the wall
334  epsilon_.boundaryField().updateCoeffs();
335 
336  // Dissipation equation
337  tmp<fvScalarMatrix> epsEqn
338  (
339  fvm::ddt(epsilon_)
340  + fvm::div(phi_, epsilon_)
341  - fvm::laplacian(DepsilonEff(), epsilon_)
342  ==
343  C1_*G*epsilon_/k_
344  - fvm::Sp(C2_*epsilon_/k_, epsilon_)
345  );
346 
347  epsEqn().relax();
348 
349  epsEqn().boundaryManipulate(epsilon_.boundaryField());
350 
351  solve(epsEqn);
352  bound(epsilon_, epsilon0_);
353 
354 
355  // Turbulent kinetic energy equation
356 
358  (
359  fvm::ddt(k_)
360  + fvm::div(phi_, k_)
361  - fvm::laplacian(DkEff(), k_)
362  ==
363  G
364  - fvm::Sp(epsilon_/k_, k_)
365  );
366 
367  kEqn().relax();
368  solve(kEqn);
369  bound(k_, k0_);
370 
371 
372  // Re-calculate viscosity
373 
374  eta_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())));
375  ksi_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())));
376  Cmu_ = 2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_));
377  fEta_ = A2_ + pow(eta_, 3.0);
378 
379  C5viscosity_ =
380  - 2.0*pow(Cmu_, 3.0)*pow(k_, 4.0)/pow(epsilon_, 3.0)
381  *(magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T()));
382 
383  nut_ = Cmu_*sqr(k_)/epsilon_ + C5viscosity_;
385 
386  nonlinearStress_ = symm
387  (
388  // quadratic terms
389  pow(k_, 3.0)/sqr(epsilon_)*
390  (
391  Ctau1_/fEta_*
392  (
393  (gradU_ & gradU_)
394  + (gradU_ & gradU_)().T()
395  )
396  + Ctau2_/fEta_*(gradU_ & gradU_.T())
397  + Ctau3_/fEta_*(gradU_.T() & gradU_)
398  )
399  // cubic term C4
400  - 20.0*pow(k_, 4.0)/pow(epsilon_, 3.0)
401  *pow(Cmu_, 3.0)
402  *(
403  ((gradU_ & gradU_) & gradU_.T())
404  + ((gradU_ & gradU_.T()) & gradU_.T())
405  - ((gradU_.T() & gradU_) & gradU_)
406  - ((gradU_.T() & gradU_.T()) & gradU_)
407  )
408  );
409 }
410 
411 
412 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
413 
414 } // End namespace RASModels
415 } // End namespace incompressible
416 } // End namespace Foam
417 
418 // ************************ vim: set sw=4 sts=4 et: ************************ //