Commit e38e5d67 authored by Kutalmis Bercin's avatar Kutalmis Bercin
Browse files

ENH: new RAS model: kEpsilonPhitF

    ENH: modify fWallFunction for kEpsilonPhitF model

    The k-epsilon-phit-f turbulence closure model for incompressible and
    compressible flows.

    The model is a three-transport-equation linear-eddy-viscosity turbulence
    closure model alongside an elliptic relaxation equation:
      - Turbulent kinetic energy, \c k,
      - Turbulent kinetic energy dissipation rate, \c epsilon,
      - Normalised wall-normal fluctuating velocity scale, \c phit,
      - Elliptic relaxation factor, \c f.

    Reference:
    \verbatim
        Standard model (Tag:LUU):
            Laurence, D. R., Uribe, J. C., & Utyuzhnikov, S. V. (2005).
            A robust formulation of the v2−f model.
            Flow, Turbulence and Combustion, 73(3-4), 169–185.
            DOI:10.1007/s10494-005-1974-8
    \endverbatim

    The default model coefficients are (LUU:Eqs. 19-20):
    \verbatim
        kEpsilonPhitFCoeffs
        {
            Cmu         0.22,    // Turbulent viscosity constant
            Ceps1a      1.4,     // Model constant for epsilon
            Ceps1b      1.0,     // Model constant for epsilon
            Ceps1c      0.05,    // Model constant for epsilon
            Ceps2       1.9,     // Model constant for epsilon
            Cf1         1.4,     // Model constant for f
            Cf2         0.3,     // Model constant for f
            CL          0.25,    // Model constant for L
            Ceta        110.0,   // Model constant for L
            CT          6.0,     // Model constant for T
            sigmaK      1.0,     // Turbulent Prandtl number for k
            sigmaEps    1.3,     // Turbulent Prandtl number for epsilon
            sigmaPhit   1.0,     // Turbulent Prandtl number for phit = sigmaK
        }
    \endverbatim

Note
    The name of the original variable replacing 'v2' is 'phi' (LUU:Eq. 14).
    However, the name 'phi' preexisted in OpenFOAM; therefore, this name was
    replaced by 'phit'
