5.4 黏性通量的离散

5.4 黏性通量的离散

计算式（5.2）的扩散通量 $ \bar{F}_{ v} $，需要知道控制体面上的流动量及其一阶导数。为方便，黏性通量采用与对流通量相同的控制体，以获得相容的空间离散并简化数据结构。

由于黏性通量具有椭圆性质，计算黏性项（2.23）、（2.24）及应力（2.15）所需的速度分量 $ (u, v, w) $、动力黏度 $ \mu $ 和热导率 k，直接在面上取平均。因此，对单元中心格式（图 5.13(a)），控制体面 IJ 上的数值为：

$$ U_{I J}=\frac{1}{2}\left(U_{I}+U_{J}\right),\tag{5.70} $$

其中，U 为上述任一流动变量。中位对偶格式的面 ij 也有类似表达式，见图 5.13(b)。

剩下的任务，是计算式（2.15）中速度分量及式（2.24）中温度的一阶导数（梯度）。可采用两种方式之一：

• 基于单元的梯度；

• 梯度的平均。

下面详细介绍这两种方法。

5.4.1 基于单元的梯度

这类梯度计算的共同特点，是必须存储网格单元信息，或与单元几何相关的系数。因此，数据结构必须在前述对流通量的面／边形式之外进一步扩充。下文讨论适用于单元中心和中位对偶离散的三种成熟方法。

面中心控制体

在控制体面上求梯度的一种方式，是定义以该面为中心的辅助控制体，再采用 Green–Gauss 定理。4.4 节已在结构有限体积离散框架内讨论这一方法。例如，对中位对偶格式，可将共享某条边的所有单元梯度按体积平均，得到边中点处的梯度 [94]。基于单元的梯度按式（5.49）计算：遍历全部网格单元，并将梯度贡献累加到边上。单元面上的 U 值由节点值平均得到，类似式（5.53）。该方法所需内存和运算次数较多，但可用于任意单元混合网格。

近似 Galerkin 有限元方法

另一种适用于中位对偶格式的方法，由 Galerkin 有限元方法导出 $ [31] $。其基本思想，是将控制体表面上的梯度积分，转化为中心节点处 Hessian 矩阵（二阶导数）的计算。此时，黏性项采用 笛卡尔 坐标下 Navier–Stokes 方程的微分形式，即式（A.4）中取 $ \xi = x $、$ \eta = y $、$ \zeta = z $、$ J^{-1} = 1 $，其中包含以下项：

$$ \partial_{x}(\mu\partial_{x}u),\text{等} $$

其中，$ \partial_{m}(\cdot)=\partial(\cdot)/\partial m $。因此，无需再在控制体各面上积分黏性通量。

原格式用于纯三角形／四面体网格，采用包含某个节点的所有单元的并集。为简化实现，对节点的动力黏度系数取平均，这是它与 Galerkin 方法的区别。于是，节点 i 处的二阶导数可按下式计算 $ [95] $：

$$ \begin{aligned}&\begin{bmatrix}\partial_{x}(\mu\partial_{x}U)&\partial_{ y}(\mu\partial_{x}U)&\partial_{z}(\mu\partial_{x}U)\\\partial_{x}(\mu\partial_{ y}U)&\partial_{ y}(\mu\partial_{ y}U)&\partial_{z}(\mu\partial_{ y}U)\\\partial_{x}(\mu\partial_{z}U)&\partial_{ y}(\mu\partial_{z}U)&\partial_{z}(\mu\partial_{z}U)\end{bmatrix}_{i}\\&=\frac{1}{\Omega^{\prime}}\sum_{j=1}^{N_{A}}\left\{\begin{bmatrix}\alpha_{xx}&\alpha_{xy}&\alpha_{xz}\\\alpha_{yx}&\alpha_{yy}&\alpha_{yz}\\\alpha_{zx}&\alpha_{zy}&\alpha_{zz}\end{bmatrix}_{ij}\frac{\mu_{i}+\mu_{j}}{2}\begin{pmatrix}U_{i}-U_{j}\end{pmatrix}\right\}.\end{aligned}\tag{5.71} $$

体积 $ \Omega^{\prime} $ 包含共享节点 i 的所有四面体。式（5.71）中的系数矩阵 $ \bar{\alpha} $ 关于对角线对称 [95]，即 $ \alpha_{xy} = \alpha_{yx}, \alpha_{xz} = \alpha_{zx} $ 及 $ \alpha_{ y z} = \alpha_{z y} $。因此，每条边仅需存储六个系数。这些系数为 [95]：

$$ \alpha_{n k}=\sum_{e}\frac{(S_{e}^{j})_{n}(S_{e}^{i})_{k}}{\Omega_{e}},\tag{5.72} $$

