Implement turbulence model

We will here investigate how turbulence models are incorporated in the solution procedure of the solvers, and how a particular turbulence model is chosen by the user and pointed at in run-time by the code. We will then investigate how the kEpsilon, kOmegaSST and kOmegaSSTSAS models are implemented. Finally we will implement a new version of the kOmegaSST model, named kOmegaSSTF, which puts a filter to the turbulent viscosity.

Initial notes

The implementation of turbulence models was changed completely in OpenFOAM-3.0. The current descriptions are based on the new implementation (in fact OpenFOAM-v1706+). If you have an old version of OpenFOAM, or if you are using FOAM-extend, you can have a look at the slides from previous years.

If you are only interested in the basic steps to copy an existing turbulence model and compile your own modified version of that model, please go to the section named "Implement the kOmegaSSTF model". In fact, it is good to anyway go there and start by following the instructions how to copy and compile the libraries, since that compilation takes some time.

If you try to follow the descriptions below by yourself it is recommended to use GDB and the Debug compilation to verify the exact paths that are taken while a simulation is executed, in this case the simpleFoam/pitzDaily tutorial. It is then recommended to set a breakpoint at the autoPtr line in createFields.H:

b createFields.H:40
run -noFunctionObjects

and then step (s) through the lines (entering the function calls). It is not possible to set a breakpoint in some of the classes, since simpleFoam does not know about all of those until it is running.

Basic usage of turbulence models

We will use the simpleFoam/pitzDaily tutorial when we examine the implementation of turbulence models. We start by going through the basic usage of turbulence models for that case.

Start by making sure that you have a fresh copy of the case in your $FOAM_RUN directory:

OF1706+ #Or OF1706+Debug, if you intend to use GDB
run
rm -r pitzDaily
cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily .
cd pitzDaily
blockMesh

The default settings for turbulence in the constant/turbulenceProperties file/dictionary are:

simulationType RAS;

RAS
{
    // Tested with kEpsilon, realizableKE, kOmega, kOmegaSST, v2f,
    // ShihQuadraticKE, LienCubicKE.
    RASModel        kEpsilon;

    turbulence      on;

    printCoeffs     on;
}

i.e. the kEpsilon RAS model will be used, it is active, and it will print its coefficients to the terminal window (commonly forwarded to a log file).

If you like you can run the case by issuing the command simpleFoam, and investigate the results in paraFoam. You will see that there are two fields relating to turbulence: k for turbulent kinetic energy, and epsilon for turbulence dissipation. Those fields are solved using discreticed equations, so there is a need for solver and scheme settings in system/{fvSolution,fvSchemes}. You can of course also see those files in the time directories, where internalField (initial conditions or solution) and boundary conditions are specified/shown. Here it should be noted that other turbulence models may be based on other variables, and it thus makes sense that it is the turbulence model that defines the variables to be used. The k and epsilon objects are not constructed in the createFields file of the simpleFoam solver. Instead they are constructed by the turbulence model that has been chosen at run-time. We will see this later, when we look at the code of turbulence models.

Some of the initial output in the terminal window when running simpleFoam is:

Selecting turbulence model type RAS
Selecting RAS turbulence model kEpsilon
kEpsilonCoeffs
{
    Cmu             0.09;
    C1              1.44;
    C2              1.92;
    C3              -0.33;
    sigmak          1;
    sigmaEps        1.3;
}

This shows that the kEpsilon model was in fact chosen for the simulation, and it also reports the model coefficients that are used in the simulation (since the RAS sub-dictionary above specified 'printCoeffs on').

Now we will figure out how the simpleFoam solver incorporates turbulence modeling in general, and the kEpsilon model in particular.

Turbulence modeling in the simpleFoam solver

The simpleFoam solver is implemented in $FOAM_SOLVERS/incompressible/simpleFoam, where the following files and directories can be found:

createFields.H overSimpleFoam porousSimpleFoam SRFSimpleFoam
Make pEqn.H simpleFoam.C UEqn.H

The directories overSimpleFoam, porousSimpleFoam and SRFSimpleFoam are other versions of the simpleFoam solver, which we do not consider here. The Make directory contains compiler instructions. The main files for the simpleFoam solver are thus:

createFields.H pEqn.H simpleFoam.C UEqn.H

Here the simpleFoam.C file is the main file, in which the other files are included by #include-statements.

The turbulence modeling declarations are added to the simpleFoam solver by including a header file in the simpleFoam.C file:

#include "turbulentTransportModel.H"

That file is located in $FOAM_SRC/TurbulenceModels/incompressible/turbulentTransportModels, and it contains:

#include "IncompressibleTurbulenceModel.H"
#include "laminarModel.H"
#include "RASModel.H"
#include "LESModel.H"
#include "incompressible/transportModel/transportModel.H"

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

namespace Foam
{
    namespace incompressible
    {
        typedef IncompressibleTurbulenceModel<transportModel> turbulenceModel;

        typedef laminarModel<turbulenceModel> laminarModel;
typedef RASModel<turbulenceModel> RASModel;
        typedef LESModel<turbulenceModel> LESModel;

        template<class BasicCompressibleTurbulenceModel>
        autoPtr<BasicCompressibleTurbulenceModel> New
        (
            const volVectorField& U,
            const surfaceScalarField& phi,
            const typename BasicCompressibleTurbulenceModel::transportModel&
                transport,
            const word& propertiesName = turbulenceModel::propertiesName
        );
    }
}

I.e. it adds the declations of some additional classes to the simpleFoam solver: IncompressibleTurbulenceModel, laminarModel, RASModel, LESModel and transportModel. Then it uses those additional classes while defining some typeDefs (we will soon refer back to the first one, so please have a look at it) and declaring a templated function named New. The effect of the New function will depend on which turbulence model we have chosen, and the New function that will eventually be executed belongs to a specific turbulence model, as we will see later. However, at this point the New function is only declared and defined according to the above pieces of code. It is not executed.

The next step in simpleFoam (with respect to turbulence modeling) is to construct a pointer for the turbulence, which is done at the end of the included file createFields.H:

autoPtr<incompressible::turbulenceModel> turbulence
(
    incompressible::turbulenceModel::New(U, phi, laminarTransport)
);

