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) 2011-2015 OpenFOAM Foundation
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"
27 #include "bound.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 namespace incompressible
35 {
36 namespace RASModels
37 {
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
41 defineTypeNameAndDebug(qZeta, 0);
42 addToRunTimeSelectionTable(RASModel, qZeta, dictionary);
43 
44 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
45 
47 {
48  const volScalarField Rt(q_*k_/(2.0*nu()*zeta_));
49 
50  if (anisotropic_)
51  {
52  return exp((-scalar(2.5) + Rt/20.0)/pow3(scalar(1) + Rt/130.0));
53  }
54  else
55  {
56  return
57  exp(-6.0/sqr(scalar(1) + Rt/50.0))
58  *(scalar(1) + 3.0*exp(-Rt/10.0));
59  }
60 }
61 
62 
64 {
65  tmp<volScalarField> Rt = q_*k_/(2.0*nu()*zeta_);
66  return scalar(1) - 0.3*exp(-sqr(Rt));
67 }
68 
69 
71 {
72  nut_ = Cmu_*fMu()*sqr(k_)/epsilon_;
74 }
75 
76 
77 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
78 
80 (
81  const geometricOneField& alpha,
82  const geometricOneField& rho,
83  const volVectorField& U,
84  const surfaceScalarField& alphaRhoPhi,
85  const surfaceScalarField& phi,
86  const transportModel& transport,
87  const word& propertiesName,
88  const word& type
89 )
90 :
92  (
93  type,
94  alpha,
95  rho,
96  U,
97  alphaRhoPhi,
98  phi,
99  transport,
100  propertiesName
101  ),
102 
103  Cmu_
104  (
106  (
107  "Cmu",
108  coeffDict_,
109  0.09
110  )
111  ),
112  C1_
113  (
115  (
116  "C1",
117  coeffDict_,
118  1.44
119  )
120  ),
121  C2_
122  (
124  (
125  "C2",
126  coeffDict_,
127  1.92
128  )
129  ),
130  sigmaZeta_
131  (
133  (
134  "sigmaZeta",
135  coeffDict_,
136  1.3
137  )
138  ),
140  (
142  (
143  "anisotropic",
144  coeffDict_,
145  false
146  )
147  ),
148 
149  qMin_("qMin", sqrt(kMin_)),
150  zetaMin_("zetaMin", epsilonMin_/(2*qMin_)),
151 
152  k_
153  (
154  IOobject
155  (
156  IOobject::groupName("k", U.group()),
157  runTime_.timeName(),
158  mesh_,
161  ),
162  mesh_
163  ),
164 
165  epsilon_
166  (
167  IOobject
168  (
169  IOobject::groupName("epsilon", U.group()),
170  runTime_.timeName(),
171  mesh_,
174  ),
175  mesh_
176  ),
177 
178  q_
179  (
180  IOobject
181  (
182  IOobject::groupName("q", U.group()),
183  runTime_.timeName(),
184  mesh_,
187  ),
188  sqrt(bound(k_, kMin_)),
189  k_.boundaryField().types()
190  ),
191 
192  zeta_
193  (
194  IOobject
195  (
196  IOobject::groupName("zeta", U.group()),
197  runTime_.timeName(),
198  mesh_,
201  ),
202  bound(epsilon_, epsilonMin_)/(2.0*q_),
203  epsilon_.boundaryField().types()
204  )
205 {
206  bound(zeta_, zetaMin_);
207 
208  if (type == typeName)
209  {
210  correctNut();
211  printCoeffs(type);
212  }
213 }
214 
215 
216 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
217 
219 {
221  {
222  Cmu_.readIfPresent(coeffDict());
223  C1_.readIfPresent(coeffDict());
224  C2_.readIfPresent(coeffDict());
225  sigmaZeta_.readIfPresent(coeffDict());
226  anisotropic_.readIfPresent("anisotropic", coeffDict());
227 
228  qMin_.readIfPresent(*this);
229  zetaMin_.readIfPresent(*this);
230 
231  return true;
232  }
233  else
234  {
235  return false;
236  }
237 }
238 
239 
241 {
242  if (!turbulence_)
243  {
244  return;
245  }
246 
248 
249  volScalarField G(GName(), nut_/(2.0*q_)*2*magSqr(symm(fvc::grad(U_))));
251 
252  // Zeta equation
253  tmp<fvScalarMatrix> zetaEqn
254  (
255  fvm::ddt(zeta_)
256  + fvm::div(phi_, zeta_)
258  ==
259  (2.0*C1_ - 1)*G*zeta_/q_
260  - fvm::SuSp((2.0*C2_*f2() - dimensionedScalar(1.0))*zeta_/q_, zeta_)
261  + E
262  );
263 
264  zetaEqn().relax();
265  solve(zetaEqn);
266  bound(zeta_, zetaMin_);
267 
268 
269  // q equation
271  (
272  fvm::ddt(q_)
273  + fvm::div(phi_, q_)
274  - fvm::laplacian(DqEff(), q_)
275  ==
276  G - fvm::Sp(zeta_/q_, q_)
277  );
278 
279  qEqn().relax();
280  solve(qEqn);
281  bound(q_, qMin_);
282 
283 
284  // Re-calculate k and epsilon
285  k_ = sqr(q_);
287 
288  epsilon_ = 2*q_*zeta_;
290 
291  correctNut();
292 }
293 
294 
295 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
296 
297 } // End namespace RASModels
298 } // End namespace incompressible
299 } // End namespace Foam
300 
301 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensionedScalar zetaMin_
Lower limit of zeta.
Definition: qZeta.H:94
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
Bound the given scalar field if it has gone unbounded.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
volScalarField & bound(volScalarField &, const dimensionedScalar &lowerBound)
Bound the given scalar field if it has gone unbounded.
Definition: bound.C:33
static Switch lookupOrAddToDict(const word &, dictionary &, const Switch &defaultValue=false)
Construct from dictionary, supplying default value so that if the.
Definition: Switch.C:107
A class for handling words, derived from string.
Definition: word.H:59
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:46
tmp< volScalarField > magSqrGradGrad(const GeometricField< Type, fvPatchField, volMesh > &vf)
virtual void correct()=0
Solve the turbulence equations and correct the turbulence viscosity.
dimensionedScalar exp(const dimensionedScalar &ds)
tmp< volScalarField > DqEff() const
Return the effective diffusivity for q.
Definition: qZeta.H:168
tmp< fvMatrix< Type > > Sp(const DimensionedField< scalar, volMesh > &, const GeometricField< Type, fvPatchField, volMesh > &)
Namespace for OpenFOAM.
static word groupName(Name name, const word &group)
tmp< volScalarField > DzetaEff() const
Return the effective diffusivity for epsilon.
Definition: qZeta.H:177
Generic dimensioned Type class.
dimensionedScalar sigmaZeta_
Definition: qZeta.H:87
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
tmp< fvMatrix< Type > > SuSp(const DimensionedField< scalar, volMesh > &, const GeometricField< Type, fvPatchField, volMesh > &)
dimensionedScalar qMin_
Lower limit of q.
Definition: qZeta.H:91
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
bool readIfPresent(const dictionary &)
Update the value of dimensioned<Type> if found in the dictionary.
Macros for easy insertion into run-time selection tables.
addToRunTimeSelectionTable(RASModel, kkLOmega, dictionary)
tmp< volScalarField > f2() const
Definition: qZeta.C:63
fileName::Type type(const fileName &)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:589
surfaceScalarField & phi
bool readIfPresent(const word &, const dictionary &)
Update the value of the Switch if it is found in the dictionary.
Definition: Switch.C:131
solverPerformance solve(fvMatrix< Type > &, const dictionary &)
Solve returning the solution statistics given convergence tolerance.
void correctBoundaryConditions()
Correct boundary field.
incompressible::RASModel::transportModel transportModel
Definition: eddyViscosity.H:75
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
Definition: qZeta.C:240
volScalarField & nu
virtual bool read()
Read RASProperties dictionary.
Definition: qZeta.C:218
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
qZeta(const geometricOneField &alpha, const geometricOneField &rho, const volVectorField &U, const surfaceScalarField &alphaRhoPhi, const surfaceScalarField &phi, const transportModel &transport, const word &propertiesName=turbulenceModel::propertiesName, const word &type=typeName)
Construct from components.
Definition: qZeta.C:80
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
RASModel< turbulenceModel > RASModel
const dimensionedScalar G
Newtonian constant of gravitation.
U
Definition: pEqn.H:82
tmp< volScalarField > fMu() const
Definition: qZeta.C:46