Merge branch 'feature-spalart-allmaras-model' into 'develop'

ENH: SpalartAllmaras: add estimation functions for k, epsilon and omega

See merge request Development/openfoam!441
This commit is contained in:
Andrew Heather 2021-05-27 08:59:05 +00:00
commit f3d715dd13
25 changed files with 231 additions and 140182 deletions

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2019-2020 OpenCFD Ltd.
Copyright (C) 2019-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -54,36 +54,40 @@ tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::fv1
) const
{
const volScalarField chi3(pow3(chi));
return chi3/(chi3 + pow3(Cv1_));
}
template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::fv2
tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::fv2
(
const volScalarField& chi,
const volScalarField& fv1
const volScalarField::Internal& chi,
const volScalarField::Internal& fv1
) const
{
return 1.0 - chi/(1.0 + chi*fv1);
return scalar(1) - chi/(scalar(1) + chi*fv1);
}
template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::Stilda
(
const volScalarField& chi,
const volScalarField& fv1
) const
tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::Stilda()
const
{
volScalarField Omega(::sqrt(2.0)*mag(skew(fvc::grad(this->U_))));
const volScalarField chi(this->chi());
const volScalarField fv1(this->fv1(chi));
const volScalarField::Internal Omega
(
::sqrt(scalar(2))*mag(skew(fvc::grad(this->U_)().v()))
);
return
(
max
(
Omega
+ fv2(chi, fv1)*nuTilda_/sqr(kappa_*y_),
Omega + fv2(chi(), fv1())*nuTilda_()/sqr(kappa_*y_),
Cs_*Omega
)
);
@ -91,12 +95,12 @@ tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::Stilda
template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::fw
tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::fw
(
const volScalarField& Stilda
const volScalarField::Internal& Stilda
) const
{
volScalarField r
const volScalarField::Internal r
(
min
(
@ -105,39 +109,33 @@ tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::fw
max
(
Stilda,
dimensionedScalar("SMALL", Stilda.dimensions(), SMALL)
dimensionedScalar(Stilda.dimensions(), SMALL)
)
*sqr(kappa_*y_)
),
scalar(10)
)
);
r.boundaryFieldRef() == 0.0;
const volScalarField g(r + Cw2_*(pow6(r) - r));
const volScalarField::Internal g(r + Cw2_*(pow6(r) - r));
return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
}
template<class BasicTurbulenceModel>
void SpalartAllmaras<BasicTurbulenceModel>::correctNut
(
const volScalarField& fv1
)
{
this->nut_ = nuTilda_*fv1;
this->nut_.correctBoundaryConditions();
fv::options::New(this->mesh_).correct(this->nut_);
BasicTurbulenceModel::correctNut();
return
g*pow
(
(scalar(1) + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)),
scalar(1)/scalar(6)
);
}
template<class BasicTurbulenceModel>
void SpalartAllmaras<BasicTurbulenceModel>::correctNut()
{
correctNut(fv1(this->chi()));
this->nut_ = nuTilda_*this->fv1(this->chi());
this->nut_.correctBoundaryConditions();
fv::options::New(this->mesh_).correct(this->nut_);
BasicTurbulenceModel::correctNut();
}
@ -174,7 +172,7 @@ SpalartAllmaras<BasicTurbulenceModel>::SpalartAllmaras
(
"sigmaNut",
this->coeffDict_,
0.66666
scalar(2)/scalar(3)
)
),
kappa_
@ -205,7 +203,7 @@ SpalartAllmaras<BasicTurbulenceModel>::SpalartAllmaras
0.622
)
),
Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
Cw1_(Cb1_/sqr(kappa_) + (scalar(1) + Cb2_)/sigmaNut_),
Cw2_
(
dimensioned<scalar>::getOrAddToDict
@ -277,7 +275,7 @@ bool SpalartAllmaras<BasicTurbulenceModel>::read()
Cb1_.readIfPresent(this->coeffDict());
Cb2_.readIfPresent(this->coeffDict());
Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
Cw1_ = Cb1_/sqr(kappa_) + (scalar(1) + Cb2_)/sigmaNut_;
Cw2_.readIfPresent(this->coeffDict());
Cw3_.readIfPresent(this->coeffDict());
Cv1_.readIfPresent(this->coeffDict());
@ -293,9 +291,10 @@ bool SpalartAllmaras<BasicTurbulenceModel>::read()
template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::DnuTildaEff() const
{
return tmp<volScalarField>
return tmp<volScalarField>::New
(
new volScalarField("DnuTildaEff", (nuTilda_ + this->nu())/sigmaNut_)
"DnuTildaEff",
(nuTilda_ + this->nu())/sigmaNut_
);
}
@ -303,22 +302,22 @@ tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::DnuTildaEff() const
template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::k() const
{
WarningInFunction
<< "Turbulence kinetic energy not defined for "
<< "Spalart-Allmaras model. Returning zero field"
<< endl;
// (B:Eq. 4.50)
const scalar Cmu = 0.09;
return tmp<volScalarField>::New
(
IOobject
(
"k",
IOobject::groupName("k", this->alphaRhoPhi_.group()),
this->runTime_.timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar(sqr(dimLength)/sqr(dimTime), Zero),
zeroGradientFvPatchField<scalar>::typeName
cbrt(this->fv1(this->chi()))
*nuTilda_
*::sqrt(scalar(2)/Cmu)
*mag(symm(fvc::grad(this->U_))),
this->nut_.boundaryField().types()
);
}
@ -326,22 +325,22 @@ tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::k() const
template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::epsilon() const
{
WarningInFunction
<< "Turbulence kinetic energy dissipation rate not defined for "
<< "Spalart-Allmaras model. Returning zero field"
<< endl;
// (B:Eq. 4.50)
const scalar Cmu = 0.09;
const dimensionedScalar nutSMALL(sqr(dimLength)/dimTime, SMALL);
return tmp<volScalarField>::New
(
IOobject
(
"epsilon",
IOobject::groupName("epsilon", this->alphaRhoPhi_.group()),
this->runTime_.timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar(sqr(dimLength)/pow3(dimTime), Zero),
zeroGradientFvPatchField<scalar>::typeName
pow(this->fv1(this->chi()), 0.5)
*pow(::sqrt(Cmu)*this->k(), 2)
/(nuTilda_ + this->nut_ + nutSMALL),
this->nut_.boundaryField().types()
);
}
@ -349,10 +348,9 @@ tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::epsilon() const
template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::omega() const
{
WarningInFunction
<< "Specific dissipation rate not defined for "
<< "Spalart-Allmaras model. Returning zero field"
<< endl;
// (P:p. 384)
const scalar betaStar = 0.09;
const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL);
return tmp<volScalarField>::New
(
@ -362,8 +360,8 @@ tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::omega() const
this->runTime_.timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar(dimless/dimTime, Zero)
this->epsilon()/(betaStar*(this->k() + k0)),
this->nut_.boundaryField().types()
);
}
@ -377,7 +375,7 @@ void SpalartAllmaras<BasicTurbulenceModel>::correct()
}
{
// Local references
// Construct local convenience references
const alphaField& alpha = this->alpha_;
const rhoField& rho = this->rho_;
const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
@ -385,10 +383,7 @@ void SpalartAllmaras<BasicTurbulenceModel>::correct()
eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
const volScalarField chi(this->chi());
const volScalarField fv1(this->fv1(chi));
const volScalarField Stilda(this->Stilda(chi, fv1));
const volScalarField::Internal Stilda(this->Stilda());
tmp<fvScalarMatrix> nuTildaEqn
(
@ -397,8 +392,8 @@ void SpalartAllmaras<BasicTurbulenceModel>::correct()
- fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
- Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
==
Cb1_*alpha*rho*Stilda*nuTilda_
- fvm::Sp(Cw1_*alpha*rho*fw(Stilda)*nuTilda_/sqr(y_), nuTilda_)
Cb1_*alpha()*rho()*Stilda*nuTilda_()
- fvm::Sp(Cw1_*alpha()*rho()*fw(Stilda)*nuTilda_()/sqr(y_), nuTilda_)
+ fvOptions(alpha, rho, nuTilda_)
);
@ -410,7 +405,7 @@ void SpalartAllmaras<BasicTurbulenceModel>::correct()
nuTilda_.correctBoundaryConditions();
}
// Update nut with latest available k,epsilon
// Update nut with latest available nuTilda
correctNut();
}

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2019-2020 OpenCFD Ltd.
Copyright (C) 2019-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -31,41 +31,82 @@ Group
grpRASTurbulence
Description
Spalart-Allmaras one-eqn mixing-length model for incompressible and
compressible external flows.
Spalart-Allmaras one-transport-equation linear-eddy-viscosity turbulence
closure model for incompressible and compressible external flows.
Reference:
Required fields
\verbatim
Spalart, P.R., & Allmaras, S.R. (1994).
A one-equation turbulence model for aerodynamic flows.
La Recherche Aerospatiale, 1, 5-21.
nuTilda | Modified kinematic viscosity [m2/s]
\endverbatim
The model is implemented without the trip-term and hence the ft2 term is
not needed.
It is necessary to limit the Stilda generation term as the model generates
unphysical results if this term becomes negative which occurs for complex
flow. Several approaches have been proposed to limit Stilda but it is not
clear which is the most appropriate. Here the limiter proposed by Spalart
is implemented in which Stilda is clipped at Cs*Omega with the default value
of Cs = 0.3.
The default model coefficients are
References:
\verbatim
Standard model:
Spalart, P.R., & Allmaras, S.R. (1994).
A one-equation turbulence model for aerodynamic flows.
La Recherche Aerospatiale, 1, 5-21.
Standard model without trip and ft2 terms (tag:R):
Rumsey, C. (2020).
The Spalart-Allmaras Turbulence Model.
Spalart-Allmaras One-Equation Model without ft2 Term (SA-noft2).
https://turbmodels.larc.nasa.gov/spalart.html#sanoft2
(Retrieved:12-01-2021).
Estimation expression for k and epsilon (tag:B), Eq. 4.50:
Bourgoin, A. (2019).
Bathymetry induced turbulence modelling the
Alderney Race site: regional approach with TELEMAC-LES.
Normandie Université.
Estimation expressions for omega (tag:P):
Pope, S. B. (2000).
Turbulent flows.
Cambridge, UK: Cambridge Univ. Press
DOI:10.1017/CBO9780511840531
\endverbatim
Usage
Example by using \c constant/turbulenceProperties:
\verbatim
RAS
{
// Mandatory entries (unmodifiable)
RASModel SpalartAllmaras;
// Optional entries (runtime modifiable)
turbulence on;
printCoeffs on;
SpalartAllmarasCoeffs
{
sigmaNut 0.66666;
kappa 0.41;
Cb1 0.1355;
Cb2 0.622;
Cw2 0.3;
Cw3 2.0;
Cv1 7.1;
Cs 0.3;
sigmaNut 0.66666;
kappa 0.41;
}
}
\endverbatim
Note
- The model is implemented without the trip-term since the model has almost
always been used in fully turbulent applications rather than those where
laminar-turbulent transition occurs.
- It has been argued that the \c ft2 term is not needed in the absence of the
trip-term, hence \c ft2 term is also not implementated.
- The \c Stilda generation term should never be allowed to be zero or negative
to avoid potential numerical issues and unphysical results for complex
flows. To this end, a limiter proposed by Spalart (R:Note-1(b)) is applied
onto \c Stilda where \c Stilda is clipped at \c Cs*Omega with the default
value of \c Cs=0.3.
- The model does not produce \c k, \c epsilon or \c omega. Nevertheless,
these quantities can be estimated by using an approximate expressions for
turbulent kinetic energy and dissipation rate reported in (B:Eq. 4.50).
SourceFiles
SpalartAllmaras.C
@ -104,13 +145,12 @@ class SpalartAllmaras
protected:
// Protected data
// Protected Data
// Model coefficients
dimensionedScalar sigmaNut_;
dimensionedScalar kappa_;
dimensionedScalar Cb1_;
dimensionedScalar Cb2_;
dimensionedScalar Cw1_;
@ -122,12 +162,13 @@ protected:
// Fields
//- Modified kinematic viscosity [m2/s]
volScalarField nuTilda_;
//- Wall distance
// Note: different to wall distance in parent RASModel
// which is for near-wall cells only
const volScalarField& y_;
const volScalarField::Internal& y_;
// Protected Member Functions
@ -136,21 +177,21 @@ protected:
tmp<volScalarField> fv1(const volScalarField& chi) const;
tmp<volScalarField> fv2
tmp<volScalarField::Internal> fv2
(
const volScalarField& chi,
const volScalarField& fv1
const volScalarField::Internal& chi,
const volScalarField::Internal& fv1
) const;
tmp<volScalarField> Stilda
tmp<volScalarField::Internal> Stilda() const;
tmp<volScalarField::Internal> fw
(
const volScalarField& chi,
const volScalarField& fv1
const volScalarField::Internal& Stilda
) const;
tmp<volScalarField> fw(const volScalarField& Stilda) const;
void correctNut(const volScalarField& fv1);
//- Update nut with the latest available nuTilda
virtual void correctNut();
@ -193,16 +234,16 @@ public:
//- Return the effective diffusivity for nuTilda
tmp<volScalarField> DnuTildaEff() const;
//- Return the turbulence kinetic energy
//- Return the (estimated) turbulent kinetic energy
virtual tmp<volScalarField> k() const;
//- Return the turbulence kinetic energy dissipation rate
//- Return the (estimated) turbulent kinetic energy dissipation rate
virtual tmp<volScalarField> epsilon() const;
//- Return the (estimated) specific dissipation rate
virtual tmp<volScalarField> omega() const;
//- Solve the turbulence equations and correct the turbulence viscosity
//- Solve the turbulence equations and correct the turbulent viscosity
virtual void correct();
};

