/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 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 .
\*---------------------------------------------------------------------------*/
#include "kOmegaSSTBase.H"
#include "fvOptions.H"
#include "bound.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template
tmp
kOmegaSST::kOmegaSST::F1
(
const volScalarField& CDkOmega
) const
{
tmp CDkOmegaPlus = max
(
CDkOmega,
dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
);
tmp 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
tmp
kOmegaSST::kOmegaSST::F2() const
{
tmp 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
tmp
kOmegaSST::kOmegaSST::F3() const
{
tmp arg3 = min
(
150*(this->mu()/this->rho_)/(omega_*sqr(y_)),
scalar(10)
);
return 1 - tanh(pow4(arg3));
}
template
tmp
kOmegaSST::kOmegaSST::F23() const
{
tmp f23(F2());
if (F3_)
{
f23.ref() *= F3();
}
return f23;
}
template
void kOmegaSST::correctNut
(
const volScalarField& S2,
const volScalarField& F2
)
{
this->nut_ = a1_*k_/max(a1_*omega_, b1_*F2*sqrt(S2));
this->nut_.correctBoundaryConditions();
fv::options::New(this->mesh_).correct(this->nut_);
BasicTurbulenceModel::correctNut();
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template
void kOmegaSST::correctNut()
{
correctNut(2*magSqr(symm(fvc::grad(this->U_))), F23());
}
template
tmp kOmegaSST::epsilonByk
(
const volScalarField& F1,
const volScalarField& F2
) const
{
return betaStar_*omega_;
}
template
tmp
kOmegaSST::kSource() const
{
return tmp
(
new fvScalarMatrix
(
k_,
dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
)
);
}
template
tmp
kOmegaSST::omegaSource() const
{
return tmp
(
new fvScalarMatrix
(
omega_,
dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
)
);
}
template
tmp kOmegaSST::Qsas
(
const volScalarField& S2,
const volScalarField& gamma,
const volScalarField& beta
) const
{
return tmp
(
new fvScalarMatrix
(
omega_,
dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
)
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template
kOmegaSST::kOmegaSST
(
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
(
type,
alpha,
rho,
U,
alphaRhoPhi,
phi,
transport,
propertiesName
),
alphaK1_
(
dimensioned::lookupOrAddToDict
(
"alphaK1",
this->coeffDict_,
0.85
)
),
alphaK2_
(
dimensioned::lookupOrAddToDict
(
"alphaK2",
this->coeffDict_,
1.0
)
),
alphaOmega1_
(
dimensioned::lookupOrAddToDict
(
"alphaOmega1",
this->coeffDict_,
0.5
)
),
alphaOmega2_
(
dimensioned::lookupOrAddToDict
(
"alphaOmega2",
this->coeffDict_,
0.856
)
),
gamma1_
(
dimensioned::lookupOrAddToDict
(
"gamma1",
this->coeffDict_,
5.0/9.0
)
),
gamma2_
(
dimensioned::lookupOrAddToDict
(
"gamma2",
this->coeffDict_,
0.44
)
),
beta1_
(
dimensioned::lookupOrAddToDict
(
"beta1",
this->coeffDict_,
0.075
)
),
beta2_
(
dimensioned::lookupOrAddToDict
(
"beta2",
this->coeffDict_,
0.0828
)
),
betaStar_
(
dimensioned::lookupOrAddToDict
(
"betaStar",
this->coeffDict_,
0.09
)
),
a1_
(
dimensioned::lookupOrAddToDict
(
"a1",
this->coeffDict_,
0.31
)
),
b1_
(
dimensioned::lookupOrAddToDict
(
"b1",
this->coeffDict_,
1.0
)
),
c1_
(
dimensioned::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
bool kOmegaSST::read()
{
if (TurbulenceModel::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
void kOmegaSST::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_));
BasicTurbulenceModel::correct();
volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
tmp 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_.boundaryFieldRef().updateCoeffs();
volScalarField CDkOmega
(
(2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
);
volScalarField F1(this->F1(CDkOmega));
volScalarField F23(this->F23());
{
volScalarField gamma(this->gamma(F1));
volScalarField beta(this->beta(F1));
// Turbulent frequency equation
tmp 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()
+ fvOptions(alpha, rho, omega_)
);
omegaEqn.ref().relax();
fvOptions.constrain(omegaEqn.ref());
omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
solve(omegaEqn);
fvOptions.correct(omega_);
bound(omega_, this->omegaMin_);
}
// Turbulent kinetic energy equation
tmp 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*epsilonByk(F1, F23), k_)
+ kSource()
+ fvOptions(alpha, rho, k_)
);
kEqn.ref().relax();
fvOptions.constrain(kEqn.ref());
solve(kEqn);
fvOptions.correct(k_);
bound(k_, this->kMin_);
correctNut(S2, F23);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //