Spalart-Shur Rotation/Curvature Correction for Spalart Allmaras model

Spalart-Shur Rotation/Curvature Correction[1] was employed for incompressible SpalartAllmaras model in OpenFOAM.

New source codes are shown below. Differences from the original ones are indicated in red.

References
1. M.L. Shur, M.K. Strelets, A.K. Travin, P.R. Spalart, Turbulence Modeling in Rotating and Curved Channels: Assessing the Spalart-Shur Correction, AIAA Journal Vol. 38, No. 5, 2000, pp. 784-792.

SpalartAllmaras.H
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 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::incompressible::RASModels::SpalartAllmaras

Group
    grpIcoRASTurbulence

Description
    Spalart-Allmaras 1-eqn mixing-length model for incompressible external
    flows.

    References:
    \verbatim
        "A One-Equation Turbulence Model for Aerodynamic Flows"
        P.R. Spalart,
        S.R. Allmaras,
        La Recherche Aerospatiale, No. 1, 1994, pp. 5-21.
    \endverbatim

    Extended according to
    \verbatim
        "An Unstructured Grid Generation and Adaptive Solution Technique
        for High Reynolds Number Compressible Flows"
        G.A. Ashford,
        Ph.D. thesis, University of Michigan, 1996.
    \endverbatim
    using the optional flag \c ashfordCorrection,

    and
    \verbatim
        "Turbulence Modeling in Rotating and Curved Channels: Assessing the
        Spalart-Shur Correction"
        M.L. Shur,
        M.K. Strelets,
        A.K. Travin,
        P.R. Spalart,
        AIAA Journal Vol. 38, No. 5, 2000, pp. 784-792.
    \endverbatim
    using the optional flag \c spalartShurCorrection

    The default model coefficients correspond to the following:
    \verbatim
        SpalartAllmarasCoeffs
        {
            Cb1         0.1355;
            Cb2         0.622;
            Cw2         0.3;
            Cw3         2.0;
            Cv1         7.1;
            Cv2         5.0;
            sigmaNut    0.66666;
            kappa       0.41;
            Cr1         1.0;
            Cr2         12.0;
            Cr3         1.0;
        }
    \endverbatim

SourceFiles
    SpalartAllmaras.C

\*---------------------------------------------------------------------------*/

#ifndef SpalartAllmaras_H
#define SpalartAllmaras_H

