我们已经在“OpenFOAM 湍流模型模板框架分析”系列博文中介绍过,OpenFOAM 将不可压缩流和可压缩流的湍流模型代码复用,放在一个统一的框架下。本文将从流体力学理论上分析 RANS 方程代码的具体实现。

动量方程

在“雷诺平均 Navier-Stokes 方程”一文中,我们提到了流体控制方程中的动量方程表示为以下形式

(1)(ρui)t时间项+(ρujui)xj对流项xj[μ(uixj+ujxi)]粘性项=pxi

将其写作矢量形式,如下

(2)(ρu)t+(ρuu){μ[u+(u)T]}=p

实际上对于可压缩流体,上式还缺少一项,完整的动量方程如下

(3)(ρu)t+(ρuu){μ[u+(u)T23(u)I]}=p

上式中的红色部分为可压缩流体多出来的项。I 为单位张量(identity tensor)

(4)I=[100010001]

进行雷诺平均后,得到如下形式的 RANS 方程

(5)(ρu¯)t+(ρu¯u¯){μ[u¯+(u¯)T23(u¯)I]+ρR}=p¯

其中,R 为雷诺应力

(6)R=νt[u¯+(u¯)T23(u¯)I]23kI

将雷诺应力和扩散项相结合,可以写作如下形式

(7)(ρu¯)t+(ρu¯u¯){μeff[u¯+(u¯)T23(u¯)I]}=(p¯+23ρk)

有效雷诺应力

定义有效雷诺应力如下:

(8)Reff=νeff[u¯+(u¯)T23(u¯)I]

那么 RANS 方程可表示为:

(9)(ρu¯)t+(ρu¯u¯)(ρReff)=(p¯+23ρk)

上式为 OpenFOAM 采用的形式。

(u¯)T 分为偏应力张量(deviatoric part)和球应力(或称静水应力)张量(spherical/hydrostatic part)两部分:

(10)(u¯)T=dev2((u¯)T)+hyd2((u¯)T)=(u¯)T23tr((u¯)T)I偏应力张量+23tr((u¯)T)I球应力张量

注意上式中球应力张量的系数为 2/3,而不是通常的 1/3。

(u¯)T 的球应力张量又有以下等式关系:

(11)tr((u¯)T)=ux+vy+wz=u¯

所以有效雷诺应力这一项可进一步化简为

(12)Reff={νeff[u¯+((u¯)T23(u¯)I)]}=(νeffu¯)+[νeffdev2((u¯)T)]

对于不可压缩流体,实际上方程右端第二项为零,但是 OpenFOAM 为了一致性(consistency)考虑,还是保留了第二项。

这正是动量方程离散代码中 divDevReffdivDevRhoReff 名称的由来。

代码实现

不可压求解器中,动量方程的代码如下(以 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)
    );
}

对应 (12) 式,其中的 fvm::laplacian 操作返回 (νeffu¯),而 fvc::div 操作返回 [νeffdev2((u¯)T)]

T(fvc::grad(U)) 表示速度梯度的转置。dev2 表示 (10) 中的球应力张量,定义如下:

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);
}