qZeta.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2020 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 // * * * * * * * * * * * Protected 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& type
88 )
89 :
91  (
92  type,
93  alpha,
94  rho,
95  U,
96  alphaRhoPhi,
97  phi,
98  transport
99  ),
100 
101  Cmu_
102  (
104  (
105  "Cmu",
106  coeffDict_,
107  0.09
108  )
109  ),
110  C1_
111  (
113  (
114  "C1",
115  coeffDict_,
116  1.44
117  )
118  ),
119  C2_
120  (
122  (
123  "C2",
124  coeffDict_,
125  1.92
126  )
127  ),
128  sigmaZeta_
129  (
131  (
132  "sigmaZeta",
133  coeffDict_,
134  1.3
135  )
136  ),
138  (
140  (
141  "anisotropic",
142  coeffDict_,
143  false
144  )
145  ),
146 
147  qMin_("qMin", sqrt(kMin_)),
148  zetaMin_("zetaMin", epsilonMin_/(2*qMin_)),
149 
150  k_
151  (
152  IOobject
153  (
154  IOobject::groupName("k", alphaRhoPhi.group()),
155  runTime_.timeName(),
156  mesh_,
159  ),
160  mesh_
161  ),
162 
163  epsilon_
164  (
165  IOobject
166  (
167  IOobject::groupName("epsilon", alphaRhoPhi.group()),
168  runTime_.timeName(),
169  mesh_,
172  ),
173  mesh_
174  ),
175 
176  q_
177  (
178  IOobject
179  (
180  IOobject::groupName("q", alphaRhoPhi.group()),
181  runTime_.timeName(),
182  mesh_,
185  ),
186  sqrt(bound(k_, kMin_)),
187  k_.boundaryField().types()
188  ),
189 
190  zeta_
191  (
192  IOobject
193  (
194  IOobject::groupName("zeta", alphaRhoPhi.group()),
195  runTime_.timeName(),
196  mesh_,
199  ),
200  bound(epsilon_, epsilonMin_)/(2.0*q_),
201  epsilon_.boundaryField().types()
202  )
203 {
204  bound(zeta_, zetaMin_);
205 
206  if (type == typeName)
207  {
208  printCoeffs(type);
209  }
210 }
211 
212 
213 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
214 
216 {
218  {
219  Cmu_.readIfPresent(coeffDict());
220  C1_.readIfPresent(coeffDict());
221  C2_.readIfPresent(coeffDict());
222  sigmaZeta_.readIfPresent(coeffDict());
223  anisotropic_.readIfPresent("anisotropic", coeffDict());
224 
225  qMin_.readIfPresent(*this);
226  zetaMin_.readIfPresent(*this);
227 
228  return true;
229  }
230  else
231  {
232  return false;
233  }
234 }
235 
236 
238 {
239  if (!turbulence_)
240  {
241  return;
242  }
243 
245 
246  volScalarField G(GName(), nut_/(2.0*q_)*2*magSqr(symm(fvc::grad(U_))));
247  const volScalarField E(nu()*nut_/q_*fvc::magSqrGradGrad(U_));
248 
249  // Zeta equation
250  tmp<fvScalarMatrix> zetaEqn
251  (
252  fvm::ddt(zeta_)
253  + fvm::div(phi_, zeta_)
255  ==
256  (2.0*C1_ - 1)*G*zeta_/q_
257  - fvm::SuSp((2.0*C2_*f2() - dimensionedScalar(1.0))*zeta_/q_, zeta_)
258  + E
259  );
260 
261  zetaEqn.ref().relax();
262  solve(zetaEqn);
263  bound(zeta_, zetaMin_);
264 
265 
266  // q equation
268  (
269  fvm::ddt(q_)
270  + fvm::div(phi_, q_)
271  - fvm::laplacian(DqEff(), q_)
272  ==
273  G - fvm::Sp(zeta_/q_, q_)
274  );
275 
276  qEqn.ref().relax();
277  solve(qEqn);
278  bound(q_, qMin_);
279 
280 
281  // Re-calculate k and epsilon
282  k_ = sqr(q_);
284 
285  epsilon_ = 2*q_*zeta_;
287 
288  correctNut();
289 }
290 
291 
292 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
293 
294 } // End namespace RASModels
295 } // End namespace incompressible
296 } // End namespace Foam
297 
298 // ************************************************************************* //
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
tmp< fvMatrix< Type > > SuSp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
bool readIfPresent(const word &, const dictionary &)
Update the value of the Switch if it is found in the dictionary.
Definition: Switch.C:135
const Boundary & boundaryField() const
Return const-reference to the boundary field.
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
dimensionedSymmTensor sqr(const dimensionedVector &dv)
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
dimensionedScalar sqrt(const dimensionedScalar &ds)
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
Definition: qZeta.C:237
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.
const dimensionedScalar G
Newtonian constant of gravitation.
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
tmp< volScalarField > DqEff() const
Return the effective diffusivity for q.
Definition: qZeta.H:164
static Switch lookupOrAddToDict(const word &, dictionary &, const Switch &defaultValue=false)
Construct from dictionary, supplying default value so that if the.
Definition: Switch.C:111
Macros for easy insertion into run-time selection tables.
tmp< volScalarField > DzetaEff() const
Return the effective diffusivity for epsilon.
Definition: qZeta.H:174
virtual void correct()=0
Solve the turbulence equations and correct the turbulence viscosity.
tmp< volScalarField > f2() const
Definition: qZeta.C:63
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
dimensionedScalar exp(const dimensionedScalar &ds)
incompressible::RASModel ::transportModel transportModel
Definition: eddyViscosity.H:72
A class for handling words, derived from string.
Definition: word.H:59
static word groupName(Name name, const word &group)
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
phi
Definition: correctPhi.H:3
RASModel< momentumTransportModel > RASModel
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)
U
Definition: pEqn.H:72
volScalarField & bound(volScalarField &, const dimensionedScalar &lowerBound)
Bound the given scalar field if it has gone unbounded.
Definition: bound.C:33
tmp< volScalarField > fMu() const
Definition: qZeta.C:46
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
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
void correctBoundaryConditions()
Correct boundary field.
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.
qZeta(const geometricOneField &alpha, const geometricOneField &rho, const volVectorField &U, const surfaceScalarField &alphaRhoPhi, const surfaceScalarField &phi, const transportModel &transport, const word &type=typeName)
Construct from components.
Definition: qZeta.C:80
A class for managing temporary objects.
Definition: PtrList.H:53
virtual bool read()
Read RASProperties dictionary.
Definition: qZeta.C:215
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
dimensionedScalar sigmaZeta_
Definition: qZeta.H:84
dimensionedScalar qMin_
Lower limit of q.
Definition: qZeta.H:88
dimensionedScalar zetaMin_
Lower limit of zeta.
Definition: qZeta.H:91
Namespace for OpenFOAM.