#include "RASModel.H"
#include "wallDist.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{
namespace incompressible
{
namespace RASModels
{

/*---------------------------------------------------------------------------*\
                       Class SpalartAllmaras Declaration
\*---------------------------------------------------------------------------*/

class SpalartAllmaras
:
    public RASModel
{

protected:

    // Protected data

        // Model coefficients

            dimensionedScalar sigmaNut_;
            dimensionedScalar kappa_;

            dimensionedScalar Cb1_;
            dimensionedScalar Cb2_;
            dimensionedScalar Cw1_;
            dimensionedScalar Cw2_;
            dimensionedScalar Cw3_;
            dimensionedScalar Cv1_;
            dimensionedScalar Cv2_;
            dimensionedScalar Cr1_;
            dimensionedScalar Cr2_;
            dimensionedScalar Cr3_;


        //- Optional flag to activate the Ashford correction
        Switch ashfordCorrection_;

        //- Optional flag to activate the Rotation/Curvature correction
        Switch spalartShurCorrection_;

        // Fields

            volScalarField nuTilda_;
            volScalarField nut_;
            volScalarField Sxx_;
            volScalarField Sxy_;
            volScalarField Sxz_;
            volScalarField Syx_;
            volScalarField Syy_;
            volScalarField Syz_;
            volScalarField Szx_;
            volScalarField Szy_;
            volScalarField Szz_;

            wallDist d_;


    // 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> fv3
        (
            const volScalarField& chi,
            const volScalarField& fv1
        ) const;

        tmp<volScalarField> fw(const volScalarField& Stilda) const;

        tmp<volScalarField> fr1(const volSymmTensorField& S, const volTensorField& W) const;


public:

    //- Runtime type information
    TypeName("SpalartAllmaras");


    // Constructors

        //- Construct from components
        SpalartAllmaras
        (
            const volVectorField& U,
            const surfaceScalarField& phi,
            transportModel& transport,
            const word& turbulenceModelName = turbulenceModel::typeName,
            const word& modelName = typeName
        );


    //- Destructor
    virtual ~SpalartAllmaras()
    {}


    // Member Functions

        //- Return the turbulence viscosity
        virtual tmp<volScalarField> nut() const
        {
            return nut_;
        }

        //- Return the effective diffusivity for nuTilda
        tmp<volScalarField> DnuTildaEff() const;

        //- Return the turbulence kinetic energy
        virtual tmp<volScalarField> k() const;

        //- Return the turbulence kinetic energy dissipation rate
        virtual tmp<volScalarField> epsilon() const;

        //- Return the Reynolds stress tensor
        virtual tmp<volSymmTensorField> R() const;

        //- Return the effective stress tensor including the laminar stress
        virtual tmp<volSymmTensorField> devReff() const;

        //- Return the source term for the momentum equation
        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;

        //- Return the source term for the momentum equation
        virtual tmp<fvVectorMatrix> divDevRhoReff
        (
            const volScalarField& rho,
            volVectorField& U
        ) const;

        //- Solve the turbulence equations and correct the turbulence viscosity
        virtual void correct();

        //- Read RASProperties dictionary
        virtual bool read();
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //


SpalartAllmaras.C
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 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 "SpalartAllmaras.H"
#include "addToRunTimeSelectionTable.H"
#include "fvcDdt.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{
namespace incompressible
{
namespace RASModels
{

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

defineTypeNameAndDebug(SpalartAllmaras, 0);
addToRunTimeSelectionTable(RASModel, SpalartAllmaras, dictionary);

// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //

tmp<volScalarField> SpalartAllmaras::chi() const
{
    return nuTilda_/nu();
}


tmp<volScalarField> SpalartAllmaras::fv1(const volScalarField& chi) const
{
    const volScalarField chi3(pow3(chi));
    return chi3/(chi3 + pow3(Cv1_));
}


tmp<volScalarField> SpalartAllmaras::fv2
(
    const volScalarField& chi,
    const volScalarField& fv1
) const
{
    if (ashfordCorrection_)
    {
        return 1.0/pow3(scalar(1) + chi/Cv2_);
    }
    else
    {
        return 1.0 - chi/(1.0 + chi*fv1);
    }
}


tmp<volScalarField> SpalartAllmaras::fv3
(
    const volScalarField& chi,
    const volScalarField& fv1
) const
{
    if (ashfordCorrection_)
    {
        const volScalarField chiByCv2((1/Cv2_)*chi);

        return
            (scalar(1) + chi*fv1)
           *(1/Cv2_)
           *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
           /pow3(scalar(1) + chiByCv2);
    }
    else
    {
        return tmp<volScalarField>
        (
            new volScalarField
            (
                IOobject
                (
                    "fv3",
                    mesh_.time().timeName(),
                    mesh_,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                mesh_,
                dimensionedScalar("fv3", dimless, 1),
                zeroGradientFvPatchScalarField::typeName
            )
        );
    }
}


tmp<volScalarField> SpalartAllmaras::fw(const volScalarField& Stilda) const
{
    volScalarField r
    (
        min
        (
            nuTilda_
           /(
               max
               (
                   Stilda,
                   dimensionedScalar("SMALL", Stilda.dimensions(), SMALL)
               )
              *sqr(kappa_*d_)
            ),
            scalar(10.0)
        )
    );
    r.boundaryField() == 0.0;

    const volScalarField g(r + Cw2_*(pow6(r) - r));

    return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
}

tmp<volScalarField> SpalartAllmaras::fr1(const volSymmTensorField& S, const volTensorField& W) const
{
    if (spalartShurCorrection_)
    {
        volScalarField sqrS(2.0*magSqr(S));
        volScalarField sqrW(2.0*magSqr(W));
        volScalarField sqrD(0.5*(sqrS + sqrW));

        volScalarField rStar(sqrt(sqrS/sqrW));

        volScalarField Wxx(W.component(tensor::XX));
        volScalarField Wxy(W.component(tensor::XY));
        volScalarField Wxz(W.component(tensor::XZ));
        volScalarField Wyx(W.component(tensor::YX));
        volScalarField Wyy(W.component(tensor::YY));
        volScalarField Wyz(W.component(tensor::YZ));
        volScalarField Wzx(W.component(tensor::ZX));
        volScalarField Wzy(W.component(tensor::ZY));
        volScalarField Wzz(W.component(tensor::ZZ));
        
        volScalarField rTilda
        (
            2.0/sqr(sqrD)
               *(
                   (Wxx*Sxx_ + Wxy*Sxy_ + Wxz*Sxz_)
                      *(fvc::ddt(Sxx_) + fvc::div(phi_, Sxx_)) //i,j=1,1
                 + (Wxx*Syx_ + Wxy*Syy_ + Wxz*Syz_)
                      *(fvc::ddt(Sxy_) + fvc::div(phi_, Sxy_)) //i,j=1,2
                 + (Wxx*Szx_ + Wxy*Szy_ + Wxz*Szz_)
                      *(fvc::ddt(Sxz_) + fvc::div(phi_, Sxz_)) //i,j=1,3
                 + (Wyx*Sxx_ + Wyy*Sxy_ + Wyz*Sxz_)
                      *(fvc::ddt(Syx_) + fvc::div(phi_, Syx_)) //i,j=2,1
                 + (Wyx*Syx_ + Wyy*Syy_ + Wyz*Syz_)
                      *(fvc::ddt(Syy_) + fvc::div(phi_, Syy_)) //i,j=2,2
                 + (Wyx*Szx_ + Wyy*Szy_ + Wyz*Szz_)
                      *(fvc::ddt(Syz_) + fvc::div(phi_, Syz_)) //i,j=2,3
                 + (Wzx*Sxx_ + Wzy*Sxy_ + Wzz*Sxz_)
                      *(fvc::ddt(Szx_) + fvc::div(phi_, Szx_)) //i,j=3,1
                 + (Wzx*Syx_ + Wzy*Syy_ + Wzz*Syz_)
                      *(fvc::ddt(Szy_) + fvc::div(phi_, Szy_)) //i,j=3,2
                 + (Wzx*Szx_ + Wzy*Szy_ + Wzz*Szz_)
                      *(fvc::ddt(Szz_) + fvc::div(phi_, Szz_)) //i,j=3,3
                )
        );

        return 
            (1 + Cr1_)*2.0*rStar/(1.0 + rStar)*(1.0 - Cr3_*atan(Cr2_*rTilda))
             - Cr1_;
    }
    else
    {
        return tmp<volScalarField>
        (
            new volScalarField
            (
                IOobject
                (
                    "fr1",
                    mesh_.time().timeName(),
                    mesh_,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                mesh_,
                dimensionedScalar("fr1", dimless, 1),
                zeroGradientFvPatchScalarField::typeName
            )
        );
    }
}

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

SpalartAllmaras::SpalartAllmaras
(
    const volVectorField& U,
    const surfaceScalarField& phi,
    transportModel& transport,
    const word& turbulenceModelName,
    const word& modelName
)
:
    RASModel(modelName, U, phi, transport, turbulenceModelName),

    sigmaNut_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "sigmaNut",
            coeffDict_,
            0.66666
        )
    ),
    kappa_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "kappa",
            coeffDict_,
            0.41
        )
    ),

    Cb1_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cb1",
            coeffDict_,
            0.1355
        )
    ),
    Cb2_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cb2",
            coeffDict_,
            0.622
        )
    ),
    Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
    Cw2_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cw2",
            coeffDict_,
            0.3
        )
    ),
    Cw3_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cw3",
            coeffDict_,
            2.0
        )
    ),
    Cv1_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cv1",
            coeffDict_,
            7.1
        )
    ),
    Cv2_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cv2",
            coeffDict_,
            5.0
        )
    ),
    Cr1_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cr1",
            coeffDict_,
            1.0
        )
    ),
    Cr2_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cr2",
            coeffDict_,
            12.0
        )
    ),
    Cr3_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cr3",
            coeffDict_,
            1.0
        )
    ),

    ashfordCorrection_(coeffDict_.lookupOrDefault("ashfordCorrection", true)),

    spalartShurCorrection_
    (
        coeffDict_.lookupOrDefault("spalartShurCorrection", true)
    ),

    nuTilda_
    (
        IOobject
        (
            "nuTilda",
            runTime_.timeName(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_
    ),

    nut_
    (
        IOobject
        (
            "nut",
            runTime_.timeName(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_
    ),

    Sxx_
    (
        IOobject
        (
            "Sxx",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Sxx", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
        zeroGradientFvPatchScalarField::typeName
    ),

    Sxy_
    (
        IOobject
        (
            "Sxy",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Sxy", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
        zeroGradientFvPatchScalarField::typeName
    ),

    Sxz_
    (
        IOobject
        (
            "Sxz",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Sxz", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
        zeroGradientFvPatchScalarField::typeName
    ),

    Syx_
    (
        IOobject
        (
            "Syx",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Syx", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
        zeroGradientFvPatchScalarField::typeName
    ),

    Syy_
    (
        IOobject
        (
            "Syy",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Syy", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
        zeroGradientFvPatchScalarField::typeName
    ),

    Syz_
    (
        IOobject
        (
            "Syz",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Syz", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
        zeroGradientFvPatchScalarField::typeName
    ),

    Szx_
    (
        IOobject
        (
            "Szx",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Szx", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
        zeroGradientFvPatchScalarField::typeName
    ),

    Szy_
    (
        IOobject
        (
            "Szy",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Szy", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
        zeroGradientFvPatchScalarField::typeName
    ),

    Szz_
    (
        IOobject
        (
            "Szz",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Szz", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
        zeroGradientFvPatchScalarField::typeName
    ),

    d_(mesh_)
{
    printCoeffs();

    if (ashfordCorrection_)
    {
        Info<< "    Employing Ashford correction" << endl;
    }
    if (spalartShurCorrection_)
    {
        Info<< "    Employing Rotation/Curvature correction" << endl;
    }
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

tmp<volScalarField> SpalartAllmaras::DnuTildaEff() const
{
    return tmp<volScalarField>
    (
        new volScalarField("DnuTildaEff", (nuTilda_ + nu())/sigmaNut_)
    );
}


tmp<volScalarField> SpalartAllmaras::k() const
{
    WarningIn("tmp<volScalarField> SpalartAllmaras::k() const")
        << "Turbulence kinetic energy not defined for Spalart-Allmaras model. "
        << "Returning zero field" << endl;

    return tmp<volScalarField>
    (
        new volScalarField
        (
            IOobject
            (
                "k",
                runTime_.timeName(),
                mesh_
            ),
            mesh_,
            dimensionedScalar("0", dimensionSet(0, 2, -2, 0, 0), 0)
        )
    );
}


tmp<volScalarField> SpalartAllmaras::epsilon() const
{
    WarningIn("tmp<volScalarField> SpalartAllmaras::epsilon() const")
        << "Turbulence kinetic energy dissipation rate not defined for "
        << "Spalart-Allmaras model. Returning zero field"
        << endl;

    return tmp<volScalarField>
    (
        new volScalarField
        (
            IOobject
            (
                "epsilon",
                runTime_.timeName(),
                mesh_
            ),
            mesh_,
            dimensionedScalar("0", dimensionSet(0, 2, -3, 0, 0), 0)
        )
    );
}


tmp<volSymmTensorField> SpalartAllmaras::R() const
{
    return tmp<volSymmTensorField>
    (
        new volSymmTensorField
        (
            IOobject
            (
                "R",
                runTime_.timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            ((2.0/3.0)*I)*k() - nut()*twoSymm(fvc::grad(U_))
        )
    );
}


tmp<volSymmTensorField> SpalartAllmaras::devReff() const
{
    return tmp<volSymmTensorField>
    (
        new volSymmTensorField
        (
            IOobject
            (
                "devRhoReff",
                runTime_.timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
           -nuEff()*dev(twoSymm(fvc::grad(U_)))
        )
    );
}


tmp<fvVectorMatrix> SpalartAllmaras::divDevReff(volVectorField& U) const
{
    const volScalarField nuEff_(nuEff());

    return
    (
      - fvm::laplacian(nuEff_, U)
      - fvc::div(nuEff_*dev(T(fvc::grad(U))))
    );
}


tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff
(
    const volScalarField& rho,
    volVectorField& U
) const
{
    volScalarField muEff("muEff", rho*nuEff());

    return
    (
      - fvm::laplacian(muEff, U)
      - fvc::div(muEff*dev(T(fvc::grad(U))))
    );
}


bool SpalartAllmaras::read()
{
    if (RASModel::read())
    {
        sigmaNut_.readIfPresent(coeffDict());
        kappa_.readIfPresent(coeffDict());

        Cb1_.readIfPresent(coeffDict());
        Cb2_.readIfPresent(coeffDict());
        Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
        Cw2_.readIfPresent(coeffDict());
        Cw3_.readIfPresent(coeffDict());
        Cv1_.readIfPresent(coeffDict());
        Cv2_.readIfPresent(coeffDict());
        Cr1_.readIfPresent(coeffDict());
        Cr2_.readIfPresent(coeffDict());
        Cr3_.readIfPresent(coeffDict());

        ashfordCorrection_.readIfPresent("ashfordCorrection", coeffDict());
        spalartShurCorrection_.readIfPresent
        (
            "spalartShurCorrection", coeffDict()
        );

        return true;
    }
    else
    {
        return false;
    }
}


void SpalartAllmaras::correct()
{
    RASModel::correct();

    if (!turbulence_)
    {
        // Re-calculate viscosity
        nut_ = nuTilda_*fv1(this->chi());
        nut_.correctBoundaryConditions();

        return;
    }

    if (mesh_.changing())
    {
        d_.correct();
    }

    const volScalarField chi(this->chi());
    const volScalarField fv1(this->fv1(chi));

    const volScalarField Stilda
    (
        fv3(chi, fv1)*::sqrt(2.0)*mag(skew(fvc::grad(U_)))
      + fv2(chi, fv1)*nuTilda_/sqr(kappa_*d_)
    );

    volTensorField gradU(fvc::grad(U_));
    volSymmTensorField S(symm(gradU));
    volTensorField W(skew(gradU));

    Sxx_ = S.component(tensor::XX);
    Sxy_ = S.component(tensor::XY);
    Sxz_ = S.component(tensor::XZ);
    Syx_ = S.component(tensor::YX);
    Syy_ = S.component(tensor::YY);
    Syz_ = S.component(tensor::YZ);
    Szx_ = S.component(tensor::ZX);
    Szy_ = S.component(tensor::ZY);
    Szz_ = S.component(tensor::ZZ);

    tmp<fvScalarMatrix> nuTildaEqn
    (
        fvm::ddt(nuTilda_)
      + fvm::div(phi_, nuTilda_)
      - fvm::laplacian(DnuTildaEff(), nuTilda_)
      - Cb2_/sigmaNut_*magSqr(fvc::grad(nuTilda_))
     ==
        Cb1_*Stilda*nuTilda_*fr1(S, W)
      - fvm::Sp(Cw1_*fw(Stilda)*nuTilda_/sqr(d_), nuTilda_)
    );

    nuTildaEqn().relax();
    solve(nuTildaEqn);
    bound(nuTilda_, dimensionedScalar("0", nuTilda_.dimensions(), 0.0));
    nuTilda_.correctBoundaryConditions();

    // Re-calculate viscosity
    nut_.internalField() = fv1*nuTilda_.internalField();
    nut_.correctBoundaryConditions();
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam

// ************************************************************************* //

Comments