kEpsilonLopesdaCosta.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) 2018-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 "kEpsilonLopesdaCosta.H"
27 #include "fvModels.H"
28 #include "fvConstraints.H"
29 #include "explicitPorositySource.H"
30 #include "bound.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 namespace RASModels
37 {
38 
39 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
40 
41 template<class BasicMomentumTransportModel>
43 (
46 )
47 {
48  if (pm.dict().found(C.name()))
49  {
50  const labelList& cellZoneIDs = pm.cellZoneIDs();
51 
52  const scalar Cpm = pm.dict().lookup<scalar>(C.name());
53 
54  forAll(cellZoneIDs, zonei)
55  {
56  const labelList& cells =
57  this->mesh_.cellZones()[cellZoneIDs[zonei]];
58 
59  forAll(cells, i)
60  {
61  const label celli = cells[i];
62  C[celli] = Cpm;
63  }
64  }
65  }
66 }
67 
68 
69 template<class BasicMomentumTransportModel>
71 (
74 )
75 {
76  if (pm.dict().found(C.name()))
77  {
78  const labelList& cellZoneIDs = pm.cellZoneIDs();
79  const scalarField& Av = pm.Av();
80 
81  const scalar Cpm = pm.dict().lookup<scalar>(C.name());
82 
83  forAll(cellZoneIDs, zonei)
84  {
85  const labelList& cells =
86  this->mesh_.cellZones()[cellZoneIDs[zonei]];
87 
88  forAll(cells, i)
89  {
90  const label celli = cells[i];
91  C[celli] = Cpm*Av[celli];
92  }
93  }
94  }
95 }
96 
97 
98 template<class BasicMomentumTransportModel>
101 {
103  (
104  Foam::fvModels::New(this->mesh_)
105  );
106 
107  forAll(fvModels, i)
108  {
109  if (isA<fv::explicitPorositySource>(fvModels[i]))
110  {
111  const fv::explicitPorositySource& eps =
112  refCast<const fv::explicitPorositySource>(fvModels[i]);
113 
114  if (isA<porosityModels::powerLawLopesdaCosta>(eps.model()))
115  {
117  refCast<const porosityModels::powerLawLopesdaCosta>
118  (
119  eps.model()
120  );
121 
122  setPorosityCoefficient(Cmu_, pm);
123  setPorosityCoefficient(C1_, pm);
124  setPorosityCoefficient(C2_, pm);
125  setPorosityCoefficient(sigmak_, pm);
126  setPorosityCoefficient(sigmaEps_, pm);
127 
128  setCdAv(CdAv_, pm);
129  setPorosityCoefficient(betap_, pm);
130  setPorosityCoefficient(betad_, pm);
131  setPorosityCoefficient(C4_, pm);
132  setPorosityCoefficient(C5_, pm);
133  }
134  }
135  }
136 }
137 
138 
139 template<class BasicMomentumTransportModel>
141 {
142  this->nut_ = Cmu_*sqr(k_)/epsilon_;
143  this->nut_.correctBoundaryConditions();
144  fvConstraints::New(this->mesh_).constrain(this->nut_);
145 }
146 
147 
148 template<class BasicMomentumTransportModel>
150 (
151  const volScalarField::Internal& magU,
152  const volScalarField::Internal& magU3
153 ) const
154 {
155  return fvm::Su(CdAv_*(betap_*magU3 - betad_*magU*k_()), k_);
156 }
157 
158 
159 template<class BasicMomentumTransportModel>
162 (
163  const volScalarField::Internal& magU,
164  const volScalarField::Internal& magU3
165 ) const
166 {
167  return fvm::Su
168  (
169  CdAv_
170  *(C4_*betap_*epsilon_()/k_()*magU3 - C5_*betad_*magU*epsilon_()),
171  epsilon_
172  );
173 }
174 
175 
176 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
177 
178 template<class BasicMomentumTransportModel>
180 (
181  const alphaField& alpha,
182  const rhoField& rho,
183  const volVectorField& U,
184  const surfaceScalarField& alphaRhoPhi,
185  const surfaceScalarField& phi,
186  const viscosity& viscosity,
187  const word& type
188 )
189 :
190  eddyViscosity<RASModel<BasicMomentumTransportModel>>
191  (
192  type,
193  alpha,
194  rho,
195  U,
196  alphaRhoPhi,
197  phi,
198  viscosity
199  ),
200 
201  Cmu_
202  (
203  IOobject
204  (
205  "Cmu",
206  this->runTime_.name(),
207  this->mesh_
208  ),
209  this->mesh_,
210  dimensioned<scalar>::lookupOrAddToDict
211  (
212  "Cmu",
213  this->coeffDict_,
214  0.09
215  )
216  ),
217  C1_
218  (
219  IOobject
220  (
221  "C1",
222  this->runTime_.name(),
223  this->mesh_
224  ),
225  this->mesh_,
226  dimensioned<scalar>::lookupOrAddToDict
227  (
228  "C1",
229  this->coeffDict_,
230  1.44
231  )
232  ),
233  C2_
234  (
235  IOobject
236  (
237  "C2",
238  this->runTime_.name(),
239  this->mesh_
240  ),
241  this->mesh_,
242  dimensioned<scalar>::lookupOrAddToDict
243  (
244  "C2",
245  this->coeffDict_,
246  1.92
247  )
248  ),
249  sigmak_
250  (
251  IOobject
252  (
253  "sigmak",
254  this->runTime_.name(),
255  this->mesh_
256  ),
257  this->mesh_,
258  dimensioned<scalar>::lookupOrAddToDict
259  (
260  "sigmak",
261  this->coeffDict_,
262  1.0
263  )
264  ),
265  sigmaEps_
266  (
267  IOobject
268  (
269  "sigmaEps",
270  this->runTime_.name(),
271  this->mesh_
272  ),
273  this->mesh_,
274  dimensioned<scalar>::lookupOrAddToDict
275  (
276  "sigmaEps",
277  this->coeffDict_,
278  1.3
279  )
280  ),
281 
282  CdAv_
283  (
284  IOobject
285  (
286  "CdAv",
287  this->runTime_.name(),
288  this->mesh_
289  ),
290  this->mesh_,
292  ),
293  betap_
294  (
295  IOobject
296  (
297  "betap",
298  this->runTime_.name(),
299  this->mesh_
300  ),
301  this->mesh_,
303  ),
304  betad_
305  (
306  IOobject
307  (
308  "betad",
309  this->runTime_.name(),
310  this->mesh_
311  ),
312  this->mesh_,
314  ),
315  C4_
316  (
317  IOobject
318  (
319  "C4",
320  this->runTime_.name(),
321  this->mesh_
322  ),
323  this->mesh_,
325  ),
326  C5_
327  (
328  IOobject
329  (
330  "C5",
331  this->runTime_.name(),
332  this->mesh_
333  ),
334  this->mesh_,
336  ),
337 
338  k_
339  (
340  IOobject
341  (
342  IOobject::groupName("k", alphaRhoPhi.group()),
343  this->runTime_.name(),
344  this->mesh_,
345  IOobject::MUST_READ,
346  IOobject::AUTO_WRITE
347  ),
348  this->mesh_
349  ),
350  epsilon_
351  (
352  IOobject
353  (
354  IOobject::groupName("epsilon", alphaRhoPhi.group()),
355  this->runTime_.name(),
356  this->mesh_,
357  IOobject::MUST_READ,
358  IOobject::AUTO_WRITE
359  ),
360  this->mesh_
361  )
362 {
363  bound(k_, this->kMin_);
364  bound(epsilon_, this->epsilonMin_);
365 
366  if (type == typeName)
367  {
368  this->printCoeffs(type);
369  }
370 
372 }
373 
374 
375 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
376 
377 template<class BasicMomentumTransportModel>
379 {
381  {
382  return true;
383  }
384  else
385  {
386  return false;
387  }
388 }
389 
390 
391 template<class BasicMomentumTransportModel>
393 {
394  if (!this->turbulence_)
395  {
396  return;
397  }
398 
399  // Local references
400  const alphaField& alpha = this->alpha_;
401  const rhoField& rho = this->rho_;
402  const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
403  const volVectorField& U = this->U_;
404  volScalarField& nut = this->nut_;
405  const Foam::fvModels& fvModels(Foam::fvModels::New(this->mesh_));
407  (
408  Foam::fvConstraints::New(this->mesh_)
409  );
410 
412 
414  (
415  fvc::div(fvc::absolute(this->phi(), U))().v()
416  );
417 
418  tmp<volTensorField> tgradU = fvc::grad(U);
420  (
421  this->GName(),
422  nut.v()*(dev(twoSymm(tgradU().v())) && tgradU().v())
423  );
424  tgradU.clear();
425 
426  // Update epsilon and G at the wall
427  epsilon_.boundaryFieldRef().updateCoeffs();
428 
430  volScalarField::Internal magU3(pow3(magU));
431 
432  // Dissipation equation
433  tmp<fvScalarMatrix> epsEqn
434  (
435  fvm::ddt(alpha, rho, epsilon_)
436  + fvm::div(alphaRhoPhi, epsilon_)
437  - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
438  ==
439  C1_*alpha()*rho()*G*epsilon_()/k_()
440  - fvm::SuSp(((2.0/3.0)*C1_)*alpha()*rho()*divU, epsilon_)
441  - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_)
442  + epsilonSource(magU, magU3)
443  + fvModels.source(alpha, rho, epsilon_)
444  );
445 
446  epsEqn.ref().relax();
447  fvConstraints.constrain(epsEqn.ref());
448  epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
449  solve(epsEqn);
450  fvConstraints.constrain(epsilon_);
451  bound(epsilon_, this->epsilonMin_);
452 
453  // Turbulent kinetic energy equation
455  (
456  fvm::ddt(alpha, rho, k_)
457  + fvm::div(alphaRhoPhi, k_)
458  - fvm::laplacian(alpha*rho*DkEff(), k_)
459  ==
460  alpha()*rho()*G
461  - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
462  - fvm::Sp(alpha()*rho()*epsilon_()/k_(), k_)
463  + kSource(magU, magU3)
464  + fvModels.source(alpha, rho, k_)
465  );
466 
467  kEqn.ref().relax();
468  fvConstraints.constrain(kEqn.ref());
469  solve(kEqn);
471  bound(k_, this->kMin_);
472 
473  correctNut();
474 }
475 
476 
477 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
478 
479 } // End namespace RASModels
480 } // End namespace Foam
481 
482 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Bound the given scalar field where it is below the specified minimum.
Graphite solid properties.
Definition: C.H:51
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
Generic GeometricField class.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
Templated abstract base class for RAS turbulence models.
Definition: RASModel.H:56
BasicMomentumTransportModel::alphaField alphaField
virtual tmp< fvScalarMatrix > kSource(const volScalarField::Internal &magU, const volScalarField::Internal &magU3) const
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
kEpsilonLopesdaCosta(const alphaField &alpha, const rhoField &rho, const volVectorField &U, const surfaceScalarField &alphaRhoPhi, const surfaceScalarField &phi, const viscosity &viscosity, const word &type=typeName)
Construct from components.
void setCdAv(volScalarField::Internal &C, const porosityModels::powerLawLopesdaCosta &pm)
virtual tmp< fvScalarMatrix > epsilonSource(const volScalarField::Internal &magU, const volScalarField::Internal &magU3) const
void setPorosityCoefficient(volScalarField::Internal &C, const porosityModels::powerLawLopesdaCosta &pm)
virtual bool read()
Re-read model coefficients if they have changed.
BasicMomentumTransportModel::rhoField rhoField
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:860
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:659
static autoPtr< dictionary > New(Istream &)
Construct top-level dictionary on freestore from Istream.
Definition: dictionaryIO.C:100
Generic dimensioned Type class.
Eddy viscosity turbulence model base class.
Definition: eddyViscosity.H:52
Finite volume constraints.
Definition: fvConstraints.H:62
bool constrain(fvMatrix< Type > &eqn) const
Apply constraints to an equation.
Finite volume models.
Definition: fvModels.H:65
const porosityModel & model() const
Return the porosity model.
const labelList & cellZoneIDs() const
Return const access to the cell zone IDs.
const dictionary & dict() const
Return dictionary used for model construction.
const scalarField & Av() const
Return the porosity surface area per unit volume zone field.
Variant of the power law porosity model with spatially varying drag coefficient.
A class for managing temporary objects.
Definition: tmp.H:55
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:237
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
Abstract base class for all fluid physical properties.
Definition: viscosity.H:50
A class for handling words, derived from string.
Definition: word.H:62
Foam::fvConstraints & fvConstraints(Foam::fvConstraints::New(mesh))
Foam::fvModels & fvModels(Foam::fvModels::New(mesh))
const scalar nut
const cellShapeList & cells
U
Definition: pEqn.H:72
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
void correct(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiCorr, const SpType &Sp, const SuType &Su)
const char *const group
Group name for atomic constants.
const dimensionedScalar G
Newtonian constant of gravitation.
tmp< VolField< typename outerProduct< vector, Type >::type > > grad(const SurfaceField< Type > &ssf)
Definition: fvcGrad.C:46
tmp< VolField< Type > > div(const SurfaceField< Type > &ssf)
Definition: fvcDiv.C:47
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
Definition: fvcMeshPhi.C:202
tmp< fvMatrix< Type > > Su(const DimensionedField< Type, volMesh > &, const VolField< Type > &)
tmp< fvMatrix< Type > > laplacian(const VolField< Type > &vf, const word &name)
Definition: fvmLaplacian.C:47
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const VolField< Type > &)
tmp< fvMatrix< Type > > SuSp(const volScalarField::Internal &, const VolField< Type > &)
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const VolField< Type > &vf, const word &name)
Definition: fvmDiv.C:48
tmp< fvMatrix< Type > > ddt(const VolField< Type > &vf)
Definition: fvmDdt.C:46
Namespace for OpenFOAM.
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
bool read(const char *, int32_t &)
Definition: int32IO.C:85
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
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
const dimensionSet dimless
const dimensionSet dimLength
dimensioned< scalar > mag(const dimensioned< Type > &)
bool bound(volScalarField &, const dimensionedScalar &min)
Bound the given scalar field where it is below the specified min value.
Definition: bound.C:31
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
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.