LamBremhorstKE.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-2023 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 "LamBremhorstKE.H"
27 #include "wallDist.H"
28 #include "bound.H"
30 
32 (
33  geometricOneField,
34  geometricOneField,
35  incompressibleMomentumTransportModel
36 )
37 
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42 namespace incompressible
43 {
44 namespace RASModels
45 {
46 
47 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
48 
49 defineTypeNameAndDebug(LamBremhorstKE, 0);
51 (
52  RASincompressibleMomentumTransportModel,
53  LamBremhorstKE,
54  dictionary
55 );
56 
57 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
58 
59 tmp<volScalarField> LamBremhorstKE::boundEpsilon()
60 {
61  tmp<volScalarField> tCmuk2(Cmu_*sqr(k_));
62  epsilon_ = max(epsilon_, tCmuk2()/(nutMaxCoeff_*nu()));
63  return tCmuk2;
64 }
65 
66 
67 tmp<volScalarField> LamBremhorstKE::Rt() const
68 {
69  return sqr(k_)/(nu()*epsilon_);
70 }
71 
72 
73 tmp<volScalarField> LamBremhorstKE::fMu(const volScalarField& Rt) const
74 {
75  tmp<volScalarField> Ry(sqrt(k_)*y()/nu());
76  return sqr(scalar(1) - exp(-0.0165*Ry))*(scalar(1) + 20.5/(Rt + small));
77 }
78 
79 
80 tmp<volScalarField> LamBremhorstKE::f1(const volScalarField& fMu) const
81 {
82  return scalar(1) + pow3(0.05/(fMu + small));
83 }
84 
85 
86 tmp<volScalarField> LamBremhorstKE::f2(const volScalarField& Rt) const
87 {
88  return scalar(1) - exp(-sqr(Rt));
89 }
90 
91 
93 {
94  nut_ = fMu*boundEpsilon()/epsilon_;
96 }
97 
98 
99 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
100 
102 {
103  correctNut(fMu(Rt()));
104 }
105 
106 
107 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
108 
110 (
111  const geometricOneField& alpha,
112  const geometricOneField& rho,
113  const volVectorField& U,
114  const surfaceScalarField& alphaRhoPhi,
115  const surfaceScalarField& phi,
116  const viscosity& viscosity,
117  const word& type
118 )
119 :
120  eddyViscosity<incompressible::RASModel>
121  (
122  type,
123  alpha,
124  rho,
125  U,
126  alphaRhoPhi,
127  phi,
128  viscosity
129  ),
130 
131  Cmu_
132  (
134  (
135  "Cmu",
136  coeffDict_,
137  0.09
138  )
139  ),
140  Ceps1_
141  (
143  (
144  "Ceps1",
145  coeffDict_,
146  1.44
147  )
148  ),
149  Ceps2_
150  (
152  (
153  "Ceps2",
154  coeffDict_,
155  1.92
156  )
157  ),
158  sigmaEps_
159  (
161  (
162  "alphaEps",
163  coeffDict_,
164  1.3
165  )
166  ),
167 
168  k_
169  (
170  IOobject
171  (
172  this->groupName("k"),
173  runTime_.name(),
174  mesh_,
177  ),
178  mesh_
179  ),
180 
181  epsilon_
182  (
183  IOobject
184  (
185  this->groupName("epsilon"),
186  runTime_.name(),
187  mesh_,
190  ),
191  mesh_
192  )
193 {
194  bound(k_, kMin_);
195  boundEpsilon();
196 
197  if (type == typeName)
198  {
199  printCoeffs(type);
200  }
201 }
202 
203 
204 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
205 
207 {
209  {
210  Cmu_.readIfPresent(coeffDict());
211  Ceps1_.readIfPresent(coeffDict());
212  Ceps2_.readIfPresent(coeffDict());
213  sigmaEps_.readIfPresent(coeffDict());
214 
215  return true;
216  }
217  else
218  {
219  return false;
220  }
221 }
222 
223 
225 {
226  if (!turbulence_)
227  {
228  return;
229  }
230 
232 
233  tmp<volTensorField> tgradU = fvc::grad(U_);
234  volScalarField G(GName(), nut_*(twoSymm(tgradU()) && tgradU()));
235  tgradU.clear();
236 
237  // Update epsilon and G at the wall
239 
240  const volScalarField Rt(this->Rt());
241  const volScalarField fMu(this->fMu(Rt));
242 
243  // Dissipation equation
244  tmp<fvScalarMatrix> epsEqn
245  (
247  + fvm::div(phi_, epsilon_)
249  ==
250  Ceps1_*f1(fMu)*G*epsilon_/k_
251  - fvm::Sp(Ceps2_*f2(Rt)*epsilon_/k_, epsilon_)
252  );
253 
254  epsEqn.ref().relax();
255  epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
256  solve(epsEqn);
257  boundEpsilon();
258 
259  // Turbulent kinetic energy equation
260  tmp<fvScalarMatrix> kEqn
261  (
262  fvm::ddt(k_)
263  + fvm::div(phi_, k_)
264  - fvm::laplacian(DkEff(), k_)
265  ==
266  G - fvm::Sp(epsilon_/k_, k_)
267  );
268 
269  kEqn.ref().relax();
270  solve(kEqn);
271  bound(k_, kMin_);
272 
273  correctNut(fMu);
274 }
275 
276 
277 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
278 
279 } // End namespace RASModels
280 } // End namespace incompressible
281 } // End namespace Foam
282 
283 // ************************************************************************* //
scalar y
makeMomentumTransportModelTypes(geometricOneField, geometricOneField, incompressibleMomentumTransportModel)
Bound the given scalar field where it is below the specified minimum.
void updateCoeffs()
Update the boundary condition coefficients.
void correctBoundaryConditions()
Correct boundary field.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
static dimensioned< scalar > lookupOrAddToDict(const word &, dictionary &, const dimensionSet &dims=dimless, const scalar &defaultValue=pTraits< scalar >::zero)
Construct from dictionary, with default value.
bool readIfPresent(const dictionary &)
Update the value of dimensioned<Type> if found in the dictionary.
virtual bool read()=0
Re-read model coefficients if they have changed.
Definition: eddyViscosity.C:74
virtual void correct()=0
Solve the turbulence equations and correct the turbulence viscosity.
virtual void correctNut()
Correct the eddy-viscosity nut.
tmp< volScalarField > DkEff() const
Return the effective diffusivity for k.
tmp< volScalarField > boundEpsilon()
Bound epsilon and return Cmu*sqr(k) for nut.
LamBremhorstKE(const geometricOneField &alpha, const geometricOneField &rho, const volVectorField &U, const surfaceScalarField &alphaRhoPhi, const surfaceScalarField &phi, const viscosity &viscosity, const word &type=typeName)
Construct from components.
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
tmp< volScalarField > DepsilonEff() const
Return the effective diffusivity for epsilon.
virtual bool read()
Read RASProperties dictionary.
U
Definition: pEqn.H:72
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
const dimensionedScalar G
Newtonian constant of gravitation.
tmp< VolField< typename outerProduct< vector, Type >::type > > grad(const SurfaceField< Type > &ssf)
Definition: fvcGrad.C:46
tmp< fvMatrix< Type > > laplacian(const VolField< Type > &vf, const word &name)
Definition: fvmLaplacian.C:47
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const VolField< Type > &vf, const word &name)
Definition: fvmDiv.C:48
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const VolField< Type > &)
tmp< fvMatrix< Type > > ddt(const VolField< Type > &vf)
Definition: fvmDdt.C:46
Namespace for OpenFOAM.
dimensionedScalar exp(const dimensionedScalar &ds)
VolField< vector > volVectorField
Definition: volFieldsFwd.H:65
addToRunTimeSelectionTable(polyPatch, mergedCyclicPolyPatch, word)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
SurfaceField< scalar > surfaceScalarField
tensor Ry(const scalar &omega)
Rotational transformation tensor about the y-axis by omega radians.
Definition: transform.H:106
dimensionedScalar sqrt(const dimensionedScalar &ds)
VolField< scalar > volScalarField
Definition: volFieldsFwd.H:64
defineTypeNameAndDebug(combustionModel, 0)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
bool bound(volScalarField &, const dimensionedScalar &min)
Bound the given scalar field where it is below the specified min value.
Definition: bound.C:31
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.
Info<< "Reading field U\n"<< endl;volVectorField U(IOobject("U", runTime.name(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Creating face flux\n"<< endl;surfaceScalarField phi(IOobject("phi", runTime.name(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), mesh, dimensionedScalar(mesh.Sf().dimensions() *U.dimensions(), 0));autoPtr< viscosityModel > viscosity(viscosityModel::New(mesh))