Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
335 lines (283 sloc) 8.09 KB
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "realizableKE.H"
#include "fvOptions.H"
#include "bound.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace RASModels
{
// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
tmp<volScalarField> realizableKE<BasicTurbulenceModel>::rCmu
(
const volTensorField& gradU,
const volScalarField& S2,
const volScalarField& magS
)
{
tmp<volSymmTensorField> tS = dev(symm(gradU));
const volSymmTensorField& S = tS();
volScalarField W
(
(2*sqrt(2.0))*((S&S)&&S)
/(
magS*S2
+ dimensionedScalar("small", dimensionSet(0, 0, -3, 0, 0), small)
)
);
tS.clear();
volScalarField phis
(
(1.0/3.0)*acos(min(max(sqrt(6.0)*W, -scalar(1)), scalar(1)))
);
volScalarField As(sqrt(6.0)*cos(phis));
volScalarField Us(sqrt(S2/2.0 + magSqr(skew(gradU))));
return 1.0/(A0_ + As*Us*k_/epsilon_);
}
template<class BasicTurbulenceModel>
void realizableKE<BasicTurbulenceModel>::correctNut
(
const volTensorField& gradU,
const volScalarField& S2,
const volScalarField& magS
)
{
this->nut_ = rCmu(gradU, S2, magS)*sqr(k_)/epsilon_;
this->nut_.correctBoundaryConditions();
fv::options::New(this->mesh_).correct(this->nut_);
BasicTurbulenceModel::correctNut();
}
template<class BasicTurbulenceModel>
void realizableKE<BasicTurbulenceModel>::correctNut()
{
tmp<volTensorField> tgradU = fvc::grad(this->U_);
volScalarField S2(2*magSqr(dev(symm(tgradU()))));
volScalarField magS(sqrt(S2));
correctNut(tgradU(), S2, magS);
}
template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> realizableKE<BasicTurbulenceModel>::kSource() const
{
return tmp<fvScalarMatrix>
(
new fvScalarMatrix
(
k_,
dimVolume*this->rho_.dimensions()*k_.dimensions()
/dimTime
)
);
}
template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> realizableKE<BasicTurbulenceModel>::epsilonSource() const
{
return tmp<fvScalarMatrix>
(
new fvScalarMatrix
(
epsilon_,
dimVolume*this->rho_.dimensions()*epsilon_.dimensions()
/dimTime
)
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
realizableKE<BasicTurbulenceModel>::realizableKE
(
const alphaField& alpha,
const rhoField& rho,
const volVectorField& U,
const surfaceScalarField& alphaRhoPhi,
const surfaceScalarField& phi,
const transportModel& transport,
const word& propertiesName,
const word& type
)
:
eddyViscosity<RASModel<BasicTurbulenceModel>>
(
type,
alpha,
rho,
U,
alphaRhoPhi,
phi,
transport,
propertiesName
),
A0_
(
dimensioned<scalar>::lookupOrAddToDict
(
"A0",
this->coeffDict_,
4.0
)
),
C2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C2",
this->coeffDict_,
1.9
)
),
sigmak_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmak",
this->coeffDict_,
1.0
)
),
sigmaEps_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmaEps",
this->coeffDict_,
1.2
)
),
k_
(
IOobject
(
"k",
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
epsilon_
(
IOobject
(
"epsilon",
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
)
{
bound(k_, this->kMin_);
bound(epsilon_, this->epsilonMin_);
if (type == typeName)
{
this->printCoeffs(type);
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
bool realizableKE<BasicTurbulenceModel>::read()
{
if (eddyViscosity<RASModel<BasicTurbulenceModel>>::read())
{
A0_.readIfPresent(this->coeffDict());
C2_.readIfPresent(this->coeffDict());
sigmak_.readIfPresent(this->coeffDict());
sigmaEps_.readIfPresent(this->coeffDict());
return true;
}
else
{
return false;
}
}
template<class BasicTurbulenceModel>
void realizableKE<BasicTurbulenceModel>::correct()
{
if (!this->turbulence_)
{
return;
}
// Local references
const alphaField& alpha = this->alpha_;
const rhoField& rho = this->rho_;
const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
const volVectorField& U = this->U_;
volScalarField& nut = this->nut_;
fv::options& fvOptions(fv::options::New(this->mesh_));
eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
tmp<volTensorField> tgradU = fvc::grad(U);
volScalarField S2(2*magSqr(dev(symm(tgradU()))));
volScalarField magS(sqrt(S2));
volScalarField eta(magS*k_/epsilon_);
volScalarField C1(max(eta/(scalar(5) + eta), scalar(0.43)));
volScalarField G(this->GName(), nut*(tgradU() && dev(twoSymm(tgradU()))));
// Update epsilon and G at the wall
epsilon_.boundaryFieldRef().updateCoeffs();
// Dissipation equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(alpha, rho, epsilon_)
+ fvm::div(alphaRhoPhi, epsilon_)
- fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
==
C1*alpha*rho*magS*epsilon_
- fvm::Sp
(
C2_*alpha*rho*epsilon_/(k_ + sqrt(this->nu()*epsilon_)),
epsilon_
)
+ epsilonSource()
+ fvOptions(alpha, rho, epsilon_)
);
epsEqn.ref().relax();
fvOptions.constrain(epsEqn.ref());
epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
solve(epsEqn);
fvOptions.correct(epsilon_);
bound(epsilon_, this->epsilonMin_);
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(alpha, rho, k_)
+ fvm::div(alphaRhoPhi, k_)
- fvm::laplacian(alpha*rho*DkEff(), k_)
==
alpha*rho*G
- fvm::SuSp(2.0/3.0*alpha*rho*divU, k_)
- fvm::Sp(alpha*rho*epsilon_/k_, k_)
+ kSource()
+ fvOptions(alpha, rho, k_)
);
kEqn.ref().relax();
fvOptions.constrain(kEqn.ref());
solve(kEqn);
fvOptions.correct(k_);
bound(k_, this->kMin_);
correctNut(tgradU(), S2, magS);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace Foam
// ************************************************************************* //
You can’t perform that action at this time.