ENH: Refactored Spalart-Allmaras turbulence models

- Added a new S-A base class: SpalartAllmarasBase
- RAS and DES models derived from new base class
- Removed code duplication

src/TurbulenceModels/turbulenceModels/Base/SpalartAllmaras/SpalartAllmarasBase.C

@@ -0,0 +1,471 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2016-2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "SpalartAllmarasBase.H"
+#include "wallDist.H"
+#include "bound.H"
+#include "fvOptions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::chi() const
+{
+    return nuTilda_/this->nu();
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::fv1
+(
+    const volScalarField& chi
+) const
+{
+    const volScalarField chi3("chi3", pow3(chi));
+    return chi3/(chi3 + pow3(Cv1_));
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::fv2
+(
+    const volScalarField& chi,
+    const volScalarField& fv1
+) const
+{
+    return scalar(1) - chi/(scalar(1) + chi*fv1);
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::ft2
+(
+    const volScalarField& chi
+) const
+{
+    return Ct3_*exp(-Ct4_*sqr(chi));
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::Omega
+(
+    const volTensorField& gradU
+) const
+{
+    return sqrt(2.0)*mag(skew(gradU));
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::r
+(
+    const volScalarField& nur,
+    const volScalarField& Stilda,
+    const volScalarField& dTilda
+) const
+{
+    const dimensionedScalar eps("SMALL", Stilda.dimensions(), SMALL);
+
+    tmp<volScalarField> tr =
+        min(nur/(max(Stilda, eps)*sqr(kappa_*dTilda)), scalar(10));
+
+    tr.ref().boundaryFieldRef() == 0;
+
+    return tr;
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField::Internal> SpalartAllmarasBase<BasicEddyViscosityModel>::fw
+(
+    const volScalarField& Stilda,
+    const volScalarField& dTilda
+) const
+{
+    const volScalarField::Internal r(this->r(nuTilda_, Stilda, dTilda)()());
+    const volScalarField::Internal g(r + Cw2_*(pow6(r) - r));
+
+    return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::Stilda
+(
+    const volScalarField& chi,
+    const volScalarField& fv1,
+    const volTensorField& gradU,
+    const volScalarField& dTilda
+) const
+{
+    const volScalarField Omega(this->Omega(gradU));
+
+    return
+        max
+        (
+            Omega + fv2(chi, fv1)*nuTilda_/sqr(kappa_*dTilda),
+            Cs_*Omega
+        );
+}
+
+
+template<class BasicEddyViscosityModel>
+void SpalartAllmarasBase<BasicEddyViscosityModel>::correctNut
+(
+    const volScalarField& fv1
+)
+{
+    this->nut_ = nuTilda_*fv1;
+    this->nut_.correctBoundaryConditions();
+    fv::options::New(this->mesh_).correct(this->nut_);
+}
+
+
+template<class BasicEddyViscosityModel>
+void SpalartAllmarasBase<BasicEddyViscosityModel>::correctNut()
+{
+    correctNut(fv1(this->chi()));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicEddyViscosityModel>
+SpalartAllmarasBase<BasicEddyViscosityModel>::SpalartAllmarasBase
+(
+    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
+    ),
+
+    sigmaNut_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "sigmaNut",
+            this->coeffDict_,
+            0.66666
+        )
+    ),
+    kappa_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "kappa",
+            this->coeffDict_,
+            0.41
+        )
+    ),
+    Cb1_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cb1",
+            this->coeffDict_,
+            0.1355
+        )
+    ),
+    Cb2_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cb2",
+            this->coeffDict_,
+            0.622
+        )
+    ),
+    Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
+    Cw2_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cw2",
+            this->coeffDict_,
+            0.3
+        )
+    ),
+    Cw3_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cw3",
+            this->coeffDict_,
+            2.0
+        )
+    ),
+    Cv1_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cv1",
+            this->coeffDict_,
+            7.1
+        )
+    ),
+    Cs_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cs",
+            this->coeffDict_,
+            0.3
+        )
+    ),
+    ck_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "ck",
+            this->coeffDict_,
+            0.07
+        )
+    ),
+    Ct3_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ct3",
+            this->coeffDict_,
+            1.2
+        )
+    ),
+    Ct4_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ct4",
+            this->coeffDict_,
+            0.5
+        )
+    ),
+
+    nuTilda_
+    (
+        IOobject
+        (
+            "nuTilda",
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+
+    y_(wallDist::New(this->mesh_).y())
+{
+    if (mag(Ct3_.value()) > SMALL)
+    {
+        Info<< "    ft2 term: active" << nl;
+    }
+    else
+    {
+        Info<< "    ft2 term: inactive" << nl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicEddyViscosityModel>
+bool SpalartAllmarasBase<BasicEddyViscosityModel>::read()
+{
+    if (BasicEddyViscosityModel::read())
+    {
+        sigmaNut_.readIfPresent(this->coeffDict());
+        kappa_.readIfPresent(*this);
+
+        Cb1_.readIfPresent(this->coeffDict());
+        Cb2_.readIfPresent(this->coeffDict());
+        Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
+        Cw2_.readIfPresent(this->coeffDict());
+        Cw3_.readIfPresent(this->coeffDict());
+        Cv1_.readIfPresent(this->coeffDict());
+        Cs_.readIfPresent(this->coeffDict());
+
+        ck_.readIfPresent(this->coeffDict());
+
+        Ct3_.readIfPresent(this->coeffDict());
+        Ct4_.readIfPresent(this->coeffDict());
+
+        if (mag(Ct3_.value()) > SMALL)
+        {
+            Info<< "    ft2 term: active" << nl;
+        }
+        else
+        {
+            Info<< "    ft2 term: inactive" << nl;
+        }
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField>
+SpalartAllmarasBase<BasicEddyViscosityModel>::DnuTildaEff() const
+{
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("DnuTildaEff", this->alphaRhoPhi_.group()),
+        (nuTilda_ + this->nu())/sigmaNut_
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::k() const
+{
+    // (B:Eq. 4.50)
+    const scalar Cmu = 0.09;
+    const auto fv1 = this->fv1(chi());
+
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("k", this->alphaRhoPhi_.group()),
+        cbrt(fv1)*nuTilda_*::sqrt(scalar(2)/Cmu)*mag(symm(fvc::grad(this->U_)))
+    );
+}
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField>
+SpalartAllmarasBase<BasicEddyViscosityModel>::epsilon() const
+{
+    // (B:Eq. 4.50)
+    const scalar Cmu = 0.09;
+    const auto fv1 = this->fv1(chi());
+    const dimensionedScalar nutSMALL(nuTilda_.dimensions(), SMALL);
+
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("epsilon", this->alphaRhoPhi_.group()),
+        sqrt(fv1)*sqr(::sqrt(Cmu)*this->k())/(nuTilda_ + this->nut_ + nutSMALL)
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::omega() const
+{
+    // (P:p. 384)
+    const scalar betaStar = 0.09;
+    const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL);
+
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("omega", this->alphaRhoPhi_.group()),
+        this->epsilon()/(betaStar*(this->k() + k0))
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+void SpalartAllmarasBase<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_;
+    fv::options& fvOptions(fv::options::New(this->mesh_));
+
+    BasicEddyViscosityModel::correct();
+
+    const volScalarField chi(this->chi());
+    const volScalarField fv1(this->fv1(chi));
+    const volScalarField ft2(this->ft2(chi));
+
+    tmp<volTensorField> tgradU = fvc::grad(U);
+    volScalarField dTilda(this->dTilda(chi, fv1, tgradU()));
+    volScalarField Stilda(this->Stilda(chi, fv1, tgradU(), dTilda));
+    tgradU.clear();
+
+    tmp<fvScalarMatrix> nuTildaEqn
+    (
+        fvm::ddt(alpha, rho, nuTilda_)
+      + fvm::div(alphaRhoPhi, nuTilda_)
+      - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
+      - Cb2_/sigmaNut_*alpha()*rho()*magSqr(fvc::grad(nuTilda_)()())
+     ==
+        Cb1_*alpha()*rho()*Stilda()*nuTilda_()*(scalar(1) - ft2())
+      - fvm::Sp
+        (
+            (Cw1_*fw(Stilda, dTilda) - Cb1_/sqr(kappa_)*ft2())
+           *alpha()*rho()*nuTilda_()/sqr(dTilda()),
+            nuTilda_
+        )
+      + fvOptions(alpha, rho, nuTilda_)
+    );
+
+    nuTildaEqn.ref().relax();
+    fvOptions.constrain(nuTildaEqn.ref());
+    solve(nuTildaEqn);
+    fvOptions.correct(nuTilda_);
+    bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), Zero));
+    nuTilda_.correctBoundaryConditions();
+
+    correctNut();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/Base/SpalartAllmaras/SpalartAllmarasBase.H

@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2015-2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::LESModels::SpalartAllmarasBase
+
+Group
+    grpDESTurbulence
+
+Description
+    SpalartAllmarasBase DES turbulence model for incompressible and
+    compressible flows
+
+    Reference:
+    \verbatim
+        Spalart, P. R., Jou, W. H., Strelets, M., & Allmaras, S. R. (1997).
+        Comments on the feasibility of LES for wings, and on a hybrid
+        RANS/LES approach.
+        Advances in DNS/LES, 1, 4-8.
+    \endverbatim
+
+SourceFiles
+    SpalartAllmarasBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_SpalartAllmarasBase_H
+#define Foam_SpalartAllmarasBase_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class SpalartAllmarasBase Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicEddyViscosityModel>
+class SpalartAllmarasBase
+:
+    public BasicEddyViscosityModel
+{
+    // Private Member Functions
+
+        //- No copy construct
+        SpalartAllmarasBase(const SpalartAllmarasBase&) = delete;
+
+        //- No copy assignment
+        void operator=(const SpalartAllmarasBase&) = delete;
+
+
+protected:
+
+    // Protected data
+
+        // Model constants
+
+            dimensionedScalar sigmaNut_;
+            dimensionedScalar kappa_;
+
+            dimensionedScalar Cb1_;
+            dimensionedScalar Cb2_;
+            dimensionedScalar Cw1_;
+            dimensionedScalar Cw2_;
+            dimensionedScalar Cw3_;
+            dimensionedScalar Cv1_;
+            dimensionedScalar Cs_;
+            dimensionedScalar ck_;
+
+
+            dimensionedScalar Ct3_;
+            dimensionedScalar Ct4_;
+
+
+        // Fields
+
+            volScalarField nuTilda_;
+
+            //- Wall distance
+            //  Note: different to wall distance in parent RASModel
+            //  which is for near-wall cells only
+            const volScalarField& y_;
+
+
+    // Protected Member Functions
+
+        tmp<volScalarField> chi() const;
+
+        tmp<volScalarField> fv1(const volScalarField& chi) const;
+
+        tmp<volScalarField> fv2
+        (
+            const volScalarField& chi,
+            const volScalarField& fv1
+        ) const;
+
+        tmp<volScalarField> ft2(const volScalarField& chi) const;
+
+        tmp<volScalarField> Omega(const volTensorField& gradU) const;
+
+        tmp<volScalarField> r
+        (
+            const volScalarField& nur,
+            const volScalarField& Stilda,
+            const volScalarField& dTilda
+        ) const;
+
+        tmp<volScalarField::Internal> fw
+        (
+            const volScalarField& Stilda,
+            const volScalarField& dTilda
+        ) const;
+
+        virtual tmp<volScalarField> Stilda
+        (
+            const volScalarField& chi,
+            const volScalarField& fv1,
+            const volTensorField& gradU,
+            const volScalarField& dTilda
+        ) const;
+
+        //- Length scale
+        virtual tmp<volScalarField> dTilda
+        (
+            const volScalarField& chi,
+            const volScalarField& fv1,
+            const volTensorField& gradU
+        ) const = 0;
+
+        void correctNut(const volScalarField& fv1);
+        virtual void correctNut();
+
+
+public:
+
+    typedef typename BasicEddyViscosityModel::alphaField alphaField;
+    typedef typename BasicEddyViscosityModel::rhoField rhoField;
+    typedef typename BasicEddyViscosityModel::transportModel transportModel;
+
+
+    // Constructors
+
+        //- Construct from components
+        SpalartAllmarasBase
+        (
+            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 ~SpalartAllmarasBase() = default;
+
+
+    // Member Functions
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return the effective diffusivity for nuTilda
+        tmp<volScalarField> DnuTildaEff() const;
+
+        //- Return the (estimated) turbulent kinetic energy
+        virtual tmp<volScalarField> k() const;
+
+        //- Return the (estimated) turbulent kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const;
+
+        //- Return the (estimated) specific dissipation rate
+        virtual tmp<volScalarField> omega() const;
+
+        tmp<volScalarField> nuTilda() const
+        {
+            return nuTilda_;
+        }
+
+        //- Correct nuTilda and related properties
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "SpalartAllmarasBase.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/DES/DESModel/DESModel.C

@@ -63,7 +63,7 @@ DESModel<BasicTurbulenceModel>::DESModel
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
 DESModel<BasicTurbulenceModel>::~DESModel()

src/TurbulenceModels/turbulenceModels/DES/DESModel/DESModel.H

@@ -92,7 +92,7 @@ public:
 
 
     //- Destructor
-    virtual ~DESModel();
+    virtual ~DESModel() = default;
 
 
     // Public Member Functions

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasDDES/SpalartAllmarasDDES.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,41 +37,15 @@ namespace LESModels
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::rd
-(
-    const volScalarField& magGradU
-) const
-{
-    tmp<volScalarField> tr
-    (
-        min
-        (
-            this->nuEff()
-           /(
-                max
-                (
-                    magGradU,
-                    dimensionedScalar("SMALL", magGradU.dimensions(), SMALL)
-                )
-               *sqr(this->kappa_*this->y_)
-            ),
-            scalar(10)
-        )
-    );
-    tr.ref().boundaryFieldRef() == 0.0;
-
-    return tr;
-}
-
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::fd
 (
     const volScalarField& magGradU
 ) const
 {
-    return 1 - tanh(pow(Cd1_*rd(magGradU), Cd2_));
+    return
+        1
+      - tanh(pow(this->Cd1_*this->r(this->nuEff(), magGradU, this->y_), Cd2_));
 }
 
 
@@ -86,7 +60,7 @@ tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::dTilda
 ) const
 {
     const volScalarField& lRAS(this->y_);
-    const volScalarField lLES(this->psi(chi, fv1)*this->CDES_*this->delta());
+    const volScalarField lLES(this->lengthScaleLES(chi, fv1));
 
     return max
     (

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasDDES/SpalartAllmarasDDES.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,10 +48,10 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SpalartAllmarasDDES_H
+#ifndef Foam_SpalartAllmarasDDES_H
-#define SpalartAllmarasDDES_H
+#define Foam_SpalartAllmarasDDES_H
 
-#include "SpalartAllmarasDES.H"
+#include "SpalartAllmarasBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -71,10 +71,9 @@ class SpalartAllmarasDDES
 {
     // Private Member Functions
 
+        //- Delay function
         tmp<volScalarField> fd(const volScalarField& magGradU) const;
 
-        tmp<volScalarField> rd(const volScalarField& magGradU) const;
-
         //- No copy construct
         SpalartAllmarasDDES(const SpalartAllmarasDDES&) = delete;
 

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasDES/SpalartAllmarasDES.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,120 +38,6 @@ namespace LESModels
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::chi() const
-{
-    return nuTilda_/this->nu();
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fv1
-(
-    const volScalarField& chi
-) const
-{
-    const volScalarField chi3("chi3", pow3(chi));
-    return chi3/(chi3 + pow3(Cv1_));
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fv2
-(
-    const volScalarField& chi,
-    const volScalarField& fv1
-) const
-{
-    return 1.0 - chi/(1.0 + chi*fv1);
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::ft2
-(
-    const volScalarField& chi
-) const
-{
-    return Ct3_*exp(-Ct4_*sqr(chi));
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Omega
-(
-    const volTensorField& gradU
-) const
-{
-    return sqrt(2.0)*mag(skew(gradU));
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Stilda
-(
-    const volScalarField& chi,
-    const volScalarField& fv1,
-    const volScalarField& Omega,
-    const volScalarField& dTilda
-) const
-{
-    return
-    (
-        max
-        (
-            Omega
-          + fv2(chi, fv1)*nuTilda_/sqr(kappa_*dTilda),
-            Cs_*Omega
-        )
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::r
-(
-    const volScalarField& nur,
-    const volScalarField& Omega,
-    const volScalarField& dTilda
-) const
-{
-    tmp<volScalarField> tr
-    (
-        min
-        (
-            nur
-           /(
-                max
-                (
-                    Omega,
-                    dimensionedScalar("SMALL", Omega.dimensions(), SMALL)
-                )
-                *sqr(kappa_*dTilda)
-            ),
-            scalar(10)
-        )
-    );
-    tr.ref().boundaryFieldRef() == 0.0;
-
-    return tr;
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fw
-(
-    const volScalarField& Omega,
-    const volScalarField& dTilda
-) const
-{
-    const volScalarField r(this->r(nuTilda_, Omega, dTilda));
-    const volScalarField g(r + Cw2_*(pow6(r) - r));
-
-    return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
-}
-
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi
 (
@@ -159,26 +45,23 @@ tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi
     const volScalarField& fv1
 ) const
 {
-    tmp<volScalarField> tpsi
+    auto tpsi = tmp<volScalarField>::New
-    (
-        new volScalarField
     (
         IOobject
         (
-                type() + ":psi",
+            IOobject::scopedName(type(), "psi"),
             this->time().timeName(),
             this->mesh(),
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         this->mesh(),
-            dimensionedScalar("one", dimless, 1)
+        dimensionedScalar(dimless, Zero)
-        )
     );
 
     if (lowReCorrection_)
     {
-        volScalarField& psi = tpsi.ref();
+        auto& psi = tpsi.ref();
 
         const volScalarField fv2(this->fv2(chi, fv1));
         const volScalarField ft2(this->ft2(chi));
@@ -189,7 +72,8 @@ tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi
                 min
                 (
                     scalar(100),
-                    (1 - Cb1_/(Cw1_*sqr(kappa_)*fwStar_)*(ft2 + (1 - ft2)*fv2))
+                    (1 - this->Cb1_/(this->Cw1_*sqr(this->kappa_)*fwStar_)
+                   *(ft2 + (1 - ft2)*fv2))
                    /max(SMALL, (fv1*max(scalar(1e-10), 1 - ft2)))
                 )
             );
@@ -199,6 +83,17 @@ tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi
 }
 
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::lengthScaleLES
+(
+    const volScalarField& chi,
+    const volScalarField& fv1
+) const
+{
+    return psi(chi, fv1)*CDES_*this->delta();
+}
+
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::dTilda
 (
@@ -207,33 +102,15 @@ tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::dTilda
     const volTensorField& gradU
 ) const
 {
-    tmp<volScalarField> tdTilda(psi(chi, fv1)*CDES_*this->delta());
+    // Initialise with LES length scale
-    min(tdTilda.ref().ref(), tdTilda(), y_);
+    tmp<volScalarField> tdTilda = lengthScaleLES(chi, fv1);
+
+    // Take min vs wall distance
+    min(tdTilda.ref(), tdTilda(), this->y_);
     return tdTilda;
 }
 
 
-template<class BasicTurbulenceModel>
-void SpalartAllmarasDES<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 SpalartAllmarasDES<BasicTurbulenceModel>::correctNut()
-{
-    correctNut(fv1(this->chi()));
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
@@ -249,7 +126,7 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
     const word& type
 )
 :
-    DESModel<BasicTurbulenceModel>
+    SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>
     (
         type,
         alpha,
@@ -261,79 +138,6 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
         propertiesName
     ),
 
-    sigmaNut_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "sigmaNut",
-            this->coeffDict_,
-            0.66666
-        )
-    ),
-    kappa_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "kappa",
-            this->coeffDict_,
-            0.41
-        )
-    ),
-    Cb1_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cb1",
-            this->coeffDict_,
-            0.1355
-        )
-    ),
-    Cb2_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cb2",
-            this->coeffDict_,
-            0.622
-        )
-    ),
-    Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
-    Cw2_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cw2",
-            this->coeffDict_,
-            0.3
-        )
-    ),
-    Cw3_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cw3",
-            this->coeffDict_,
-            2.0
-        )
-    ),
-    Cv1_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cv1",
-            this->coeffDict_,
-            7.1
-        )
-    ),
-    Cs_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cs",
-            this->coeffDict_,
-            0.3
-        )
-    ),
     CDES_
     (
         dimensioned<scalar>::getOrAddToDict
@@ -343,15 +147,6 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
             0.65
         )
     ),
