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

ENH: SpalartAllmaras: reduce peak-memory usage

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<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::Internal& chi,
-    const volScalarField& fv1
+    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
+tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::Stilda()
-(
+const
-    const volScalarField& chi,
-    const volScalarField& fv1
-) 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
+            Omega + fv2(chi(), fv1())*nuTilda_()/sqr(kappa_*y_),
-          + 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);
+    return
-}
+        g*pow
-
+        (
-
+            (scalar(1) + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)),
-template<class BasicTurbulenceModel>
+            scalar(1)/scalar(6)
-void SpalartAllmaras<BasicTurbulenceModel>::correctNut
+        );
-(
-    const volScalarField& fv1
-)
-{
-    this->nut_ = nuTilda_*fv1;
-    this->nut_.correctBoundaryConditions();
-    fv::options::New(this->mesh_).correct(this->nut_);
-
-    BasicTurbulenceModel::correctNut();
 }
 
 
 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
+    // (B:Eq. 4.50)
-        << "Turbulence kinetic energy not defined for "
+    const scalar Cmu = 0.09;
-        << "Spalart-Allmaras model. Returning zero field"
-        << endl;
 
     return tmp<volScalarField>::New
     (
         IOobject
         (
-            "k",
+            IOobject::groupName("k", this->alphaRhoPhi_.group()),
             this->runTime_.timeName(),
             this->mesh_
         ),
-        this->mesh_,
+        cbrt(this->fv1(this->chi()))
-        dimensionedScalar(sqr(dimLength)/sqr(dimTime), Zero),
+        *nuTilda_
-        zeroGradientFvPatchField<scalar>::typeName
+        *::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
+    // (B:Eq. 4.50)
-        << "Turbulence kinetic energy dissipation rate not defined for "
+    const scalar Cmu = 0.09;
-        << "Spalart-Allmaras model. Returning zero field"
+    const dimensionedScalar nutSMALL(sqr(dimLength)/dimTime, SMALL);
-        << endl;
 
     return tmp<volScalarField>::New
     (
         IOobject
         (
-            "epsilon",
+            IOobject::groupName("epsilon", this->alphaRhoPhi_.group()),
             this->runTime_.timeName(),
             this->mesh_
         ),
-        this->mesh_,
+        pow(this->fv1(this->chi()), 0.5)
-        dimensionedScalar(sqr(dimLength)/pow3(dimTime), Zero),
+        *pow(::sqrt(Cmu)*this->k(), 2)
-        zeroGradientFvPatchField<scalar>::typeName
+        /(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
+    // (P:p. 384)
-        << "Specific dissipation rate not defined for "
+    const scalar betaStar = 0.09;
-        << "Spalart-Allmaras model. Returning zero field"
+    const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL);
-        << endl;
 
     return tmp<volScalarField>::New
     (
@@ -362,8 +360,8 @@ tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::omega() const
             this->runTime_.timeName(),
             this->mesh_
         ),
-        this->mesh_,
+        this->epsilon()/(betaStar*(this->k() + k0)),
-        dimensionedScalar(dimless/dimTime, Zero)
+        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::Internal Stilda(this->Stilda());
-        const volScalarField fv1(this->fv1(chi));
-
-        const volScalarField Stilda(this->Stilda(chi, fv1));
 
         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_
+            Cb1_*alpha()*rho()*Stilda*nuTilda_()
-          - fvm::Sp(Cw1_*alpha*rho*fw(Stilda)*nuTilda_/sqr(y_), 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();
 }
 

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
+    Spalart-Allmaras one-transport-equation linear-eddy-viscosity turbulence
-    compressible external flows.
+    closure model for incompressible and compressible external flows.
 
-    Reference:
+    Required fields
     \verbatim
-        Spalart, P.R., & Allmaras, S.R. (1994).
+        nuTilda   | Modified kinematic viscosity [m2/s]
-        A one-equation turbulence model for aerodynamic flows.
-        La Recherche Aerospatiale, 1, 5-21.
     \endverbatim
 
-    The model is implemented without the trip-term and hence the ft2 term is
+    References:
-    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
     \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::Internal& chi,
-            const volScalarField& fv1
+            const volScalarField::Internal& fv1
         ) const;
 
-        tmp<volScalarField> Stilda
+        tmp<volScalarField::Internal> Stilda() const;
+
+        tmp<volScalarField::Internal> fw
         (
-            const volScalarField& chi,
+            const volScalarField::Internal& Stilda
-            const volScalarField& fv1
         ) 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();
 };