I.e., the function New that is called (see that single line with the word "New" above) belongs to namespace incompressible (see the turbulentTransportModel.H file) and to the class turbulenceModel (which is a typeDef of IncompressibleTurbulenceModel<transportModel> in turbulentTransportModel.H - I asked you to look specifically at it before). The New function in turbulentTransportModel.H returns an autoPtr<BasicCompressibleTurbulenceModel>, and at the first line in the above piece of code we make sure that the template parameter BasicCompressibleTurbulenceModel in turbulentTransportModel.H should be interpreted as incompressible::turbulenceModel, where "incompressible" and "turbulenceModel" were described just above, ending up with a BasicCompressibleTurbulenceModel template parameter as IncompressibleTurbulenceModel<transportModel>, and a fact that the New function of that class will be called. We thus need to dig further into the class IncompressibleTurbulenceModel<transportModel> later (in the following section, since we here have to stay at a somewhat high level in the code). The arguments that are passed to the function New in the above piece of code are U, phi, and laminarTransport. The default value is used for the fourth argument propertiesName (see the first declaration above), i.e. turbulenceModel::propertiesName. It is set in $FOAM_SRC/TurbulenceModels/turbulenceModels/turbulenceModel.C (we are explicitly naming that class in turbulenceModel::propertiesName) as:

const Foam::word Foam::turbulenceModel::propertiesName("turbulenceProperties");

Here we see that the default fourth argument for the New function is the file name of the turbulence properties dictionary, which we can thus change if we like in createFields.H (you can try that as an exercise if you like). The fields U and phi are the velocity and face flux fields that are part of the solution in simpleFoam, and those fields are constructed in createFields.H. The object laminarTransport is also constructed in createFields.H by:

singlePhaseTransportModel laminarTransport(U, phi);

Note that for that line to work it was necessary to add the declaration of that class in simpleFoam.C:

#include "singlePhaseTransportModel.H"

This constructs the object named laminarTransport as a singlePhaseTransportModel. That class is implemented in $FOAM_SRC/transportModels/incompressible/singlePhaseTransportModel, and it makes sure that the transportProperties dictionary is read and that the viscosity model is set (Newtonian or one of the non-Newtonian models). The laminarTransport object is supplied as one of the parameters to the function named New, since eddy-viscosity turbulence models also affect the flow through the viscosity, by computing an efficient viscosity, and the turbulence model must thus also take into account the effects of laminar viscosity and non-Newtonian behaviour on the laminar viscosity. We will see this later.

The pointer named turbulence and its member functions are then used in the simpleFoam solver (in simpleFoam.C and UEqn.H) to add the solution and effects of the turbulence model of choise, mainly:

turbulence->validate();
+ turbulence->divDevReff(U)
turbulence->correct();

We realize that the object named turbulence is a pointer, since the functions of the object are called using '->'.The operation of those member functions depend on which turbulence model we are using, which is specified when we run the case. We must thus later figure out where those functions are defined for each turbulence model.

It is worth to repeat that is not obvious, when looking only at the files in the simpleFoam directory, how the fields for the turbulence model are constructed, and why they are read from the start time directory and written to the following time directories. Different turbulence models use different fields, so it is each turbulence model that constructs the necessary fields for that particular turbulence model, and makes sure that they are written to the time directories when the line runTime.write() in simpleFoam.C is executed.

In the next section we will dig deeper into the code to figure out what the New function is doing, and how the turbulence pointer ends up pointing at a particular turbulence model. The section can be skipped if you just want to accept that the turbulence pointer points at the turbulence model class that you have specified in the constant/turbulenceProperties dictionary.

The function named New

This section is highly complicated, and not 100% clear. Those who are VERY interested can go through it and get back to me with suggested improvements.

It was discussed in the previous section that the function named New was called when constructing the pointer named turbulence, and that it belonged to the IncompressibleTurbulenceModel<transportModel> class. We will now have a look at the implementation of that class and we find the declaration of the function New in $FOAM_SRC/TurbulenceModels/incompressible/IncompressibleTurbulenceModel/IncompressibleTurbulenceModel.H:

        //- Return a reference to the selected turbulence model
        static autoPtr<IncompressibleTurbulenceModel> New
        (
            const volVectorField& U,
            const surfaceScalarField& phi,
            const TransportModel& trasportModel,
            const word& propertiesName = turbulenceModel::propertiesName
        );

The definition of the same function is found in the correcponding IncompressibleTurbulenceModel.C file:

template<class TransportModel>
Foam::autoPtr<Foam::IncompressibleTurbulenceModel<TransportModel>>
Foam::IncompressibleTurbulenceModel<TransportModel>::New
(
const volVectorField& U,
const surfaceScalarField& phi,
const TransportModel& transport,
const word& propertiesName
)
{
return autoPtr<IncompressibleTurbulenceModel>
(
  static_cast<IncompressibleTurbulenceModel*>(
  TurbulenceModel
  <
  geometricOneField,
geometricOneField,
  incompressibleTurbulenceModel,
  TransportModel
  >::New
  (
  geometricOneField(),
  geometricOneField(),
U,
phi,
phi,
transport,
  propertiesName
).ptr())
);
}

We see that the function returns the output of the New function (taking 7 arguments) of the base class TurbulenceModel, which IncompressibleTurbulenceModel Inherits from. That function is defined in $FOAM_SRC/TurbulenceModels/ turbulenceModels/TurbulenceModel.C, as:

template
<
    class Alpha,
    class Rho,
    class BasicTurbulenceModel,
    class TransportModel
>
Foam::autoPtr
<
    Foam::TurbulenceModel<Alpha, Rho, BasicTurbulenceModel, TransportModel>
>
Foam::TurbulenceModel<Alpha, Rho, BasicTurbulenceModel, TransportModel>::New
(
    const alphaField& alpha,
    const rhoField& rho,
    const volVectorField& U,
    const surfaceScalarField& alphaRhoPhi,
    const surfaceScalarField& phi,
    const transportModel& transport,
    const word& propertiesName
)
{
    // get model name, but do not register the dictionary
    // otherwise it is registered in the database twice
    const word modelType
    (
        IOdictionary
        (
            IOobject
            (
                IOobject::groupName(propertiesName, U.group()),
                U.time().constant(),
                U.db(),
                IOobject::MUST_READ_IF_MODIFIED,
                IOobject::NO_WRITE,
                false
            )
        ).lookup("simulationType")
    );

    Info<< "Selecting turbulence model type " << modelType << endl;

    typename dictionaryConstructorTable::iterator cstrIter =
        dictionaryConstructorTablePtr_->find(modelType);

    if (!cstrIter.found())
    {
        FatalErrorInFunction
            << "Unknown TurbulenceModel type "
            << modelType << nl << nl
            << "Valid TurbulenceModel types:" << endl
            << dictionaryConstructorTablePtr_->sortedToc()
            << exit(FatalError);
    }

    return autoPtr<TurbulenceModel>
    (
        cstrIter()(alpha, rho, U, alphaRhoPhi, phi, transport, propertiesName)
    );
}

