我们已经在“OpenFOAM 湍流模型模板框架分析”系列博文中介绍过,OpenFOAM 将不可压缩流和可压缩流的湍流模型代码复用,放在一个统一的框架下。本文将从流体力学理论上分析 RANS 方程代码的具体实现。
动量方程
在“雷诺平均 Navier-Stokes 方程”一文中,我们提到了流体控制方程中的动量方程表示为以下形式
将其写作矢量形式,如下
实际上对于可压缩流体,上式还缺少一项,完整的动量方程如下
上式中的红色部分为可压缩流体多出来的项。 为单位张量(identity tensor)
进行雷诺平均后,得到如下形式的 RANS 方程
其中, 为雷诺应力
将雷诺应力和扩散项相结合,可以写作如下形式
有效雷诺应力
定义有效雷诺应力如下:
那么 RANS 方程可表示为:
上式为 OpenFOAM 采用的形式。
将 分为偏应力张量(deviatoric part)和球应力(或称静水应力)张量(spherical/hydrostatic part)两部分:
注意上式中球应力张量的系数为 2/3,而不是通常的 1/3。
的球应力张量又有以下等式关系:
所以有效雷诺应力这一项可进一步化简为
对于不可压缩流体,实际上方程右端第二项为零,但是 OpenFOAM 为了一致性(consistency)考虑,还是保留了第二项。
这正是动量方程离散代码中 divDevReff
和 divDevRhoReff
名称的由来。
代码实现
不可压求解器中,动量方程的代码如下(以 pimpleFoam 求解器为例):
1
2
3
4
5
6
7
8
|
tmp<fvVectorMatrix> tUEqn
(
fvm::ddt(U) + fvm::div(phi, U)
+ MRF.DDt(U)
+ turbulence->divDevReff(U)
==
fvOptions(U)
);
|
可压缩求解器中,动量方程的代码如下(以 rhoPimpleFoam 求解器为例):
1
2
3
4
5
6
7
8
|
tmp<fvVectorMatrix> tUEqn
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(U)
==
fvOptions(rho, U)
);
|
代码中的 fvm::ddt(U)
和 fvm::ddt(rho, U)
为时间项,fvm::div(phi, U)
为对流项。我们重点看 turbulence->divDevReff(U)
和 turbulence->divDevRhoReff(U)
。
不可压缩版本 divDevReff(U)
的具体实现在 IncompressibleTurbulenceModel
这个类中:
1
2
3
4
5
6
7
8
9
|
template<class TransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::IncompressibleTurbulenceModel<TransportModel>::divDevReff
(
volVectorField& U
) const
{
return divDevRhoReff(U);
}
|
若采用的湍流模型是线性涡粘模型(即派生自 linearViscousStress
),那么其调用的是 linearViscousStress
中的 divDevRhoReff(U)
:
1
2
3
4
5
6
7
8
9
10
11
12
13
|
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)
);
}
|
可压缩版本的 divDevRhoReff(rho, U)
调用的也是 linearViscousStress
中的实现:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
|
template<class BasicTurbulenceModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::linearViscousStress<BasicTurbulenceModel>::divDevRhoReff
(
const volScalarField& rho,
volVectorField& U
) const
{
return
(
- fvc::div((this->alpha_*rho*this->nuEff())*dev2(T(fvc::grad(U))))
- fvm::laplacian(this->alpha_*rho*this->nuEff(), U)
);
}
|
对应 式,其中的 fvm::laplacian
操作返回 ,而 fvc::div
操作返回 。
T(fvc::grad(U))
表示速度梯度的转置。dev2
表示 中的球应力张量,定义如下:
1
2
3
4
5
6
|
//- Return the deviatoric part of a tensor
template<class Cmpt>
inline Tensor<Cmpt> dev2(const Tensor<Cmpt>& t)
{
return t - SphericalTensor<Cmpt>::twoThirdsI*tr(t);
}
|