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

动量方程

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

\[ \begin{equation} \underbrace{\frac {\partial (\rho {u}_i)}{\partial t}}_{时间项} + \underbrace{\frac{\partial (\rho u_j u_i)}{\partial x_j}}_{对流项} - \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]}_{粘性项} = - \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\{\mu \left[ \nabla \mathbf u + (\nabla \mathbf u)^T \right] \right\} = - \nabla p \end{equation} \]

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

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

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

\[ \begin{equation} \mathbf I = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \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\{ \mu \left[ \nabla \overline{\mathbf u} + (\nabla \overline{\mathbf u})^T - \frac{2}{3} \nabla \cdot \overline{\mathbf u} \right] + \rho \mathbf R \right\} = - \nabla \overline{p} \end{equation} \]

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

\[ \begin{equation} \mathbf R = \nu_t \left[ \nabla \overline{\mathbf u} + (\nabla \overline{\mathbf u})^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\{ \mu_\mathit{eff} \left[ \nabla \overline{\mathbf u} + (\nabla \overline{\mathbf u})^T - \frac{2}{3} (\nabla \cdot \overline{\mathbf u}) \mathbf I \right] \right\} = - \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})^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) \end{equation} \label{eq:ranse} \]

上式为 OpenFOAM 采用的形式。

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

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

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

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

\[ \begin{equation} \mathrm{tr} \left((\nabla \overline{\mathbf u})^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 \nu_\mathit{eff} \left[ \nabla \overline{\mathbf u} + \left( (\nabla \overline{\mathbf u})^T - \frac{2}{3} (\nabla \cdot \overline{\mathbf u}) \mathbf I \right) \right] \\ =& \nu_\mathit{eff} \nabla^2 \overline{\mathbf u} + \nabla \cdot \left [ \nu_\mathit{eff} \mathbf{dev2} \left( (\nabla \overline{\mathbf u})^T \right) \right] \end{aligned} \end{equation} \label{eq:divdevreff} \]

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

代码实现

不可压求解器中,动量方程的代码如下(以 pimpleFoam 求解器为例):

1
2
3
4
5
6
7
8
9
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
9
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
10
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
14
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 操作返回 $\nu_\mathit{eff} \nabla^2 \overline{\mathbf u}​$,而 fvc::div 操作返回 $\nabla \cdot \left [ \nu_\mathit{eff} \mathbf{dev2} \left( (\nabla \overline{\mathbf u})^T \right) \right]​$ 。

T 表示转置。dev2 表示 $\eqref{eq:stress_tensor}$ 中的球应力张量,定义如下:

1
2
3
4
5
6
7
//- 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);
}