| muycfd |
June 28, 2021 15:17 |
Underprediction Issues (Turbulence Model Implementation)
1 Attachment(s)
Dear Foamers,
I am a fairly beginner user, trying to implement a RANS turbulence model (k-e-t) that I found in the literature, which has multiple time scales. As far as I know, this model and its coefficients are tuned for plane jet case (for shear-free flow). It is an incompressible flow case. Up to this point, I have tried many things but couldn’t match it neither with the experimental results nor with published simulation results of it. My simulations show slight underprediction in turbulence quantities (i.e. k/Uc^2, uv/Uc^2). However, velocity distribution is almost the same. Any suggestion or corrections would be invaluable. Model equations is attached as an image. And the code can be found below:
Code:
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 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 "kET.H"
#include "fvOptions.H"
#include "bound.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace RASModels
{
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template<class BasicMomentumTransportModel>
void kET<BasicMomentumTransportModel>::correctNut()
{
this->nut_ = Cmu_*sqr(k_)/epsilon_;
this->nut_.correctBoundaryConditions();
fv::options::New(this->mesh_).correct(this->nut_);
}
template<class BasicMomentumTransportModel>
tmp<fvScalarMatrix> kET<BasicMomentumTransportModel>::kSource() const
{
return tmp<fvScalarMatrix>
(
new fvScalarMatrix
(
k_,
dimVolume*this->rho_.dimensions()*k_.dimensions()
/dimTime
)
);
}
template<class BasicMomentumTransportModel>
tmp<fvScalarMatrix> kET<BasicMomentumTransportModel>::epsilonSource() const
{
return tmp<fvScalarMatrix>
(
new fvScalarMatrix
(
epsilon_,
dimVolume*this->rho_.dimensions()*epsilon_.dimensions()
/dimTime
)
);
}
template<class BasicMomentumTransportModel>
tmp<fvScalarMatrix> kET<BasicMomentumTransportModel>::taoSource() const
{
return tmp<fvScalarMatrix>
(
new fvScalarMatrix
(
tao_,
dimVolume*this->rho_.dimensions()*tao_.dimensions()
/dimTime
)
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasicMomentumTransportModel>
kET<BasicMomentumTransportModel>::kET
(
const alphaField& alpha,
const rhoField& rho,
const volVectorField& U,
const surfaceScalarField& alphaRhoPhi,
const surfaceScalarField& phi,
const transportModel& transport,
const word& type
)
:
eddyViscosity<RASModel<BasicMomentumTransportModel>>
(
type,
alpha,
rho,
U,
alphaRhoPhi,
phi,
transport
),
Cmu_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cmu",
this->coeffDict_,
0.09
)
),
Ceps1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ceps1",
this->coeffDict_,
0.75
)
),
Ceps2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ceps2",
this->coeffDict_,
1.0
)
),
Ceps3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ceps3",
this->coeffDict_,
0.67
)
),
Ctao0_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ctao0",
this->coeffDict_,
1.054
)
),
Ctao1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ctao1",
this->coeffDict_,
1.1
)
),
Ctao2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ctao2",
this->coeffDict_,
0.59
)
),
Ctao3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ctao3",
this->coeffDict_,
0.83
)
),
sigmak_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmak",
this->coeffDict_,
1.0
)
),
sigmaEps_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmaEps",
this->coeffDict_,
1.2
)
),
sigmaTao_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmaTao",
this->coeffDict_,
1.1
)
),
k_
(
IOobject
(
IOobject::groupName("k", alphaRhoPhi.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
epsilon_
(
IOobject
(
IOobject::groupName("epsilon", alphaRhoPhi.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
tao_
(
IOobject
(
IOobject::groupName("tao", alphaRhoPhi.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
)
{
bound(k_, this->kMin_);
bound(epsilon_, this->epsilonMin_);
bound(tao_, this->taoMin_);
if (type == typeName)
{
this->printCoeffs(type);
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasicMomentumTransportModel>
bool kET<BasicMomentumTransportModel>::read()
{
if (eddyViscosity<RASModel<BasicMomentumTransportModel>>::read())
{
Cmu_.readIfPresent(this->coeffDict());
Ceps1_.readIfPresent(this->coeffDict());
Ceps2_.readIfPresent(this->coeffDict());
Ceps3_.readIfPresent(this->coeffDict());
Ctao0_.readIfPresent(this->coeffDict());
Ctao1_.readIfPresent(this->coeffDict());
Ctao2_.readIfPresent(this->coeffDict());
Ctao3_.readIfPresent(this->coeffDict());
sigmak_.readIfPresent(this->coeffDict());
sigmaEps_.readIfPresent(this->coeffDict());
sigmaTao_.readIfPresent(this->coeffDict());
return true;
}
else
{
return false;
}
}
template<class BasicMomentumTransportModel>
void kET<BasicMomentumTransportModel>::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<BasicMomentumTransportModel>>::correct();
volScalarField::Internal divU
(
fvc::div(fvc::absolute(this->phi(), U))()
);
tmp<volTensorField> tgradU = fvc::grad(U);
volScalarField Str(2*magSqr(dev(symm(tgradU()))));
volScalarField magStr(sqrt(Str));
volScalarField::Internal G
(
this->GName(),
nut()*(dev(twoSymm(tgradU().v())) && tgradU().v())
);
tgradU.clear();
// 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_)
==
Ceps1_*alpha()*rho()*G/tao_()
- fvm::SuSp(((2.0/3.0)*Ceps1_)*alpha()*rho()*k_()*divU/tao_()/epsilon_(), epsilon_)
- fvm::Sp(Ceps2_*alpha()*rho()/tao_(), epsilon_)
- fvm::SuSp(Ceps3_*alpha()*rho()*magStr(), 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_);
// Tao equation
tmp<fvScalarMatrix> taoEqn
(
fvm::ddt(alpha, rho, tao_)
+ fvm::div(alphaRhoPhi, tao_)
- fvm::laplacian(alpha*rho*DtaoEff(), tao_)
==
Ctao0_*alpha()*rho()
- fvm::Sp(Ctao1_*alpha()*rho()*epsilon_()/k_(),tao_)
- fvm::Sp(Ctao2_*alpha()*rho()*G/k_(),tao_)
+ fvm::SuSp(((2.0/3.0)*Ctao2_)*alpha()*rho()*divU, tao_)
- fvm::SuSp(Ctao3_*alpha()*rho()*magStr(), tao_)
+ taoSource()
+ fvOptions(alpha, rho, tao_)
);
taoEqn.ref().relax();
fvOptions.constrain(taoEqn.ref());
solve(taoEqn);
fvOptions.correct(tao_);
bound(tao_, this->taoMin_);
correctNut();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace Foam
// ************************************************************************* //
Thanks,
P.S.: I just used tao instead of tau to prevent any unforeseen match of variables. |