其中，n、k 表示 x、$ y $、z 下标；$ (S_{e}^{i})_{k} $ 和 $ (S_{e}^{j})_{n} $ 为图 5.18 所示外侧面矢量 $ \vec{S}_{e}^{i} $ 与 $ \vec{S}_{e}^{j} $ 的分量。求和遍及共享边 ij 的所有四面体，各四面体的体积为 $ \Omega_{e} $。若网格静止，可在预处理阶段计算系数。该方法的一个优点，是能采用与对流通量相同的基于边的数据结构。

不过，其缺点是只有在三角形或四面体网格上才保留完整黏性项。对棱柱或六面体等单元，该技术退化为在三个坐标方向上均类似于薄剪切层（TSL）的 Navier–Stokes 方程近似 $ [8] $。文献 $ [96] $ 将其扩展到非单纯形单元，并保留完整黏性项。但由于扩展模板包含未通过边与点 i 相连的节点，便不能再使用高效的基于边的数据结构。

节点值的平均

该格式用于单元中心控制体和纯四面体网格，采用文献 $ [4] $ 对文献 $ [62] $ 所提出模板的修正形式，计算单元面上的梯度。其基础是将定义该面的三个节点值取平均，并与单元形心处的数值结合。单元面的一阶导数由求解线性方程组 $ [4] $ 得到。

$$ \begin{bmatrix}x_{J}-x_{I}& y_{J}- y_{I}&z_{J}-z_{I}\\ \frac{1}{2}(x_{2}+x_{3})-x_{1}&\frac{1}{2}( y_{2}+ y_{3})- y_{1}&\frac{1}{2}(z_{2}+z_{3})-z_{1}\\ \frac{1}{2}(x_{1}+x_{3})-x_{2}&\frac{1}{2}( y_{1}+ y_{3})- y_{2}&\frac{1}{2}(z_{1}+z_{3})-z_{2}\end{bmatrix}\begin{bmatrix}\partial_{x}U\\ \partial_{ y}U\\ \partial_{z}U\end{bmatrix}=\begin{bmatrix}U_{J}-U_{I}\\ \frac{1}{2}(U_{2}+U_{3})-U_{1}\\ \frac{1}{2}(U_{1}+U_{3})-U_{2}\end{bmatrix},\tag{5.73} $$

对图 5.19，下标 I、J 表示单元形心，下标 1、2、3 分别表示节点 $ P_{1} $、$ P_{2} $、$ P_{3} $。节点处的流动变量，可用距离倒数加权（式（5.45））求得，也可用伪 Laplace 加权 [2] 求得，后者类似式（5.24）。

5.4.2 梯度的平均

既然已经计算了每个控制体内部的梯度，例如采用式（5.41）或式（5.42）的分片线性重构，自然会考虑直接用简单平均求面中点的梯度 $ [97] $：

$$ \overline{{\nabla U}}_{I J}=\frac{1}{2}\left[\nabla U_{I}+\nabla U_{J}\right].\tag{5.74} $$

该方法很有吸引力，因为它只需要基本的面／边数据结构，不需增加存储量。不过，如文献 $ ^{[95]} $ 等指出，它产生的模板较宽，且权重分布不利 $ [50] $。此外，文献 $ [50] $ 证明，这一模板允许四边形或六面体网格上的解发生解耦。

对单元中心格式，利用连接两个单元形心的方向上的方向导数，可以改善这一方法的性质，尤其能防止解耦，即：

$$ \left(\frac{\partial U}{\partial\ell}\right)_{I J}\approx U_{J}-U_{I}\ell_{I J},\tag{5.75} $$

校注：原书式（5.75）未显示分式结构，译文保留原式。按方向导数定义和量纲，应核对为 $(\partial U/\partial\ell)_{IJ}\approx(U_J-U_I)/\ell_{IJ}$。

其中，$ \ell_{IJ} $ 为单元形心 I、J 的距离（图 5.19 的虚线）。中位对偶格式也有类似表达式，其中 $ \vec{r}_{ij} $ 按式（5.43）给出。定义沿 I、J 连线的单位矢量 $ \vec{t}_{IJ} $：

$$ \vec{t}_{I J}=\frac{\vec{r}_{I J}}{\ell_{I J}},\tag{5.76} $$

修正后的平均可写为 $ [98, 99] $：

$$ \nabla U_{I J}=\overline{{\nabla U}}_{I J}-\left[\overline{{\nabla U}}_{I J}\cdot\vec{t}_{I J}-\left(\frac{\partial U}{\partial\ell}\right)_{I J}\right]\vec{t}_{I J},\tag{5.77} $$

其中，$ \overline{V U}_{I j} $ 由式（5.74）给出。该修正能在四面体、棱柱和六面体网格上形成强耦合模板 [50]，仍与基于面／边的数据结构相容，且无需增加存储量。因此，只要控制体内部梯度本来就用于计算对流通量，这一方法便比基于单元的方法更有吸引力。