kOmega2006: New RAS turbulence model

This is the 2006 version of Wilcox's k-omega RAS turbulence model which has some
similarities in formulation and behaviour to the k-omega-SST model but is much
simpler and cleaner.  This model is likely to perform just as well as the
k-omega-SST model for a wide range of engineering cases.

Description
    Standard (2006) high Reynolds-number k-omega turbulence model for
    incompressible and compressible flows.

    References:
    \verbatim
        Wilcox, D. C. (2006).
        Turbulence modeling for CFD, 3rd edition
        La Canada, CA: DCW industries, Inc.

        Wilcox, D. C. (2008).
        Formulation of the kw turbulence model revisited.
        AIAA journal, 46(11), 2823-2838.
    \endverbatim

    The default model coefficients are
    \verbatim
        kOmega2006Coeffs
        {
            Cmu         0.09;
            beta0       0.0708;
            gamma       0.52;
            Clim        0.875;
            alphak      0.6;
            alphaOmega  0.5;
        }
    \endverbatim

src/MomentumTransportModels/compressible/dynamicMomentumTransportModels/dynamicMomentumTransportModels.C

@@ -85,6 +85,9 @@ makeRASModel(LaunderSharmaKE);
 #include "kOmega.H"
 makeRASModel(kOmega);
 
+#include "kOmega2006.H"
+makeRASModel(kOmega2006);
+
 #include "kOmegaSST.H"
 makeRASModel(kOmegaSST);
 

src/MomentumTransportModels/incompressible/kinematicMomentumTransportModels/kinematicMomentumTransportModels.C

@@ -82,6 +82,9 @@ makeRASModel(LaunderSharmaKE);
 #include "kOmega.H"
 makeRASModel(kOmega);
 
+#include "kOmega2006.H"
+makeRASModel(kOmega2006);
+
 #include "kOmegaSST.H"
 makeRASModel(kOmegaSST);
 

src/MomentumTransportModels/momentumTransportModels/RAS/kOmega2006/kOmega2006.C

