5.1 控制体的几何量

第 5 章 非结构有限体积格式

如第 3 章引言所述，绝大多数求解 Euler 方程和 Navier–Stokes 方程的数值格式采用线方法，即分别进行空间离散和时间离散。其主要优点是可以为空间导数与时间导数选择不同精度的数值近似。因此，与采用时空耦合离散的 Lax–Wendroff 类方法相比，这种方法具有明显更大的灵活性。后者包括显式 MacCormack 预估—校正格式、隐式 Lerat 格式等，详细介绍可参见文献 $ ^{[1]} $。由于线方法应用广泛，本章也采用这一方法。

本章讨论的有限体积格式以守恒定律为基础，其数学表达为 Navier–Stokes 方程（2.19）或 Euler 方程（2.44）。在预处理阶段，首先将物理域划分为若干单元（网格单元）。二维单元通常为三角形，有时也与四边形组合使用；三维则最常采用四面体 $ [2-7] $。不过，在高 Reynolds 数黏性流动模拟中，越来越多的流动求解器采用由四面体、棱柱、棱锥以及某些情况下的六面体组成的网格（图 5.1）$ [8-16] $。包含多种单元类型的非结构网格称为混合单元网格（mixed grid），实例见图 3.3 和图 5.2。这一名称应与混合结构网格（hybrid grid）区分；后者指结构网格与非结构网格的组合，例如文献 $ [17-19] $ 所用的网格。

网格生成必须保证控制方程的守恒性质，即：

• 网格必须完全覆盖物理域；

• 单元之间不得留有空隙；

• 单元不得彼此重叠。

除满足上述要求外，网格还应保持光滑，即相邻网格单元的体积或拉伸比不应有很大差别。此外，单元应尽可能规整，尤其要避免薄片状四面体等退化单元（见图 11.24）。否则，数值误差可能使解的精度完全丧失 $ [22, 23] $。这类单元还会造成严重的收敛问题。

在网格基础上定义适当的控制体，以计算对流通量、黏性通量及源项的积分。为简单起见，假设某一控制体不随时间变化；控制体随时间变化的情况见 A.5 节。于是，守恒变量 $ \vec{W} $ 的时间导数可写为：

$$ \frac{\partial}{\partial t}\int_{\Omega}\vec{W}\mathrm{d}\Omega=\Omega\frac{\partial\vec{W}}{\partial t}. $$

因此，式（2.19）成为：

$$ \frac{\partial\vec{W}}{\partial t}=-\frac{1}{\Omega}\left[\oint_{\partial\Omega}(\vec{F}_{c}-\vec{F}_{v})\mathrm{d}S-\int_{\Omega}\vec{Q}\mathrm{d}\Omega\right].\tag{5.1} $$

式（5.1）右端的面积分，用穿过控制体各面的通量之和近似。这种近似称为空间离散。通常假设通量在每个面上为常量，并在面的中点计算。对二阶精度格式而言，这样处理已足够。一般还假设源项在控制体内为常量。不过，若源项占主导地位，则宜将 $ \vec{Q} $ 按相邻控制体中数值的加权和计算，见文献 [24] 及其所引文献。对某一控制体 $ \Omega_{I} $，由式（5.1）可得：

$$ \frac{\mathrm{d}\vec{W}_{I}}{\mathrm{d}t}=-\frac{1}{\Omega_{I}}\left[\sum_{m=1}^{N_{\mathrm{F}}}(\vec{F}_{c}-\vec{F}_{ v})_{m}\Delta S_{m}-(\vec{Q}\Omega)_{I}\right].\tag{5.2} $$

上式中，大写下标 I 表示控制体，因为控制体未必与网格单元重合，后文将对此作说明。此外，$ N_{F} $ 表示控制体 $ \Omega_{I} $ 的面数，变量 $ \Delta S_{m} $ 表示第 m 个面的面积。面数 $ N_{F} $ 既取决于单元类型，也取决于控制体类型。一般而言，不同控制体的面数也不相同，这是非结构网格与结构网格的主要区别之一。不过，已有数值过程和数据结构可以避免预先知道 $ N_{F} $ 的具体值。后续各节将再讨论这一点。

式（5.2）右端方括号内的项通常称为残差。因此，式（5.2）可简写为：