Following the code above, it reads the turbulenceProperties dictionary, determines the simulationType/modelType, and checks that the simulationType/modelType is one of the available types. Without describing the details related to the object cstrIter and its class (I don't know all the details, so please help me sort this out), the return statement at the end returns a pointer to the particular turbulence model type (laminar/RAS/LES) that is specified in the turbulenceProperties dictionary. In the case of the default simpleFoam/pitzDaily tutorial the turbulence model type is set to RAS. We end up in RASModel.C in the corresponding function named New:

template<class BasicTurbulenceModel>
Foam::autoPtr<Foam::RASModel<BasicTurbulenceModel>>
Foam::RASModel<BasicTurbulenceModel>::New
(
    const alphaField& alpha,
    const rhoField& rho,
    const volVectorField& U,
    const surfaceScalarField& alphaRhoPhi,
    const surfaceScalarField& phi,
    const transportModel& transport,
    const word& propertiesName
)
{
    // get model name, but do not register the dictionary
    // otherwise it is registered in the database twice
    const word modelType
    (
        IOdictionary
        (
            IOobject
            (
                IOobject::groupName(propertiesName, U.group()),
                U.time().constant(),
                U.db(),
                IOobject::MUST_READ_IF_MODIFIED,
                IOobject::NO_WRITE,
                false
            )
        ).subDict("RAS").lookup("RASModel")
    );
    
    Info<< "Selecting RAS turbulence model " << modelType << endl;
    
    typename dictionaryConstructorTable::iterator cstrIter =
        dictionaryConstructorTablePtr_->find(modelType);
    
    if (cstrIter == dictionaryConstructorTablePtr_->end())
    {
        FatalErrorInFunction
            << "Unknown RASModel type "
            << modelType << nl << nl
            << "Valid RASModel types:" << endl
            << dictionaryConstructorTablePtr_->sortedToc()
            << exit(FatalError);
    }
    
    return autoPtr<RASModel>
    (
        cstrIter()(alpha, rho, U, alphaRhoPhi, phi, transport, propertiesName)
    );
}

We see that the particular RAS turbulence model is selected based on the setting in the sub-dictionary RASModel, in the turbulenceProperties dictionary. The same procedure as before is again used to point at the particular RAS turbulence model, and we end up in the constructor of the kEpsilon model (according to the standard settings of the simpleFoam/pitzDaily tutorial), in $FOAM_SRC/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C, which initializes the un-named object that the object named turbulence is pointing at, and writes out the model coefficients (as we will see later). The object named turbulence can now be used to call the member fucntions in the kEpsilon class, such as validate(), divDevRef() and correct(), which we will have a look at soon.

The kEpsilon Model

The kEpsilon model is implemented in $FOAM_SRC/TurbulenceModels/turbulenceModels/RAS/kEpsilon. The class inheritance is given by:

template<class BasicTurbulenceModel>
class kEpsilon
:
public eddyViscosity<RASModel<BasicTurbulenceModel>>
{

I.e., the kEpsilon class is an eddyViscosity model, with the RASModel specialization using the template parameter BasicTurbulenceModel (in this case interpreted as IncompressibleTurbulenceModel<transportModel>.

The class declares and defines some member data and member functions as well as some typeDefs and a TypeName.

The contructor is defined in kEpsilon.C:

template<class BasicTurbulenceModel>
kEpsilon<BasicTurbulenceModel>::kEpsilon
(
    const alphaField& alpha,
    const rhoField& rho,
    const volVectorField& U,
    const surfaceScalarField& alphaRhoPhi,
    const surfaceScalarField& phi,
    const transportModel& transport,
    const word& propertiesName,
    const word& type
)
:
    eddyViscosity<RASModel<BasicTurbulenceModel>>
    (
        type,
        alpha,
        rho,
        U,
        alphaRhoPhi,
        phi,
        transport,
        propertiesName
    ),

    Cmu_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cmu",
            this->coeffDict_,
            0.09
        )
    ),
    C1_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "C1",
            this->coeffDict_,
            1.44
        )
    ),
    C2_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "C2",
            this->coeffDict_,
            1.92
        )
    ),
    C3_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "C3",
            this->coeffDict_,
            -0.33
        )
    ),
    sigmak_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "sigmak",
            this->coeffDict_,
            1.0
        )
    ),
    sigmaEps_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "sigmaEps",
            this->coeffDict_,
            1.3
        )
    ),

    k_
    (
        IOobject
        (
            IOobject::groupName("k", U.group()),
            this->runTime_.timeName(),
            this->mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        this->mesh_
    ),
    epsilon_
    (
        IOobject
        (
            IOobject::groupName("epsilon", U.group()),
            this->runTime_.timeName(),
            this->mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        this->mesh_
    )
{
    bound(k_, this->kMin_);
    bound(epsilon_, this->epsilonMin_);

    if (type == typeName)
    {
        this->printCoeffs(type);
    }
}

The constructor first uses the constructor of the base class eddyViscosity to initialize, followed by initialization of the properties that are specific to kEpsilon, finishing with a bounding of k and epsilon and writing the model coefficients to the standard output. We see that the fields k_ and epsilon_ are constructed similar to the p field in createFields.H in the top-level solver. This will make sure that those fields are available for the turbulence model, and that they are written to disk the same way as the p field.

The three functions that are called in the top-level solver simpleFoam are validate(), divDevReff(U) and correct(). We find only the correct() function in kEpsilon.C, while we have to look for the other functions in the base class(es). The correct() function is implemented as:

template<class BasicTurbulenceModel>
void kEpsilon<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_;
volScalarField& nut = this->nut_;
fv::options& fvOptions(fv::options::New(this->mesh_));

eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();

volScalarField::Internal divU
(
fvc::div(fvc::absolute(this->phi(), U))().v()
);

tmp<volTensorField> tgradU = fvc::grad(U);
volScalarField::Internal G
(
this->GName(),
nut.v()*(dev(twoSymm(tgradU().v())) && tgradU().v())
);
tgradU.clear();

// Update epsilon and G at the wall
epsilon_.boundaryFieldRef().updateCoeffs();

// Dissipation equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(alpha, rho, epsilon_)
+ fvm::div(alphaRhoPhi, epsilon_)
- fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
==
C1_*alpha()*rho()*G*epsilon_()/k_()
- fvm::SuSp(((2.0/3.0)*C1_ + C3_)*alpha()*rho()*divU, epsilon_)
- fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_)
+ epsilonSource()
+ fvOptions(alpha, rho, epsilon_)
);