@@ -0,0 +1,363 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "kOmega2006.H"
+#include "fvModels.H"
+#include "fvConstraints.H"
+#include "bound.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+template<class BasicMomentumTransportModel>
+void kOmega2006<BasicMomentumTransportModel>::correctNut
+(
+    const volTensorField& gradU
+)
+{
+    this->nut_ = k_/max(omega_, Clim_*sqrt(2/betaStar_)*mag(dev(symm(gradU))));
+    this->nut_.correctBoundaryConditions();
+    fvConstraints::New(this->mesh_).constrain(this->nut_);
+}
+
+
+template<class BasicMomentumTransportModel>
+tmp<volScalarField::Internal> kOmega2006<BasicMomentumTransportModel>::beta
+(
+    const volTensorField& gradU
+) const
+{
+    const volSymmTensorField::Internal S(symm(gradU()));
+    const volSymmTensorField::Internal Shat(S - 0.5*tr(S)*I);
+    const volTensorField::Internal Omega(skew(gradU.v()));
+
+    const volScalarField::Internal ChiOmega
+    (
+        modelName("ChiOmega"),
+        mag((Omega & Omega) && Shat)/pow3(betaStar_*omega_.v())
+    );
+
+    const volScalarField::Internal fBeta
+    (
+        modelName("fBeta"),
+        (1 + 85*ChiOmega)/(1 + 100*ChiOmega)
+    );
+
+    return beta0_*fBeta;
+}
+
+
+template<class BasicMomentumTransportModel>
+tmp<volScalarField::Internal>
+kOmega2006<BasicMomentumTransportModel>::CDkOmega() const
+{
+    return max
+    (
+        sigmaDo_*(fvc::grad(k_)().v() & fvc::grad(omega_)().v())/omega_(),
+        dimensionedScalar(dimless/sqr(dimTime), 0)
+    );
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicMomentumTransportModel>
+void kOmega2006<BasicMomentumTransportModel>::correctNut()
+{
+    correctNut(fvc::grad(this->U_));
+}
+
+
+template<class BasicMomentumTransportModel>
+tmp<fvScalarMatrix> kOmega2006<BasicMomentumTransportModel>::kSource() const
+{
+    return tmp<fvScalarMatrix>
+    (
+        new fvScalarMatrix
+        (
+            k_,
+            dimVolume*this->rho_.dimensions()*k_.dimensions()
+            /dimTime
+        )
+    );
+}
+
+
+template<class BasicMomentumTransportModel>
+tmp<fvScalarMatrix> kOmega2006<BasicMomentumTransportModel>::omegaSource() const
+{
+    return tmp<fvScalarMatrix>
+    (
+        new fvScalarMatrix
+        (
+            omega_,
+            dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicMomentumTransportModel>
+kOmega2006<BasicMomentumTransportModel>::kOmega2006
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& type
+)
+:
+    eddyViscosity<RASModel<BasicMomentumTransportModel>>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport
+    ),
+
+    betaStar_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "betaStar",
+            this->coeffDict_,
+            0.09
+        )
+    ),
+    beta0_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "beta0",
+            this->coeffDict_,
+            0.0708
+        )
+    ),
+    gamma_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "gamma",
+            this->coeffDict_,
+            0.52
+        )
+    ),
+    Clim_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Clim",
+            this->coeffDict_,
+            0.875
+        )
+    ),
+    sigmaDo_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "sigmaDo",
+            this->coeffDict_,
+            0.125
+        )
+    ),
+    alphaK_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "alphaK",
+            this->coeffDict_,
+            0.6
+        )
+    ),
+    alphaOmega_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "alphaOmega",
+            this->coeffDict_,
+            0.5
+        )
+    ),
+
+    k_
+    (
+        IOobject
+        (
+            IOobject::groupName("k", alphaRhoPhi.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+    omega_
+    (
+        IOobject
+        (
+            IOobject::groupName("omega", alphaRhoPhi.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    )
+{
+    bound(k_, this->kMin_);
+    bound(omega_, this->omegaMin_);
+
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicMomentumTransportModel>
+bool kOmega2006<BasicMomentumTransportModel>::read()
+{
+    if (eddyViscosity<RASModel<BasicMomentumTransportModel>>::read())
+    {
+        betaStar_.readIfPresent(this->coeffDict());
+        beta0_.readIfPresent(this->coeffDict());
+        gamma_.readIfPresent(this->coeffDict());
+        sigmaDo_.readIfPresent(this->coeffDict());
+        alphaK_.readIfPresent(this->coeffDict());
+        alphaOmega_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicMomentumTransportModel>
+void kOmega2006<BasicMomentumTransportModel>::correct()
+{
+    if (!this->turbulence_)
+    {
+        return;
+    }
+
+    // Local references
+    const alphaField& alpha = this->alpha_;
+    const rhoField& rho = this->rho_;
+    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
+    const volVectorField& U = this->U_;
+    volScalarField& nut = this->nut_;
+    const Foam::fvModels& fvModels(Foam::fvModels::New(this->mesh_));
+    const Foam::fvConstraints& fvConstraints
+    (
+        Foam::fvConstraints::New(this->mesh_)
+    );
+
+    eddyViscosity<RASModel<BasicMomentumTransportModel>>::correct();
+
+    volScalarField::Internal divU
+    (
+        modelName("divU"),
+        fvc::div(fvc::absolute(this->phi(), U))().v()
+    );
+
+    const volTensorField gradU(fvc::grad(U));
+    volScalarField::Internal G
+    (
+        this->GName(),
+        nut.v()*(dev(twoSymm(gradU.v())) && gradU.v())
+    );
+
+    // Update omega and G at the wall
+    omega_.boundaryFieldRef().updateCoeffs();
+
+    // Turbulence specific dissipation rate equation
+    tmp<fvScalarMatrix> omegaEqn
+    (
+        fvm::ddt(alpha, rho, omega_)
+      + fvm::div(alphaRhoPhi, omega_)
+      - fvm::laplacian(alpha*rho*DomegaEff(), omega_)
+     ==
+        gamma_*alpha()*rho()*G*omega_()/k_()
+      - fvm::SuSp(((2.0/3.0)*gamma_)*alpha()*rho()*divU, omega_)
+      - fvm::Sp(beta(gradU)*alpha()*rho()*omega_(), omega_)
+      + alpha()*rho()*CDkOmega()
+      + omegaSource()
+      + fvModels.source(alpha, rho, omega_)
+    );
+
+    omegaEqn.ref().relax();
+    fvConstraints.constrain(omegaEqn.ref());
+    omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
+    solve(omegaEqn);
+    fvConstraints.constrain(omega_);
+    bound(omega_, this->omegaMin_);
+
+
+    // Turbulent kinetic energy equation
+    tmp<fvScalarMatrix> kEqn
+    (
+        fvm::ddt(alpha, rho, k_)
+      + fvm::div(alphaRhoPhi, k_)
+      - fvm::laplacian(alpha*rho*DkEff(), k_)
+     ==
+        alpha()*rho()*G
+      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
+      - fvm::Sp(betaStar_*alpha()*rho()*omega_(), k_)
+      + kSource()
+      + fvModels.source(alpha, rho, k_)
+    );
+
+    kEqn.ref().relax();
+    fvConstraints.constrain(kEqn.ref());
+    solve(kEqn);
+    fvConstraints.constrain(k_);
+    bound(k_, this->kMin_);
+
+    correctNut(gradU);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/MomentumTransportModels/momentumTransportModels/RAS/kOmega2006/kOmega2006.H

@@ -0,0 +1,219 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::RASModels::kOmega2006
+
+Description
+    Standard (2006) high Reynolds-number k-omega turbulence model for
+    incompressible and compressible flows.
+
+    References:
+    \verbatim
+        Wilcox, D. C. (2006).
+        Turbulence modeling for CFD, 3rd edition
+        La Canada, CA: DCW industries, Inc.
+
+        Wilcox, D. C. (2008).
+        Formulation of the kw turbulence model revisited.
+        AIAA journal, 46(11), 2823-2838.
+    \endverbatim
+
+    The default model coefficients are
+    \verbatim
+        kOmega2006Coeffs
+        {
+            Cmu         0.09;
+            beta0       0.0708;
+            gamma       0.52;
+            Clim        0.875;
+            alphak      0.6;
+            alphaOmega  0.5;
+        }
+    \endverbatim
+
+SourceFiles
+    kOmega2006.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kOmega2006_H
+#define kOmega2006_H
+
+#include "RASModel.H"
+#include "eddyViscosity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class kOmega2006 Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicMomentumTransportModel>
+class kOmega2006
+:
+    public eddyViscosity<RASModel<BasicMomentumTransportModel>>
+{
+    // Private Member Functions
+
+        using IOobject::modelName;
+
+        void correctNut(const volTensorField& gradU);
+
+        tmp<volScalarField::Internal> beta(const volTensorField& gradU) const;
+
+        tmp<volScalarField::Internal> CDkOmega() const;
+
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar betaStar_;
+            dimensionedScalar beta0_;
+            dimensionedScalar gamma_;
+            dimensionedScalar Clim_;
+            dimensionedScalar sigmaDo_;
+            dimensionedScalar alphaK_;
+            dimensionedScalar alphaOmega_;
+
+
+        // Fields
+
+            volScalarField k_;
+            volScalarField omega_;
+
+
+    // Protected Member Functions
+
+        virtual void correctNut();
+        virtual tmp<fvScalarMatrix> kSource() const;
+        virtual tmp<fvScalarMatrix> omegaSource() const;
+
+
+public:
+
+    typedef typename BasicMomentumTransportModel::alphaField alphaField;
+    typedef typename BasicMomentumTransportModel::rhoField rhoField;
+    typedef typename BasicMomentumTransportModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("kOmega2006");
+
+
+    // Constructors
+
+        //- Construct from components
+        kOmega2006
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~kOmega2006()
+    {}
+
+
+    // Member Functions
+
+        //- Read RASProperties dictionary
+        virtual bool read();
+
+        //- Return the effective diffusivity for k
+        tmp<volScalarField> DkEff() const
+        {
+            return volScalarField::New
+            (
+                "DkEff",
+                alphaK_*this->nut_ + this->nu()
+            );
+        }
+
+        //- Return the effective diffusivity for omega
+        tmp<volScalarField> DomegaEff() const
+        {
+            return volScalarField::New
+            (
+                "DomegaEff",
+                alphaOmega_*this->nut_ + this->nu()
+            );
+        }
+
+        //- Return the turbulence kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return k_;
+        }
+
+        //- Return the turbulence specific dissipation rate
+        virtual tmp<volScalarField> omega() const
+        {
+            return omega_;
+        }
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const
+        {
+            return volScalarField::New
+            (
+                "epsilon",
+                betaStar_*k_*omega_,
+                omega_.boundaryField().types()
+            );
+        }
+
+        //- Solve the turbulence equations and correct the turbulence viscosity
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#ifdef NoRepository
+    #include "kOmega2006.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/MomentumTransportModels/momentumTransportModels/RAS/realizableKE/realizableKE.C

@@ -48,7 +48,7 @@ tmp<volScalarField> realizableKE<BasicMomentumTransportModel>::rCmu
     tmp<volSymmTensorField> tS = dev(symm(gradU));
     const volSymmTensorField& S = tS();
 
-    volScalarField W
+    const volScalarField W
     (
         (2*sqrt(2.0))*((S&S)&&S)
        /(
@@ -59,12 +59,12 @@ tmp<volScalarField> realizableKE<BasicMomentumTransportModel>::rCmu
 
     tS.clear();
 
-    volScalarField phis
+    const volScalarField phis
     (
         (1.0/3.0)*acos(min(max(sqrt(6.0)*W, -scalar(1)), scalar(1)))
     );
-    volScalarField As(sqrt(6.0)*cos(phis));
+    const volScalarField As(sqrt(6.0)*cos(phis));
-    volScalarField Us(sqrt(S2/2.0 + magSqr(skew(gradU))));
+    const volScalarField Us(sqrt(S2/2.0 + magSqr(skew(gradU))));
 
     return 1.0/(A0_ + As*Us*k_/epsilon_);
 }
@@ -87,11 +87,11 @@ void realizableKE<BasicMomentumTransportModel>::correctNut
 template<class BasicMomentumTransportModel>
 void realizableKE<BasicMomentumTransportModel>::correctNut()
 {
-    tmp<volTensorField> tgradU = fvc::grad(this->U_);
+    const volTensorField gradU(fvc::grad(this->U_));
-    volScalarField S2(modelName("S2"), 2*magSqr(dev(symm(tgradU()))));
+    const volScalarField S2(modelName("S2"), 2*magSqr(dev(symm(gradU))));
-    volScalarField magS(modelName("magS"), sqrt(S2));
+    const volScalarField magS(modelName("magS"), sqrt(S2));
 
-    correctNut(tgradU(), S2, magS);
+    correctNut(gradU, S2, magS);
 }
 
 
@@ -271,21 +271,24 @@ void realizableKE<BasicMomentumTransportModel>::correct()
         fvc::div(fvc::absolute(this->phi(), U))()
     );
 
-    tmp<volTensorField> tgradU = fvc::grad(U);
+    const volTensorField gradU(fvc::grad(U));
-    volScalarField S2(modelName("S2"), 2*magSqr(dev(symm(tgradU()))));
+    const volScalarField S2(modelName("S2"), 2*magSqr(dev(symm(gradU))));
-    volScalarField magS(modelName("magS"), sqrt(S2));
+    const volScalarField magS(modelName("magS"), sqrt(S2));
 
-    volScalarField::Internal eta(modelName("eta"), magS()*k_()/epsilon_());
+    const volScalarField::Internal eta
-    volScalarField::Internal C1
+    (
+        modelName("eta"), magS()*k_()/epsilon_()
+    );
+    const volScalarField::Internal C1
     (
         modelName("C1"),
         max(eta/(scalar(5) + eta), scalar(0.43))
     );
 
-    volScalarField::Internal G
+    const volScalarField::Internal G
     (
         this->GName(),
-        nut*(tgradU().v() && dev(twoSymm(tgradU().v())))
+        nut*(gradU.v() && dev(twoSymm(gradU.v())))
     );
 
     // Update epsilon and G at the wall
@@ -337,7 +340,7 @@ void realizableKE<BasicMomentumTransportModel>::correct()
     fvConstraints.constrain(k_);
     bound(k_, this->kMin_);
 
-    correctNut(tgradU(), S2, magS);
+    correctNut(gradU, S2, magS);
 }
 
 

tutorials/incompressible/simpleFoam/pitzDaily/constant/momentumTransport

@@ -18,7 +18,7 @@ simulationType RAS;
 
 RAS
 {
-    // Tested with kEpsilon, realizableKE, kOmega, kOmegaSST, v2f,
+    // Tested with kEpsilon, realizableKE, kOmega, kOmega2006, kOmegaSST, v2f,
     // ShihQuadraticKE, LienCubicKE.
     model           kEpsilon;