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

动量方程

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

$$ \begin{equation} \underbrace{\frac {\partial (\rho u_i)}{\partial t}}_ % workaround for hugo {\text{时间项}} + \underbrace{\frac{\partial (\rho u_j u_i)}{\partial x_j}}_ % workaround for hugo {\text{对流项}} - \underbrace{\frac{\partial}{\partial x_j} \left[ \mu \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \right]}_ % workaround for hugo {\text{粘性项}} = - \frac{\partial p}{\partial x_i} \end{equation} $$

将其写作矢量形式，如下

$$ \begin{equation} \frac {\partial (\rho \mathbf u)}{\partial t} + \nabla \cdot (\rho \mathbf u \mathbf u) - \nabla \cdot \left\lbrace\mu \left[ \nabla \mathbf u + (\nabla \mathbf u)^\mathsf{T} \right] \right\rbrace = - \nabla p \end{equation} $$

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

$$ \begin{equation} \frac {\partial (\rho \mathbf u)}{\partial t} + \nabla \cdot (\rho \mathbf u \mathbf u) - \nabla \cdot \left\lbrace \mu \left[ \nabla \mathbf u + (\nabla \mathbf u)^\mathsf{T} \color{red}{- \frac{2}{3} (\nabla \cdot \mathbf u}) \mathbf I \right] \right\rbrace = - \nabla p \end{equation} $$

上式中的红色部分为可压缩流体多出来的项。$\mathbf I$ 为单位张量（identity tensor）

$$ \begin{equation} \mathbf I = \left[ \begin{array}{ccc} 1 & 0 & 0 \newline 0 & 1 & 0 \newline 0 & 0 & 1 \newline \end{array} \right] \end{equation} $$

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

$$ \begin{equation} \frac {\partial (\rho \overline{\mathbf u})}{\partial t} + \nabla \cdot (\rho \overline{\mathbf u} \overline{\mathbf u}) - \nabla \cdot \left\lbrace \mu \left[ \nabla \overline{\mathbf u} + (\nabla \overline{\mathbf u})^\mathsf{T} - \frac{2}{3} \left( \nabla \cdot \overline{\mathbf u} \right) \mathbf I \right] + \rho \mathbf R \right\rbrace = - \nabla \overline{p} \end{equation} $$

其中，$\mathbf R$ 为雷诺应力

$$ \begin{equation} \mathbf R = \nu_t \left[ \nabla \overline{\mathbf u} + (\nabla \overline{\mathbf u})^\mathsf{T} - \frac{2}{3} (\nabla \cdot \overline{\mathbf u}) \mathbf I \right] - \frac{2}{3} k \mathbf I \end{equation} $$

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

$$ \begin{equation} \frac {\partial (\rho \overline{\mathbf u})}{\partial t} + \nabla \cdot (\rho \overline{\mathbf u} \overline{\mathbf u}) - \nabla \cdot \left\lbrace \mu_\mathit{eff} \left[ \nabla \overline{\mathbf u} + (\nabla \overline{\mathbf u})^\mathsf{T} - \frac{2}{3} (\nabla \cdot \overline{\mathbf u}) \mathbf I \right] \right\rbrace = - \nabla \left( \overline{p} + \frac{2}{3} \rho k\right) \end{equation} $$

有效雷诺应力

定义有效雷诺应力如下：

$$ \begin{equation} \mathbf R_\mathrm{eff} = \nu_\mathit{eff} \left[ \nabla \overline{\mathbf u} + (\nabla \overline{\mathbf u})^\mathsf{T} - \frac{2}{3} (\nabla \cdot \overline{\mathbf u}) \mathbf I \right] \end{equation} $$

那么 RANS 方程可表示为：

$$ \begin{equation} \frac {\partial (\rho \overline{\mathbf u})}{\partial t} + \nabla \cdot (\rho \overline{\mathbf u} \overline{\mathbf u}) - \nabla \cdot (\rho \mathbf R_\mathrm{eff}) = - \nabla \left( \overline{p} + \frac{2}{3} \rho k\right) \label{eq:ranse} \end{equation} $$

上式为 OpenFOAM 采用的形式。

将 $(\nabla \overline{\mathbf u})^T$ 分为偏应力张量（deviatoric part）和球应力（或称静水应力）张量（spherical/hydrostatic part）两部分：

$$ \begin{equation} \begin{aligned} (\nabla \overline{\mathbf u})^\mathsf{T} =& \mathbf{dev2} \left( (\nabla \overline{\mathbf u})^\mathsf{T} \right) + \mathbf{hyd2} \left( (\nabla \overline{\mathbf u})^\mathsf{T} \right) \newline = & \underbrace{(\nabla \overline{\mathbf u})^\mathsf{T} - \frac{2}{3}\mathrm{tr} \left((\nabla \overline{\mathbf u})^\mathsf{T} \right) \mathbf I }_ % workaround for hugo {\text{偏应力张量}} + \underbrace{\frac{2}{3}\mathrm{tr} \left((\nabla \overline{\mathbf u})^\mathsf{T} \right)\mathbf I}_ % workaround for hugo {\text{球应力张量}} \end{aligned} \label{eq:stress_tensor} \end{equation} $$

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

$(\nabla \overline{\mathbf u})^\mathsf{T}$ 的球应力张量又有以下等式关系：

$$ \begin{equation} \mathrm{tr} \left((\nabla \overline{\mathbf u})^\mathsf{T} \right) = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = \nabla \cdot \overline{\mathbf u} \end{equation} $$

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

$$ \begin{equation} \begin{aligned} \nabla \cdot \mathbf R_\mathrm{eff} =& \nabla \cdot \left\lbrace \nu_\mathit{eff} \left[ \nabla \overline{\mathbf u} + \left( (\nabla \overline{\mathbf u})^\mathsf{T} - \frac{2}{3} (\nabla \cdot \overline{\mathbf u}) \mathbf I \right) \right] \right\rbrace \newline =& \nabla \cdot (\nu_\mathit{eff} \nabla \overline{\mathbf u}) + \nabla \cdot \left[ \nu_\mathit{eff} \mathbf{dev2} \left( (\nabla \overline{\mathbf u})^\mathsf{T} \right) \right] \end{aligned} \label{eq:divdevreff} \end{equation} $$

对于不可压缩流体，实际上方程右端第二项为零，但是 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)
    );
}

对应 $\eqref{eq:divdevreff}$ 式，其中的 fvm::laplacian 操作返回 $\nabla \cdot (\nu_\mathit{eff} \nabla \overline{\mathbf u})$，而 fvc::div 操作返回 $\nabla \cdot \left [ \nu_\mathit{eff} \mathbf{dev2} \left( (\nabla \overline{\mathbf u})^\mathsf{T} \right) \right]$ 。

T(fvc::grad(U)) 表示速度梯度的转置。dev2 表示 $\eqref{eq:stress_tensor}$ 中的球应力张量，定义如下：

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