songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

ENH: Turbulence - refactored kOmegaSST models to enable easier addition of derived types

Browse Source
This commit is contained in:
Andrew Heather 2015-11-30 16:17:21 +00:00
parent 981b6303a4
commit 3e783accd2
5 changed files with 819 additions and 551 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

405
src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C
View File

@ -24,8 +24,6 @@ License
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "kOmegaSST.H" #include "kOmegaSST.H"
#include "bound.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -34,91 +32,22 @@ namespace Foam
namespace RASModels namespace RASModels
{ {
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template<class BasicTurbulenceModel>
tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F1
(
const volScalarField& CDkOmega
) const
{
tmp<volScalarField> CDkOmegaPlus = max
(
CDkOmega,
dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
);
tmp<volScalarField> arg1 = min
(
min
(
max
(
(scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_)
),
(4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
),
scalar(10)
);
return tanh(pow4(arg1));
}
template<class BasicTurbulenceModel>
tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F2() const
{
tmp<volScalarField> arg2 = min
(
max
(
(scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_)
),
scalar(100)
);
return tanh(sqr(arg2));
}
template<class BasicTurbulenceModel>
tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F3() const
{
tmp<volScalarField> arg3 = min
(
150*(this->mu()/this->rho_)/(omega_*sqr(y_)),
scalar(10)
);
return 1 - tanh(pow4(arg3));
}
template<class BasicTurbulenceModel>
tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F23() const
{
tmp<volScalarField> f23(F2());
if (F3_)
{
f23() *= F3();
}
return f23;
}
template<class BasicTurbulenceModel> template<class BasicTurbulenceModel>
void kOmegaSST<BasicTurbulenceModel>::correctNut(const volScalarField& S2) void kOmegaSST<BasicTurbulenceModel>::correctNut(const volScalarField& S2)
{ {
this->nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2)); // Correct the turbulence viscosity
this->nut_.correctBoundaryConditions(); kOmegaSSTBase<eddyViscosity<RASModel<BasicTurbulenceModel> > >::correctNut
(
S2
);
// Correct the turbulence thermal diffusivity
BasicTurbulenceModel::correctNut(); BasicTurbulenceModel::correctNut();
} }
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template<class BasicTurbulenceModel> template<class BasicTurbulenceModel>
void kOmegaSST<BasicTurbulenceModel>::correctNut() void kOmegaSST<BasicTurbulenceModel>::correctNut()
{ {
@ -126,53 +55,6 @@ void kOmegaSST<BasicTurbulenceModel>::correctNut()
} }
template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> kOmegaSST<BasicTurbulenceModel>::kSource() const
{
return tmp<fvScalarMatrix>
(
new fvScalarMatrix
(
k_,
dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
)
);
}
template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> kOmegaSST<BasicTurbulenceModel>::omegaSource() const
{
return tmp<fvScalarMatrix>
(
new fvScalarMatrix
(
omega_,
dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
)
);
}
template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> kOmegaSST<BasicTurbulenceModel>::Qsas
(
const volScalarField& S2,
const volScalarField& gamma,
const volScalarField& beta
) const
{
return tmp<fvScalarMatrix>
(
new fvScalarMatrix
(
omega_,
dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
)
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasicTurbulenceModel> template<class BasicTurbulenceModel>
@ -188,7 +70,7 @@ kOmegaSST<BasicTurbulenceModel>::kOmegaSST
const word& type const word& type
) )
: :
eddyViscosity<RASModel<BasicTurbulenceModel> > kOmegaSSTBase<eddyViscosity<RASModel<BasicTurbulenceModel> > >
( (
type, type,
alpha, alpha,
@ -198,286 +80,17 @@ kOmegaSST<BasicTurbulenceModel>::kOmegaSST
phi, phi,
transport, transport,
propertiesName propertiesName
),
alphaK1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaK1",
this->coeffDict_,
0.85
)
),
alphaK2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaK2",
this->coeffDict_,
1.0
)
),
alphaOmega1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaOmega1",
this->coeffDict_,
0.5
)
),
alphaOmega2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaOmega2",
this->coeffDict_,
0.856
)
),
gamma1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"gamma1",
this->coeffDict_,
5.0/9.0
)
),
gamma2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"gamma2",
this->coeffDict_,
0.44
)
),
beta1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"beta1",
this->coeffDict_,
0.075
)
),
beta2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"beta2",
this->coeffDict_,
0.0828
)
),
betaStar_
(
dimensioned<scalar>::lookupOrAddToDict
(
"betaStar",
this->coeffDict_,
0.09
)
),
a1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"a1",
this->coeffDict_,
0.31
)
),
b1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"b1",
this->coeffDict_,
1.0
)
),
c1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"c1",
this->coeffDict_,
10.0
)
),
F3_
(
Switch::lookupOrAddToDict
(
"F3",
this->coeffDict_,
false
)
),
y_(wallDist::New(this->mesh_).y()),
k_
(
IOobject
(
IOobject::groupName("k", U.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
omega_
(
IOobject
(
IOobject::groupName("omega", U.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
) )
{ {
bound(k_, this->kMin_); correctNut();
bound(omega_, this->omegaMin_);
if (type == typeName) if (type == typeName)
{ {
correctNut();
this->printCoeffs(type); this->printCoeffs(type);
} }
} }
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
bool kOmegaSST<BasicTurbulenceModel>::read()
{
if (eddyViscosity<RASModel<BasicTurbulenceModel> >::read())
{
alphaK1_.readIfPresent(this->coeffDict());
alphaK2_.readIfPresent(this->coeffDict());
alphaOmega1_.readIfPresent(this->coeffDict());
alphaOmega2_.readIfPresent(this->coeffDict());
gamma1_.readIfPresent(this->coeffDict());
gamma2_.readIfPresent(this->coeffDict());
beta1_.readIfPresent(this->coeffDict());
beta2_.readIfPresent(this->coeffDict());
betaStar_.readIfPresent(this->coeffDict());
a1_.readIfPresent(this->coeffDict());
b1_.readIfPresent(this->coeffDict());
c1_.readIfPresent(this->coeffDict());
F3_.readIfPresent("F3", this->coeffDict());
return true;
}
else
{
return false;
}
}
template<class BasicTurbulenceModel>
void kOmegaSST<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_;
eddyViscosity<RASModel<BasicTurbulenceModel> >::correct();
volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
tmp<volTensorField> tgradU = fvc::grad(U);
volScalarField S2(2*magSqr(symm(tgradU())));
volScalarField GbyNu((tgradU() && dev(twoSymm(tgradU()))));
volScalarField G(this->GName(), nut*GbyNu);
tgradU.clear();
// Update omega and G at the wall
omega_.boundaryField().updateCoeffs();
volScalarField CDkOmega
(
(2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
);
volScalarField F1(this->F1(CDkOmega));
{
volScalarField gamma(this->gamma(F1));
volScalarField beta(this->beta(F1));
// Turbulent frequency equation
tmp<fvScalarMatrix> omegaEqn
(
fvm::ddt(alpha, rho, omega_)
+ fvm::div(alphaRhoPhi, omega_)
- fvm::laplacian(alpha*rho*DomegaEff(F1), omega_)
==
alpha*rho*gamma
*min
(
GbyNu,
(c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2))
)
- fvm::SuSp((2.0/3.0)*alpha*rho*gamma*divU, omega_)
- fvm::Sp(alpha*rho*beta*omega_, omega_)
- fvm::SuSp
(
alpha*rho*(F1 - scalar(1))*CDkOmega/omega_,
omega_
)
+ Qsas(S2, gamma, beta)
+ omegaSource()
);
omegaEqn().relax();
omegaEqn().boundaryManipulate(omega_.boundaryField());
solve(omegaEqn);
bound(omega_, this->omegaMin_);
}
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(alpha, rho, k_)
+ fvm::div(alphaRhoPhi, k_)
- fvm::laplacian(alpha*rho*DkEff(F1), k_)
==
min(alpha*rho*G, (c1_*betaStar_)*alpha*rho*k_*omega_)
- fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
- fvm::Sp(alpha*rho*betaStar_*omega_, k_)
+ kSource()
);
kEqn().relax();
solve(kEqn);
bound(k_, this->kMin_);
correctNut(S2);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels } // End namespace RASModels

163
src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H
View File

@ -31,6 +31,8 @@ Description
Implementation of the k-omega-SST turbulence model for Implementation of the k-omega-SST turbulence model for
incompressible and compressible flows. incompressible and compressible flows.
Light wrapper around base class.
Turbulence model described in: Turbulence model described in:
\verbatim \verbatim
Menter, F. R. & Esch, T. (2001). Menter, F. R. & Esch, T. (2001).
@ -98,11 +100,15 @@ Description
SourceFiles SourceFiles
kOmegaSST.C kOmegaSST.C
SeeAlso
kOmegaSSTBase.H
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#ifndef kOmegaSST_H #ifndef kOmegaSST_H
#define kOmegaSST_H #define kOmegaSST_H
#include "kOmegaSSTBase.H"
#include "RASModel.H" #include "RASModel.H"
#include "eddyViscosity.H" #include "eddyViscosity.H"
@ -120,104 +126,15 @@ namespace RASModels
template<class BasicTurbulenceModel> template<class BasicTurbulenceModel>
class kOmegaSST class kOmegaSST
: :
public eddyViscosity<RASModel<BasicTurbulenceModel> > public kOmegaSSTBase<eddyViscosity<RASModel<BasicTurbulenceModel> > >
{ {
// Private Member Functions
// Disallow default bitwise copy construct and assignment
kOmegaSST(const kOmegaSST&);
kOmegaSST& operator=(const kOmegaSST&);
protected: protected:
// Protected data // Protecetd Member Functions
// Model coefficients
dimensionedScalar alphaK1_;
dimensionedScalar alphaK2_;
dimensionedScalar alphaOmega1_;
dimensionedScalar alphaOmega2_;
dimensionedScalar gamma1_;
dimensionedScalar gamma2_;
dimensionedScalar beta1_;
dimensionedScalar beta2_;
dimensionedScalar betaStar_;
dimensionedScalar a1_;
dimensionedScalar b1_;
dimensionedScalar c1_;
Switch F3_;
// Fields
//- Wall distance
// Note: different to wall distance in parent RASModel
// which is for near-wall cells only
const volScalarField& y_;
volScalarField k_;
volScalarField omega_;
// Private Member Functions
tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
tmp<volScalarField> F2() const;
tmp<volScalarField> F3() const;
tmp<volScalarField> F23() const;
tmp<volScalarField> blend
(
const volScalarField& F1,
const dimensionedScalar& psi1,
const dimensionedScalar& psi2
) const
{
return F1*(psi1 - psi2) + psi2;
}
tmp<volScalarField> alphaK(const volScalarField& F1) const
{
return blend(F1, alphaK1_, alphaK2_);
}
tmp<volScalarField> alphaOmega(const volScalarField& F1) const
{
return blend(F1, alphaOmega1_, alphaOmega2_);
}
tmp<volScalarField> beta(const volScalarField& F1) const
{
return blend(F1, beta1_, beta2_);
}
tmp<volScalarField> gamma(const volScalarField& F1) const
{
return blend(F1, gamma1_, gamma2_);
}
void correctNut(const volScalarField& S2);
// Protected Member Functions
virtual void correctNut(const volScalarField& S2);
virtual void correctNut(); virtual void correctNut();
virtual tmp<fvScalarMatrix> kSource() const;
virtual tmp<fvScalarMatrix> omegaSource() const;
virtual tmp<fvScalarMatrix> Qsas
(
const volScalarField& S2,
const volScalarField& gamma,
const volScalarField& beta
) const;
public: public:
@ -250,68 +167,6 @@ public:
//- Destructor //- Destructor
virtual ~kOmegaSST() virtual ~kOmegaSST()
{} {}
// Member Functions
//- Re-read model coefficients if they have changed
virtual bool read();
//- Return the effective diffusivity for k
tmp<volScalarField> DkEff(const volScalarField& F1) const
{
return tmp<volScalarField>
(
new volScalarField("DkEff", alphaK(F1)*this->nut_ + this->nu())
);
}
//- Return the effective diffusivity for omega
tmp<volScalarField> DomegaEff(const volScalarField& F1) const
{
return tmp<volScalarField>
(
new volScalarField
(
"DomegaEff",
alphaOmega(F1)*this->nut_ + this->nu()
)
);
}
//- Return the turbulence kinetic energy
virtual tmp<volScalarField> k() const
{
return k_;
}
//- Return the turbulence kinetic energy dissipation rate
virtual tmp<volScalarField> epsilon() const
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"epsilon",
this->mesh_.time().timeName(),
this->mesh_
),
betaStar_*k_*omega_,
omega_.boundaryField().types()
)
);
}
//- Return the turbulence kinetic energy dissipation rate
virtual tmp<volScalarField> omega() const
{
return omega_;
}
//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct();
}; };

