LienLeschziner.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-2024 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 "LienLeschziner.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(LienLeschziner, 0);
51 (
52  RASincompressibleMomentumTransportModel,
53  LienLeschziner,
54  dictionary
55 );
56 
57 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
58 
59 tmp<volScalarField> LienLeschziner::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> LienLeschziner::fMu() const
68 {
69  const volScalarField yStar(sqrt(k_)*y()/nu());
70 
71  return
72  (scalar(1) - exp(-Anu_*yStar))
73  /((scalar(1) + small) - exp(-Aeps_*yStar));
74 }
75 
76 
77 tmp<volScalarField> LienLeschziner::f2() const
78 {
79  tmp<volScalarField> Rt = sqr(k_)/(nu()*epsilon_);
80 
81  return scalar(1) - 0.3*exp(-sqr(Rt));
82 }
83 
84 
85 tmp<volScalarField> LienLeschziner::E(const volScalarField& f2) const
86 {
87  const volScalarField yStar(sqrt(k_)*y()/nu());
88  const volScalarField le
89  (
90  kappa_*y()*((scalar(1) + small) - exp(-Aeps_*yStar))
91  );
92 
93  return
94  (Ceps2_*pow(Cmu_, 0.75))
95  *(f2*sqrt(k_)*epsilon_/le)*exp(-AE_*sqr(yStar));
96 }
97 
98 
100 {
103 }
104 
105 
106 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
107 
109 (
110  const geometricOneField& alpha,
111  const geometricOneField& rho,
112  const volVectorField& U,
113  const surfaceScalarField& alphaRhoPhi,
114  const surfaceScalarField& phi,
115  const viscosity& viscosity,
116  const word& type
117 )
118 :
119  eddyViscosity<incompressible::RASModel>
120  (
121  type,
122  alpha,
123  rho,
124  U,
125  alphaRhoPhi,
126  phi,
127  viscosity
128  ),
129 
130  Ceps1_("Ceps1", coeffDict(), 1.44),
131  Ceps2_("Ceps2", coeffDict(), 1.92),
132  sigmak_("sigmak", coeffDict(), 1.0),
133  sigmaEps_("sigmaEps", coeffDict(), 1.3),
134  Cmu_("Cmu", coeffDict(), 0.09),
135  kappa_("kappa", coeffDict(), 0.41),
136  Anu_("Anu", coeffDict(), 0.016),
137  Aeps_("Aeps", coeffDict(), 0.263),
138  AE_("AE", coeffDict(), 0.00222),
139 
140  k_
141  (
142  IOobject
143  (
144  this->groupName("k"),
145  runTime_.name(),
146  mesh_,
149  ),
150  mesh_
151  ),
152 
153  epsilon_
154  (
155  IOobject
156  (
157  this->groupName("epsilon"),
158  runTime_.name(),
159  mesh_,
162  ),
163  mesh_
164  )
165 {
166  bound(k_, kMin_);
167  boundEpsilon();
168 }
169 
170 
171 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
172 
174 {
176  {
177  Ceps1_.readIfPresent(coeffDict());
178  Ceps2_.readIfPresent(coeffDict());
179  sigmak_.readIfPresent(coeffDict());
180  sigmaEps_.readIfPresent(coeffDict());
181  Cmu_.readIfPresent(coeffDict());
182  kappa_.readIfPresent(coeffDict());
183  Anu_.readIfPresent(coeffDict());
184  Aeps_.readIfPresent(coeffDict());
185  AE_.readIfPresent(coeffDict());
186 
187  return true;
188  }
189  else
190  {
191  return false;
192  }
193 }
194 
195 
197 {
198  if (!turbulence_)
199  {
200  return;
201  }
202 
204 
205  tmp<volTensorField> tgradU = fvc::grad(U_);
207  (
208  GName(),
209  nut_*(tgradU() && twoSymm(tgradU()))
210  );
211  tgradU.clear();
212 
213  // Update epsilon and G at the wall
215 
216  const volScalarField f2(this->f2());
217 
218  // Dissipation equation
219  tmp<fvScalarMatrix> epsEqn
220  (
222  + fvm::div(phi_, epsilon_)
224  ==
227  + E(f2)
228  );
229 
230  epsEqn.ref().relax();
231  epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
232  solve(epsEqn);
233  boundEpsilon();
234 
235 
236  // Turbulent kinetic energy equation
237  tmp<fvScalarMatrix> kEqn
238  (
239  fvm::ddt(k_)
240  + fvm::div(phi_, k_)
241  - fvm::laplacian(DkEff(), k_)
242  ==
243  G
244  - fvm::Sp(epsilon_/k_, k_)
245  );
246 
247  kEqn.ref().relax();
248  solve(kEqn);
249  bound(k_, kMin_);
250 
251  correctNut();
252 }
253 
254 
255 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
256 
257 } // End namespace RASModels
258 } // End namespace incompressible
259 } // End namespace Foam
260 
261 // ************************************************************************* //
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.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
void correctBoundaryConditions()
Correct boundary field.
bool readIfPresent(const dictionary &, const unitConversion &defaultUnits=NullObjectRef< unitConversion >())
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.
tmp< volScalarField > E(const volScalarField &f2) const
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
tmp< volScalarField > DepsilonEff() const
Return the effective diffusivity for epsilon.
LienLeschziner(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 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:63
addToRunTimeSelectionTable(polyPatch, mergedCyclicPolyPatch, word)
void twoSymm(LagrangianPatchField< tensor > &f, const LagrangianPatchField< tensor > &f1)
SurfaceField< scalar > surfaceScalarField
void pow(LagrangianPatchField< typename powProduct< Type, r >::type > &f, const LagrangianPatchField< Type > &f1)
VolField< scalar > volScalarField
Definition: volFieldsFwd.H:62
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
void sqr(LagrangianPatchField< typename outerProduct< Type, Type >::type > &f, const LagrangianPatchField< Type > &f1)
void sqrt(LagrangianPatchField< scalar > &f, const LagrangianPatchField< scalar > &f1)
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))