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 incompressible
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)*nu()/(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)*nu()/(sqr(y_)*omega_)
80  ),
81  scalar(100)
82  );
83 
84  return tanh(sqr(arg2));
85 }
86 
87 
88 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
89 
91 (
92  const volVectorField& U,
93  const surfaceScalarField& phi,
94  transportModel& lamTransportModel
95 )
96 :
97  RASModel(typeName, U, phi, lamTransportModel),
98 
99  alphaK1_
100  (
102  (
103  "alphaK1",
104  coeffDict_,
105  0.85034
106  )
107  ),
108  alphaK2_
109  (
111  (
112  "alphaK2",
113  coeffDict_,
114  1.0
115  )
116  ),
117  alphaOmega1_
118  (
120  (
121  "alphaOmega1",
122  coeffDict_,
123  0.5
124  )
125  ),
126  alphaOmega2_
127  (
129  (
130  "alphaOmega2",
131  coeffDict_,
132  0.85616
133  )
134  ),
135  gamma1_
136  (
138  (
139  "gamma1",
140  coeffDict_,
141  0.5532
142  )
143  ),
144  gamma2_
145  (
147  (
148  "gamma2",
149  coeffDict_,
150  0.4403
151  )
152  ),
153  beta1_
154  (
156  (
157  "beta1",
158  coeffDict_,
159  0.075
160  )
161  ),
162  beta2_
163  (
165  (
166  "beta2",
167  coeffDict_,
168  0.0828
169  )
170  ),
171  betaStar_
172  (
174  (
175  "betaStar",
176  coeffDict_,
177  0.09
178  )
179  ),
180  a1_
181  (
183  (
184  "a1",
185  coeffDict_,
186  0.31
187  )
188  ),
189  c1_
190  (
192  (
193  "c1",
194  coeffDict_,
195  10.0
196  )
197  ),
198 
199  y_(mesh_),
200 
201  k_
202  (
203  IOobject
204  (
205  "k",
206  runTime_.timeName(),
207  mesh_,
210  ),
211  autoCreateK("k", mesh_)
212  ),
213  omega_
214  (
215  IOobject
216  (
217  "omega",
218  runTime_.timeName(),
219  mesh_,
222  ),
223  autoCreateOmega("omega", mesh_)
224  ),
225  nut_
226  (
227  IOobject
228  (
229  "nut",
230  runTime_.timeName(),
231  mesh_,
234  ),
235  autoCreateNut("nut", mesh_)
236  )
237 {
238  bound(omega_, omega0_);
239 
240  nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(2.0)*mag(symm(fvc::grad(U_))));
241  nut_.correctBoundaryConditions();
242 
243  printCoeffs();
244 }
245 
246 
247 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
248 
250 {
252  (
254  (
255  IOobject
256  (
257  "R",
258  runTime_.timeName(),
259  mesh_,
262  ),
263  ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
264  k_.boundaryField().types()
265  )
266  );
267 }
268 
269 
271 {
273  (
275  (
276  IOobject
277  (
278  "devRhoReff",
279  runTime_.timeName(),
280  mesh_,
283  ),
285  )
286  );
287 }
288 
289 
291 {
292  return
293  (
294  - fvm::laplacian(nuEff(), U)
295  - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
296  );
297 }
298 
299 
301 {
302  if (RASModel::read())
303  {
304  alphaK1_.readIfPresent(coeffDict());
305  alphaK2_.readIfPresent(coeffDict());
306  alphaOmega1_.readIfPresent(coeffDict());
307  alphaOmega2_.readIfPresent(coeffDict());
308  gamma1_.readIfPresent(coeffDict());
309  gamma2_.readIfPresent(coeffDict());
310  beta1_.readIfPresent(coeffDict());
311  beta2_.readIfPresent(coeffDict());
312  betaStar_.readIfPresent(coeffDict());
313  a1_.readIfPresent(coeffDict());
314  c1_.readIfPresent(coeffDict());
315 
316  return true;
317  }
318  else
319  {
320  return false;
321  }
322 }
323 
324 
326 {
328 
329  if (!turbulence_)
330  {
331  return;
332  }
333 
334  if (mesh_.changing())
335  {
336  y_.correct();
337  }
338 
340  volScalarField G("RASModel::G", nut_*2*S2);
341 
342  // Update omega and G at the wall
343  omega_.boundaryField().updateCoeffs();
344 
345  volScalarField CDkOmega =
346  (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
347 
348  volScalarField F1 = this->F1(CDkOmega);
349 
350  // Turbulent frequency equation
351  tmp<fvScalarMatrix> omegaEqn
352  (
353  fvm::ddt(omega_)
354  + fvm::div(phi_, omega_)
355  - fvm::Sp(fvc::div(phi_), omega_)
356  - fvm::laplacian(DomegaEff(F1), omega_)
357  ==
358  gamma(F1)*2*S2
359  - fvm::Sp(beta(F1)*omega_, omega_)
360  - fvm::SuSp
361  (
362  (F1 - scalar(1))*CDkOmega/omega_,
363  omega_
364  )
365  );
366 
367  omegaEqn().relax();
368 
369  omegaEqn().boundaryManipulate(omega_.boundaryField());
370 
371  solve(omegaEqn);
372  bound(omega_, omega0_);
373 
374  // Turbulent kinetic energy equation
376  (
377  fvm::ddt(k_)
378  + fvm::div(phi_, k_)
379  - fvm::Sp(fvc::div(phi_), k_)
380  - fvm::laplacian(DkEff(F1), k_)
381  ==
382  min(G, c1_*betaStar_*k_*omega_)
383  - fvm::Sp(betaStar_*omega_, k_)
384  );
385 
386  kEqn().relax();
387  solve(kEqn);
388  bound(k_, k0_);
389 
390 
391  // Re-calculate viscosity
392  nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(2*S2));
394 }
395 
396 
397 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
398 
399 } // End namespace RASModels
400 } // End namespace incompressible
401 } // End namespace Foam
402 
403 // ************************ vim: set sw=4 sts=4 et: ************************ //