LRR.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-2021 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 "LRR.H"
27 #include "fvModels.H"
28 #include "fvConstraints.H"
29 #include "wallDist.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 namespace RASModels
36 {
37 
38 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
39 
40 template<class BasicMomentumTransportModel>
42 {
43  this->nut_ = this->Cmu_*sqr(k_)/epsilon_;
44  this->nut_.correctBoundaryConditions();
45  fvConstraints::New(this->mesh_).constrain(this->nut_);
46 }
47 
48 
49 template<class BasicMomentumTransportModel>
51 {
52  return tmp<fvScalarMatrix>
53  (
54  new fvScalarMatrix
55  (
56  epsilon_,
57  dimVolume*this->rho_.dimensions()*epsilon_.dimensions()/dimTime
58  )
59  );
60 }
61 
62 
63 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
64 
65 template<class BasicMomentumTransportModel>
67 (
68  const alphaField& alpha,
69  const rhoField& rho,
70  const volVectorField& U,
71  const surfaceScalarField& alphaRhoPhi,
72  const surfaceScalarField& phi,
73  const transportModel& transport,
74  const word& type
75 )
76 :
78  (
79  type,
80  alpha,
81  rho,
82  U,
83  alphaRhoPhi,
84  phi,
85  transport
86  ),
87 
88  Cmu_
89  (
91  (
92  "Cmu",
93  this->coeffDict_,
94  0.09
95  )
96  ),
97  C1_
98  (
100  (
101  "C1",
102  this->coeffDict_,
103  1.8
104  )
105  ),
106  C2_
107  (
109  (
110  "C2",
111  this->coeffDict_,
112  0.6
113  )
114  ),
115  Ceps1_
116  (
118  (
119  "Ceps1",
120  this->coeffDict_,
121  1.44
122  )
123  ),
124  Ceps2_
125  (
127  (
128  "Ceps2",
129  this->coeffDict_,
130  1.92
131  )
132  ),
133  Cs_
134  (
136  (
137  "Cs",
138  this->coeffDict_,
139  0.25
140  )
141  ),
142  Ceps_
143  (
145  (
146  "Ceps",
147  this->coeffDict_,
148  0.15
149  )
150  ),
151 
152  wallReflection_
153  (
155  (
156  "wallReflection",
157  this->coeffDict_,
158  true
159  )
160  ),
161  kappa_
162  (
164  (
165  "kappa",
166  this->coeffDict_,
167  0.41
168  )
169  ),
170  Cref1_
171  (
173  (
174  "Cref1",
175  this->coeffDict_,
176  0.5
177  )
178  ),
179  Cref2_
180  (
182  (
183  "Cref2",
184  this->coeffDict_,
185  0.3
186  )
187  ),
188 
189  k_
190  (
191  IOobject
192  (
193  IOobject::groupName("k", alphaRhoPhi.group()),
194  this->runTime_.timeName(),
195  this->mesh_,
198  ),
199  0.5*tr(this->R_)
200  ),
201  epsilon_
202  (
203  IOobject
204  (
205  IOobject::groupName("epsilon", alphaRhoPhi.group()),
206  this->runTime_.timeName(),
207  this->mesh_,
210  ),
211  this->mesh_
212  )
213 {
214  if (type == typeName)
215  {
216  this->printCoeffs(type);
217 
218  this->boundNormalStress(this->R_);
219  bound(epsilon_, this->epsilonMin_);
220  k_ = 0.5*tr(this->R_);
221  }
222 }
223 
224 
225 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
226 
227 template<class BasicMomentumTransportModel>
229 {
231  {
232  Cmu_.readIfPresent(this->coeffDict());
233  C1_.readIfPresent(this->coeffDict());
234  C2_.readIfPresent(this->coeffDict());
235  Ceps1_.readIfPresent(this->coeffDict());
236  Ceps2_.readIfPresent(this->coeffDict());
237  Cs_.readIfPresent(this->coeffDict());
238  Ceps_.readIfPresent(this->coeffDict());
239 
240  wallReflection_.readIfPresent("wallReflection", this->coeffDict());
241  kappa_.readIfPresent(this->coeffDict());
242  Cref1_.readIfPresent(this->coeffDict());
243  Cref2_.readIfPresent(this->coeffDict());
244 
245  return true;
246  }
247  else
248  {
249  return false;
250  }
251 }
252 
253 
254 template<class BasicMomentumTransportModel>
256 {
258  (
259  "DREff",
260  (Cs_*(this->k_/this->epsilon_))*this->R_ + I*this->nu()
261  );
262 }
263 
264 
265 template<class BasicMomentumTransportModel>
267 {
269  (
270  "DepsilonEff",
271  (Ceps_*(this->k_/this->epsilon_))*this->R_ + I*this->nu()
272  );
273 }
274 
275 
276 template<class BasicMomentumTransportModel>
278 {
279  if (!this->turbulence_)
280  {
281  return;
282  }
283 
284  // Local references
285  const alphaField& alpha = this->alpha_;
286  const rhoField& rho = this->rho_;
287  const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
288  const volVectorField& U = this->U_;
289  volSymmTensorField& R = this->R_;
290  const Foam::fvModels& fvModels(Foam::fvModels::New(this->mesh_));
292  (
293  Foam::fvConstraints::New(this->mesh_)
294  );
295 
297 
298  tmp<volTensorField> tgradU(fvc::grad(U));
299  const volTensorField& gradU = tgradU();
300 
301  volSymmTensorField P(-twoSymm(R & gradU));
302  volScalarField G(this->GName(), 0.5*mag(tr(P)));
303 
304  // Update epsilon and G at the wall
305  epsilon_.boundaryFieldRef().updateCoeffs();
306 
307  // Dissipation equation
308  tmp<fvScalarMatrix> epsEqn
309  (
310  fvm::ddt(alpha, rho, epsilon_)
311  + fvm::div(alphaRhoPhi, epsilon_)
312  - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
313  ==
314  Ceps1_*alpha*rho*G*epsilon_/k_
315  - fvm::Sp(Ceps2_*alpha*rho*epsilon_/k_, epsilon_)
316  + epsilonSource()
317  + fvModels.source(alpha, rho, epsilon_)
318  );
319 
320  epsEqn.ref().relax();
321  fvConstraints.constrain(epsEqn.ref());
322  epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
323  solve(epsEqn);
324  fvConstraints.constrain(epsilon_);
325  bound(epsilon_, this->epsilonMin_);
326 
327 
328  // Correct the trace of the tensorial production to be consistent
329  // with the near-wall generation from the wall-functions
330  const fvPatchList& patches = this->mesh_.boundary();
331 
332  forAll(patches, patchi)
333  {
334  const fvPatch& curPatch = patches[patchi];
335 
336  if (isA<wallFvPatch>(curPatch))
337  {
338  forAll(curPatch, facei)
339  {
340  label celli = curPatch.faceCells()[facei];
341  P[celli] *= min
342  (
343  G[celli]/(0.5*mag(tr(P[celli])) + small),
344  1.0
345  );
346  }
347  }
348  }
349 
350  // Reynolds stress equation
352  (
353  fvm::ddt(alpha, rho, R)
354  + fvm::div(alphaRhoPhi, R)
355  - fvm::laplacian(alpha*rho*DREff(), R)
356  + fvm::Sp(C1_*alpha*rho*epsilon_/k_, R)
357  ==
358  alpha*rho*P
359  - (2.0/3.0*(1 - C1_)*I)*alpha*rho*epsilon_
360  - C2_*alpha*rho*dev(P)
361  + this->RSource()
362  + fvModels.source(alpha, rho, R)
363  );
364 
365  // Optionally add wall-refection term
366  if (wallReflection_)
367  {
368  const volVectorField& n_(wallDist::New(this->mesh_).n());
369  const volScalarField& y_(wallDist::New(this->mesh_).y());
370 
371  const volSymmTensorField reflect
372  (
373  Cref1_*R - ((Cref2_*C2_)*(k_/epsilon_))*dev(P)
374  );
375 
376  REqn.ref() +=
377  ((3*pow(Cmu_, 0.75)/kappa_)*(alpha*rho*sqrt(k_)/y_))
378  *dev(symm((n_ & reflect)*n_));
379  }
380 
381  REqn.ref().relax();
382  fvConstraints.constrain(REqn.ref());
383  solve(REqn);
384  fvConstraints.constrain(R);
385 
386  this->boundNormalStress(R);
387 
388  k_ = 0.5*tr(R);
389 
390  correctNut();
391 
392  // Correct wall shear-stresses when applying wall-functions
393  this->correctWallShearStress(R);
394 }
395 
396 
397 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
398 
399 } // End namespace RASModels
400 } // End namespace Foam
401 
402 // ************************************************************************* //
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
virtual tmp< fvScalarMatrix > epsilonSource() const
Source term for the epsilon equation.
Definition: LRR.C:50
tmp< volSymmTensorField > DREff() const
Return the effective diffusivity for R.
Definition: LRR.C:255
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))
static tmp< GeometricField< symmTensor, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< symmTensor >> &)
Return a temporary field constructed from name,.
dimensionedScalar sqrt(const dimensionedScalar &ds)
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:62
Generic dimensioned Type class.
patches[0]
const dimensionedScalar G
Newtonian constant of gravitation.
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
virtual bool read()
Read model coefficients if they have changed.
Definition: LRR.C:228
BasicMomentumTransportModel::alphaField alphaField
Definition: LRR.H:143
static Switch lookupOrAddToDict(const word &, dictionary &, const Switch &defaultValue=false)
Construct from dictionary, supplying default value so that if the.
Definition: Switch.C:111
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
tmp< volSymmTensorField > DepsilonEff() const
Return the effective diffusivity for epsilon.
Definition: LRR.C:266
Templated abstract base class for RAS turbulence models.
Definition: RASModel.H:49
const dimensionSet dimTime
scalar y
bool read(const char *, int32_t &)
Definition: int32IO.C:85
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
static const Identity< scalar > I
Definition: Identity.H:93
Foam::fvConstraints & fvConstraints
BasicMomentumTransportModel::transportModel transportModel
Definition: LRR.H:145
virtual void correctNut()
Update the eddy-viscosity.
Definition: LRR.C:41
A class for handling words, derived from string.
Definition: word.H:59
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-9/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
Definition: pEqn.H:68
virtual const labelUList & faceCells() const
Return faceCells.
Definition: fvPatch.C:99
static word groupName(Name name, const word &group)
Finite volume models.
Definition: fvModels.H:60
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
phi
Definition: correctPhi.H:3
dimensionedScalar tr(const dimensionedSphericalTensor &dt)
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
static autoPtr< dictionary > New(Istream &)
Construct top-level dictionary on freestore from Istream.
Definition: dictionaryIO.C:96
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
bool constrain(fvMatrix< Type > &eqn) const
Apply constraints to an equation.
BasicMomentumTransportModel::rhoField rhoField
Definition: LRR.H:144
Foam::fvModels & fvModels
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
virtual void correct()
Solve the turbulence equations and correct eddy-Viscosity and.
Definition: LRR.C:277
label patchi
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
#define R(A, B, C, D, E, F, K, M)
Finite volume constraints.
Definition: fvConstraints.H:57
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:46
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
const dimensionSet dimVolume
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
LRR(const alphaField &alpha, const rhoField &rho, const volVectorField &U, const surfaceScalarField &alphaRhoPhi, const surfaceScalarField &phi, const transportModel &transport, const word &type=typeName)
Construct from components.
Definition: LRR.C:67
A class for managing temporary objects.
Definition: PtrList.H:53
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Namespace for OpenFOAM.
Reynolds-stress turbulence model base class.