View File

@ -43,4 +43,5 @@ boundaryField
}
}
// ************************************************************************* //

View File

@ -44,4 +44,5 @@ boundaryField
}
}
// ************************************************************************* //

View File

@ -44,4 +44,5 @@ boundaryField
}
}
// ************************************************************************* //

View File

@ -43,4 +43,5 @@ boundaryField
}
}
// ************************************************************************* //

View File

@ -1,9 +1,10 @@
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions # Tutorial clean functions
#------------------------------------------------------------------------------
# Clean time directories only
rm -rf *[1-9]*
rm -f log.*
cleanCase0
rm -rf constant/polyMesh
#------------------------------------------------------------------------------

View File

@ -3,6 +3,32 @@ cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
runApplication $(getApplication)
# settings
# flag to enable computations in parallel mode
parallel=true
#------------------------------------------------------------------------------
restore0Dir
[ -d constant/polyMesh ] || cp -rf constant/polyMesh.orig constant/polyMesh
if [ "$parallel" = true ]
then
runApplication decomposePar
runParallel $(getApplication)
runApplication reconstructPar
else
runApplication $(getApplication)
fi
#------------------------------------------------------------------------------

View File

@ -10,15 +10,13 @@ FoamFile
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
transportModel Newtonian;
rho 1;
nu 1e-05;
// ************************************************************************* //