$$ \frac{\mathrm{d}\vec{W}_{I}}{\mathrm{d}t}=-\frac{1}{\Omega_{I}}\vec{R}_{I}.\tag{5.3} $$

对所有控制体 $ \Omega_{I} $ 写出式（5.3），便得到一组一阶常微分方程。这些方程在时间上具有双曲性质，因此必须从已知初始解出发沿时间推进。同时，还必须为黏性通量和无黏通量给定适当的边界条件，具体见第 8 章。

数值求解离散后的控制方程组 $ (5.3) $ 时，首先要确定如何定义控制体，以及流动变量相对于网格点应放在何处。

在有限体积格式中，可以采用三种基本策略：

• 单元中心格式 $ [2, 16, 25, 26] $：控制体与网格单元重合，流动变量位于单元形心（图 5.6）。

• 采用重叠控制体的单元顶点格式 $ [27, 28] $：流动量存放在网格顶点处，控制体由共享该顶点的所有网格单元的并集构成。因此，与两个相邻顶点对应的控制体彼此重叠。

• 采用中位对偶控制体的单元顶点格式 $ [13, 15, 29-33] $：流动变量仍存放在网格顶点，但控制体由周围单元的形心、面的形心及边的中点连接而成（图 5.7 和图 5.8）。这样，每个网格点被其对应的控制体包围。这些控制体构成对偶网格，彼此不重叠。

采用重叠控制体的单元顶点格式如今已不再使用，因此本章重点讨论单元中心格式和中位对偶格式。5.2 节将详细介绍这两种方法。

应注意，这里所有流动变量，即守恒变量 $ (\rho, \rho u, \rho v, \rho w, \rho E) $ 和依赖变量 $ (p, T, c, \operatorname{etc.}) $，都位于同一位置：单元中心或网格点。这种方法称为同位网格格式。相比之下，许多较早的基于压力的结构网格方法（参见 3.1 节）采用所谓交错网格格式，将压力和速度分量存放在不同位置，以抑制空间中心差分引起的解的振荡。

对流通量的计算有多种选择。基本困难在于：必须知道控制体全部 $ N_{F} $ 个面上的通量，但这些面上并无直接给定的流动变量。因此，必须将通量或流动变量插值到控制体面上。将控制体内部的流动变量插值至面上的过程称为解的重构（见 5.3.3 节）。原则上，插值可采用以下两种方式之一：

• 像中心离散格式那样进行算术平均；

• 像迎风离散格式那样采用某种偏向插值，使其体现流动方程的特征性质。

5.3 节除介绍最常用的对流通量离散格式外，还将讨论它们的精度、适用范围和计算量。

计算控制体面上黏性通量的一种常用方法，是对流动量作算术平均。式（2.15）和式（2.24）中的速度梯度与温度梯度则较难计算，混合单元网格尤其如此。5.4 节将给出完整的处理过程。

5.1 控制体的几何量

在讨论对流通量与黏性通量的离散方法之前，先考察控制体 $ \Omega_{I} $ 几何量的计算。这些量包括控制体积、外指单位法向量 $ \vec{n}_{m} $、第 m 个面的面积 $ \Delta S_{m} $，以及单元的形心。法向量与面面积也称为控制体的度量。下文分别讨论二维和三维情况。

5.1.1 二维情况

通常，将平面流动视为三维问题的一种特殊情形，即解对某一坐标方向（例如 z 方向）具有对称性。考虑这种对称性，并为体积、压力等量获得正确的物理单位，将所有网格单元和控制体的厚度设为同一常数 b。于是，二维控制体的体积为其面积与厚度 b 的乘积。由于 b 可任意选取，为方便可令 b = 1。以下只讨论三角形和四边形单元。中位对偶格式的控制体虽然可能具有相当复杂的形状，却总能分解为三角形和／或四边形。

三角形单元

一般三角形的面积可用 Gauss 公式方便而精确地计算。按图 5.3(a) 的节点编号，体积为：

$$ \begin{aligned}\Omega=&\frac{b}{2}\big[(x_{1}-x_{2})( y_{1}+ y_{2})\\&+(x_{2}-x_{3})( y_{2}+ y_{3})\\&+(x_{3}-x_{1})( y_{3}+ y_{1})\big].\end{aligned}\tag{5.4} $$

