FreeFOAM The Cross-Platform CFD Toolkit
qZeta.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 "qZeta.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(qZeta, 0);
43 addToRunTimeSelectionTable(RASModel, qZeta, dictionary);
44 
45 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
46 
47 tmp<volScalarField> qZeta::fMu() const
48 {
49  volScalarField Rt = q_*k_/(2.0*nu()*zeta_);
50 
51  if (anisotropic_)
52  {
53  return exp((-scalar(2.5) + Rt/20.0)/pow(scalar(1) + Rt/130.0, 3.0));
54  }
55  else
56  {
57  return
58  exp(-6.0/sqr(scalar(1) + Rt/50.0))
59  *(scalar(1) + 3.0*exp(-Rt/10.0));
60  }
61 }
62 
63 
64 tmp<volScalarField> qZeta::f2() const
65 {
66  volScalarField Rt = q_*k_/(2.0*nu()*zeta_);
67  return scalar(1) - 0.3*exp(-sqr(Rt));
68 }
69 
70 
71 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
72 
74 (
75  const volVectorField& U,
76  const surfaceScalarField& phi,
77  transportModel& lamTransportModel
78 )
79 :
80  RASModel(typeName, U, phi, lamTransportModel),
81 
82  Cmu_
83  (
85  (
86  "Cmu",
87  coeffDict_,
88  0.09
89  )
90  ),
91  C1_
92  (
94  (
95  "C1",
96  coeffDict_,
97  1.44
98  )
99  ),
100  C2_
101  (
103  (
104  "C2",
105  coeffDict_,
106  1.92
107  )
108  ),
109  sigmaZeta_
110  (
112  (
113  "sigmaZeta",
114  coeffDict_,
115  1.3
116  )
117  ),
118  anisotropic_
119  (
121  (
122  "anisotropic",
123  coeffDict_,
124  false
125  )
126  ),
127 
128  k_
129  (
130  IOobject
131  (
132  "k",
133  runTime_.timeName(),
134  mesh_,
137  ),
138  mesh_
139  ),
140 
141  epsilon_
142  (
143  IOobject
144  (
145  "epsilon",
146  runTime_.timeName(),
147  mesh_,
150  ),
151  mesh_
152  ),
153 
154  q_
155  (
156  IOobject
157  (
158  "q",
159  runTime_.timeName(),
160  mesh_,
163  ),
164  sqrt(k_),
165  k_.boundaryField().types()
166  ),
167 
168  zeta_
169  (
170  IOobject
171  (
172  "zeta",
173  runTime_.timeName(),
174  mesh_,
177  ),
178  epsilon_/(2.0*q_),
179  epsilon_.boundaryField().types()
180  ),
181 
182  nut_
183  (
184  IOobject
185  (
186  "nut",
187  runTime_.timeName(),
188  mesh_,
191  ),
192  autoCreateNut("nut", mesh_)
193  )
194 {
195  nut_ = Cmu_*fMu()*sqr(k_)/(epsilon_ + epsilonSmall_);
196  nut_.correctBoundaryConditions();
197 
198  printCoeffs();
199 }
200 
201 
202 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
203 
205 {
207  (
209  (
210  IOobject
211  (
212  "R",
213  runTime_.timeName(),
214  mesh_,
217  ),
218  ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
219  k_.boundaryField().types()
220  )
221  );
222 }
223 
224 
226 {
228  (
230  (
231  IOobject
232  (
233  "devRhoReff",
234  runTime_.timeName(),
235  mesh_,
238  ),
240  )
241  );
242 }
243 
244 
246 {
247  return
248  (
249  - fvm::laplacian(nuEff(), U)
250  - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
251  );
252 }
253 
254 
256 {
257  if (RASModel::read())
258  {
259  Cmu_.readIfPresent(coeffDict());
260  C1_.readIfPresent(coeffDict());
261  C2_.readIfPresent(coeffDict());
262  sigmaZeta_.readIfPresent(coeffDict());
263  anisotropic_.readIfPresent("anisotropic", coeffDict());
264 
265  return true;
266  }
267  else
268  {
269  return false;
270  }
271 }
272 
273 
275 {
277 
278  if (!turbulence_)
279  {
280  return;
281  }
282 
284 
285  volScalarField G("RASModel::G", nut_/(2.0*q_)*S2);
287 
288 
289  // Zeta equation
290 
291  tmp<fvScalarMatrix> zetaEqn
292  (
293  fvm::ddt(zeta_)
294  + fvm::div(phi_, zeta_)
295  - fvm::laplacian(DzetaEff(), zeta_)
296  ==
297  (2.0*C1_ - 1)*G*zeta_/q_
298  - fvm::Sp((2.0*C2_ - dimensionedScalar(1.0))*f2()*zeta_/q_, zeta_)
299  + E
300  );
301 
302  zetaEqn().relax();
303  solve(zetaEqn);
304  bound(zeta_, epsilon0_/(2*sqrt(k0_)));
305 
306 
307  // q equation
308 
310  (
311  fvm::ddt(q_)
312  + fvm::div(phi_, q_)
313  - fvm::laplacian(DqEff(), q_)
314  ==
315  G - fvm::Sp(zeta_/q_, q_)
316  );
317 
318  qEqn().relax();
319  solve(qEqn);
320  bound(q_, sqrt(k0_));
321 
322 
323  // Re-calculate k and epsilon
324  k_ = sqr(q_);
326 
327  epsilon_ = 2*q_*zeta_;
328  epsilon_.correctBoundaryConditions();
329 
330 
331  // Re-calculate viscosity
332  nut_ = Cmu_*fMu()*sqr(k_)/epsilon_;
334 }
335 
336 
337 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
338 
339 } // End namespace RASModels
340 } // End namespace incompressible
341 } // End namespace Foam
342 
343 // ************************ vim: set sw=4 sts=4 et: ************************ //