对OpenFOAM中 turbulenceFields 和湍流模型如何从模板类具体实现的笔记

广义胖鱼

我在查阅function Object turbulenceFields 的源代码时一直在思考一个问题：实现对湍流场的访问如R,k等是如何实现的。因为在源码turbulenceFields.C中没有定义任何显式的公式，经过很长时间的查阅之后我终于大概明白了，这涉及到OpenFOAM中通过复杂的宏实现湍流模型的具体化，所以就在这里记录一下。（在查阅时为了方便，用英语写的注解就懒得改了）

Turbulence model is defined by looking up in obr_(Object registry) which stored all accessible filed and models (in turbulenceFields.C)

in turbulenceFields.C/turbulenceFields::execute()

const auto& model =           obr_.lookupObject<incompressible::turbulenceModel>(modelName_);
/***/
processField<symmTensor>(f, model.R());

in turbulentTransportModel.H

this file defines name space incompressible and its type turbulenceModel. which gives incompressible::turbulenceModel is equal to type IncompressibleTurbulenceModel<transportModel>.
```
namespace Foam
{
    namespace incompressible
    {
        typedef IncompressibleTurbulenceModel<transportModel> turbulenceModel;
        /***/
    }
}
```
in IncompressibleTurbulenceModel.H

include src/TurbulenceModels/turbulenceModels/TurbulenceModel/TurbulenceModel.H
```
#include "TurbulenceModel.H"
```
then in TurbulenceModel.H, include src/TurbulenceModels/turbulenceModels/turbulenceModel.H
```
#include "turbulenceModel.H"
```
where in turbulenceModel.H, turbulence fields i. e. R and k, are defined as pure virtual functions: (has to be realized by subclass)
```
virtual tmp<volScalarField> k() const = 0;
virtual tmp<volSymmTensorField> R() const = 0;
```

LES Model is then derived from template parameters BasicTurbulenceModels (indicates the compressibility of the flow), in turbulentTransportModels.H, macro is defined to make LES models:

#define makeLESModel(Type)\
	makeTemplatedTurbulenceModel \
	(transportModelIncompressibleTurbulenceModel, LES, Type)

this macro is based in makeTurbulenceModel.H:

 #define makeTemplatedTurbulenceModel(BaseModel, SType, Type) \
 /***/

which gives template parameters:

BaseModel: transportModelIncompressibleTurbulenceModel
Stype: LES

in turbulentTransportModels.C, BaseModel is made by macro realization:

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

this macro is defined in makeTurbulenceModel.H as well:

#define defineTurbulenceModelTypes( \
    Alpha, Rho, baseModel, BaseModel, Transport) \
    namespace Foam \
    {               \
        typedef TurbulenceModel \
        < \
            Alpha, \
            Rho,    \
            baseModel,  \
            Transport  \
        > Transport##baseModel;                  
        typedef BaseModel<Transport>  \
            Transport##BaseModel;      \
        typedef laminarModel<Transport##BaseModel> \
            laminar##Transport##BaseModel; \
        typedef RASModel<Transport##BaseModel> \
            RAS##Transport##BaseModel;    \
        typedef LESModel<Transport##BaseModel> \
            LES##Transport##BaseModel;             
    }

