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-2016 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  printCoeffs(type);
211  }
212 }
213 
214 
215 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
216 
218 {
220  {
221  Cmu_.readIfPresent(coeffDict());
222  C1_.readIfPresent(coeffDict());
223  C2_.readIfPresent(coeffDict());
224  sigmaZeta_.readIfPresent(coeffDict());
225  anisotropic_.readIfPresent("anisotropic", coeffDict());
226 
227  qMin_.readIfPresent(*this);
228  zetaMin_.readIfPresent(*this);
229 
230  return true;
231  }
232  else
233  {
234  return false;
235  }
236 }
237 
238 
240 {
241  if (!turbulence_)
242  {
243  return;
244  }
245 
247 
248  volScalarField G(GName(), nut_/(2.0*q_)*2*magSqr(symm(fvc::grad(U_))));
250 
251  // Zeta equation
252  tmp<fvScalarMatrix> zetaEqn
253  (
254  fvm::ddt(zeta_)
255  + fvm::div(phi_, zeta_)
257  ==
258  (2.0*C1_ - 1)*G*zeta_/q_
259  - fvm::SuSp((2.0*C2_*f2() - dimensionedScalar(1.0))*zeta_/q_, zeta_)
260  + E
261  );
262 
263  zetaEqn.ref().relax();
264  solve(zetaEqn);
265  bound(zeta_, zetaMin_);
266 
267 
268  // q equation
270  (
271  fvm::ddt(q_)
272  + fvm::div(phi_, q_)
273  - fvm::laplacian(DqEff(), q_)
274  ==
275  G - fvm::Sp(zeta_/q_, q_)
276  );
277 
278  qEqn.ref().relax();
279  solve(qEqn);
280  bound(q_, qMin_);
281 
282 
283  // Re-calculate k and epsilon
284  k_ = sqr(q_);
286 
287  epsilon_ = 2*q_*zeta_;
289 
290  correctNut();
291 }
292 
293 
294 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
295 
296 } // End namespace RASModels
297 } // End namespace incompressible
298 } // End namespace Foam
299 
300 // ************************************************************************* //
surfaceScalarField & phi
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
U
Definition: pEqn.H:83
bool readIfPresent(const word &, const dictionary &)
Update the value of the Switch if it is found in the dictionary.
Definition: Switch.C:129
tmp< fvMatrix< Type > > Sp(const DimensionedField< scalar, volMesh > &, const GeometricField< Type, fvPatchField, volMesh > &)
const dimensionedScalar G
Newtonian constant of gravitation.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
RASModel< turbulenceModel > RASModel
dimensionedScalar sqrt(const dimensionedScalar &ds)
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
Definition: qZeta.C:239
tmp< volScalarField > fMu() const
Definition: qZeta.C:46
tmp< volScalarField > magSqrGradGrad(const GeometricField< Type, fvPatchField, volMesh > &vf)
bool readIfPresent(const dictionary &)
Update the value of dimensioned<Type> if found in the dictionary.
Generic dimensioned Type class.
static Switch lookupOrAddToDict(const word &, dictionary &, const Switch &defaultValue=false)
Construct from dictionary, supplying default value so that if the.
Definition: Switch.C:105
Macros for easy insertion into run-time selection tables.
tmp< volScalarField > f2() const
Definition: qZeta.C:63
incompressible::RASModel::transportModel transportModel
Definition: eddyViscosity.H:75
virtual void correct()=0
Solve the turbulence equations and correct the turbulence viscosity.
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
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
dimensionedScalar exp(const dimensionedScalar &ds)
const Boundary & boundaryField() const
Return const-reference to the boundary field.
tmp< fvMatrix< Type > > SuSp(const DimensionedField< scalar, volMesh > &, const GeometricField< Type, fvPatchField, volMesh > &)
A class for handling words, derived from string.
Definition: word.H:59
SolverPerformance< Type > solve(fvMatrix< Type > &, const dictionary &)
Solve returning the solution statistics given convergence tolerance.
static word groupName(Name name, const word &group)
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
dimensioned< scalar > magSqr(const dimensioned< Type > &)
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
Bound the given scalar field if it has gone unbounded.
dimensionedScalar pow3(const dimensionedScalar &ds)
tmp< volScalarField > DqEff() const
Return the effective diffusivity for q.
Definition: qZeta.H:168
volScalarField & bound(volScalarField &, const dimensionedScalar &lowerBound)
Bound the given scalar field if it has gone unbounded.
Definition: bound.C:33
fileName::Type type(const fileName &)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:461
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
addToRunTimeSelectionTable(RASModel, kkLOmega, dictionary)
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:46
void correctBoundaryConditions()
Correct boundary field.
tmp< volScalarField > DzetaEff() const
Return the effective diffusivity for epsilon.
Definition: qZeta.H:177
A class for managing temporary objects.
Definition: PtrList.H:54
virtual bool read()
Read RASProperties dictionary.
Definition: qZeta.C:217
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
dimensionedScalar sigmaZeta_
Definition: qZeta.H:87
dimensionedScalar qMin_
Lower limit of q.
Definition: qZeta.H:91
volScalarField & nu
dimensionedScalar zetaMin_
Lower limit of zeta.
Definition: qZeta.H:94
Namespace for OpenFOAM.