OpenFOAM 中的单相流湍流模型之一

相信有不少 OpenFOAM 用户有添加湍流模型的需求,我自己最早用 OpenFOAM 完成的一项工作就是在其中添加了一些单相流的湍流模型,并进行了一些计算。这里将我对单相湍流模型代码框架的理解记录下来,供大家参考。本系列将包含三篇,第一篇介绍湍流模型类的继承派生关系,第二篇具体分析几个 OpenFOAM 中带的湍流模型,并给出修改或增加新模型的方法,第三篇分析湍流模型的运行时选择机制(Run Time Selection)的原理。

1. 湍流模型类的继承派生关系

这一部分是最简单的,只要有一点C++的知识,看一下湍流模型的代码头文件的类声明部分,就能理解。OpenFOAM 里的单相湍流模型包含两大类,RASLES,下面将分别分析。

OpenFOAM 单相流湍流模型的代码在 src/turbulenceModels 目录下,目录结构如下:

1
Allwmake  compressible  derivedFvPatchFields  incompressible  LES

其中, compressibleincompressible 分别是可压缩和不可压缩湍流模型的代码, derivedFvPatchFields 是两个湍流相关的边界条件的代码, LES 是大涡模拟的两个相关的类( LESdeltasLESfilters ),具体在后面会介绍。这里我主要关心不可压缩湍流模型。
进入incompressible,目录结构为:

1
2
3
$  cd  incompressible
$ ls
Allwmake LES RAS turbulenceModel

这里, 目录turbulenceModel 里是基类 turbulenceModel 相关的代码, RASLES ,顾名思义,分别是雷诺时均和大涡模拟湍流模型的代码。

1.1 基类 turbulenceModel

首先看基类 turbulenceModel ,这里我挑着我觉得重要的部分代码列出来:

先看头文件 turbulenceModel.H:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
namespace Foam
{

// incompressible 命名空间,注意这个是很重要的,作用是将类隔离开。比如,可压和不可压都有kEpsilon模型,这两个模型的类名是一样的。要区分这两个类,靠的就是这个命名空间。可压的是 compressible::kEpsilon, 而不可压的则是 incompressible::kEpsilon
namespace incompressible
{

class turbulenceModel //定义一个 turbulenceModel 类,继承自 regIOobject 类
:
public regIOobject
{

protected:

// 数据成员
const Time& runTime_;
const fvMesh& mesh_;

const volVectorField& U_;
const surfaceScalarField& phi_;

transportModel& transportModel_; // 输运模型,涉及到分子粘度的设置

//- Near wall distance boundary field
nearWallDist y_;

public:

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


// Declare run-time New selection table 这个跟运行时选择有关,具体以后会涉及

declareRunTimeNewSelectionTable
(
autoPtr,
turbulenceModel,
turbulenceModel,
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName
),
(U, phi, transport, turbulenceModelName)
);


// Constructors 构造函数

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


// Selectors // 这个也是跟运行时选择机制有关,涉及到具体湍流模型的选择过程
//- Return a reference to the selected turbulence model
static autoPtr<turbulenceModel> New
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName = typeName
)
;



// Member Functions

//- Const access to the coefficients dictionary
virtual const dictionary& coeffDict() const = 0;

//- Helper function to return the nam eof the turbulence G field
inline word GName() const
{

return word(type() + ":G");
}

//- Access function to velocity field
inline const volVectorField& U() const
{
return U_;
}

//- Access function to flux field
inline const surfaceScalarField& phi() const
{
return phi_;
}

//- Access function to incompressible transport model
inline transportModel& transport() const
{
return transportModel_;
}

//- Return the near wall distances
const nearWallDist& y() const
{
return y_;
}

//- Return the laminar viscosity
// 分子粘度,注意这里的返回值是输运模型类对象的 nu 函数的返回值。
// 具体来说,如果是牛顿流体,那么这个返回值就是我们在transportProperties里设置的粘度;
// 如果是非牛顿流体,那么粘度是根据具体的非牛顿流体模型计算得到的。
inline tmp<volScalarField> nu() const
{

return transportModel_.nu();
}

// 以下形如 "virtual xxxx () const = 0" 的函数,都是纯虚函数,这是很重要的一部分。
// 当一个类中有纯虚函数时,这个类就被称作抽象类,抽象类本身不能建立对象,一般都是作为接口类来使用。
// 这里的turbulenceModel类就是接口类,“接口类”三个字的具体含义后面会解释。
//- Return the turbulence viscosity
virtual tmp<volScalarField> nut() const = 0;
//- Return the effective viscosity
virtual tmp<volScalarField> nuEff() const = 0;
//- Return the turbulence kinetic energy
virtual tmp<volScalarField> k() const = 0;
//- Return the turbulence kinetic energy dissipation rate
virtual tmp<volScalarField> epsilon() const = 0;
//- Return the Reynolds stress tensor
virtual tmp<volSymmTensorField> R() const = 0;
//- Return the effective stress tensor including the laminar stress
virtual tmp<volSymmTensorField> devReff() const = 0;
//- Return the source term for the momentum equation
virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const = 0;
//- Return the source term for the momentum equation
virtual tmp<fvVectorMatrix> divDevRhoReff
(
const volScalarField& rho,
volVectorField& U
) const
= 0;

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

//- Read LESProperties or RASProperties dictionary
virtual bool read() = 0;
};

再来看 turbulenceModel.C,重点关注构造函数和 New 函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
namespace Foam
{
namespace incompressible
{

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

//这句与运行时选择机制有关,后面再说
defineRunTimeSelectionTable(turbulenceModel, turbulenceModel);

//构造函数的定义。构造函数包括四个参数,其中最后一个"turbulenceModelName"是带缺省参数的,所以,只需要提供三个参数。
turbulenceModel::turbulenceModel
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName
)
// 这里往下一段叫做成员初始化列表,用于对当前类以及其基类成员进行初始化
:
regIOobject
(
IOobject
(
turbulenceModelName,
U.time().constant(),
U.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
)
),
runTime_(U.time()),
mesh_(U.mesh()),

U_(U),
phi_(phi),
transportModel_(transport), // 输运模型从构造函数的参数中读取
y_(mesh_)
{}


// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //

这个函数,起着选择具体湍流模型的作用,后面我会结合求解器代码仔细说说这个函数。更详细的机制,将在结合运行时选择机制来解释。
autoPtr<turbulenceModel> turbulenceModel::New
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName
)
{

// 注意这里了,算例中的"turbulenceProperties" 文件即由这段代码来读取。
// 需要从"turbulenceProperties" 文件中查找一个关键字"simulationType"
// 然后将"simulationType"对应的值赋值给变量"modelType"(对于单相流,modelType 只可能是 "RAS" 或者 "LES")。
const word modelType
(
IOdictionary
(
IOobject
(
"turbulenceProperties",
U.time().constant(),
U.db(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
).lookup("simulationType")
);

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

turbulenceModelConstructorTable::iterator cstrIter =
turbulenceModelConstructorTablePtr_->find(modelType);


if (cstrIter == turbulenceModelConstructorTablePtr_->end())
{
FatalErrorIn
(
"turbulenceModel::New(const volVectorField&, "
"const surfaceScalarField&, transportModel&, const word&)"
) << "Unknown turbulenceModel type "
<< modelType << nl << nl
<< "Valid turbulenceModel types:" << endl
<< turbulenceModelConstructorTablePtr_->sortedToc()

<< exit(FatalError);
}

return autoPtr<turbulenceModel>

(
cstrIter()(U, phi, transport, turbulenceModelName)
);
}

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

void turbulenceModel::correct()
{
transportModel_.correct();

if (mesh_.changing())
{
y_.correct();
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

从基类 turbulenceModel 以下,就要花开两朵,各表一枝了。先来看 RAS 类的湍流模型。

1.2 RAS 模型

我们从头文件 RASModel.H 看起

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
namespace Foam
{
namespace incompressible
{

/*---------------------------------------------------------------------------*\
Class RASModel Declaration
\*---------------------------------------------------------------------------*/


class RASModel
:
public turbulenceModel, // RASModel 类是前面讲的turbulenceModel类的派生类
public IOdictionary // 同时也继承 IOdictionary 类
{

protected:

// 两个开关变量,类似于 C++ 中的bool变量,其值要么是true,要么是false。
// 只是这里将 true 换成 on,false 换成 off 也是一样的。
Switch turbulence_;
Switch printCoeffs_;

//- Model coefficients dictionary
dictionary coeffDict_;

//- Lower limit of k
dimensionedScalar kMin_;

//- Lower limit of epsilon
dimensionedScalar epsilonMin_;

//- Lower limit for omega
dimensionedScalar omegaMin_;

// Protected Member Functions

//- Print model coefficients
virtual void printCoeffs();

private:
// Private Member Functions
//- Disallow default bitwise copy construct
RASModel(const RASModel&);
//- Disallow default bitwise assignment
void operator=(const RASModel&);

public:

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

// 这里也涉及到 运行时选择机制,以后一起讲。
declareRunTimeSelectionTable
(
autoPtr,
RASModel,
dictionary,
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName
),
(U, phi, transport, turbulenceModelName)
);

// 构造函数
RASModel
(
const word& type,
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName = turbulenceModel::typeName
);

//- Return a reference to the selected RAS model 这个函数作为运行时选择机制里的选择器。
static autoPtr<RASModel> New
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName = turbulenceModel::typeName
);


// Member Functions

// Access

// 注意,一下这些成员函数,是将来 RASModel 类的派生类们可能会用到的,这里一起定义了,
// 后面的派生类将可以直接继承这些函数,这样达到了代码重复利用的作用。
//- Return the lower allowable limit for k (default: SMALL)
const dimensionedScalar& kMin() const
{
return kMin_;
}

//- Return the lower allowable limit for epsilon (default: SMALL)
const dimensionedScalar& epsilonMin() const
{
return epsilonMin_;
}

//- Return the lower allowable limit for omega (default: SMALL)
const dimensionedScalar& omegaMin() const
{
return omegaMin_;
}

//- Allow kMin to be changed
dimensionedScalar& kMin()
{
return kMin_;
}

//- Allow epsilonMin to be changed
dimensionedScalar& epsilonMin()
{
return epsilonMin_;
}

//- Allow omegaMin to be changed
dimensionedScalar& omegaMin()
{
return omegaMin_;
}

//- Const access to the coefficients dictionary
virtual const dictionary& coeffDict() const
{
return coeffDict_;
}


//- Return the effective viscosity
// 这个函数很重要,默认情况下,雷诺时均湍流模型中,有效粘度等于湍流粘度加层流粘度,即nut + nu。
// 这里的 nut 和 nu 两个函数在基类 turbulenceModel 声明了,请往前翻以证实这一点。
virtual tmp<volScalarField> nuEff() const
{
return tmp<volScalarField>
(
new volScalarField("nuEff", nut() + nu())
);
}

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

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

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

再来看 RASModel.C,这里跟 turbulenceModel.C 类似,重点关注构造函数和 选择器(Selector)函数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
namespace Foam
{
namespace incompressible
{

defineTypeNameAndDebug(RASModel, 0);
defineRunTimeSelectionTable(RASModel, dictionary);
addToRunTimeSelectionTable(turbulenceModel, RASModel, turbulenceModel);

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

void RASModel::printCoeffs()
{
if (printCoeffs_)
{
Info<< type() << "Coeffs" << coeffDict_ << endl;
}
}

// 构造函数
RASModel::RASModel
(
const word& type,
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName
)
:
类似的,这里也是使用成员初始化列表来初始化当前类及其父类的数据成员

// 这一句是传递参数给父类 turbulenceModel,注意这里传给父类构造函数的参数要与父类中的构造函数的参数表一致哦
turbulenceModel(U, phi, transport, turbulenceModelName),

//这里是建立一个IOobject来初始化另一个父类IOdictionary,注意这里建立的是一个名为“RASProperties”的IO对象,这个文件相信用过 OpenFOAM 湍流模拟的一定很熟悉吧
IOdictionary
(
IOobject
(
"RASProperties",
U.time().constant(), // 这里表示文件的位置是在constant文件夹下
U.db(),
IOobject::MUST_READ_IF_MODIFIED, //这里的意思是,如果修改了,则需要重新读取,所以,如果你在算例运行时修改了这个文件,你的修改会即时生效的
IOobject::NO_WRITE
)
),

// 查找“turbulence” 关键字,并用其对应的值来初始化变量“turbulnce_”。
turbulence_(lookup("turbulence")),
// 同上,区别是这里是带缺省参数的,也就是说如果找不到“printCoeffs” 就用缺省值“false
printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),

// 这里是读取模型参数字典来初始化变量“coeffDict_”。注意,这里的 type 是啥?这个要等到进一步看一个具体的湍流模型类时才能明了。
// 另外,注意这里的 `lookup` , `lookupOrDefault` 和 `subOrEmptyDict` 三个函数都是继承自 IOdictionary 类的函数。
coeffDict_(subOrEmptyDict(type + "Coeffs")),

kMin_("kMin", sqr(dimVelocity), SMALL),
epsilonMin_("epsilonMin", kMin_.dimensions()/dimTime, SMALL),
omegaMin_("omegaMin", dimless/dimTime, SMALL)
{

// 这里是从 RASProperties 里读取kMin等的值,注意如果 RASProperties 里不设置那就使用上面初始化的值“SMALL”。
// 注意这里的“*this”代表的是 RASModel 类本身,但是readIfPresent函数的参数应该是 dictionary 类,由于RASModel 类继承自 IODictionary 类,所以,这里其实隐含了一个派生类指针向基类指针的转换
kMin_.readIfPresent(*this);
epsilonMin_.readIfPresent(*this);
omegaMin_.readIfPresent(*this);

// Force the construction of the mesh deltaCoeffs which may be needed
// for the construction of the derived models and BCs
mesh_.deltaCoeffs();
}


//这个也是跟运行时选择有关,以后会细说
autoPtr<RASModel> RASModel::New
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName
)
{
// get model name, but do not register the dictionary
// otherwise it is registered in the database twice
const word modelType
(
IOdictionary
(
IOobject
(
"RASProperties",
U.time().constant(),
U.db(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
).lookup("RASModel")
);

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

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

if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorIn
(
"RASModel::New"
"("
"const volVectorField&, "
"const surfaceScalarField&, "
"transportModel&, "
"const word&"
")"
) << "Unknown RASModel type "
<< modelType << nl << nl
<< "Valid RASModel types:" << endl
<< dictionaryConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}

return autoPtr<RASModel>
(
cstrIter()(U, phi, transport, turbulenceModelName)
);
}


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

void RASModel::correct()
{
turbulenceModel::correct();
}

bool RASModel::read()
{
//if (regIOobject::read())

// Bit of trickery : we are both IOdictionary ('RASProperties') and
// an regIOobject from the turbulenceModel level. Problem is to distinguish
// between the two - we only want to reread the IOdictionary.

bool ok = IOdictionary::readData
(
IOdictionary::readStream
(
IOdictionary::type()
)
);
IOdictionary::close();

if (ok)
{
lookup("turbulence") >> turbulence_;

if (const dictionary* dictPtr = subDictPtr(type() + "Coeffs"))
{
coeffDict_ <<= *dictPtr;
}
kMin_.readIfPresent(*this);
epsilonMin_.readIfPresent(*this);
omegaMin_.readIfPresent(*this);
return true;
}
else
{
return false;
}
}

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

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

1.3 LES 模型

LESModel 类的结构与 RASModel 非常接近,所以这里就简单提一下区别之处。

  • delta_ 成员
    1
    autoPtr<Foam::LESdelta> delta_;

这个成员是 LESdelta 类的对象,定义滤波尺度。这个类的定义在 src/turbulenceModels/LES/LESdeltas ,当中定义了几个可选的 delta 模型

  • 亚格子粘度和应力
1
2
3
4
5
6
7
8
9
10
11
12
13
14
//- Return the SGS viscosity. 亚格子粘度
virtual tmp<volScalarField> nuSgs() const = 0;

//- Return the effective viscosity
virtual tmp<volScalarField> nuEff() const //有效粘度等于亚格子粘度与分子粘度之和
{

return tmp<volScalarField>
(
new volScalarField("nuEff", nuSgs() + nu())
);
}

//- Return the sub-grid stress tensor. // 亚格子应力
virtual tmp<volSymmTensorField> B() const = 0;

LESModel.C 的结构与 RASModel.C 几乎一样,所以这里就不重复了。

前面提到,基类 turbulenceModel 里声明了很多纯虚函数,所以,turbulenceModel 类是抽象类,不能直接创建对象。注意这里的 RASModel 类和 LESModel 类由于继承了 turbulenceModel 类的纯虚函数,所以这两个依然是抽象类。这一点在后面解析具体湍流模型类的时候还会提到,这里先提个醒。

1.4 继承派生关系

前面看完了基类 turbulenceModel,RASModel 以及 LESModel,可以看出这三个类的继承派生关系为:

RAS 和 LES 与 turbulenceModel 的继承关系

下面继续看具体湍流模型类与基类的继承关系。
先看雷诺时均类湍流模型:
位于 src/turbulenceModels/incompressible/RAS 目录下的都是雷诺时均类的湍流模型,这个类型的继承关系很简单:所有具体湍流模型类都继承自 RASModel 类,关系示意图如下:

RAS 类湍流模型的继承关系

位于 src/turbulenceModels/incompressible/LES 目录下的是大涡模拟类型的湍流模型,这类湍流模型的继承关系比雷诺时均类的要复杂一点,但是也不难捋清,过程我就不详述了,下面给出我整理的继承关系图:

LES 类湍流模型的继承关系

上图中,虚线框里的是抽象类,实线框里的是具体的可以调用的湍流模型类。

1.5 求解器中湍流模型的调用

最后简单提一下求解器中调用湍流模型的接口。以 pisoFoam 求解器为例:
纵观 pisoFoam 的代码,跟湍流模型有关的共有三处,第一处是在 createField.H 的最后面,创建一个智能指针

1
2
3
4
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
)
;

注意,这个指针的类型是 incompressible::turbulenceModel ,也就是说,创建的是基类的指针。
第二处位于 UEqn 中,

1
+ turbulence->divDevReff(U)

从这里可以看出,这是用指针 turbulence 调用成员函数 divDevReff
第三处在压力修正之后

1
turbulence->correct();

这里是调用成员函数 correct
这简单的几行代码,就完成了湍流模型的调用。这里先大致说一下调用的原理和过程,详细的留待后面跟运行时选择机制一起说。
首先,注意刚才提到的智能指针 turbulence 的类型是 incompressible::turbulenceModel ,是基类类型的。这里就不得不说一下 C++ 虚函数的作用了。还记得上面提到的基类 turbulenceModel 中声明的那些纯虚函数吧,如果你往上翻翻,你会发现, divDevReffcorrect 在基类 turbulenceModel 中都被声明为纯虚函数。这里只要把握两点,就能理解湍流模型的调用原理:

  1. 基类类型的指针可以指向派生类的对象;
  2. 在基类中声明的纯虚函数可以在派生类中进行定义,当基类类型指针指向派生类对象以后,用这个指针调用成员函数时,实际调用的是指针指向的派生类中定义的函数。

把握这两点,然后再去理解湍流模型的调用过程:
首先是调用 incompressible::turbulenceModel::New 函数来初始化指针 turbulence ,查看上面 turbulenceModel 类中 New 函数的定义,可以知道,函数要先从“turbulenceProperties” 文件里 simulationType 关键字,从而决定是调用 RAS 模型还是 LES 模型。如果用户设定的是 RAS ,那么 turbulenceModel 类的 New 函数将返回一个指向 RASModel 类的指针,这个指针将继续调用 RASModel 类的 New 函数,并在这个函数中读取 RASProperties 文件,查找关键字 RASModel ,从而决定具体调用的湍流模型。假设用户指定的是 kEpsilon 那么最终 createField.H 中定义的指针 turbulence 将指向一个 kEpsilon 类的对象,由此, turbulence->divDevReff(U)turbulence->correct() 都将分别调用定义在 kEpsilon 类中的成员函数 divDevReffcorrect
一个简单的湍流模型调用示意图如下:

湍流模型调用示意图

以上是 RAS 类型模型的调用,LES 类型的基本上差不多,但是,从上面的继承关系图也能看出,LES 模型类的结构更复杂一点。根本原因在于,LES 模型不像 RAS 模型那样都是同一的对湍流粘度 nut 建模,而是有一部分是对亚格子粘度 nuSgs 建模(这一部分湍流模型均继承自 GenEddyVisc ),还有另一部分是直接对亚格子应力 B 建模(这一部分湍流模型均继承自 GenSGSStress ),此外,分离涡模型(Detached eddy model, DESModel)也放在这个目录下,而且还有一个直接继承自 LESModel 的模型 kOmegaSSTSAS (这个模型与 LESModel 的关系就跟 kEpsilonRASModel 的关系一样)。虽然有这么多种类型,但是调用过程其实跟 RAS 类型的是差不多的。如果从字典文件“turbulenceProperties” 里读到的 simulationTypeLES 的话,那么将继续从字典文件 LESProperties 里读取具体的 LES 模型。上面总结的那张继承关系图中,所有的实线框中的模型都可以选择。但是,具体的模型需要具体的设置,比如,需要设置滤波尺度 delta 模型,还可能需要设置 filter 模型,具体的要求可以去具体的湍流模型类的代码中去看。如果你了解一点你需要使用的湍流模型的基本理论,能写出模型的方程,那要看懂这个湍流模型在 OpenFOAM 中实现的代码是很容易的。

上述关于调用过程的叙述,只是我的理解,其实不严谨,但大致原理应该是这样。计划中这个系列将写三篇,我将在第三篇中叙述运行时选择机制,到时候还会深入说说这个调用过程。

P.S:本系列在筹划时,OpenFOAM-3.0 版本还没发布。随着 3.0 版本的发布,本系列里对湍流模型的描述已经“过时”了,因为在 3.0 版中,湍流模型类进行了重新模板化,将单相湍流和多相湍流模型整合在一起了,所以这里的描述只适用于 3.0 以下的版本。