480
src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.C Normal file
View File

@ -0,0 +1,480 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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 "kOmegaSSTBase.H"
#include "bound.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template<class BasicEddyViscosityModel>
tmp<volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::F1
(
const volScalarField& CDkOmega
) const
{
tmp<volScalarField> CDkOmegaPlus = max
(
CDkOmega,
dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
);
tmp<volScalarField> arg1 = min
(
min
(
max
(
(scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_)
),
(4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
),
scalar(10)
);
return tanh(pow4(arg1));
}
template<class BasicEddyViscosityModel>
tmp<volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::F2() const
{
tmp<volScalarField> arg2 = min
(
max
(
(scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_)
),
scalar(100)
);
return tanh(sqr(arg2));
}
template<class BasicEddyViscosityModel>
tmp<volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::F3() const
{
tmp<volScalarField> arg3 = min
(
150*(this->mu()/this->rho_)/(omega_*sqr(y_)),
scalar(10)
);
return 1 - tanh(pow4(arg3));
}
template<class BasicEddyViscosityModel>
tmp<volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::F23() const
{
tmp<volScalarField> f23(F2());
if (F3_)
{
f23() *= F3();
}
return f23;
}
template<class BasicEddyViscosityModel>
void kOmegaSSTBase<BasicEddyViscosityModel>::correctNut
(
const volScalarField& S2
)
{
// Correct the turbulence viscosity
this->nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
this->nut_.correctBoundaryConditions();
}
template<class BasicEddyViscosityModel>
void kOmegaSSTBase<BasicEddyViscosityModel>::correctNut()
{
correctNut(2*magSqr(symm(fvc::grad(this->U_))));
}
template<class BasicEddyViscosityModel>
tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::kSource() const
{
return tmp<fvScalarMatrix>
(
new fvScalarMatrix
(
k_,
dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
)
);
}
template<class BasicEddyViscosityModel>
tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::omegaSource() const
{
return tmp<fvScalarMatrix>
(
new fvScalarMatrix
(
omega_,
dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
)
);
}
template<class BasicEddyViscosityModel>
tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::Qsas
(
const volScalarField& S2,
const volScalarField& gamma,
const volScalarField& beta
) const
{
return tmp<fvScalarMatrix>
(
new fvScalarMatrix
(
omega_,
dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
)
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasicEddyViscosityModel>
kOmegaSSTBase<BasicEddyViscosityModel>::kOmegaSSTBase
(
const word& type,
const alphaField& alpha,
const rhoField& rho,
const volVectorField& U,
const surfaceScalarField& alphaRhoPhi,
const surfaceScalarField& phi,
const transportModel& transport,
const word& propertiesName
)
:
BasicEddyViscosityModel
(
type,
alpha,
rho,
U,
alphaRhoPhi,
phi,
transport,
propertiesName
),
alphaK1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaK1",
this->coeffDict_,
0.85
)
),
alphaK2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaK2",
this->coeffDict_,
1.0
)
),
alphaOmega1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaOmega1",
this->coeffDict_,
0.5
)
),
alphaOmega2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaOmega2",
this->coeffDict_,
0.856
)
),
gamma1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"gamma1",
this->coeffDict_,
5.0/9.0
)
),
gamma2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"gamma2",
this->coeffDict_,
0.44
)
),
beta1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"beta1",
this->coeffDict_,
0.075
)
),
beta2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"beta2",
this->coeffDict_,
0.0828
)
),
betaStar_
(
dimensioned<scalar>::lookupOrAddToDict
(
"betaStar",
this->coeffDict_,
0.09
)
),
a1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"a1",
this->coeffDict_,
0.31
)
),
b1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"b1",
this->coeffDict_,
1.0
)
),
c1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"c1",
this->coeffDict_,
10.0
)
),
F3_
(
Switch::lookupOrAddToDict
(
"F3",
this->coeffDict_,
false
)
),
y_(wallDist::New(this->mesh_).y()),
k_
(
IOobject
(
IOobject::groupName("k", U.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
omega_
(
IOobject
(
IOobject::groupName("omega", U.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
)
{
bound(k_, this->kMin_);
bound(omega_, this->omegaMin_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasicEddyViscosityModel>
bool kOmegaSSTBase<BasicEddyViscosityModel>::read()
{
if (BasicEddyViscosityModel::read())
{
alphaK1_.readIfPresent(this->coeffDict());
alphaK2_.readIfPresent(this->coeffDict());
alphaOmega1_.readIfPresent(this->coeffDict());
alphaOmega2_.readIfPresent(this->coeffDict());
gamma1_.readIfPresent(this->coeffDict());
gamma2_.readIfPresent(this->coeffDict());
beta1_.readIfPresent(this->coeffDict());
beta2_.readIfPresent(this->coeffDict());
betaStar_.readIfPresent(this->coeffDict());
a1_.readIfPresent(this->coeffDict());
b1_.readIfPresent(this->coeffDict());
c1_.readIfPresent(this->coeffDict());
F3_.readIfPresent("F3", this->coeffDict());
return true;
}
else
{
return false;
}
}
template<class BasicEddyViscosityModel>
void kOmegaSSTBase<BasicEddyViscosityModel>::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_;
BasicEddyViscosityModel::correct();
volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
tmp<volTensorField> tgradU = fvc::grad(U);
volScalarField S2(2*magSqr(symm(tgradU())));
volScalarField GbyNu((tgradU() && dev(twoSymm(tgradU()))));
volScalarField G(this->GName(), nut*GbyNu);
tgradU.clear();
// Update omega and G at the wall
omega_.boundaryField().updateCoeffs();
volScalarField CDkOmega
(
(2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
);
volScalarField F1(this->F1(CDkOmega));
{
volScalarField gamma(this->gamma(F1));
volScalarField beta(this->beta(F1));
// Turbulent frequency equation
tmp<fvScalarMatrix> omegaEqn
(
fvm::ddt(alpha, rho, omega_)
+ fvm::div(alphaRhoPhi, omega_)
- fvm::laplacian(alpha*rho*DomegaEff(F1), omega_)
==
alpha*rho*gamma
*min
(
GbyNu,
(c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2))
)
- fvm::SuSp((2.0/3.0)*alpha*rho*gamma*divU, omega_)
- fvm::Sp(alpha*rho*beta*omega_, omega_)
- fvm::SuSp
(
alpha*rho*(F1 - scalar(1))*CDkOmega/omega_,
omega_
)
+ Qsas(S2, gamma, beta)
+ omegaSource()
);
omegaEqn().relax();
omegaEqn().boundaryManipulate(omega_.boundaryField());
solve(omegaEqn);
bound(omega_, this->omegaMin_);
}
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(alpha, rho, k_)
+ fvm::div(alphaRhoPhi, k_)
- fvm::laplacian(alpha*rho*DkEff(F1), k_)
==
min(alpha*rho*G, (c1_*betaStar_)*alpha*rho*k_*omega_)
- fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
- fvm::Sp(alpha*rho*betaStar_*omega_, k_)
+ kSource()
);
kEqn().relax();
solve(kEqn);
bound(k_, this->kMin_);
correctNut(S2);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

320
src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.H Normal file
View File

@ -0,0 +1,320 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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/>.
Class
Foam::kOmegaSSTBase
Description
Base class implementation of the k-omega-SST turbulence model for
incompressible and compressible flows.
Turbulence model described in:
\verbatim
Menter, F. R. & Esch, T. (2001).
Elements of Industrial Heat Transfer Prediction.
16th Brazilian Congress of Mechanical Engineering (COBEM).
\endverbatim
with updated coefficients from
\verbatim
Menter, F. R., Kuntz, M., and Langtry, R. (2003).
Ten Years of Industrial Experience with the SST Turbulence Model.
Turbulence, Heat and Mass Transfer 4, ed: K. Hanjalic, Y. Nagano,
& M. Tummers, Begell House, Inc., 625 - 632.
\endverbatim
but with the consistent production terms from the 2001 paper as form in the
2003 paper is a typo, see
\verbatim
http://turbmodels.larc.nasa.gov/sst.html
\endverbatim
and the addition of the optional F3 term for rough walls from
\verbatim
Hellsten, A. (1998).
"Some Improvements in Menter’s k-omega-SST turbulence model"
29th AIAA Fluid Dynamics Conference, AIAA-98-2554.
\endverbatim
Note that this implementation is written in terms of alpha diffusion
coefficients rather than the more traditional sigma (alpha = 1/sigma) so
that the blending can be applied to all coefficuients in a consistent
manner. The paper suggests that sigma is blended but this would not be
consistent with the blending of the k-epsilon and k-omega models.
Also note that the error in the last term of equation (2) relating to
sigma has been corrected.
Wall-functions are applied in this implementation by using equations (14)
to specify the near-wall omega as appropriate.
The blending functions (15) and (16) are not currently used because of the
uncertainty in their origin, range of applicability and that if y+ becomes
sufficiently small blending u_tau in this manner clearly becomes nonsense.
The default model coefficients are
\verbatim
kOmegaSSTBaseCoeffs
{
alphaK1 0.85;
alphaK2 1.0;
alphaOmega1 0.5;
alphaOmega2 0.856;
beta1 0.075;
beta2 0.0828;
betaStar 0.09;
gamma1 5/9;
gamma2 0.44;
a1 0.31;
b1 1.0;
c1 10.0;
F3 no;
}
\endverbatim
SourceFiles
kOmegaSSTBase.C
\*---------------------------------------------------------------------------*/
#ifndef kOmegaSSTBase_H
#define kOmegaSSTBase_H
#include "RASModel.H"
#include "eddyViscosity.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class kOmegaSSTBase Declaration
\*---------------------------------------------------------------------------*/
template<class BasicEddyViscosityModel>
class kOmegaSSTBase
:
public BasicEddyViscosityModel
{
// Private Member Functions
// Disallow default bitwise copy construct and assignment
kOmegaSSTBase(const kOmegaSSTBase&);
kOmegaSSTBase& operator=(const kOmegaSSTBase&);
protected:
// Protected data
// Model coefficients
dimensionedScalar alphaK1_;
dimensionedScalar alphaK2_;
dimensionedScalar alphaOmega1_;
dimensionedScalar alphaOmega2_;
dimensionedScalar gamma1_;
dimensionedScalar gamma2_;
dimensionedScalar beta1_;
dimensionedScalar beta2_;
dimensionedScalar betaStar_;
dimensionedScalar a1_;
dimensionedScalar b1_;
dimensionedScalar c1_;
Switch F3_;
// Fields
//- Wall distance
// Note: different to wall distance in parent RASModel
// which is for near-wall cells only
const volScalarField& y_;
volScalarField k_;
volScalarField omega_;
// Protected Member Functions
tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
tmp<volScalarField> F2() const;
tmp<volScalarField> F3() const;
tmp<volScalarField> F23() const;
tmp<volScalarField> blend
(
const volScalarField& F1,
const dimensionedScalar& psi1,
const dimensionedScalar& psi2
) const
{
return F1*(psi1 - psi2) + psi2;
}
tmp<volScalarField> alphaK(const volScalarField& F1) const
{
return blend(F1, alphaK1_, alphaK2_);
}
tmp<volScalarField> alphaOmega(const volScalarField& F1) const
{
return blend(F1, alphaOmega1_, alphaOmega2_);
}
tmp<volScalarField> beta(const volScalarField& F1) const
{
return blend(F1, beta1_, beta2_);
}
tmp<volScalarField> gamma(const volScalarField& F1) const
{
return blend(F1, gamma1_, gamma2_);
}
void correctNut(const volScalarField& S2);
virtual void correctNut();
virtual tmp<fvScalarMatrix> kSource() const;
virtual tmp<fvScalarMatrix> omegaSource() const;
virtual tmp<fvScalarMatrix> Qsas
(
const volScalarField& S2,
const volScalarField& gamma,
const volScalarField& beta
) const;
public:
typedef typename BasicEddyViscosityModel::alphaField alphaField;
typedef typename BasicEddyViscosityModel::rhoField rhoField;
typedef typename BasicEddyViscosityModel::transportModel transportModel;
// Constructors
//- Construct from components
kOmegaSSTBase
(
const word& type,
const alphaField& alpha,
const rhoField& rho,
const volVectorField& U,
const surfaceScalarField& alphaRhoPhi,
const surfaceScalarField& phi,
const transportModel& transport,
const word& propertiesName = turbulenceModel::propertiesName
);
//- Destructor
virtual ~kOmegaSSTBase()
{}
// Member Functions
//- Re-read model coefficients if they have changed
virtual bool read();
//- Return the effective diffusivity for k
tmp<volScalarField> DkEff(const volScalarField& F1) const
{
return tmp<volScalarField>
(
new volScalarField("DkEff", alphaK(F1)*this->nut_ + this->nu())
);
}
//- Return the effective diffusivity for omega
tmp<volScalarField> DomegaEff(const volScalarField& F1) const
{
return tmp<volScalarField>
(
new volScalarField
(
"DomegaEff",
alphaOmega(F1)*this->nut_ + this->nu()
)
);
}
//- Return the turbulence kinetic energy
virtual tmp<volScalarField> k() const
{
return k_;
}
//- Return the turbulence kinetic energy dissipation rate
virtual tmp<volScalarField> epsilon() const
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"epsilon",
this->mesh_.time().timeName(),
this->mesh_
),
betaStar_*k_*omega_,
omega_.boundaryField().types()
)
);
}
//- Return the turbulence kinetic energy dissipation rate
virtual tmp<volScalarField> omega() const
{
return omega_;
}
//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "kOmegaSSTBase.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

0
src/TurbulenceModels/turbulenceModels/turbulenceModelDoc.H → src/TurbulenceModels/turbulenceModels/doc/turbulenceModelDoc.H
View File