diff --git a/src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.C b/src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.C index 0c9ecfd1a8..3e38df486d 100644 --- a/src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.C +++ b/src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.C @@ -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 SpalartAllmaras::fv1 ) const { const volScalarField chi3(pow3(chi)); + return chi3/(chi3 + pow3(Cv1_)); } template -tmp SpalartAllmaras::fv2 +tmp SpalartAllmaras::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 -tmp SpalartAllmaras::Stilda -( - const volScalarField& chi, - const volScalarField& fv1 -) const +tmp SpalartAllmaras::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 SpalartAllmaras::Stilda template -tmp SpalartAllmaras::fw +tmp SpalartAllmaras::fw ( - const volScalarField& Stilda + const volScalarField::Internal& Stilda ) const { - volScalarField r + const volScalarField::Internal r ( min ( @@ -105,39 +109,33 @@ tmp SpalartAllmaras::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 -void SpalartAllmaras::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 void SpalartAllmaras::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::SpalartAllmaras ( "sigmaNut", this->coeffDict_, - 0.66666 + scalar(2)/scalar(3) ) ), kappa_ @@ -205,7 +203,7 @@ SpalartAllmaras::SpalartAllmaras 0.622 ) ), - Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_), + Cw1_(Cb1_/sqr(kappa_) + (scalar(1) + Cb2_)/sigmaNut_), Cw2_ ( dimensioned::getOrAddToDict @@ -277,7 +275,7 @@ bool SpalartAllmaras::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::read() template tmp SpalartAllmaras::DnuTildaEff() const { - return tmp + return tmp::New ( - new volScalarField("DnuTildaEff", (nuTilda_ + this->nu())/sigmaNut_) + "DnuTildaEff", + (nuTilda_ + this->nu())/sigmaNut_ ); } @@ -303,22 +302,22 @@ tmp SpalartAllmaras::DnuTildaEff() const template tmp SpalartAllmaras::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::New ( IOobject ( - "k", + IOobject::groupName("k", this->alphaRhoPhi_.group()), this->runTime_.timeName(), this->mesh_ ), - this->mesh_, - dimensionedScalar(sqr(dimLength)/sqr(dimTime), Zero), - zeroGradientFvPatchField::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 SpalartAllmaras::k() const template tmp SpalartAllmaras::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::New ( IOobject ( - "epsilon", + IOobject::groupName("epsilon", this->alphaRhoPhi_.group()), this->runTime_.timeName(), this->mesh_ ), - this->mesh_, - dimensionedScalar(sqr(dimLength)/pow3(dimTime), Zero), - zeroGradientFvPatchField::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 SpalartAllmaras::epsilon() const template tmp SpalartAllmaras::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::New ( @@ -362,8 +360,8 @@ tmp SpalartAllmaras::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::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::correct() eddyViscosity>::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 nuTildaEqn ( @@ -397,8 +392,8 @@ void SpalartAllmaras::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::correct() nuTilda_.correctBoundaryConditions(); } - // Update nut with latest available k,epsilon + // Update nut with latest available nuTilda correctNut(); } diff --git a/src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.H b/src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.H index 853bf82c3d..5391c7d12c 100644 --- a/src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.H +++ b/src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.H @@ -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 fv1(const volScalarField& chi) const; - tmp fv2 + tmp fv2 ( - const volScalarField& chi, - const volScalarField& fv1 + const volScalarField::Internal& chi, + const volScalarField::Internal& fv1 ) const; - tmp Stilda + tmp Stilda() const; + + tmp fw ( - const volScalarField& chi, - const volScalarField& fv1 + const volScalarField::Internal& Stilda ) const; - tmp 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 DnuTildaEff() const; - //- Return the turbulence kinetic energy + //- Return the (estimated) turbulent kinetic energy virtual tmp k() const; - //- Return the turbulence kinetic energy dissipation rate + //- Return the (estimated) turbulent kinetic energy dissipation rate virtual tmp epsilon() const; //- Return the (estimated) specific dissipation rate virtual tmp omega() const; - //- Solve the turbulence equations and correct the turbulence viscosity + //- Solve the turbulence equations and correct the turbulent viscosity virtual void correct(); };