节点必须按逆时针方向编号，才能得到正体积。

四边形单元

一般四边形的面积也可由 Gauss 公式精确计算。经代数整理，得到：

$$ \Omega=\frac{b}{2}\left[(x_{1}-x_{3})( y_{2}- y_{4})+(x_{4}-x_{2})( y_{1}- y_{3})\right],\tag{5.5} $$

其中，节点按图 5.3(b) 以逆时针方向编号。这里假设控制体位于 x–y 平面，z 坐标轴为对称轴。

二维控制体的边为直线，因此单位法向量沿边保持不变。按式（5.2）的近似计算通量积分时，需要计算面的面积 $ \Delta S $ 与相应单位法向量 $ \vec{n} $ 的乘积，即面矢量 $ \vec{S} $。以图 5.3 中的边 2–3 为例，其外指面矢量为：

$$ \vec{S}_{23}=\vec{n}_{23}\Delta S_{23}=b\begin{bmatrix} y_{3}- y_{2}\\ x_{2}-x_{3}\end{bmatrix}.\tag{5.6} $$

由于对称性，面矢量及单位法向量的 z 分量为零，所以式（5.6）略去该分量。由式（5.6）可得单位法向量：

$$ \Delta S=\mid\vec{S}\mid=\sqrt{S_{x}^{2}+S_{ y}^{2}},\tag{5.7} $$

其中，$ S_{x} $ 和 $ S_{ y} $ 为面矢量的 笛卡尔 分量。

单元中心

图 5.3(a) 中三角形的中心定义为：

$$ \vec{r}_{c}=\frac{1}{3}\left(\vec{r}_{1}+\vec{r}_{2}+\vec{r}_{3}\right)\tag{5.8} $$

其中，$ \vec{r}_{1/2/3} $ 表示节点的 笛卡尔 坐标。四边形单元的中心可由文献 $ [34] $ 给出的公式计算。为此，将四边形分成两个共享两个节点的三角形。若采用图 5.3(b) 的节点编号，且公共节点为 1 和 3，则有：

$$ \vec{r}_{c}=\Omega_{123}\vec{r}_{c,123}+\Omega_{134}\vec{r}_{c,134}\Omega_{123}+\Omega_{134}.\tag{5.9} $$

校注：原书式（5.9）未显示分式结构，译文保留原式。由前文的体积加权形心定义，应核对为 $\vec r_c=(\Omega_{123}\vec r_{c,123}+\Omega_{134}\vec r_{c,134})/(\Omega_{123}+\Omega_{134})$。

两个三角形的体积 $ \Omega_{123} $ 与 $ \Omega_{134} $ 由式（5.4）计算，其形心 $ \vec{r}_{c} $ 由式（5.8）求得。一般多边形的相应公式见文献 [35]。

5.1.2 三维情况

与二维情况不同，对具有四边形面的三维单元或控制体，面矢量与体积的计算存在一些困难。主要原因是：控制体四边形面的四个顶点一般未必共面，此时法向量在该面上不再为常量（见图 4.2）。为解决这一问题，可以将每个四边形面分成两个甚至更多三角形。不过，对光滑网格上的二阶格式，这样获得的精度提升几乎难以察觉。只有对三阶及更高阶的空间离散，增加的计算量才有意义，而且实际上是必需的。因此，下文采用基于平均法向量的简化方法处理四边形面。

三角形面

三角形面的面矢量 S 可由 Gauss 公式精确计算。按图 5.4(a) 定义节点，三角形 1–2–3 的边差矢量为：

$$ \begin{aligned}\Delta x y_{A}&=(x_{1}-x_{2})( y_{1}+ y_{2}),\\\Delta y z_{A}&=( y_{1}- y_{2})(z_{1}+z_{2}),\\\Delta x y_{B}&=(x_{2}-x_{3})( y_{2}+ y_{3}),\\\Delta y z_{B}&=( y_{2}- y_{3})(z_{2}+z_{3}),\\\Delta x y_{C}&=(x_{3}-x_{1})( y_{3}+ y_{1}),\\\Delta y z_{C}&=( y_{3}- y_{1})(z_{3}+z_{1}),\\\Delta z x_{A}&=(z_{1}-z_{2})(x_{1}+x_{2}),\\\Delta z x_{B}&=(z_{2}-z_{3})(x_{2}+x_{3}),\\\Delta z x_{C}&=(z_{3}-z_{1})(x_{3}+x_{1}).\end{aligned}\tag{5.10} $$