#define makeBaseTurbulenceModel(Alpha, Rho, baseModel, BaseModel, Transport)   \
    /* Turbulence typedefs */\
    defineTurbulenceModelTypes(Alpha, Rho, baseModel, BaseModel, Transport)    \
    namespace Foam       \
    {                    \
        /* Turbulence selection table */\
        defineTemplateRunTimeSelectionTable \
        (                                    \
            Transport##baseModel,            \
            dictionary                        \
        );                                                  /***/
        /* LES models */  \
defineNamedTemplateTypeNameAndDebug(LES##Transport##BaseModel, 0);   \
	defineTemplateRunTimeSelectionTable \
        (LES##Transport##BaseModel, dictionary); \
        addToRunTimeSelectionTable               \
        (                                       \
            Transport##baseModel,               \
            LES##Transport##BaseModel,        \
            dictionary                         \
        );                                    \
    }

by macro, transportModel(Transport) is attached to IncompressibleTurbulenceModel(BaseModel) making a BaseModel of type IncompressibleTurbulenceModel<transportModel>. As noted above, this is also aliased as incompressible::turbulenceModel

the WALE model is made by macro makeLESModel():

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

which gives template parameters:

BaseModel: transportModelIncompressibleTurbulenceModel
Stype: LES
Type: WALE

Finally, this leads to the conclusion that the function object turbulenceFields is linked to turbulence models original definition. It uses processField() to look up registered fields in turbulence model defined in constant/turbulenceProperties. i. e. turbulence fields for WALE model is defined in WALE.C. As for k, it is realized directly in WALE.C:

template<class BasicTurbulenceModel>
 tmp<volSymmTensorField> WALE<BasicTurbulenceModel>::Sd
 (
     const volTensorField& gradU
 ) const
 {
     return dev(symm(gradU & gradU));
 }

$$
S^d_{ij}=\frac12(\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_k}+\frac{\partial u_k}{\partial x_j}\frac{\partial u_j}{\partial x_i})
$$

template<class BasicTurbulenceModel>
tmp<volScalarField> WALE<BasicTurbulenceModel>::k
(
   const volTensorField& gradU
) const
{
volScalarField magSqrSd(magSqr(Sd(gradU)));
return tmp<volScalarField>
(
new volScalarField
(
   IOobject(IOobject::groupName("k", this-> alphaRhoPhi_.group()), this->runTime_.timeName(),this->mesh_),
      
sqr(sqr(Cw_)*this->delta()/Ck_)*(pow3(magSqrSd) / (
sqr(pow(magSqr(symm(gradU)),5.0/2.0)
+pow(magSqrSd,5.0/4.0))
+dimensionedScalar("SMALL",dimensionSet(0, 0, -10, 0, 0),SMALL))
)
));}

$$
\text{magSqr}S^d=(S^d_{ij}S^d_{ij})^{\frac12},
S_{ij}=\frac12(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i})
$$
$$\
k=\left(\frac {C_w^2\Delta}{C_k}\right)^2\frac{(S^d_{ij}S^d_{ij})^{\frac32}}{(S_{ij}S_{ij})^{\frac52}+(S^d_{ij}S^d_{ij})^{\frac54}}
$$

While R is inherited from eddyViscosity::R(), because in WALE.H:

 #include "LESeddyViscosity.H"

then in LESeddyViscosity.H:

#include "eddyViscosity.H"

finally, in eddyViscosity.H:

template<class BasicTurbulenceModel>
 Foam::tmp<Foam::volSymmTensorField>
 Foam::eddyViscosity<BasicTurbulenceModel>::R() const
 {
tmp<volScalarField> tk(k());
// Get list of patchField type names from k
wordList patchFieldTypes(tk().boundaryField().types());
// For k patchField types which do not have an equivalent for symmTensor
// set to calculated
forAll(patchFieldTypes, i){       if(!fvPatchField<symmTensor>::patchConstructorTablePtr_->found(patchFieldTypes[i])){
patchFieldTypes[i] = calculatedFvPatchField<symmTensor>::typeName;}
}
return tmp<volSymmTensorField>
(
new volSymmTensorField(
	IOobject(IOobject::groupName("R", this->alphaRhoPhi_.group()),this->runTime_.timeName(),this->mesh_,IOobject::NO_READ,IOobject::NO_WRITE,false),
	((2.0/3.0)*I)*tk() - (nut_)*dev(twoSymm(fvc::grad(this->U_))),
   patchFieldTypes)
);}

$$
\tau_{ij}=\frac23k\delta_{ij}-2\nu_tS_{ij}
$$

While nut_ is updated in WALE.C:

template<class BasicTurbulenceModel>
 void WALE<BasicTurbulenceModel>::correctNut()
 {
     this->nut_ = Ck_*this->delta()*sqrt(this->k(fvc::grad(this->U_)));
     this->nut_.correctBoundaryConditions();
     fv::options::New(this->mesh_).correct(this->nut_);
  
     BasicTurbulenceModel::correctNut();
 }

$$
\nu_t=C_k\delta\sqrt k
$$

李东岳

感谢分享！！

CFD中文网

对OpenFOAM中 turbulenceFields 和湍流模型如何从模板类具体实现的笔记