epsEqn.ref().relax();
fvOptions.constrain(epsEqn.ref());
epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
solve(epsEqn);
fvOptions.correct(epsilon_);
bound(epsilon_, this->epsilonMin_);

// 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(alpha()*rho()*epsilon_()/k_(), k_)
+ kSource()
+ fvOptions(alpha, rho, k_)
);

kEqn.ref().relax();
fvOptions.constrain(kEqn.ref());
solve(kEqn);
fvOptions.correct(k_);
bound(k_, this->kMin_);

correctNut();
}

We see that it first calls the correct() function of the base class, in case there are some general things to be done. Then it discretizes and solves the epsilon_ and k_ fields, followed by relaxation and bounding of the fields. There is also a possibility to manipulate the model through the fvOptions functionality. The final line calls the function correctNut(), which is implemented in the same file and calculates the new value for the turbulent viscosity:

template<class BasicTurbulenceModel>
void kEpsilon<BasicTurbulenceModel>::correctNut()
{
this->nut_ = Cmu_*sqr(k_)/epsilon_;
this->nut_.correctBoundaryConditions();
fv::options::New(this->mesh_).correct(this->nut_);

BasicTurbulenceModel::correctNut();
}

The validate() function is implemented in the base class eddyViscosity ($FOAM_SRC/TurbulenceModels/turbulenceModels/eddyViscosity/eddyViscosity.C):

template<class BasicTurbulenceModel>
void Foam::eddyViscosity<BasicTurbulenceModel>::validate()
{
correctNut();
}

This just makes sure that the turbulent viscosity is computed correctly before starting the simulation.

The divDevReff(U) function is neither implemented in the eddyViscosity class, so we have a look in the base class linearViscousStress, located in $FOAM_SRC/TurbulenceModels/turbulenceModels/linearViscousStress. We can't find it there either, so we look in the base class, which is referred to as BasicTurbulenceModel (which we know is IncompressibleTurbulenceModel<transportModel> in our case). We find in IncompressibleTurbulenceModel.C:

template<class TransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::IncompressibleTurbulenceModel<TransportModel>::divDevReff
(
volVectorField& U
) const
{
return divDevRhoReff(U);
}

This in turn calls the divDevRhoReff(U) function, assuming that we are currently in the kEpsilon class. That function can't be found there, and neither in the eddyViscosity class, but in the linearViscousStress class ($FOAM_SRC//TurbulenceModels/turbulenceModels/linearViscousStress/linearViscousStress.C):

template<class BasicTurbulenceModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::linearViscousStress<BasicTurbulenceModel>::divDevRhoReff
(
volVectorField& U
) const
{
return
(
- fvc::div((this->alpha_*this->rho_*this->nuEff())*dev2(T(fvc::grad(U))))
- fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U)
);
}

We remember that for the incompressible case the density (rho_) was set to constant 1. The nuEff() function can NOT be found by the same search pattern as before: kEpsilon - eddyViscosity - linearViscousStress, but we recall that kEpsilon is a RASModel, and find the function in $FOAM_SRC/TurbulenceModels/turbulenceModels/RAS/RASModel/RASModel.H:

 //- Return the effective viscosity
virtual tmp<volScalarField> nuEff() const
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject::groupName("nuEff", this->U_.group()),
this->nut() + this->nu()
)
);
}

The nut() function is implemented in eddyViscosity.H:

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

We can imagine that the nu() function returns the laminar kinematic viscosity, and that the nuEff() function thus returns the sum of the turbulent and laminar kinematic viscosities.

We see that only the things that are special for the kEpsilon class are implemented in that class. The rest is implemented in the base classes. If we are to implement a special version of the kEpsilon model we need to focus on the kEpsilon.H and kEpsilon.C files.

The kOmegaSST and kOmegaSSTSAS models

The kOmegaSST model is implemented in $FOAM_SRC/TurbulenceModels/turbulenceModels/RAS/kOmegaSST. If we look in those files we can't find any definition of the member functions. We realize that the kOmegaSST class inherits from another class that is also named kOmegaSST:

