Permalink
Find file
Fetching contributors…
Cannot retrieve contributors at this time
375 lines (320 sloc) 8.31 KB
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 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 "RNGkEpsilon.H"
#include "addToRunTimeSelectionTable.H"
#include "backwardsCompatibilityWallFunctions.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace compressible
{
namespace RASModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(RNGkEpsilon, 0);
addToRunTimeSelectionTable(RASModel, RNGkEpsilon, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
RNGkEpsilon::RNGkEpsilon
(
const volScalarField& rho,
const volVectorField& U,
const surfaceScalarField& phi,
const basicThermo& thermophysicalModel,
const word& turbulenceModelName,
const word& modelName
)
:
RASModel(modelName, rho, U, phi, thermophysicalModel, turbulenceModelName),
Cmu_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cmu",
coeffDict_,
0.0845
)
),
C1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C1",
coeffDict_,
1.42
)
),
C2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C2",
coeffDict_,
1.68
)
),
C3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C3",
coeffDict_,
-0.33
)
),
sigmak_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmak",
coeffDict_,
0.71942
)
),
sigmaEps_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmaEps",
coeffDict_,
0.71942
)
),
Prt_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Prt",
coeffDict_,
1.0
)
),
eta0_
(
dimensioned<scalar>::lookupOrAddToDict
(
"eta0",
coeffDict_,
4.38
)
),
beta_
(
dimensioned<scalar>::lookupOrAddToDict
(
"beta",
coeffDict_,
0.012
)
),
k_
(
IOobject
(
"k",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateK("k", mesh_)
),
epsilon_
(
IOobject
(
"epsilon",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateEpsilon("epsilon", mesh_)
),
mut_
(
IOobject
(
"mut",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateMut("mut", mesh_)
),
alphat_
(
IOobject
(
"alphat",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateAlphat("alphat", mesh_)
)
{
bound(k_, kMin_);
bound(epsilon_, epsilonMin_);
mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
mut_.correctBoundaryConditions();
alphat_ = mut_/Prt_;
alphat_.correctBoundaryConditions();
printCoeffs();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
tmp<volSymmTensorField> RNGkEpsilon::R() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"R",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
k_.boundaryField().types()
)
);
}
tmp<volSymmTensorField> RNGkEpsilon::devRhoReff() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"devRhoReff",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
-muEff()*dev(twoSymm(fvc::grad(U_)))
)
);
}
tmp<fvVectorMatrix> RNGkEpsilon::divDevRhoReff(volVectorField& U) const
{
return
(
- fvm::laplacian(muEff(), U) - fvc::div(muEff()*dev2(T(fvc::grad(U))))
);
}
bool RNGkEpsilon::read()
{
if (RASModel::read())
{
Cmu_.readIfPresent(coeffDict());
C1_.readIfPresent(coeffDict());
C2_.readIfPresent(coeffDict());
C3_.readIfPresent(coeffDict());
Prt_.readIfPresent(coeffDict());
sigmaEps_.readIfPresent(coeffDict());
Prt_.readIfPresent(coeffDict());
eta0_.readIfPresent(coeffDict());
beta_.readIfPresent(coeffDict());
return true;
}
else
{
return false;
}
}
void RNGkEpsilon::correct()
{
if (!turbulence_)
{
// Re-calculate viscosity
mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
mut_.correctBoundaryConditions();
// Re-calculate thermal diffusivity
alphat_ = mut_/Prt_;
alphat_.correctBoundaryConditions();
return;
}
RASModel::correct();
volScalarField divU(fvc::div(phi_/fvc::interpolate(rho_)));
if (mesh_.moving())
{
divU += fvc::div(mesh_.phi());
}
tmp<volTensorField> tgradU = fvc::grad(U_);
volScalarField S2((tgradU() && dev(twoSymm(tgradU()))));
tgradU.clear();
volScalarField G("RASModel::G", mut_*S2);
volScalarField eta(sqrt(mag(S2))*k_/epsilon_);
volScalarField eta3(eta*sqr(eta));
volScalarField R
(
((eta*(-eta/eta0_ + scalar(1)))/(beta_*eta3 + scalar(1)))
);
// Update epsilon and G at the wall
epsilon_.boundaryField().updateCoeffs();
// Dissipation equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(rho_, epsilon_)
+ fvm::div(phi_, epsilon_)
- fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
- fvm::laplacian(DepsilonEff(), epsilon_)
==
(C1_ - R)*G*epsilon_/k_
- fvm::SuSp(((2.0/3.0)*C1_ + C3_)*rho_*divU, epsilon_)
- fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
);
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
solve(epsEqn);
bound(epsilon_, epsilonMin_);
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(rho_, k_)
+ fvm::div(phi_, k_)
- fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
- fvm::laplacian(DkEff(), k_)
==
G - fvm::SuSp(2.0/3.0*rho_*divU, k_)
- fvm::Sp(rho_*epsilon_/k_, k_)
);
kEqn().relax();
solve(kEqn);
bound(k_, kMin_);
// Re-calculate viscosity
mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
mut_.correctBoundaryConditions();
// Re-calculate thermal diffusivity
alphat_ = mut_/Prt_;
alphat_.correctBoundaryConditions();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace compressible
} // End namespace Foam
// ************************************************************************* //