View File

@ -10,20 +10,20 @@ FoamFile
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType RAS;
simulationType RAS;
RAS
{
RASModel SpalartAllmaras;
RASModel SpalartAllmaras;
turbulence on;
turbulence on;
printCoeffs on;
printCoeffs on;
}
// ************************************************************************* //

View File

@ -10,7 +10,6 @@ FoamFile
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -49,21 +48,52 @@ functions
{
momErr
{
type momentumError;
libs (fieldFunctionObjects);
type momentumError;
libs (fieldFunctionObjects);
executeControl writeTime;
writeControl writeTime;
}
contErr
{
type div;
libs (fieldFunctionObjects);
field phi;
type div;
libs (fieldFunctionObjects);
field phi;
executeControl writeTime;
writeControl writeTime;
}
turbulenceFields1
{
type turbulenceFields;
libs (fieldFunctionObjects);
fields
(
k
epsilon
nut
nuEff
R
devReff
L
I
);
executeControl writeTime;
writeControl writeTime;
}
mag1
{
type mag;
libs (fieldFunctionObjects);
field turbulenceProperties:R;
result magR;
executeControl writeTime;
writeControl writeTime;
}
}

View File

@ -9,43 +9,19 @@ FoamFile
{
version 2.0;
format ascii;
class polyBoundaryMesh;
object boundary;
class dictionary;
object decomposeParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
4
(
inlet
numberOfSubdomains 8;
method hierarchical;
coeffs
{
type patch;
physicalType inlet;
nFaces 134;
startFace 21254;
n (4 2 1);
}
outlet
{
type patch;
physicalType outlet;
nFaces 160;
startFace 21388;
}
walls
{
type wall;
physicalType wall;
nFaces 78;
startFace 21548;
}
frontAndBack
{
type empty;
physicalType empty;
nFaces 21440;
startFace 21626;
}
)
// ************************************************************************* //

View File

@ -10,7 +10,6 @@ FoamFile
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -28,8 +27,11 @@ gradSchemes
divSchemes
{
default none;
div(phi,U) bounded Gauss linearUpwind grad(U);
div(phi,nuTilda) bounded Gauss linearUpwind grad(nuTilda);
div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
@ -50,7 +52,7 @@ snGradSchemes
wallDist
{
method meshWave;
method meshWave;
}

View File

@ -10,7 +10,6 @@ FoamFile
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //