OpenFOAM 中双流体模型的 kineticTheoryModel 是以 “ Derivation, implementation, and validation of computer simulation models for gas-solid fluidized beds, B.G.M. van Wachem, Ph.D. Thesis, Delft University of Technology, Amsterdam, 2000. “ 为蓝本来设计的，下面分析这个类的代码。需要注意的是，这个类需要调用一些别的类（如 viscosityModel 等，后面会一一分析）来完成其功能。
1. 头文件 kineticTheoryModel.H
头文件中要注意的是 kineticTheoryModel 类的继承关系，以及六个子类的智能指针作为 kineticTheoryModel 类的数据成员。
1 | class kineticTheoryModel |
2. 构造函数
注意这里向基类传递的参数，这里的湍流类的继承关系比较复杂，比单相湍流复杂很多，后面会有具体的分析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
64Foam::RASModels::kineticTheoryModel::kineticTheoryModel
(
const volScalarField& alpha,
const volScalarField& rho,
const volVectorField& U,
const surfaceScalarField& alphaRhoPhi,
const surfaceScalarField& phi,
const transportModel& phase,
const word& propertiesName,
const word& type
)
:
eddyViscosity<RASModel<PhaseCompressibleTurbulenceModel<phaseModel> > >
(
type,
alpha,
rho,
U,
alphaRhoPhi,
phi,
phase,
propertiesName
),
phase_(phase),
viscosityModel_
(
kineticTheoryModels::viscosityModel::New
(
this->coeffDict_ // coeffDict_ 是 RASModel的成员
)
),
conductivityModel_
(
kineticTheoryModels::conductivityModel::New
(
this->coeffDict_
)
),
radialModel_
(
kineticTheoryModels::radialModel::New
(
this->coeffDict_
)
),
granularPressureModel_
(
kineticTheoryModels::granularPressureModel::New
(
this->coeffDict_
)
),
frictionalStressModel_
(
kineticTheoryModels::frictionalStressModel::New
(
this->coeffDict_
)
),
......
......
3. 主要成员函数
1 | // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // |
k 和 epsilon 是两个从基类继承下来的函数，这两个函数在基类中是纯虚函数，所以虽然 kineticTheoryModel 用不到它们，但是还是需要进行定义，否则 kineticTheoryModel 类就将是一个虚基类，从而无法创建对象了。
1 | Foam::tmp<Foam::volSymmTensorField> |
R 函数返回所谓的雷诺应力，其实跟单相湍流类的 R 函数返回值有点差别：
$$
-\nu_t[\nabla U + \nabla U^T-\frac{2}{3}(\nabla \cdot U)\cdot \mathrm{I}] - \lambda(\nabla \cdot U)\cdot \mathrm{I}
$$
这里的 R 与下面的 devRoReff 函数返回值只相差一个 rho。
而单相湍流模型的 R 则为：
$$
\frac{2}{3}k\cdot \mathrm{I}- \nu_t (\nabla U + \nabla U^T)
$$
1 |
|
pPrime 和 pPrimef 两个函数，返回的是固相压力对固相孔隙率的导数（$\partial P_s/\partial \varepsilon_s$）。
两个在”UEqn.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
40Foam::tmp<Foam::volSymmTensorField>
Foam::RASModels::kineticTheoryModel::devRhoReff() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
IOobject::groupName("devRhoReff", this->U_.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
- (this->rho_*this->nut_)
*dev(twoSymm(fvc::grad(this->U_)))
- ((this->rho_*lambda_)*fvc::div(this->phi_))*symmTensor::I
)
);
}
Foam::tmp<Foam::fvVectorMatrix>
Foam::RASModels::kineticTheoryModel::divDevRhoReff
(
volVectorField& U
) const
{
return
(
- fvm::laplacian(this->rho_*this->nut_, U)
- fvc::div
(
(this->rho_*this->nut_)*dev2(T(fvc::grad(U)))
+ ((this->rho_*lambda_)*fvc::div(this->phi_))
*dimensioned<symmTensor>("I", dimless, symmTensor::I)
)
);
}
对应的公式分别为：
devRhoReff
$$
- \rho \nu_t \cdot dev(\nabla U + \nabla U^T) - \rho \lambda(\nabla \cdot U)\cdot \mathrm{I} \\
= - \rho \nu_t(\nabla U + \nabla U^T) + \rho \nu_t \frac{2}{3}(\nabla \cdot U)\cdot \mathrm{I} - \rho \lambda (\nabla \cdot U)\cdot \mathrm{I}
$$
divDevRhoReff
$$
\ - \nabla \cdot (\rho \nu_t \nabla U) - \nabla \cdot \left [\rho \nu_t \nabla U^T\ - \rho \nu_t \frac{2}{3}(\nabla \cdot U)\cdot \mathrm{I} \right ] + \rho \lambda (\nabla \cdot U) \cdot \mathrm{I} )
$$
correct是计算颗粒温度 Theta 的函数，这是 kineticTheoryModel 类最重要的部分。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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204void Foam::RASModels::kineticTheoryModel::correct()
{
// Local references
volScalarField alpha(max(this->alpha_, scalar(0)));
const volScalarField& rho = phase_.rho();
const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
const volVectorField& U = this->U_; // 当前相的速度
const volVectorField& Uc_ = phase_.fluid().otherPhase(phase_).U(); // 另一相的速度
const scalar sqrtPi = sqrt(constant::mathematical::pi);
dimensionedScalar ThetaSmall("ThetaSmall", Theta_.dimensions(), 1.0e-6);
dimensionedScalar ThetaSmallSqrt(sqrt(ThetaSmall));
tmp<volScalarField> tda(phase_.d()); // 颗粒直径
const volScalarField& da = tda();
tmp<volTensorField> tgradU(fvc::grad(this->U_));
const volTensorField& gradU(tgradU());
volSymmTensorField D(symm(gradU));
// Calculating the radial distribution function
// 调用 radialModel 类 来计算径向分布函数
gs0_ = radialModel_->g0(alpha, alphaMinFriction_, alphaMax_);
if (!equilibrium_) // 如果 equilibrium_ = off ，那么将用这个偏微分方程来计算颗粒温度，否则将改用一个代数方程，见下文
{
// Particle viscosity (Table 3.2, p.47)
// 调用viscosityModel 来更新颗粒相的粘度
nut_ = viscosityModel_->nu(alpha, Theta_, gs0_, rho, da, e_);
volScalarField ThetaSqrt(sqrt(Theta_));
// Bulk viscosity p. 45 (Lun et al. 1984).
lambda_ = (4.0/3.0)*sqr(alpha)*da*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi;
// Stress tensor, Definitions, Table 3.1, p. 43
volSymmTensorField tau // 颗粒应力
(
rho*(2.0*nut_*D + (lambda_ - (2.0/3.0)*nut_)*tr(D)*I)
);
// Dissipation (Eq. 3.24, p.50)
volScalarField gammaCoeff //颗粒动能耗散项
(
12.0*(1.0 - sqr(e_))
*max(sqr(alpha), residualAlpha_)
*rho*gs0_*(1.0/da)*ThetaSqrt/sqrtPi
);
// Drag // 调用曳力模型计算曳力系数
volScalarField beta(phase_.fluid().drag(phase_).K());
// Eq. 3.25, p. 50 Js = J1 - J2
volScalarField J1(3.0*beta);
volScalarField J2
(
0.25*sqr(beta)*da*magSqr(U - Uc_)
/(
max(alpha, residualAlpha_)*rho
*sqrtPi*(ThetaSqrt + ThetaSmallSqrt)
)
);
// particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45)
volScalarField PsCoeff
(
granularPressureModel_->granularPressureCoeff
(
alpha,
gs0_,
rho,
e_
)
);
// 'thermal' conductivity (Table 3.3, p. 49)
// 调用 conductivityModel 计算颗粒脉动能量的传导系数
kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_);
// Construct the granular temperature equation (Eq. 3.20, p. 44)
// NB. note that there are two typos in Eq. 3.20:
// Ps should be without grad
// the laplacian has the wrong sign
// 构建颗粒温度方程，注意，开头提到的文献里的颗粒温度方程有两处 typo，
//下面的代码修复了这两处错误，后文会给出正确的公式。
fvScalarMatrix ThetaEqn
(
1.5*
(
fvm::ddt(alpha, rho, Theta_)
+ fvm::div(alphaRhoPhi, Theta_)
- fvc::Sp(fvc::ddt(alpha, rho) + fvc::div(alphaRhoPhi), Theta_)
)
- fvm::laplacian(kappa_, Theta_, "laplacian(kappa,Theta)")
==
fvm::SuSp(-((PsCoeff*I) && gradU), Theta_)
+ (tau && gradU)
+ fvm::Sp(-gammaCoeff, Theta_)
+ fvm::Sp(-J1, Theta_)
+ fvm::Sp(J2/(Theta_ + ThetaSmall), Theta_)
);
ThetaEqn.relax();
ThetaEqn.solve();
}
else // 如果 equilibrium = on， 将使用一个代数方程来计算颗粒温度。
{
// Equilibrium => dissipation == production
// Eq. 4.14, p.82
volScalarField K1(2.0*(1.0 + e_)*rho*gs0_);
volScalarField K3
(
0.5*da*rho*
(
(sqrtPi/(3.0*(3.0 - e_)))
*(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha*gs0_)
+1.6*alpha*gs0_*(1.0 + e_)/sqrtPi
)
);
volScalarField K2
(
4.0*da*rho*(1.0 + e_)*alpha*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0
);
volScalarField K4(12.0*(1.0 - sqr(e_))*rho*gs0_/(da*sqrtPi));
volScalarField trD
(
alpha/(alpha + residualAlpha_)
*fvc::div(this->phi_)
);
volScalarField tr2D(sqr(trD));
volScalarField trD2(tr(D & D));
volScalarField t1(K1*alpha + rho);
volScalarField l1(-t1*trD);
volScalarField l2(sqr(t1)*tr2D);
volScalarField l3
(
4.0
*K4
*alpha
*(2.0*K3*trD2 + K2*tr2D)
);
Theta_ = sqr
(
(l1 + sqrt(l2 + l3))
/(2.0*max(alpha, residualAlpha_)*K4)
);
kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_);
}
// 限定 颗粒温度的上下限。
// max 和 min 函数的定义，没有找到。经验证， max 的作用是让小于0的归零， min 是让大于100的等于100。
Theta_.max(0);
Theta_.min(100);
{
//利用先得到的颗粒温度更新颗粒相的粘度
// particle viscosity (Table 3.2, p.47)
nut_ = viscosityModel_->nu(alpha, Theta_, gs0_, rho, da, e_);
volScalarField ThetaSqrt(sqrt(Theta_));
// Bulk viscosity p. 45 (Lun et al. 1984).
lambda_ = (4.0/3.0)*sqr(alpha)*da*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi;
// Frictional pressure // 计算由于颗粒之间的摩擦作用产生的一个等效的颗粒相压力作用。
volScalarField pf
(
frictionalStressModel_->frictionalPressure
(
alpha,
alphaMinFriction_,
alphaMax_
)
);
// Add frictional shear viscosity, Eq. 3.30, p. 52
// 将颗粒摩擦产生的颗粒相粘度加到由颗粒温度计算得到的颗粒相粘度中，作为总的颗粒相粘度
nut_ += frictionalStressModel_->nu
(
alpha,
alphaMax_,
pf/rho,
D
);
// Limit viscosity
// 限定颗粒相粘度的上限
nut_.min(100);
}
if (debug)
{
Info<< typeName << ':' << nl
<< " max(Theta) = " << max(Theta_).value() << nl
<< " max(nut) = " << max(nut_).value() << endl;
}
}
这里重点关注一下 ThetaEqn 的写法。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16fvScalarMatrix ThetaEqn
(
1.5*
(
fvm::ddt(alpha, rho, Theta_)
+ fvm::div(alphaRhoPhi, Theta_)
- fvc::Sp(fvc::ddt(alpha, rho) + fvc::div(alphaRhoPhi), Theta_)
)
- fvm::laplacian(kappa_, Theta_, "laplacian(kappa,Theta)")
==
fvm::SuSp(-((PsCoeff*I) && gradU), Theta_)
+ (tau && gradU)
+ fvm::Sp(-gammaCoeff, Theta_)
+ fvm::Sp(-J1, Theta_)
+ fvm::Sp(J2/(Theta_ + ThetaSmall), Theta_)
);
对应的偏微分方程是，
$$
\frac{3}{2}\left [ \frac{\partial }{\partial t}(\varepsilon_s \rho_s \Theta)+ \nabla \cdot (\varepsilon_s \rho_s \Theta U_s) \right ] = (-P_s \mathrm{I} + \tau_s):\nabla U_s + \nabla \cdot (\kappa_s \nabla \Theta) - \gamma_s - J_s
$$
下面将公式与代码一一对应。代码里跟公式相比，多了一项
$$
-\Theta\frac{\partial }{\partial t}(\varepsilon_s \rho_s) - \Theta \nabla \cdot (\varepsilon_s \rho_s U_s)
$$
这样，相当于代码对应的前两项是
$$
\varepsilon_s \rho_s \frac{\partial \Theta}{\partial t} + \varepsilon_s \rho_s U_s \cdot \nabla \Theta
$$
这样做的目的仍不是很明确。
Laplacian 项不需多言，剩下的几项，全都当作了源项来处理。
$(-P_s \mathrm{I} + \tau_s)\:\nabla U_s$ 拆开成了两项，分别对应 fvm::SuSp(-((PsCoeff*I) && gradU), Theta_) 和 (tau && gradU) 。第一项，由于$P_s$是$\Theta$ 的函数，所以，在固相压力类中，返回的值是固相压力系数（固相压力除以颗粒温度），在这里将$\Theta$进行了隐式处理。而由于$\tau$与$\nabla U_s$与$\Theta$无关，所以就只当成一般的源项了。
$\gamma$项的处理与颗粒压力类似，也是只定义 gammaCoeff ，然后将$\Theta$作隐式处理。由于实际上$\gamma$的表达式
$$
\gamma = \frac{12(1-e^2)\varepsilon_s^2\rho g_0}{d_p \sqrt{\pi}} \Theta^{3/2}
$$
这里有关于$\Theta$的非线性项，程序里实际上是对$\Theta$作了部分隐式处理，从 gammaCoeff的定义可以看出来：
$$
\gamma_{Coeff} = \frac{12(1-e^2)\varepsilon_s^2\rho g_0}{d_p } \sqrt{\frac{\Theta}{\pi}}
$$
颗粒温度方程中的$J_s$ 也分成了两项来处理，$J_s = J_1 - J_2$
$J_1$的处理很简单，在公式里，$J_1 = 3\beta \Theta$，而在代码当中，变量 J1 定义为 3*beta，而 Theta 则作隐式处理。
$J_2$也人为作了隐式处理，文献中$J_2$的表达式为
$$
J_2 = \frac{\beta^2 d_p |U_g-U_s|^2}{4\varepsilon_s\rho\sqrt{\pi \Theta}}
$$
程序里则将$J_2$除以$\Theta$以后作为系数，以实现对$\Theta$进行隐式处理。
注意这里要说明一下”Sp”和”SuSp”的区别，这个在”Programmer’s Guide”里有说明，此外这个网页里也有说明。可是，为什么固相压力项要用”SuSp”，这个也暂时不明白。
计算颗粒温度的代数方程的具体公式这里不写了，可以参考文献。
最后，上面提到的由于颗粒摩擦作用产生的颗粒压力和颗粒粘度，只在颗粒体积分率很大的区域才需要启用，算例中需要设定一个值，对于颗粒相体积分率大于一个设定值时， alphaMinFriction，只有颗粒相体积分率超过这个值时，才会启用由于颗粒相摩擦而产生的压力和粘度。