于是，外指面矢量 $ \vec{S}=\vec{n}\Delta S $ 为：

$$ \vec{S}=\frac{1}{2}\begin{bmatrix}\Delta y z_{A}+\Delta y z_{B}+\Delta y z_{C}\\ \Delta z x_{A}+\Delta z x_{B}+\Delta z x_{C}\\ \Delta x y_{A}+\Delta x y_{B}+\Delta x y_{C}\end{bmatrix}.\tag{5.11} $$

四边形面

对图 5.4(b) 所示的四边形面，平均面矢量 $ \vec{S} $ 最方便的计算方法，是采用二维四边形面积计算中使用的同一 Gauss 公式。

因此，对图 5.4(b) 中由节点 5、6、7、8 确定的面，先定义以下差矢量：

$$ \begin{aligned}\Delta x_{A}&=x_{8}-x_{6},\quad&\Delta x_{B}&=x_{7}-x_{5},\\\Delta y_{A}&= y_{8}- y_{6},\quad&\Delta y_{B}&= y_{7}- y_{5},\\\Delta z_{A}&=z_{8}-z_{6},\quad&\Delta z_{B}&=z_{7}-z_{5}.\end{aligned}\tag{5.12} $$

随后由下式求外指面矢量 $ \vec{S}=\vec{n}\Delta S $：

$$ \vec{S}=\frac{1}{2}\begin{bmatrix}\Delta z_{A}\Delta y_{B}-\Delta y_{A}\Delta z_{B}\\ \Delta x_{A}\Delta z_{B}-\Delta z_{A}\Delta x_{B}\\ \Delta y_{A}\Delta x_{B}-\Delta x_{A}\Delta y_{B}\end{bmatrix}.\tag{5.13} $$

当面趋近于平行四边形，即面的全部顶点共面时，这一近似成为精确结果。

在上述两种情况下，单位法向量均由 $ \vec{n}=\vec{S}/\Delta S $ 按下式求得：

$$ \Delta S=\sqrt{S_{x}^{2}+S_{ y}^{2}+S_{z}^{2}},\tag{5.14} $$

其中，$ S_{x} $、$ S_{ y} $ 和 $ S_{z} $ 分别为式（5.11）或式（5.13）所给面矢量的 笛卡尔 分量。

体积

讨论三维结构有限体积格式时已经指出，计算体积的一种方便方法基于散度定理 $ [36] $。4.1.2 节的讨论最终得到以下体积表达式：

$$ \Omega=\frac{1}{3}\sum_{m=1}^{N_{\mathrm{F}}}(\vec{r}_{c}\cdot\vec{S})_{m}\tag{5.15} $$

其中，$ N_{F} $ 为控制体的面数，$ (\vec{r}_{c})_{m} $ 为控制体第 m 个面的中心，$ \vec{S}_{m} $ 为该面的外指面矢量。式（5.15）可直接用于非结构网格。对全部由三角形面构成的体积，或全部四边形面均为平面的体积，该式是精确的。

单元形心

构造前述中位对偶控制体，需要知道网格单元的形心。一般体积的形心定义为：

$$ \vec{r}_{c}=\frac{1}{\Omega}\int_{\Omega}\vec{r}\mathrm{d}\Omega.\tag{5.16} $$

按文献 $ [34] $ 的推导，式（5.16）可离散为：

$$ \vec{r}_{c}=\frac{3\sum_{m=1}^{N_{\mathrm{F}}}(\vec{r}_{c}\cdot\vec{n})_{m}(\vec{r}_{c})_{m}\Delta S_{m}}{4\sum_{m=1}^{N_{\mathrm{F}}}(\vec{r}_{c}\cdot\vec{n})_{m}\Delta S_{m}},\tag{5.17} $$

其中，第 m 个面的中心 $ (\vec{r}_{c})_{m} $ 对三角形面由式（5.8）计算，对四边形面由式（5.9）计算。可以看出，式（5.17）的分母正是式（5.15）中的 $ 12\Omega $。由三角形面构成的多面体所用的类似公式见文献 [35]。