FreeFOAM The Cross-Platform CFD Toolkit
kOmegaSST.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 "kOmegaSST.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(kOmegaSST, 0);
43 addToRunTimeSelectionTable(RASModel, kOmegaSST, dictionary);
44 
45 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
46 
47 tmp<volScalarField> kOmegaSST::F1(const volScalarField& CDkOmega) const
48 {
49  volScalarField CDkOmegaPlus = max
50  (
51  CDkOmega,
52  dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
53  );
54 
55  volScalarField arg1 = min
56  (
57  min
58  (
59  max
60  (
61  (scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
62  scalar(500)*(mu()/rho_)/(sqr(y_)*omega_)
63  ),
64  (4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
65  ),
66  scalar(10)
67  );
68 
69  return tanh(pow4(arg1));
70 }
71 
72 tmp<volScalarField> kOmegaSST::F2() const
73 {
74  volScalarField arg2 = min
75  (
76  max
77  (
78  (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
79  scalar(500)*(mu()/rho_)/(sqr(y_)*omega_)
80  ),
81  scalar(100)
82  );
83 
84  return tanh(sqr(arg2));
85 }
86 
87 
88 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
89 
90 kOmegaSST::kOmegaSST
91 (
92  const volScalarField& rho,
93  const volVectorField& U,
94  const surfaceScalarField& phi,
95  const basicThermo& thermophysicalModel
96 )
97 :
98  RASModel(typeName, rho, U, phi, thermophysicalModel),
99 
100  alphaK1_
101  (
103  (
104  "alphaK1",
105  coeffDict_,
106  0.85034
107  )
108  ),
109  alphaK2_
110  (
112  (
113  "alphaK2",
114  coeffDict_,
115  1.0
116  )
117  ),
118  alphaOmega1_
119  (
121  (
122  "alphaOmega1",
123  coeffDict_,
124  0.5
125  )
126  ),
127  alphaOmega2_
128  (
130  (
131  "alphaOmega2",
132  coeffDict_,
133  0.85616
134  )
135  ),
136  Prt_
137  (
139  (
140  "Prt",
141  coeffDict_,
142  1.0
143  )
144  ),
145  gamma1_
146  (
148  (
149  "gamma1",
150  coeffDict_,
151  0.5532
152  )
153  ),
154  gamma2_
155  (
157  (
158  "gamma2",
159  coeffDict_,
160  0.4403
161  )
162  ),
163  beta1_
164  (
166  (
167  "beta1",
168  coeffDict_,
169  0.075
170  )
171  ),
172  beta2_
173  (
175  (
176  "beta2",
177  coeffDict_,
178  0.0828
179  )
180  ),
181  betaStar_
182  (
184  (
185  "betaStar",
186  coeffDict_,
187  0.09
188  )
189  ),
190  a1_
191  (
193  (
194  "a1",
195  coeffDict_,
196  0.31
197  )
198  ),
199  c1_
200  (
202  (
203  "c1",
204  coeffDict_,
205  10.0
206  )
207  ),
208 
209  y_(mesh_),
210 
211  k_
212  (
213  IOobject
214  (
215  "k",
216  runTime_.timeName(),
217  mesh_,
220  ),
221  autoCreateK("k", mesh_)
222  ),
223  omega_
224  (
225  IOobject
226  (
227  "omega",
228  runTime_.timeName(),
229  mesh_,
232  ),
233  autoCreateOmega("omega", mesh_)
234  ),
235  mut_
236  (
237  IOobject
238  (
239  "mut",
240  runTime_.timeName(),
241  mesh_,
244  ),
245  autoCreateMut("mut", mesh_)
246  ),
247  alphat_
248  (
249  IOobject
250  (
251  "alphat",
252  runTime_.timeName(),
253  mesh_,
256  ),
257  autoCreateAlphat("alphat", mesh_)
258  )
259 {
260  bound(omega_, omega0_);
261 
262  mut_ =
263  (
264  a1_*rho_*k_
265  / max
266  (
267  a1_*omega_,
268  F2()*sqrt(2.0*magSqr(symm(fvc::grad(U_))))
269  )
270  );
271  mut_.correctBoundaryConditions();
272 
273  alphat_ = mut_/Prt_;
274  alphat_.correctBoundaryConditions();
275 
276  printCoeffs();
277 }
278 
279 
280 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
281 
282 tmp<volSymmTensorField> kOmegaSST::R() const
283 {
285  (
287  (
288  IOobject
289  (
290  "R",
291  runTime_.timeName(),
292  mesh_,
295  ),
296  ((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
297  k_.boundaryField().types()
298  )
299  );
300 }
301 
302 
303 tmp<volSymmTensorField> kOmegaSST::devRhoReff() const
304 {
306  (
308  (
309  IOobject
310  (
311  "devRhoReff",
312  runTime_.timeName(),
313  mesh_,
316  ),
317  -muEff()*dev(twoSymm(fvc::grad(U_)))
318  )
319  );
320 }
321 
322 
323 tmp<fvVectorMatrix> kOmegaSST::divDevRhoReff(volVectorField& U) const
324 {
325  return
326  (
327  - fvm::laplacian(muEff(), U) - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
328  );
329 }
330 
331 
332 bool kOmegaSST::read()
333 {
334  if (RASModel::read())
335  {
336  alphaK1_.readIfPresent(coeffDict());
337  alphaK2_.readIfPresent(coeffDict());
338  alphaOmega1_.readIfPresent(coeffDict());
339  alphaOmega2_.readIfPresent(coeffDict());
340  Prt_.readIfPresent(coeffDict());
341  gamma1_.readIfPresent(coeffDict());
342  gamma2_.readIfPresent(coeffDict());
343  beta1_.readIfPresent(coeffDict());
344  beta2_.readIfPresent(coeffDict());
345  betaStar_.readIfPresent(coeffDict());
346  a1_.readIfPresent(coeffDict());
347  c1_.readIfPresent(coeffDict());
348 
349  return true;
350  }
351  else
352  {
353  return false;
354  }
355 }
356 
357 
359 {
360  if (!turbulence_)
361  {
362  // Re-calculate viscosity
363  mut_ =
364  a1_*rho_*k_
365  /max(a1_*omega_, F2()*sqrt(2.0*magSqr(symm(fvc::grad(U_)))));
366  mut_.correctBoundaryConditions();
367 
368  // Re-calculate thermal diffusivity
369  alphat_ = mut_/Prt_;
370  alphat_.correctBoundaryConditions();
371 
372  return;
373  }
374 
376 
377  volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));
378 
379  if (mesh_.changing())
380  {
381  y_.correct();
382  }
383 
384  if (mesh_.moving())
385  {
386  divU += fvc::div(mesh_.phi());
387  }
388 
389  tmp<volTensorField> tgradU = fvc::grad(U_);
390  volScalarField S2 = magSqr(symm(tgradU()));
391  volScalarField GbyMu = (tgradU() && dev(twoSymm(tgradU())));
392  volScalarField G("RASModel::G", mut_*GbyMu);
393  tgradU.clear();
394 
395  // Update omega and G at the wall
396  omega_.boundaryField().updateCoeffs();
397 
398  volScalarField CDkOmega =
399  (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
400 
401  volScalarField F1 = this->F1(CDkOmega);
402  volScalarField rhoGammaF1 = rho_*gamma(F1);
403 
404  // Turbulent frequency equation
405  tmp<fvScalarMatrix> omegaEqn
406  (
407  fvm::ddt(rho_, omega_)
408  + fvm::div(phi_, omega_)
409  - fvm::laplacian(DomegaEff(F1), omega_)
410  ==
411  rhoGammaF1*GbyMu
412  - fvm::SuSp((2.0/3.0)*rhoGammaF1*divU, omega_)
413  - fvm::Sp(rho_*beta(F1)*omega_, omega_)
414  - fvm::SuSp
415  (
416  rho_*(F1 - scalar(1))*CDkOmega/omega_,
417  omega_
418  )
419  );
420 
421  omegaEqn().relax();
422 
423  omegaEqn().boundaryManipulate(omega_.boundaryField());
424 
425  solve(omegaEqn);
426  bound(omega_, omega0_);
427 
428  // Turbulent kinetic energy equation
430  (
431  fvm::ddt(rho_, k_)
432  + fvm::div(phi_, k_)
433  - fvm::laplacian(DkEff(F1), k_)
434  ==
435  min(G, (c1_*betaStar_)*rho_*k_*omega_)
436  - fvm::SuSp(2.0/3.0*rho_*divU, k_)
437  - fvm::Sp(rho_*betaStar_*omega_, k_)
438  );
439 
440  kEqn().relax();
441  solve(kEqn);
442  bound(k_, k0_);
443 
444 
445  // Re-calculate viscosity
446  mut_ = a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(2.0*S2));
447  mut_.correctBoundaryConditions();
448 
449  // Re-calculate thermal diffusivity
450  alphat_ = mut_/Prt_;
451  alphat_.correctBoundaryConditions();
452 }
453 
454 
455 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
456 
457 } // End namespace RASModels
458 } // End namespace compressible
459 } // End namespace Foam
460 
461 // ************************ vim: set sw=4 sts=4 et: ************************ //