#include "kOmegaSSTBase.H"
...
template<class BasicTurbulenceModel>
class kOmegaSST
:
    public Foam::kOmegaSST
    <
        eddyViscosity<RASModel<BasicTurbulenceModel>>,
        BasicTurbulenceModel
    >
{

The file kOmegaSSTBase.H is located in $FOAM_SRC/TurbulenceModels/turbulenceModels/Base/kOmegaSST, and we see that the files in that directory include a full description and implementation of the kOmegaSST model. The implementation can be compared with  the mathematical desrription (written for version 1.6, so it is not exact). In particular we find a source term named Qsas. It is a preparation for the kOmegaSSTSAS model, and the preparation for the kOmegaSSTSAS model is also the reason to make a separate Base version of the kOmegaSST model.

The kOmegaSSTSAS model is located in $FOAM_SRC/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTSAS. We see that the kOmegaSSTSAS class inherits from the kOmegaSST class:

#include "kOmegaSST.H"
...
template<class BasicTurbulenceModel>
class kOmegaSSTSAS
:
    public kOmegaSST<BasicTurbulenceModel>
{

The kOmegaSSTSAS sub-class adds some coefficients, adds a run-time selectable delta model, and remimplements the Qsas function. We will not go into the details at this point.

Implement the kOmegaSSTF model

We will in the rest of the sections go through the basic steps of implementing a new kOmegaSST model, named kOmegaSSTF. It will not be the most ideal way of doing it, but it will show one way of doing it. An old description of how to implement the kOmegaSSTF model is found  here. It sets a limit to the turbulent viscosity, which has a similar effect as the added source term in the kOmegaSSTSAS model. We will simply make a copy of the original kOmegaSST model and reimplement the correctNut(S2) function that computes the turbulent viscosity. We will as well see how the entire correct() function can be re-implemented (although we will not do any changes to it here).

As noticed above, the implementation is now templated. That makes it difficult to only copy the directory of the turbulence model we are modifying. Instead we copy the entire TurbulenceModels directory (the fact that the directory name starts with a capital letter tells us that it is templated), and we have the possibility to modify whatever we like in that entire directory structure:

OF1706+
foam
cp -r --parents src/TurbulenceModels $WM_PROJECT_USER_DIR
cd $WM_PROJECT_USER_DIR/src/TurbulenceModels

Find all the Make directories:

$ find . -name Make
./incompressible/Make
./compressible/Make
./turbulenceModels/Make
./schemes/Make

Change the location of the compiled files in all the Make/files files:

sed -i s/FOAM_LIBBIN/FOAM_USER_LIBBIN/g ./*/Make/files

Compile (this will take 10min):

./Allwmake

When that is finished (don't proceed until then unless you know what you are doing) there will be three new shared-object files in $WM_PROJECT_USER_DIR/platforms/linux64GccDPInt32Opt/lib:

libcompressibleTurbulenceModels.so
libincompressibleTurbulenceModels.so
libturbulenceModels.so

Those will be used instead of the ones in the main installation, since the environment variable LD_LIBRARY_PATH points at that directory before pointing at the original directory ($WM_PROJECT_DIR/platforms/linux64GccDPInt32Opt/lib). Have a look at the order of the directories listed in the LD_LIBRARY_PATH environment variable to see this. We can check which libraries named something with "urbulenceModels" will be used by the simpleFoam solver by:

ldd `which simpleFoam` | grep urbulenceModels #Note that the quotation marks must be the backquote

which returns pointers to the files in your working directory (here oscfd-plus, if your user name is oscfd - note that this output requires that the compilation above has finished) :

libturbulenceModels.so => /home/oscfd/OpenFOAM/oscfd-plus/platforms/linux64GccDPInt32Opt/lib/libturbulenceModels.so (0x00007f4267931000)
libincompressibleTurbulenceModels.so => /home/oscfd/OpenFOAM/oscfd-plus/platforms/linux64GccDPInt32Opt/lib/libincompressibleTurbulenceModels.so (0x00007f426742b000)
libcompressibleTurbulenceModels.so => /home/oscfd/OpenFOAM/oscfd-plus/platforms/linux64GccDPInt32Opt/lib/libcompressibleTurbulenceModels.so (0x00007f42615d9000)
Thus, we know that the simpleFoam solver will use the libraries that we have in our working directory.

We can now copy the kOmegaSST model and do our modifications (first rename files and class name):
cp -r turbulenceModels/RAS/kOmegaSST turbulenceModels/RAS/kOmegaSSTF
mv turbulenceModels/RAS/kOmegaSSTF/kOmegaSST.H turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.H
mv turbulenceModels/RAS/kOmegaSSTF/kOmegaSST.C turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C

Investigating the two files we see that we want to change the class name from kOmegaSST to kOmegaSSTF. That can be done by a simple sed command. However, there is then a risk that also lines containing the word kOmegaSSTBase will be changed, and they should not be changed. Here is a sed command that omits lines that contain the word kOmegaSSTBase (and also gives an example of how to also omit lines with the words kOmegaSSTF and dummy):

sed -Ei '/(kOmegaSSTBase|kOmegaSSTF|dummy)/!s/kOmegaSST/kOmegaSSTF/g' turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.*

Two important notes about the above command: 1: I experienced that it is important to have the flags in the order Ei and not iE, since otherwise the omission will not work. 2: According to comments in forums the command will not omit consecutive lines with those words. However, in this particular case it works.

Now we need to make sure that our class is compiled. The usual way of doing that is to modify the Make/files file of the library, in this case $WM_PROJECT_USER_DIR/src/TurbulenceModels/turbulenceModels/Make/files. However, we can't find the original kOmegaSST.C file listed there. That is because the class is templated, and the compilation is instead taken care of by macros. We look for the string kOmegaSST except in the explicit turbulence model classes and in the binary files:

grep -r kOmegaSST . | grep -v Base | grep -v RAS | grep -v LES | grep -v DES | grep -v lnInclude | grep -v linux

We find that for incompressible turbulence modeling there is a file incompressible/turbulentTransportModels/turbulentTransportModels.C that includes the lines:

#include "kOmegaSST.H"
makeRASModel(kOmegaSST);

The second line is a macro that makes sure that the kOmegaSST model is compiled as one of the alternatives of the templating. We simply add in that file, under the above lines:

#include "kOmegaSSTF.H"
makeRASModel(kOmegaSSTF);

Try typing ./Allwmake and it will complain that it can't find the file named kOmegaSSTF.H. The reason is that the linking in the lnInclude directory was done without that file being present. We can update the linking by:

wmakeLnInclude -u turbulenceModels

Try compiling again:

./Allwmake

Notice that it will not be shown in the output of the compilation that the kOmegaSSTF model is compiled, but it is. The compilation will take more time than expected when we only added one additonal turbulence model, and that is because all the templated alternatives mentioned in the above file must be recompiled.

Make sure that it works by running the a fresh simpleFoam/pitzDaily case:

rm -r $FOAM_RUN/pitzDaily
cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily $FOAM_RUN
sed -i s/kEpsilon/kOmegaSSTF/g $FOAM_RUN/pitzDaily/constant/turbulenceProperties
blockMesh -case $FOAM_RUN/pitzDaily
simpleFoam -noFunctionObjects -case $FOAM_RUN/pitzDaily >& $FOAM_RUN/pitzDaily/log&

Look at the start of that log file and read "Selecting RAS turbulence model kOmegaSSTF", which shows that the kOmegaSSTF model is used.

Now we will do some modifications to our new model. We will reimplement both the correct() and the correctNut(S2) functions. It should be noted that we do not have to reimplement the correct() function in this case, but it is done here as an example of how to do it if there are some moficiations that need to be done in that function. In the correctNut(S2) function we will introduce the limit of the turbulent viscosity, which is the idea of the kOmegaSSTF model (see the slides referred to earlier). The correct() function is inherited to the kOmegaSSTF sub-class from the kOmegaSSTBase sub-class (located in turbulenceModels/Base/kOmegaSST). In order to reimplement it in the sub-class kOmegaSSTF we first need to add the declaration of the function in the .H file, and then we need to add the new definition of the function in the .C file.

In turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.H, add at the end of the class declaration (after the destructor and before the curly bracket that is followed by a semicolon):

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

In turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C, add at the end of the namespace RASModels (before the two curly brackets at the end), the entire correct() function of the kOmegaSSTBase.C file but change the template parameter name from BasicEddyViscosityModel to BasicTurbulenceModel (three instances), to adapt to the template parameter names in the rest of the kOmegaSTF.C file. Also change the class name of the correct() function from kOmegaSSTBase to kOmegaSST (one instance). Add as well an Info statement at the beginning so that we can check that our new class it is in fact used. This is what the start of the correct() function in kOmegaSSTF.C should now look like (above modifications marked in yellow):

template<class BasicTurbulenceModel>
void kOmegaSSTF<BasicTurbulenceModel>::correct()
{
    if (!this->turbulence_)
    {
       return;
    }

    Info << "This is kOmegaSSTF" << endl;

    // 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_;
    fv::options& fvOptions(fv::options::New(this->mesh_));

    BasicTurbulenceModel::correct();
...

Note in the above piece of code that you can also manipulate turbulence models using fvOptions, but we are not doing that here.

Try compiling with ./Allwmake, and it will not recognize that we did any modifications. That is because we modified a file that is not listed in any Make/files file. The way we added the compilation of the kOmegaSSTF model was to modify the file incompressible/turbulentTransportModels/turbulentTransportModels.C. The model will thus be recompiled if that file is recompiled, i.e. type:

touch incompressible/turbulentTransportModels/turbulentTransportModels.C
./Allwmake

You will get lots of error messages. Don't panic! The error messages are of the following types:

‘omega_’ was not declared in this scope
there are no arguments to ‘DomegaEff’ that depend on a template parameter, so a declaration of ‘DomegaEff’ must be available
‘DomegaEff’ was not declared in this scope, and no declarations were found by argument-dependent lookup at the point of instantiatio

The reason we get those messages is that those objects and functions are not declared in our sub-class, but in some base-class. The current class does not know about those objects and functions, but the object belonging to the class does. We can therefore use the object belonging to the class to reach those objects and functions. This is what is done in the above code, where e.g. a local reference object alphaField is constructed as:

    const alphaField& alpha = this->alpha_;

The word "this" is a pointer to the object that belongs to our sub-class, and we here make alphaField refer to the object alpha_ of the object that belongs to our sub-class. We simply have to do the same for all the objects and fucntions that were mentioned in the error message. Make sure to run the following commands only once:

sed -i s/"omega_"/"this->omega_"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C
sed -i s/"alphaOmega2_"/"this->alphaOmega2_"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C
sed -i s/"k_"/"this->k_"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C
sed -i s/"DomegaEff(F1)"/"this->DomegaEff(F1)"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C
sed -i s/"GbyNu(GbyNu0, F23(), S2())"/"this->GbyNu(GbyNu0, F23(), S2())"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C
sed -i s/"Qsas(S2(), gamma, beta)"/"this->Qsas(S2(), gamma, beta)"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C
sed -i s/"omegaSource()"/"this->omegaSource()"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C
sed -i s/"DkEff(F1)"/"this->DkEff(F1)"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C
sed -i s/"Pk(G)"/"this->Pk(G)"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C
sed -i s/"epsilonByk(F1, tgradU())"/"this->epsilonByk(F1, tgradU())"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C
sed -i s/"kSource()"/"this->kSource()"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C
sed -i s/"DomegaEff(F1), omega_)"/"this->DomegaEff(F1), omega_)"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C
sed -i s/"epsilonByk(F1, tgradU()), k_)"/"this->epsilonByk(F1, tgradU()), k_)"/g turbulenceModels/RAS/kOmegaSSTF/kOmegaSSTF.C

Try compiling again, and it should work:

touch incompressible/turbulentTransportModels/turbulentTransportModels.C
./Allwmake

Test using the simpleFoam/pitzDaily case:

simpleFoam -noFunctionObjects -case $FOAM_RUN/pitzDaily >& $FOAM_RUN/pitzDaily/log&

Look before the solution of omega at each time step, and make sure that the Info statement is written.

Now we want to change how the turbulent viscosity is corrected. We see at the end of the correct() function that the turbulent viscosity is corrected by a call to a function named correctNut(S2). That function is implemented at thebeginning of the kOmegaSSTF.C file, as:

template<class BasicTurbulenceModel>
void kOmegaSSTF<BasicTurbulenceModel>::correctNut(const volScalarField& S2)
{
// Correct the turbulence viscosity
  kOmegaSSTBase<eddyViscosity<RASModel<BasicTurbulenceModel>>>::correctNut
(
  S2
  );

  // Correct the turbulence thermal diffusivity
BasicTurbulenceModel::correctNut();
}

There are calls to two functions here. The second one is to a function implemented in incompressible/incompressibleTurbulenceModel.H:

        //- ***HGW Temporary function to be removed when the run-time selectable
// thermal transport layer is complete
virtual void correctNut()
{}

The function does not do anything, so we can leave that call as it is.

The first function call in correctNut(S2), above, is to a function with the same name in the base class kOmegaSSTBase. It is implemented in turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C:

template<class BasicEddyViscosityModel>
void kOmegaSSTBase<BasicEddyViscosityModel>::correctNut
(
const volScalarField& S2
)
{
// Correct the turbulence viscosity
this->nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
this->nut_.correctBoundaryConditions();
fv::options::New(this->mesh_).correct(this->nut_);
}

In order to reimplement this in the sub-class we can simply avoid calling the base class function and instead do those operations in the sub-class function. The correctNut(S2) function should then read (where we have added "this->" as discussed above):

template<class BasicTurbulenceModel>
void kOmegaSSTF<BasicTurbulenceModel>::correctNut(const volScalarField& S2)
{
// Correct the turbulence viscosity
this->nut_ = this->a1_*this->k_/max(this->a1_*this->omega_, this->b1_*this->F23()*sqrt(S2));
this->nut_.correctBoundaryConditions();
fv::options::New(this->mesh_).correct(this->nut_);

// Correct the turbulence thermal diffusivity
BasicTurbulenceModel::correctNut();
}

 

Compile and test:

touch incompressible/turbulentTransportModels/turbulentTransportModels.C
./Allwmake
simpleFoam -noFunctionObjects -case $FOAM_RUN/pitzDaily >& $FOAM_RUN/pitzDaily/log&

It runs, but the implementation is still the same as before. Now we will reimplement the first line that we added above, i.e.:

    this->nut_ = this->a1_*this->k_/max(this->a1_*this->omega_, this->b1_*this->F23()*sqrt(S2));

to (See the slides regarding the kOmegaSSTF model linked to above. The implementation below could be improved!!!):

    // Compute Filter
    scalar alph = 3.0; // Should be in a dictionary
    scalarField Lt = sqrt(this->k_)/(this->betaStar_*this->omega_);
    scalarField lt = alph*Foam::max(Foam::pow(this->mesh_.V().field(), 1.0/3.0),
                     (mag(this->U_)*this->runTime_.deltaT())->internalField());
    // Recalculate viscosity
    this->nut_.internalField() == Foam::min(Foam::pow(lt/Lt, 4.0/3.0), 1.0)*
       (this->a1_*this->k_/max(this->a1_*this->omega_, this->b1_*this->F23()*sqrt(S2)))
       ->internalField();

Compile and check that it goes through:

touch incompressible/turbulentTransportModels/turbulentTransportModels.C
./Allwmake

Run the kOmegaSSTF model with pimpleFoam for the pitzDaily case

The kOmegaSSTF model is supposed to run in unsteady mode, so that it can resolve as much as possible of the turbulence. We therefore have to switch from the steady-state simpleFoam solver to the unsteady pimpleFoam solver. For that we need to manipulate the pitzDaily case:

run
cd pitzDaily
foamListTimes -rm

Modify the files in the system directory, for use with pimpleFoam:

cp $FOAM_TUTORIALS/incompressible/pimpleFoam/RAS/TJunction/system/{fvSolution,fvSchemes} system
sed -i s/epsilon/omega/g system/fvSchemes
echo "wallDist { method meshWave; }" >> system/fvSchemes
sed -i s/epsilon/omega/g system/fvSolution
sed -i s/simpleFoam/pimpleFoam/g system/controlDict
sed -i s/2000/0.3/g system/controlDict;
sed -i s/"1;"/"0.0001;"/g system/controlDict
sed -i s/"writeCompression off"/"writeCompression on"/g system/controlDict
echo "adjustTimeStep no;" >> system/controlDict
echo "maxCo 5;" >> system/controlDict

Run the case and watch in paraFoam how the flow becomes highly unsteady:

pimpleFoam -noFunctionObjects -case $FOAM_RUN/pitzDaily >& $FOAM_RUN/pitzDaily/log&
paraFoam

Movie:

Velocity at final time step with kOmegaSSTF:

 

Velocity at final time step with kOmegaSST:

Final file, kOmegaSSTF.H (indentation broken when pasting)

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 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::kOmegaSSTF

Group
grpRASTurbulence

Description
Implementation of the k-omega-SST turbulence model for
incompressible and compressible flows.

Light wrapper around base class.

Turbulence model described in:
\verbatim
Menter, F. R. & Esch, T. (2001).
Elements of Industrial Heat Transfer Prediction.
16th Brazilian Congress of Mechanical Engineering (COBEM).
\endverbatim

with updated coefficients from
\verbatim
Menter, F. R., Kuntz, M., and Langtry, R. (2003).
Ten Years of Industrial Experience with the SST Turbulence Model.
Turbulence, Heat and Mass Transfer 4, ed: K. Hanjalic, Y. Nagano,
& M. Tummers, Begell House, Inc., 625 - 632.
\endverbatim

but with the consistent production terms from the 2001 paper as form in the
2003 paper is a typo, see
\verbatim
http://turbmodels.larc.nasa.gov/sst.html
\endverbatim

and the addition of the optional F3 term for rough walls from
\verbatim
Hellsten, A. (1998).
"Some Improvements in Menter’s k-omega-SST turbulence model"
29th AIAA Fluid Dynamics Conference, AIAA-98-2554.
\endverbatim

Note that this implementation is written in terms of alpha diffusion
coefficients rather than the more traditional sigma (alpha = 1/sigma) so
that the blending can be applied to all coefficuients in a consistent
manner. The paper suggests that sigma is blended but this would not be
consistent with the blending of the k-epsilon and k-omega models.

Also note that the error in the last term of equation (2) relating to
sigma has been corrected.

Wall-functions are applied in this implementation by using equations (14)
to specify the near-wall omega as appropriate.

The blending functions (15) and (16) are not currently used because of the
uncertainty in their origin, range of applicability and that if y+ becomes
sufficiently small blending u_tau in this manner clearly becomes nonsense.

The default model coefficients are
\verbatim
kOmegaSSTFCoeffs
{
alphaK1 0.85;
alphaK2 1.0;
alphaOmega1 0.5;
alphaOmega2 0.856;
beta1 0.075;
beta2 0.0828;
betaStar 0.09;
gamma1 5/9;
gamma2 0.44;
a1 0.31;
b1 1.0;
c1 10.0;
F3 no;
}
\endverbatim

SourceFiles
kOmegaSSTF.C

SeeAlso
kOmegaSSTBase.H

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

#ifndef kOmegaSSTF_H
#define kOmegaSSTF_H

#include "kOmegaSSTBase.H"
#include "RASModel.H"
#include "eddyViscosity.H"

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

namespace Foam
{
namespace RASModels
{

/*---------------------------------------------------------------------------*\
Class kOmegaSSTF Declaration
\*---------------------------------------------------------------------------*/

template<class BasicTurbulenceModel>
class kOmegaSSTF
:
public kOmegaSSTBase<eddyViscosity<RASModel<BasicTurbulenceModel>>>
{
// Private Member Functions

// Disallow default bitwise copy construct and assignment
kOmegaSSTF(const kOmegaSSTF&);
void operator=(const kOmegaSSTF&);


protected:

// Protecetd Member Functions

virtual void correctNut(const volScalarField& S2);
virtual void correctNut();


public:

typedef typename BasicTurbulenceModel::alphaField alphaField;
typedef typename BasicTurbulenceModel::rhoField rhoField;
typedef typename BasicTurbulenceModel::transportModel transportModel;


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


// Constructors

//- Construct from components
kOmegaSSTF
(
const alphaField& alpha,
const rhoField& rho,
const volVectorField& U,
const surfaceScalarField& alphaRhoPhi,
const surfaceScalarField& phi,
const transportModel& transport,
const word& propertiesName = turbulenceModel::propertiesName,
const word& type = typeName
);


//- Destructor
virtual ~kOmegaSSTF()
{}

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

};


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

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

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "kOmegaSSTF.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

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

Final file, kOmegaSSTF.C (indentation broken when pasting)

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

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

namespace Foam
{
namespace RASModels
{

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

template<class BasicTurbulenceModel>
void kOmegaSSTF<BasicTurbulenceModel>::correctNut(const volScalarField& S2)
{
// Correct the turbulence viscosity
// Compute Filter
scalar alph = 3.0; // Should be in a dictionary
scalarField Lt = sqrt(this->k_)/(this->betaStar_*this->omega_);
scalarField lt = alph*Foam::max(Foam::pow(this->mesh_.V().field(), 1.0/3.0),
(mag(this->U_)*this->runTime_.deltaT())->internalField());
// Recalculate viscosity
this->nut_.internalField() == Foam::min(Foam::pow(lt/Lt, 4.0/3.0), 1.0)*
(this->a1_*this->k_/max(this->a1_*this->omega_, this->b1_*this->F23()*sqrt(S2)))
->internalField();
this->nut_.correctBoundaryConditions();
fv::options::New(this->mesh_).correct(this->nut_);

// Correct the turbulence thermal diffusivity
BasicTurbulenceModel::correctNut();
}


template<class BasicTurbulenceModel>
void kOmegaSSTF<BasicTurbulenceModel>::correctNut()
{
correctNut(2*magSqr(symm(fvc::grad(this->U_))));
}


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

template<class BasicTurbulenceModel>
kOmegaSSTF<BasicTurbulenceModel>::kOmegaSSTF
(
const alphaField& alpha,
const rhoField& rho,
const volVectorField& U,
const surfaceScalarField& alphaRhoPhi,
const surfaceScalarField& phi,
const transportModel& transport,
const word& propertiesName,
const word& type
)
:
kOmegaSSTBase<eddyViscosity<RASModel<BasicTurbulenceModel>>>
(
type,
alpha,
rho,
U,
alphaRhoPhi,
phi,
transport,
propertiesName
)
{
if (type == typeName)
{
this->printCoeffs(type);
}
}


template<class BasicTurbulenceModel>
void kOmegaSSTF<BasicTurbulenceModel>::correct()
{
if (!this->turbulence_)
{
return;
}

Info << "This is kOmegaSSTF" << endl;

// 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_;
fv::options& fvOptions(fv::options::New(this->mesh_));

BasicTurbulenceModel::correct();

volScalarField::Internal divU(fvc::div(fvc::absolute(this->phi(), U)));

tmp<volTensorField> tgradU = fvc::grad(U);
volScalarField S2(2*magSqr(symm(tgradU())));
volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
volScalarField::Internal G(this->GName(), nut*GbyNu0);

// Update omega and G at the wall
this->omega_.boundaryFieldRef().updateCoeffs();

volScalarField CDkOmega
(
(2*this->alphaOmega2_)*(fvc::grad(this->k_) & fvc::grad(this->omega_))/this->omega_
);

volScalarField F1(this->F1(CDkOmega));
volScalarField F23(this->F23());

{
volScalarField::Internal gamma(this->gamma(F1));
volScalarField::Internal beta(this->beta(F1));

// Turbulent frequency equation
tmp<fvScalarMatrix> omegaEqn
(
fvm::ddt(alpha, rho, this->omega_)
+ fvm::div(alphaRhoPhi, this->omega_)
- fvm::laplacian(alpha*rho*this->DomegaEff(F1), this->omega_)
==
alpha()*rho()*gamma*this->GbyNu(GbyNu0, F23(), S2())
- fvm::SuSp((2.0/3.0)*alpha()*rho()*gamma*divU, this->omega_)
- fvm::Sp(alpha()*rho()*beta*this->omega_(), this->omega_)
- fvm::SuSp
(
alpha()*rho()*(F1() - scalar(1))*CDkOmega()/this->omega_(),
this->omega_
)
+ this->Qsas(S2(), gamma, beta)
+ this->omegaSource()
+ fvOptions(alpha, rho, this->omega_)
);

omegaEqn.ref().relax();
fvOptions.constrain(omegaEqn.ref());
omegaEqn.ref().boundaryManipulate(this->omega_.boundaryFieldRef());
solve(omegaEqn);
fvOptions.correct(this->omega_);
bound(this->omega_, this->omegaMin_);
}

// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(alpha, rho, this->k_)
+ fvm::div(alphaRhoPhi, this->k_)
- fvm::laplacian(alpha*rho*this->DkEff(F1), this->k_)
==
alpha()*rho()*this->Pk(G)
- fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, this->k_)
- fvm::Sp(alpha()*rho()*this->epsilonByk(F1, tgradU()), this->k_)
+ this->kSource()
+ fvOptions(alpha, rho, this->k_)
);

tgradU.clear();

kEqn.ref().relax();
fvOptions.constrain(kEqn.ref());
solve(kEqn);
fvOptions.correct(this->k_);
bound(this->k_, this->kMin_);

correctNut(S2);
}

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

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

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

Final file, turbulentTransportModels.C (indentation broken when pasting)

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

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

makeBaseTurbulenceModel
(
geometricOneField,
geometricOneField,
incompressibleTurbulenceModel,
IncompressibleTurbulenceModel,
transportModel
);


// -------------------------------------------------------------------------- //
// Laminar models
// -------------------------------------------------------------------------- //

#include "Stokes.H"
makeLaminarModel(Stokes);

#include "Maxwell.H"
makeLaminarModel(Maxwell);


// -------------------------------------------------------------------------- //
// RAS models
// -------------------------------------------------------------------------- //

#include "SpalartAllmaras.H"
makeRASModel(SpalartAllmaras);

#include "kEpsilon.H"
makeRASModel(kEpsilon);

#include "RNGkEpsilon.H"
makeRASModel(RNGkEpsilon);

#include "realizableKE.H"
makeRASModel(realizableKE);

#include "LaunderSharmaKE.H"
makeRASModel(LaunderSharmaKE);

#include "kOmega.H"
makeRASModel(kOmega);

#include "kOmegaSST.H"
makeRASModel(kOmegaSST);

#include "kOmegaSSTF.H"
makeRASModel(kOmegaSSTF);

#include "kOmegaSSTSAS.H"
makeRASModel(kOmegaSSTSAS);

#include "kOmegaSSTLM.H"
makeRASModel(kOmegaSSTLM);

#include "v2f.H"
makeRASModel(v2f);

#include "LRR.H"
makeRASModel(LRR);

#include "SSG.H"
makeRASModel(SSG);


// -------------------------------------------------------------------------- //
// LES models
// -------------------------------------------------------------------------- //

#include "Smagorinsky.H"
makeLESModel(Smagorinsky);

#include "WALE.H"
makeLESModel(WALE);

#include "kEqn.H"
makeLESModel(kEqn);

#include "dynamicKEqn.H"
makeLESModel(dynamicKEqn);

#include "dynamicLagrangian.H"
makeLESModel(dynamicLagrangian);

#include "SpalartAllmarasDES.H"
makeLESModel(SpalartAllmarasDES);

#include "SpalartAllmarasDDES.H"
makeLESModel(SpalartAllmarasDDES);

#include "SpalartAllmarasIDDES.H"
makeLESModel(SpalartAllmarasIDDES);

#include "DeardorffDiffStress.H"
makeLESModel(DeardorffDiffStress);

#include "kOmegaSSTDES.H"
makeLESModel(kOmegaSSTDES);

#include "kOmegaSSTDDES.H"
makeLESModel(kOmegaSSTDDES);

#include "kOmegaSSTIDDES.H"
makeLESModel(kOmegaSSTIDDES);


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