OpenFOAM 中的不可压和可压 RANS 方程

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

动量方程

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

$$\begin{array}{}\text{(1)}& \underset{\text{时间项}}{\underbrace{\frac{\mathrm{\partial }(\rho u_i)}{\mathrm{\partial }t}}}+\underset{\text{对流项}}{\underbrace{\frac{\mathrm{\partial }(\rho u_j{u}_i)}{\mathrm{\partial }{x}_j}}}-\underset{\text{粘性项}}{\underbrace{\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[\mu (\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_i})\right]}}=-\frac{\mathrm{\partial }p}{\mathrm{\partial }{x}_i}\end{array}$$

将其写作矢量形式，如下

$$\begin{array}{}\text{(2)}& \frac{\mathrm{\partial }(\rho \mathbf{u})}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\rho \mathbf{u}\mathbf{u})-\mathrm{\nabla }\cdot \left\{\mu [\mathrm{\nabla }\mathbf{u}+(\mathrm{\nabla }\mathbf{u}{)}^{\mathsf{T}}]\right\}=-\mathrm{\nabla }p\end{array}$$

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

$$\begin{array}{}\text{(3)}& \frac{\mathrm{\partial }(\rho \mathbf{u})}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\rho \mathbf{u}\mathbf{u})-\mathrm{\nabla }\cdot \left\{\mu [\mathrm{\nabla }\mathbf{u}+(\mathrm{\nabla }\mathbf{u}{)}^{\mathsf{T}}-\frac{2}{3}(\mathrm{\nabla }\cdot \mathbf{u})\mathbf{I}]\right\}=-\mathrm{\nabla }p\end{array}$$

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

$$\begin{array}{}\text{(4)}& \mathbf{I}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\end{array}$$

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

$$\begin{array}{}\text{(5)}& \frac{\mathrm{\partial }(\rho \overline{\mathbf{u}})}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\rho \overline{\mathbf{u}}\overline{\mathbf{u}})-\mathrm{\nabla }\cdot \{\mu [\mathrm{\nabla }\overline{\mathbf{u}}+(\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}}-\frac{2}{3}(\mathrm{\nabla }\cdot \overline{\mathbf{u}})\mathbf{I}]+\rho \mathbf{R}\}=-\mathrm{\nabla }\overline{p}\end{array}$$

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

$$\begin{array}{}\text{(6)}& \mathbf{R}={\nu }_t[\mathrm{\nabla }\overline{\mathbf{u}}+(\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}}-\frac{2}{3}(\mathrm{\nabla }\cdot \overline{\mathbf{u}})\mathbf{I}]-\frac{2}{3}k\mathbf{I}\end{array}$$

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

$$\begin{array}{}\text{(7)}& \frac{\mathrm{\partial }(\rho \overline{\mathbf{u}})}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\rho \overline{\mathbf{u}}\overline{\mathbf{u}})-\mathrm{\nabla }\cdot \left\{{\mu }_{\mathit{e}\mathit{f}\mathit{f}}[\mathrm{\nabla }\overline{\mathbf{u}}+(\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}}-\frac{2}{3}(\mathrm{\nabla }\cdot \overline{\mathbf{u}})\mathbf{I}]\right\}=-\mathrm{\nabla }(\overline{p}+\frac{2}{3}\rho k)\end{array}$$

有效雷诺应力

定义有效雷诺应力如下：

$$\begin{array}{}\text{(8)}& {\mathbf{R}}_{\mathrm{e}\mathrm{f}\mathrm{f}}={\nu }_{\mathit{e}\mathit{f}\mathit{f}}[\mathrm{\nabla }\overline{\mathbf{u}}+(\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}}-\frac{2}{3}(\mathrm{\nabla }\cdot \overline{\mathbf{u}})\mathbf{I}]\end{array}$$

那么 RANS 方程可表示为：

$$\begin{array}{}\text{(9)}& \frac{\mathrm{\partial }(\rho \overline{\mathbf{u}})}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\rho \overline{\mathbf{u}}\overline{\mathbf{u}})-\mathrm{\nabla }\cdot (\rho {\mathbf{R}}_{\mathrm{e}\mathrm{f}\mathrm{f}})=-\mathrm{\nabla }(\overline{p}+\frac{2}{3}\rho k)\end{array}$$

上式为 OpenFOAM 采用的形式。

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

$$\begin{array}{}\text{(10)}& \begin{array}{rl}(\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}}=& \mathbf{d}\mathbf{e}\mathbf{v}\mathbf{2}((\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}})+\mathbf{h}\mathbf{y}\mathbf{d}\mathbf{2}((\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}})\\ =& \underset{\text{偏应力张量}}{\underbrace{(\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}}-\frac{2}{3}\mathrm{t}\mathrm{r}((\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}})\mathbf{I}}}+\underset{\text{球应力张量}}{\underbrace{\frac{2}{3}\mathrm{t}\mathrm{r}((\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}})\mathbf{I}}}\end{array}\end{array}$$

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

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

$$\begin{array}{}\text{(11)}& \mathrm{t}\mathrm{r}((\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}})=\frac{\mathrm{\partial }u}{\mathrm{\partial }x}+\frac{\mathrm{\partial }v}{\mathrm{\partial }y}+\frac{\mathrm{\partial }w}{\mathrm{\partial }z}=\mathrm{\nabla }\cdot \overline{\mathbf{u}}\end{array}$$

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

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathrm{\nabla }\cdot {\mathbf{R}}_{\mathrm{e}\mathrm{f}\mathrm{f}}=& \mathrm{\nabla }\cdot \left\{{\nu }_{\mathit{e}\mathit{f}\mathit{f}}[\mathrm{\nabla }\overline{\mathbf{u}}+((\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}}-\frac{2}{3}(\mathrm{\nabla }\cdot \overline{\mathbf{u}})\mathbf{I})]\right\}\\ =& \mathrm{\nabla }\cdot ({\nu }_{\mathit{e}\mathit{f}\mathit{f}}\mathrm{\nabla }\overline{\mathbf{u}})+\mathrm{\nabla }\cdot \left[{\nu }_{\mathit{e}\mathit{f}\mathit{f}}\mathbf{d}\mathbf{e}\mathbf{v}\mathbf{2}((\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}})\right]\end{array}\end{array}$$

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

对应 $\text{(12)}$ 式，其中的 fvm::laplacian 操作返回 $\mathrm{\nabla }\cdot ({\nu }_{\mathit{e}\mathit{f}\mathit{f}}\mathrm{\nabla }\overline{\mathbf{u}})$，而 fvc::div 操作返回 $\mathrm{\nabla }\cdot \left[{\nu }_{\mathit{e}\mathit{f}\mathit{f}}\mathbf{d}\mathbf{e}\mathbf{v}\mathbf{2}((\mathrm{\nabla }\overline{\mathbf{u}}{)}^{\mathsf{T}})\right]$ 。

T(fvc::grad(U)) 表示速度梯度的转置。dev2 表示 $\text{(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);
}