parent 427f9221
......@@ -94,6 +94,9 @@ makeRASModel(LRR);
#include "SSG.H"
makeRASModel(SSG);
#include "kEpsilonPhitF.H"
makeRASModel(kEpsilonPhitF);
// -------------------------------------------------------------------------- //
// LES models
......
......@@ -90,6 +90,9 @@ makeRASModel(LRR);
#include "SSG.H"
makeRASModel(SSG);
#include "kEpsilonPhitF.H"
makeRASModel(kEpsilonPhitF);
// -------------------------------------------------------------------------- //
// LES models
......
turbulenceModel.C
RAS/v2f/v2fBase.C
RAS/kEpsilonPhitF/kEpsilonPhitFBase.C
LESdelta = LES/LESdeltas
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "kEpsilonPhitF.H"
#include "fvOptions.H"
#include "bound.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace RASModels
{
// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
void kEpsilonPhitF<BasicTurbulenceModel>::correctNut()
{
// (LUU:p. 173)
this->nut_ = Cmu_*phit_*k_*T_;
this->nut_.correctBoundaryConditions();
fv::options::New(this->mesh_).correct(this->nut_);
BasicTurbulenceModel::correctNut();
}
template<class BasicTurbulenceModel>
volScalarField kEpsilonPhitF<BasicTurbulenceModel>::Ts() const
{
// (LUU:Eq. 7)
return
max
(
k_/epsilon_,
CT_*sqrt
(
max
(
this->nu(),
dimensionedScalar(this->nu()().dimensions(), Zero)
)/epsilon_
)
);
}
template<class BasicTurbulenceModel>
volScalarField kEpsilonPhitF<BasicTurbulenceModel>::Ls() const
{
// (LUU:Eq. 7)
return
CL_*max
(
pow(k_, 1.5)/epsilon_,
Ceta_*pow025
(
pow3
(
max
(
this->nu(),
dimensionedScalar(this->nu()().dimensions(), Zero)
)
)/epsilon_
)
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
kEpsilonPhitF<BasicTurbulenceModel>::kEpsilonPhitF
(
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
),
kEpsilonPhitFBase(),
Cmu_
(
dimensioned<scalar>::getOrAddToDict
(
"Cmu",
this->coeffDict_,
0.22
)
),
Ceps1a_
(
dimensioned<scalar>::getOrAddToDict
(
"Ceps1a",
this->coeffDict_,
1.4
)
),
Ceps1b_
(
dimensioned<scalar>::getOrAddToDict
(
"Ceps1b",
this->coeffDict_,
1.0
)
),
Ceps1c_
(
dimensioned<scalar>::getOrAddToDict
(
"Ceps1c",
this->coeffDict_,
0.05
)
),
Ceps2_
(
dimensioned<scalar>::getOrAddToDict
(
"Ceps2",
this->coeffDict_,
1.9
)
),
Cf1_
(
dimensioned<scalar>::getOrAddToDict
(
"Cf1",
this->coeffDict_,
1.4
)
),
Cf2_
(
dimensioned<scalar>::getOrAddToDict
(
"Cf2",
this->coeffDict_,
0.3
)
),
CL_
(
dimensioned<scalar>::getOrAddToDict
(
"CL",
this->coeffDict_,
0.25
)
),
Ceta_
(
dimensioned<scalar>::getOrAddToDict
(
"Ceta",
this->coeffDict_,
110.0
)
),
CT_
(
dimensioned<scalar>::getOrAddToDict
(
"CT",
this->coeffDict_,
6.0
)
),
sigmaK_
(
dimensioned<scalar>::getOrAddToDict
(
"sigmaK",
this->coeffDict_,
1.0
)
),
sigmaEps_
(
dimensioned<scalar>::getOrAddToDict
(
"sigmaEps",
this->coeffDict_,
1.3
)
),
sigmaPhit_
(
dimensioned<scalar>::getOrAddToDict
(
"sigmaPhit",
this->coeffDict_,
1.0
)
),
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_
),
phit_
(
IOobject
(
IOobject::groupName("phit", alphaRhoPhi.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
f_
(
IOobject
(
IOobject::groupName("f", alphaRhoPhi.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
T_
(
IOobject
(
"T",
this->runTime_.timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh_,
dimensionedScalar(dimTime, Zero)
),
phitMin_(dimensionedScalar("phitMin", phit_.dimensions(), SMALL)),
fMin_(dimensionedScalar("fMin", f_.dimensions(), SMALL)),
TMin_(dimensionedScalar("TMin", dimTime, SMALL)),
L2Min_(dimensionedScalar("L2Min", sqr(dimLength), SMALL))
{
bound(k_, this->kMin_);
bound(epsilon_, this->epsilonMin_);
bound(phit_, phitMin_);
bound(f_, fMin_);
if (type == typeName)
{
this->printCoeffs(type);
}
if
(
mag(sigmaK_.value()) < VSMALL
|| mag(sigmaEps_.value()) < VSMALL
|| mag(sigmaPhit_.value()) < VSMALL
)
{
FatalErrorInFunction
<< "Non-zero values are required for the model constants:" << nl
<< "sigmaK = " << sigmaK_ << nl
<< "sigmaEps = " << sigmaEps_ << nl
<< "sigmaPhit = " << sigmaPhit_ << nl
<< exit(FatalError);
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
bool kEpsilonPhitF<BasicTurbulenceModel>::read()
{
if (eddyViscosity<RASModel<BasicTurbulenceModel>>::read())
{
Cmu_.readIfPresent(this->coeffDict());
Ceps1a_.readIfPresent(this->coeffDict());
Ceps1b_.readIfPresent(this->coeffDict());
Ceps1c_.readIfPresent(this->coeffDict());
Ceps2_.readIfPresent(this->coeffDict());
Cf1_.readIfPresent(this->coeffDict());
Cf2_.readIfPresent(this->coeffDict());
CL_.readIfPresent(this->coeffDict());
Ceta_.readIfPresent(this->coeffDict());
CT_.readIfPresent(this->coeffDict());
sigmaK_.readIfPresent(this->coeffDict());
sigmaEps_.readIfPresent(this->coeffDict());
sigmaPhit_.readIfPresent(this->coeffDict());
return true;
}
return false;
}
template<class BasicTurbulenceModel>
void kEpsilonPhitF<BasicTurbulenceModel>::correct()
{
if (!this->turbulence_)
{
return;
}
// Construct local convenience 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();
const volScalarField::Internal divU
(
fvc::div(fvc::absolute(this->phi(), U))
);
tmp<volSymmTensorField> tS(symm(fvc::grad(U)));
volScalarField G(this->GName(), nut*(2.0*(dev(tS()) && tS())));
tS.clear();
T_ = Ts();
bound(T_, TMin_);
const volScalarField L2(type() + "L2", sqr(Ls()) + L2Min_);
const volScalarField::Internal Ceps1
(
"Ceps1",
Ceps1a_*(Ceps1b_ + Ceps1c_*sqrt(1.0/phit_()))
);
// Update epsilon and G at the wall
epsilon_.boundaryFieldRef().updateCoeffs();
// Turbulent kinetic energy dissipation rate equation (LUU:Eq. 4)
// k/T ~ epsilon
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(alpha, rho, epsilon_)
+ fvm::div(alphaRhoPhi, epsilon_)
- fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
==
alpha()*rho()*Ceps1*G()/T_()
- fvm::SuSp
(
(2.0/3.0*Ceps1)*(alpha()*rho()*divU),
epsilon_
)
- fvm::Sp(alpha()*rho()*Ceps2_/T_(), epsilon_)
+ 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 (LUU:Eq. 3)
// epsilon/k ~ 1/Ts
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()/T_(), k_)
+ fvOptions(alpha, rho, k_)
);
kEqn.ref().relax();
fvOptions.constrain(kEqn.ref());
solve(kEqn);
fvOptions.correct(k_);
bound(k_, this->kMin_);
// Elliptic relaxation function equation (LUU:Eq. 18)
// All source terms are non-negative functions (LUU:p. 176)
tmp<fvScalarMatrix> fEqn
(
- fvm::laplacian(f_)
==
- fvm::Sp(1.0/L2(), f_)
- (
(Cf1_ - 1.0)*(phit_() - 2.0/3.0)/T_()
-(Cf2_*G())/k_()
+(Cf2_*(2.0/3.0)*divU)
-(2.0*this->nu()*(fvc::grad(phit_) & fvc::grad(k_)))()/k_()
-(this->nu()*fvc::laplacian(phit_))()
)/L2()
);
fEqn.ref().relax();
solve(fEqn);
bound(f_, fMin_);
// Normalised wall-normal fluctuating velocity scale equation (LUU:Eq. 17)
// All source terms are non-negative functions (LUU:p. 176)
tmp<fvScalarMatrix> phitEqn
(
fvm::ddt(alpha, rho, phit_)
+ fvm::div(alphaRhoPhi, phit_)
- fvm::laplacian(alpha*rho*DphitEff(), phit_)
==
alpha()*rho()*f_()
- fvm::SuSp
(
alpha()*rho()*
(
G()/k_()
- (2.0/3.0)*divU
- (2.0*nut*(fvc::grad(phit_) & fvc::grad(k_)))()
/(k_()*sigmaPhit_*phit_())
)
, phit_
)
+ fvOptions(alpha, rho, phit_)
);
phitEqn.ref().relax();
fvOptions.constrain(phitEqn.ref());
solve(phitEqn);
fvOptions.correct(phit_);
bound(phit_, phitMin_);
correctNut();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace Foam
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
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/>.
Class
Foam::RASModels::kEpsilonPhitF
Group
grpRASTurbulence
Description
The k-epsilon-phit-f turbulence closure model for incompressible and
compressible flows.
The model is a three-transport-equation linear-eddy-viscosity turbulence
closure model alongside an elliptic relaxation equation:
- Turbulent kinetic energy, \c k,
- Turbulent kinetic energy dissipation rate, \c epsilon,
- Normalised wall-normal fluctuating velocity scale, \c phit,
- Elliptic relaxation factor, \c f.
Reference:
\verbatim
Standard model (Tag:LUU):
Laurence, D. R., Uribe, J. C., & Utyuzhnikov, S. V. (2005).
A robust formulation of the v2−f model.
Flow, Turbulence and Combustion, 73(3-4), 169–185.
DOI:10.1007/s10494-005-1974-8
\endverbatim