kEpsilonLopesdaCosta.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2018 OpenFOAM Foundation
9 Copyright (C) 2018-2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
30#include "fvOptions.H"
32#include "bound.H"
33
34// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35
36namespace Foam
37{
38namespace RASModels
39{
40
41// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
42
43template<class BasicTurbulenceModel>
45(
48)
49{
50 if (pm.dict().found(C.name()))
51 {
52 const labelList& cellZoneIDs = pm.cellZoneIDs();
53
54 const scalar Cpm = pm.dict().get<scalar>(C.name());
55
56 for (const label zonei : cellZoneIDs)
57 {
58 const labelList& cells = this->mesh_.cellZones()[zonei];
59
60 for (const label celli : cells)
61 {
62 C[celli] = Cpm;
63 }
64 }
65 }
66}
67
68
69template<class BasicTurbulenceModel>
71(
74)
75{
76 if (pm.dict().found(C.name()))
77 {
78 const labelList& cellZoneIDs = pm.cellZoneIDs();
79 const scalarField& Sigma = pm.Sigma();
80
81 const scalar Cpm = pm.dict().get<scalar>(C.name());
82
83 for (const label zonei : cellZoneIDs)
84 {
85 const labelList& cells = this->mesh_.cellZones()[zonei];
86
87 forAll(cells, i)
88 {
89 const label celli = cells[i];
90 C[celli] = Cpm*Sigma[i];
91 }
92 }
93 }
94}
95
96
97template<class BasicTurbulenceModel>
99{
101
102 forAll(fvOptions, i)
103 {
104 if (isA<fv::explicitPorositySource>(fvOptions[i]))
105 {
106 const fv::explicitPorositySource& eps =
107 refCast<const fv::explicitPorositySource>(fvOptions[i]);
108
109 if (isA<porosityModels::powerLawLopesdaCosta>(eps.model()))
110 {
112 refCast<const porosityModels::powerLawLopesdaCosta>
113 (
114 eps.model()
115 );
116
117 setPorosityCoefficient(Cmu_, pm);
118 Cmu_.correctBoundaryConditions();
119 setPorosityCoefficient(C1_, pm);
120 setPorosityCoefficient(C2_, pm);
121 setPorosityCoefficient(sigmak_, pm);
122 setPorosityCoefficient(sigmaEps_, pm);
123 sigmak_.correctBoundaryConditions();
124 sigmaEps_.correctBoundaryConditions();
125
126 setCdSigma(CdSigma_, pm);
127 setPorosityCoefficient(betap_, pm);
128 setPorosityCoefficient(betad_, pm);
129 setPorosityCoefficient(C4_, pm);
130 setPorosityCoefficient(C5_, pm);
131 }
132 }
133 }
134}
135
136
137template<class BasicTurbulenceModel>
139{
140 this->nut_ = Cmu_*sqr(k_)/epsilon_;
141 this->nut_.correctBoundaryConditions();
142 fv::options::New(this->mesh_).correct(this->nut_);
143
144 BasicTurbulenceModel::correctNut();
145}
146
147
148template<class BasicTurbulenceModel>
150(
151 const volScalarField::Internal& magU,
152 const volScalarField::Internal& magU3
153) const
154{
155 return fvm::Su(CdSigma_*(betap_*magU3 - betad_*magU*k_()), k_);
156}
157
158
159template<class BasicTurbulenceModel>
162(
163 const volScalarField::Internal& magU,
164 const volScalarField::Internal& magU3
165) const
166{
167 return fvm::Su
168 (
169 CdSigma_
170 *(C4_*betap_*epsilon_()/k_()*magU3 - C5_*betad_*magU*epsilon_()),
171 epsilon_
172 );
173}
174
175
176// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
177
178template<class BasicTurbulenceModel>
180(
181 const alphaField& alpha,
182 const rhoField& rho,
183 const volVectorField& U,
184 const surfaceScalarField& alphaRhoPhi,
185 const surfaceScalarField& phi,
186 const transportModel& transport,
187 const word& propertiesName,
188 const word& type
189)
190:
191 eddyViscosity<RASModel<BasicTurbulenceModel>>
192 (
193 type,
194 alpha,
195 rho,
196 U,
197 alphaRhoPhi,
198 phi,
199 transport,
200 propertiesName
201 ),
202
203 Cmu_
204 (
206 (
207 "Cmu",
208 this->runTime_.timeName(),
209 this->mesh_
210 ),
211 this->mesh_,
212 dimensioned<scalar>::getOrAddToDict
213 (
214 "Cmu",
215 this->coeffDict_,
216 0.09
217 )
218 ),
219 C1_
220 (
222 (
223 "C1",
224 this->runTime_.timeName(),
225 this->mesh_
226 ),
227 this->mesh_,
228 dimensioned<scalar>::getOrAddToDict
229 (
230 "C1",
231 this->coeffDict_,
232 1.44
233 )
234 ),
235 C2_
236 (
238 (
239 "C2",
240 this->runTime_.timeName(),
241 this->mesh_
242 ),
243 this->mesh_,
244 dimensioned<scalar>::getOrAddToDict
245 (
246 "C2",
247 this->coeffDict_,
248 1.92
249 )
250 ),
251 sigmak_
252 (
254 (
255 "sigmak",
256 this->runTime_.timeName(),
257 this->mesh_
258 ),
259 this->mesh_,
260 dimensioned<scalar>::getOrAddToDict
261 (
262 "sigmak",
263 this->coeffDict_,
264 1.0
265 )
266 ),
267 sigmaEps_
268 (
270 (
271 "sigmaEps",
272 this->runTime_.timeName(),
273 this->mesh_
274 ),
275 this->mesh_,
276 dimensioned<scalar>::getOrAddToDict
277 (
278 "sigmaEps",
279 this->coeffDict_,
280 1.3
281 )
282 ),
283
284 CdSigma_
285 (
287 (
288 "CdSigma",
289 this->runTime_.timeName(),
290 this->mesh_
291 ),
292 this->mesh_,
293 dimensionedScalar("CdSigma", dimless/dimLength, 0)
294 ),
295 betap_
296 (
298 (
299 "betap",
300 this->runTime_.timeName(),
301 this->mesh_
302 ),
303 this->mesh_,
304 dimensionedScalar("betap", dimless, 0)
305 ),
306 betad_
307 (
309 (
310 "betad",
311 this->runTime_.timeName(),
312 this->mesh_
313 ),
314 this->mesh_,
315 dimensionedScalar("betad", dimless, 0)
316 ),
317 C4_
318 (
320 (
321 "C4",
322 this->runTime_.timeName(),
323 this->mesh_
324 ),
325 this->mesh_,
326 dimensionedScalar("C4", dimless, 0)
327 ),
328 C5_
329 (
331 (
332 "C5",
333 this->runTime_.timeName(),
334 this->mesh_
335 ),
336 this->mesh_,
337 dimensionedScalar("C5", dimless, 0)
338 ),
339
340 k_
341 (
343 (
344 IOobject::groupName("k", alphaRhoPhi.group()),
345 this->runTime_.timeName(),
346 this->mesh_,
347 IOobject::MUST_READ,
348 IOobject::AUTO_WRITE
349 ),
350 this->mesh_
351 ),
352 epsilon_
353 (
355 (
356 IOobject::groupName("epsilon", alphaRhoPhi.group()),
357 this->runTime_.timeName(),
358 this->mesh_,
359 IOobject::MUST_READ,
360 IOobject::AUTO_WRITE
361 ),
362 this->mesh_
363 )
364{
365 bound(k_, this->kMin_);
366 bound(epsilon_, this->epsilonMin_);
367
368 if (type == typeName)
369 {
370 this->printCoeffs(type);
371 }
372
374}
375
376
377// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
378
379template<class BasicTurbulenceModel>
381{
383 {
384 return true;
385 }
386
387 return false;
388}
389
390
391template<class BasicTurbulenceModel>
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_;
406
408
410 (
411 fvc::div(fvc::absolute(this->phi(), U))().v()
412 );
413
416 (
417 this->GName(),
418 nut.v()*(dev(twoSymm(tgradU().v())) && tgradU().v())
419 );
420 tgradU.clear();
421
422 // Update epsilon and G at the wall
423 epsilon_.boundaryFieldRef().updateCoeffs();
424
426 volScalarField::Internal magU3(pow3(magU));
427
428 // Dissipation equation
430 (
431 fvm::ddt(alpha, rho, epsilon_)
432 + fvm::div(alphaRhoPhi, epsilon_)
433 - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
434 ==
435 C1_*alpha()*rho()*G*epsilon_()/k_()
436 - fvm::SuSp(((2.0/3.0)*C1_)*alpha()*rho()*divU, epsilon_)
437 - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_)
438 + epsilonSource(magU, magU3)
439 + fvOptions(alpha, rho, epsilon_)
440 );
441
442 epsEqn.ref().relax();
443 fvOptions.constrain(epsEqn.ref());
444 epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
445 solve(epsEqn);
446 fvOptions.correct(epsilon_);
447 bound(epsilon_, this->epsilonMin_);
448
449 // Turbulent kinetic energy equation
451 (
452 fvm::ddt(alpha, rho, k_)
453 + fvm::div(alphaRhoPhi, k_)
454 - fvm::laplacian(alpha*rho*DkEff(), k_)
455 ==
456 alpha()*rho()*G
457 - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
458 - fvm::Sp(alpha()*rho()*epsilon_()/k_(), k_)
459 + kSource(magU, magU3)
460 + fvOptions(alpha, rho, k_)
461 );
462
463 kEqn.ref().relax();
464 fvOptions.constrain(kEqn.ref());
465 solve(kEqn);
466 fvOptions.correct(k_);
467 bound(k_, this->kMin_);
468
469 correctNut();
470}
471
472
473// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
474
475} // End namespace RASModels
476} // End namespace Foam
477
478// ************************************************************************* //
Bound the given scalar field if it has gone unbounded.
fv::options & fvOptions
surfaceScalarField & phi
Graphite solid properties.
Definition: C.H:53
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
Templated abstract base class for RAS turbulence models.
Definition: RASModel.H:55
virtual tmp< fvScalarMatrix > epsilonSource() const
Variant of the standard k-epsilon turbulence model with additional source terms to handle the changes...
BasicTurbulenceModel::alphaField alphaField
BasicTurbulenceModel::rhoField rhoField
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
void setCdSigma(volScalarField::Internal &C, const porosityModels::powerLawLopesdaCosta &pm)
BasicTurbulenceModel::transportModel transportModel
void setPorosityCoefficient(volScalarField::Internal &C, const porosityModels::powerLawLopesdaCosta &pm)
virtual bool read()
Re-read model coefficients if they have changed.
virtual tmp< fvScalarMatrix > kSource() const
Add explicit source for turbulent kinetic energy.
Definition: kL.C:236
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
bool found(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Search for an entry (const access) with the given keyword.
Definition: dictionaryI.H:87
Generic dimensioned Type class.
Eddy viscosity turbulence model base class.
Definition: eddyViscosity.H:58
List of finite volume options.
Definition: faOptionList.H:70
Applies an explicit porosity source to the momentum equation within a specified region.
const porosityModel & model() const noexcept
Access to the porosityModel.
Finite-volume options.
Definition: fvOptions.H:59
const labelList & cellZoneIDs() const
Return const access to the cell zone IDs.
const dictionary & dict() const
Return dictionary used for model construction.
const scalarField & Sigma() 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:65
void clear() const noexcept
Definition: tmpI.H:287
T & ref() const
Definition: tmpI.H:227
A class for handling words, derived from Foam::string.
Definition: word.H:68
U
Definition: pEqn.H:72
thermo correct()
scalar nut
const cellShapeList & cells
zeroField divU
Definition: alphaSuSp.H:3
word timeName
Definition: getTimeIndex.H:3
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:54
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:49
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
Definition: fvcMeshPhi.C:190
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:48
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:48
zeroField Su(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:48
zeroField SuSp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
zeroField Sp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
Namespace for OpenFOAM.
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
const dimensionSet dimless
Dimensionless.
bool read(const char *buf, int32_t &val)
Same as readInt32.
Definition: int32.H:108
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
dimensionedScalar pow3(const dimensionedScalar &ds)
volScalarField & bound(volScalarField &, const dimensionedScalar &lowerBound)
Bound the given scalar field if it has gone unbounded.
Definition: bound.C:35
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:598
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
volScalarField & alpha
CEqn solve()
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333