-    ck_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "ck",
-            this->coeffDict_,
-            0.07
-        )
-    ),
     lowReCorrection_
     (
         Switch::getOrAddToDict
@@ -361,24 +156,6 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
             true
         )
     ),
-    Ct3_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Ct3",
-            this->coeffDict_,
-            1.2
-        )
-    ),
-    Ct4_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Ct4",
-            this->coeffDict_,
-            0.5
-        )
-    ),
     fwStar_
     (
         dimensioned<scalar>::getOrAddToDict
@@ -387,22 +164,7 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
             this->coeffDict_,
             0.424
         )
-    ),
+    )
-
-    nuTilda_
-    (
-        IOobject
-        (
-            "nuTilda",
-            this->runTime_.timeName(),
-            this->mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        this->mesh_
-    ),
-
-    y_(wallDist::New(this->mesh_).y())
 {
     if (type == typeName)
     {
@@ -416,25 +178,10 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
 template<class BasicTurbulenceModel>
 bool SpalartAllmarasDES<BasicTurbulenceModel>::read()
 {
-    if (DESModel<BasicTurbulenceModel>::read())
+    if (SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>::read())
     {
-        sigmaNut_.readIfPresent(this->coeffDict());
-        kappa_.readIfPresent(*this);
-
-        Cb1_.readIfPresent(this->coeffDict());
-        Cb2_.readIfPresent(this->coeffDict());
-        Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
-        Cw2_.readIfPresent(this->coeffDict());
-        Cw3_.readIfPresent(this->coeffDict());
-        Cv1_.readIfPresent(this->coeffDict());
-        Cs_.readIfPresent(this->coeffDict());
-
         CDES_.readIfPresent(this->coeffDict());
-        ck_.readIfPresent(this->coeffDict());
-
         lowReCorrection_.readIfPresent("lowReCorrection", this->coeffDict());
-        Ct3_.readIfPresent(this->coeffDict());
-        Ct4_.readIfPresent(this->coeffDict());
         fwStar_.readIfPresent(this->coeffDict());
 
         return true;
@@ -445,124 +192,24 @@ bool SpalartAllmarasDES<BasicTurbulenceModel>::read()
 
 
 template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::
+tmp<volScalarField>
-DnuTildaEff() const
+SpalartAllmarasDES<BasicTurbulenceModel>::LESRegion() const
-{
-    return tmp<volScalarField>
-    (
-        new volScalarField("DnuTildaEff", (nuTilda_ + this->nu())/sigmaNut_)
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::k() const
 {
     const volScalarField chi(this->chi());
     const volScalarField fv1(this->fv1(chi));
 
-    tmp<volScalarField> tdTilda
+    return tmp<volScalarField>::New
-    (
-        new volScalarField
     (
         IOobject
         (
-                "dTilda",
+            IOobject::scopedName("DES", "LESRegion"),
             this->mesh_.time().timeName(),
             this->mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-            this->mesh_,
+        neg(dTilda(chi, fv1, fvc::grad(this->U_)) - this->y_)
-            dimensionedScalar(dimLength, Zero),
-            zeroGradientFvPatchField<vector>::typeName
-        )
     );
-    volScalarField& dTildaF = tdTilda.ref();
-    dTildaF = dTilda(chi, fv1, fvc::grad(this->U_));
-    dTildaF.correctBoundaryConditions();
-
-    return sqr(this->nut()/ck_/dTildaF);
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::LESRegion() const
-{
-    const volScalarField chi(this->chi());
-    const volScalarField fv1(this->fv1(chi));
-
-    tmp<volScalarField> tLESRegion
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "DES::LESRegion",
-                this->mesh_.time().timeName(),
-                this->mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            neg(dTilda(chi, fv1, fvc::grad(this->U_)) - y_)
-        )
-    );
-
-    return tLESRegion;
-}
-
-
-template<class BasicTurbulenceModel>
-void SpalartAllmarasDES<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_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
-
-    DESModel<BasicTurbulenceModel>::correct();
-
-    const volScalarField chi(this->chi());
-    const volScalarField fv1(this->fv1(chi));
-    const volScalarField ft2(this->ft2(chi));
-
-    tmp<volTensorField> tgradU = fvc::grad(U);
-    const volScalarField Omega(this->Omega(tgradU()));
-    const volScalarField dTilda(this->dTilda(chi, fv1, tgradU()));
-    const volScalarField Stilda(this->Stilda(chi, fv1, Omega, dTilda));
-
-    tmp<fvScalarMatrix> nuTildaEqn
-    (
-        fvm::ddt(alpha, rho, nuTilda_)
-      + fvm::div(alphaRhoPhi, nuTilda_)
-      - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
-      - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
-     ==
-        Cb1_*alpha*rho*Stilda*nuTilda_*(scalar(1) - ft2)
-      - fvm::Sp
-        (
-            (Cw1_*fw(Stilda, dTilda) - Cb1_/sqr(kappa_)*ft2)
-           *alpha*rho*nuTilda_/sqr(dTilda),
-            nuTilda_
-        )
-      + fvOptions(alpha, rho, nuTilda_)
-    );
-
-    nuTildaEqn.ref().relax();
-    fvOptions.constrain(nuTildaEqn.ref());
-    solve(nuTildaEqn);
-    fvOptions.correct(nuTilda_);
-    bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), Zero));
-    nuTilda_.correctBoundaryConditions();
-
-    correctNut();
 }
 
 

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasDES/SpalartAllmarasDES.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2019 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,9 +61,10 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SpalartAllmarasDES_H
+#ifndef Foam_SpalartAllmarasDES_H
-#define SpalartAllmarasDES_H
+#define Foam_SpalartAllmarasDES_H
 
+#include "SpalartAllmarasBase.H"
 #include "DESModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -80,7 +81,7 @@ namespace LESModels
 template<class BasicTurbulenceModel>
 class SpalartAllmarasDES
 :
-    public DESModel<BasicTurbulenceModel>
+    public SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>
 {
     // Private Member Functions
 
@@ -97,76 +98,25 @@ protected:
 
         // Model constants
 
-            dimensionedScalar sigmaNut_;
+            // DES coefficient
-            dimensionedScalar kappa_;
-
-            dimensionedScalar Cb1_;
-            dimensionedScalar Cb2_;
-            dimensionedScalar Cw1_;
-            dimensionedScalar Cw2_;
-            dimensionedScalar Cw3_;
-            dimensionedScalar Cv1_;
-            dimensionedScalar Cs_;
             dimensionedScalar CDES_;
-            dimensionedScalar ck_;
-
 
             // Low Reynolds number correction
-
             Switch lowReCorrection_;
-                dimensionedScalar Ct3_;
-                dimensionedScalar Ct4_;
             dimensionedScalar fwStar_;
 
 
-        // Fields
-
-            volScalarField nuTilda_;
-
-            //- Wall distance
-            //  Note: different to wall distance in parent RASModel
-            //  which is for near-wall cells only
-            const volScalarField& y_;
-
-
     // Protected Member Functions
 
-        tmp<volScalarField> chi() const;
+        //- Shielding function
-
+        virtual tmp<volScalarField> psi
-        tmp<volScalarField> fv1(const volScalarField& chi) const;
-
-        tmp<volScalarField> fv2
         (
             const volScalarField& chi,
             const volScalarField& fv1
         ) const;
 
-        tmp<volScalarField> ft2(const volScalarField& chi) const;
+        //- LES length scale
-
+        virtual tmp<volScalarField> lengthScaleLES
-        tmp<volScalarField> Omega(const volTensorField& gradU) const;
-
-        tmp<volScalarField> Stilda
-        (
-            const volScalarField& chi,
-            const volScalarField& fv1,
-            const volScalarField& Omega,
-            const volScalarField& dTilda
-        ) const;
-
-        tmp<volScalarField> r
-        (
-            const volScalarField& nur,
-            const volScalarField& Omega,
-            const volScalarField& dTilda
-        ) const;
-
-        tmp<volScalarField> fw
-        (
-            const volScalarField& Omega,
-            const volScalarField& dTilda
-        ) const;
-
-        tmp<volScalarField> psi
         (
             const volScalarField& chi,
             const volScalarField& fv1
@@ -180,9 +130,6 @@ protected:
             const volTensorField& gradU
         ) const;
 
-        void correctNut(const volScalarField& fv1);
-        virtual void correctNut();
-
 
 public:
 
@@ -220,22 +167,8 @@ public:
         //- Re-read model coefficients if they have changed
         virtual bool read();
 
-        //- Return the effective diffusivity for nuTilda
-        tmp<volScalarField> DnuTildaEff() const;
-
-        //- Return SGS kinetic energy
-        virtual tmp<volScalarField> k() const;
-
-        tmp<volScalarField> nuTilda() const
-        {
-            return nuTilda_;
-        }
-
         //- Return the LES field indicator
         virtual tmp<volScalarField> LESRegion() const;
-
-        //- Correct nuTilda and related properties
-        virtual void correct();
 };
 
 

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -64,7 +64,7 @@ tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::ft
     const volScalarField& magGradU
 ) const
 {
-    return tanh(pow3(sqr(Ct_)*rd(this->nut_, magGradU)));
+    return tanh(pow3(sqr(Ct_)*this->r(this->nut_, magGradU, this->y_)));
 }
 
 
@@ -74,36 +74,7 @@ tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::fl
     const volScalarField& magGradU
 ) const
 {
-    return tanh(pow(sqr(Cl_)*rd(this->nu(), magGradU), 10));
+    return tanh(pow(sqr(Cl_)*this->r(this->nu(), magGradU, this->y_), 10));
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::rd
-(
-    const volScalarField& nur,
-    const volScalarField& magGradU
-) const
-{
-    tmp<volScalarField> tr
-    (
-        min
-        (
-            nur
-           /(
-                max
-                (
-                    magGradU,
-                    dimensionedScalar("SMALL", magGradU.dimensions(), SMALL)
-                )
-               *sqr(this->kappa_*this->y_)
-            ),
-            scalar(10)
-        )
-    );
-    tr.ref().boundaryFieldRef() == 0.0;
-
-    return tr;
 }
 
 
@@ -113,7 +84,7 @@ tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::fdt
     const volScalarField& magGradU
 ) const
 {
-    return 1 - tanh(pow(Cdt1_*rd(this->nut_, magGradU), Cdt2_));
+    return 1 - tanh(pow(Cdt1_*this->r(this->nut_, magGradU, this->y_), Cdt2_));
 }
 
 
@@ -130,7 +101,7 @@ tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::dTilda
     const volScalarField magGradU(mag(gradU));
     const volScalarField psi(this->psi(chi, fv1));
 
-    const volScalarField& lRAS(this->y_);
+    const volScalarField& lRAS = this->y_;
     const volScalarField lLES(psi*this->CDES_*this->delta());
 
     const volScalarField alpha(this->alpha());

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,8 +47,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SpalartAllmarasIDDES_H
+#ifndef Foam_SpalartAllmarasIDDES_H
-#define SpalartAllmarasIDDES_H
+#define Foam_SpalartAllmarasIDDES_H
 
 #include "SpalartAllmarasDES.H"
 #include "IDDESDelta.H"
@@ -75,14 +75,10 @@ class SpalartAllmarasIDDES
         const IDDESDelta& setDelta() const;
 
         tmp<volScalarField> alpha() const;
-        tmp<volScalarField> ft(const volScalarField& magGradU) const;
-        tmp<volScalarField> fl(const volScalarField& magGradU) const;
 
-        tmp<volScalarField> rd
+        tmp<volScalarField> ft(const volScalarField& magGradU) const;
-        (
+
-            const volScalarField& nur,
+        tmp<volScalarField> fl(const volScalarField& magGradU) const;
-            const volScalarField& magGradU
-        ) const;
 
         //- Delay function
         tmp<volScalarField> fdt(const volScalarField& magGradU) const;
@@ -105,10 +101,11 @@ protected:
             dimensionedScalar Cl_;
             dimensionedScalar Ct_;
 
-        // Fields
 
+        //- IDDES delta
         const IDDESDelta& IDDESDelta_;
 
+
     // Protected Member Functions
 
         //- Length scale

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTIDDES/kOmegaSSTIDDES.C

@@ -107,8 +107,6 @@ tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::rd
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::fdt
 (
@@ -119,6 +117,8 @@ tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::fdt
 }
 
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::dTilda
 (

src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.C

@@ -27,9 +27,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "SpalartAllmaras.H"
-#include "fvOptions.H"
-#include "bound.H"
-#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,101 +38,14 @@ namespace RASModels
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::chi() const
+Foam::tmp<Foam::volScalarField> SpalartAllmaras<BasicTurbulenceModel>::dTilda
-{
-    return nuTilda_/this->nu();
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::fv1
 (
-    const volScalarField& chi
+    const volScalarField& chi,
+    const volScalarField& fv1,
+    const volTensorField& gradU
 ) const
 {
-    const volScalarField chi3(pow3(chi));
+    return this->y_;
-
-    return chi3/(chi3 + pow3(Cv1_));
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::fv2
-(
-    const volScalarField::Internal& chi,
-    const volScalarField::Internal& fv1
-) const
-{
-    return scalar(1) - chi/(scalar(1) + chi*fv1);
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::Stilda()
-const
-{
-    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_),
-            Cs_*Omega
-        )
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::fw
-(
-    const volScalarField::Internal& Stilda
-) const
-{
-    const volScalarField::Internal r
-    (
-        min
-        (
-            nuTilda_
-           /(
-               max
-               (
-                   Stilda,
-                   dimensionedScalar(Stilda.dimensions(), SMALL)
-               )
-              *sqr(kappa_*y_)
-            ),
-            scalar(10)
-        )
-    );
-
-    const volScalarField::Internal g(r + Cw2_*(pow6(r) - r));
-
-    return
-        g*pow
-        (
-            (scalar(1) + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)),
-            scalar(1)/scalar(6)
-        );
-}
-
-
-template<class BasicTurbulenceModel>
-void SpalartAllmaras<BasicTurbulenceModel>::correctNut()
-{
-    this->nut_ = nuTilda_*this->fv1(this->chi());
-    this->nut_.correctBoundaryConditions();
-    fv::options::New(this->mesh_).correct(this->nut_);
-
-    BasicTurbulenceModel::correctNut();
 }
 
 
@@ -154,7 +64,7 @@ SpalartAllmaras<BasicTurbulenceModel>::SpalartAllmaras
     const word& type
 )
 :
-    eddyViscosity<RASModel<BasicTurbulenceModel>>
+    SpalartAllmarasBase<eddyViscosity<RASModel<BasicTurbulenceModel>>>
     (
         type,
         alpha,
@@ -164,97 +74,7 @@ SpalartAllmaras<BasicTurbulenceModel>::SpalartAllmaras
         phi,
         transport,
         propertiesName
-    ),
-
-    sigmaNut_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "sigmaNut",
-            this->coeffDict_,
-            scalar(2)/scalar(3)
     )
-    ),
-    kappa_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "kappa",
-            this->coeffDict_,
-            0.41
-        )
-    ),
-
-    Cb1_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cb1",
-            this->coeffDict_,
-            0.1355
-        )
-    ),
-    Cb2_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cb2",
-            this->coeffDict_,
-            0.622
-        )
-    ),
-    Cw1_(Cb1_/sqr(kappa_) + (scalar(1) + Cb2_)/sigmaNut_),
-    Cw2_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cw2",
-            this->coeffDict_,
-            0.3
-        )
-    ),
-    Cw3_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cw3",
-            this->coeffDict_,
-            2.0
-        )
-    ),
-    Cv1_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cv1",
-            this->coeffDict_,
-            7.1
-        )
-    ),
-    Cs_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cs",
-            this->coeffDict_,
-            0.3
-        )
-    ),
-
-    nuTilda_
-    (
-        IOobject
-        (
-            "nuTilda",
-            this->runTime_.timeName(),
-            this->mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        this->mesh_
-    ),
-
-    y_(wallDist::New(this->mesh_).y())
 {
     if (type == typeName)
     {
@@ -263,153 +83,6 @@ SpalartAllmaras<BasicTurbulenceModel>::SpalartAllmaras
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class BasicTurbulenceModel>
-bool SpalartAllmaras<BasicTurbulenceModel>::read()
-{
-    if (eddyViscosity<RASModel<BasicTurbulenceModel>>::read())
-    {
-        sigmaNut_.readIfPresent(this->coeffDict());
-        kappa_.readIfPresent(this->coeffDict());
-
-        Cb1_.readIfPresent(this->coeffDict());
-        Cb2_.readIfPresent(this->coeffDict());
-        Cw1_ = Cb1_/sqr(kappa_) + (scalar(1) + Cb2_)/sigmaNut_;
-        Cw2_.readIfPresent(this->coeffDict());
-        Cw3_.readIfPresent(this->coeffDict());
-        Cv1_.readIfPresent(this->coeffDict());
-        Cs_.readIfPresent(this->coeffDict());
-
-        return true;
-    }
-
-    return false;
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::DnuTildaEff() const
-{
-    return tmp<volScalarField>::New
-    (
-        "DnuTildaEff",
-        (nuTilda_ + this->nu())/sigmaNut_
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::k() const
-{
-    // (B:Eq. 4.50)
-    const scalar Cmu = 0.09;
-
-    return tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            IOobject::groupName("k", this->alphaRhoPhi_.group()),
-            this->runTime_.timeName(),
-            this->mesh_
-        ),
-        cbrt(this->fv1(this->chi()))
-        *nuTilda_
-        *::sqrt(scalar(2)/Cmu)
-        *mag(symm(fvc::grad(this->U_))),
-        this->nut_.boundaryField().types()
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::epsilon() const
-{
-    // (B:Eq. 4.50)
-    const scalar Cmu = 0.09;
-    const dimensionedScalar nutSMALL(sqr(dimLength)/dimTime, SMALL);
-
-    return tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            IOobject::groupName("epsilon", this->alphaRhoPhi_.group()),
-            this->runTime_.timeName(),
-            this->mesh_
-        ),
-        pow(this->fv1(this->chi()), 0.5)
-        *pow(::sqrt(Cmu)*this->k(), 2)
-        /(nuTilda_ + this->nut_ + nutSMALL),
-        this->nut_.boundaryField().types()
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::omega() const
-{
-    // (P:p. 384)
-    const scalar betaStar = 0.09;
-    const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL);
-
-    return tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            IOobject::groupName("omega", this->alphaRhoPhi_.group()),
-            this->runTime_.timeName(),
-            this->mesh_
-        ),
-        this->epsilon()/(betaStar*(this->k() + k0)),
-        this->nut_.boundaryField().types()
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-void SpalartAllmaras<BasicTurbulenceModel>::correct()
-{
-    if (!this->turbulence_)
-    {
-        return;
-    }
-
-    {
-        // Construct local convenience references
-        const alphaField& alpha = this->alpha_;
-        const rhoField& rho = this->rho_;
-        const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
-        fv::options& fvOptions(fv::options::New(this->mesh_));
-
-        eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
-
-        const volScalarField::Internal Stilda(this->Stilda());
-
-        tmp<fvScalarMatrix> nuTildaEqn
-        (
-            fvm::ddt(alpha, rho, nuTilda_)
-          + fvm::div(alphaRhoPhi, nuTilda_)
-          - 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_)
-          + fvOptions(alpha, rho, nuTilda_)
-        );
-
-        nuTildaEqn.ref().relax();
-        fvOptions.constrain(nuTildaEqn.ref());
-        solve(nuTildaEqn);
-        fvOptions.correct(nuTilda_);
-        bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), Zero));
-        nuTilda_.correctBoundaryConditions();
-    }
-
-    // Update nut with latest available nuTilda
-    correctNut();
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace RASModels

src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.H

@@ -96,8 +96,6 @@ 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
@@ -112,11 +110,12 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SpalartAllmaras_H
+#ifndef Foam_SpalartAllmaras_H
-#define SpalartAllmaras_H
+#define Foam_SpalartAllmaras_H
 
 #include "RASModel.H"
 #include "eddyViscosity.H"
+#include "SpalartAllmarasBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -132,7 +131,7 @@ namespace RASModels
 template<class BasicTurbulenceModel>
 class SpalartAllmaras
 :
-    public eddyViscosity<RASModel<BasicTurbulenceModel>>
+    public SpalartAllmarasBase<eddyViscosity<RASModel<BasicTurbulenceModel>>>
 {
     // Private Member Functions
 
@@ -145,55 +144,16 @@ class SpalartAllmaras
 
 protected:
 
-    // Protected Data
-
-        // Model coefficients
-
-            dimensionedScalar sigmaNut_;
-            dimensionedScalar kappa_;
-            dimensionedScalar Cb1_;
-            dimensionedScalar Cb2_;
-            dimensionedScalar Cw1_;
-            dimensionedScalar Cw2_;
-            dimensionedScalar Cw3_;
-            dimensionedScalar Cv1_;
-            dimensionedScalar Cs_;
-
-
-        // 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::Internal& y_;
-
-
     // Protected Member Functions
 
-        tmp<volScalarField> chi() const;
+        //- Length scale
-
+        virtual tmp<volScalarField> dTilda
-        tmp<volScalarField> fv1(const volScalarField& chi) const;
-
-        tmp<volScalarField::Internal> fv2
         (
-            const volScalarField::Internal& chi,
+            const volScalarField& chi,
-            const volScalarField::Internal& fv1
+            const volScalarField& fv1,
+            const volTensorField& gradU
         ) const;
 
-        tmp<volScalarField::Internal> Stilda() const;
-
-        tmp<volScalarField::Internal> fw
-        (
-            const volScalarField::Internal& Stilda
-        ) const;
-
-
-        //- Update nut with the latest available nuTilda
-        virtual void correctNut();
-
 
 public:
 
@@ -224,27 +184,6 @@ public:
 
     //- Destructor
     virtual ~SpalartAllmaras() = default;
-
-
-    // Member Functions
-
-        //- Re-read model coefficients if they have changed
-        virtual bool read();
-
-        //- Return the effective diffusivity for nuTilda
-        tmp<volScalarField> DnuTildaEff() const;
-
-        //- Return the (estimated) turbulent kinetic energy
-        virtual tmp<volScalarField> k() const;
-
-        //- 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 turbulent viscosity
-        virtual void correct();
 };