4.1 控制体的几何量

第 4 章 结构有限体积格式

如第 3 章引言所述，求解 Euler 和 Navier–Stokes 方程的绝大多数数值格式采用线方法，即分别进行空间和时间离散。因此，可以根据待求问题的需要，对空间导数和时间导数采用不同精度的数值近似，从而获得很大的灵活性。本章也将采用线方法。关于空间与时间耦合离散的数值方法，例如 Lax–Wendroff 格式族（包括显式 MacCormack 预估—校正格式、隐式 Lerat 格式等），可参阅文献 $ ^{[1]} $ 中的详细讨论。

一般的结构有限体积格式自然以守恒定律为基础，即 Navier–Stokes 方程（2.19）或 Euler 方程（2.44）。在预处理阶段，将物理空间划分为网格单元：二维为四边形，三维为六面体。网格生成应满足：

• 网格完全覆盖计算域；

• 网格单元之间不留空隙；

• 网格单元互不重叠。

所得结构网格由网格点（网格单元角点）的坐标 x、y、z，以及计算空间中的索引（见图 3.2，记为 i、j、k）唯一确定。随后，以网格为基础定义控制体，以计算对流通量、黏性通量和源项的积分。为简单起见，假定某一控制体不随时间变化（否则参见第 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}_{\mathrm{v}})\mathrm{d}S-\int_{\Omega}\vec{Q}\mathrm{d}\Omega\right].\tag{4.1} $$

式（4.1）右侧的面积分，近似为穿过控制体各面的通量之和，这一近似称为空间离散。通常假定通量在每个面上恒定，并在面的中点计算；源项通常也假定在控制体内部恒定。然而，当源项占主导时，宜以相邻控制体内数值的加权和计算 $ \vec{Q} $（见 [2] 及其所引文献）。考虑图 4.1b 所示的特定控制体 $ \Omega_{I,J,K} $，由式（4.1）可得：

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

上式中的大写索引 $ (I,J,K) $ 标识控制体在计算空间中的位置（参见第 3.1 节）。后面将看到，控制体并不一定与网格单元重合。此外，$ N_{F} $ 表示控制体的面数。

面数在二维为 $ N_{F}=4 $，三维为 $ N_{F}=6 $。变量 $ \Delta S_{m} $ 表示面 m 的面积。式（4.2）右侧方括号中的项通常也称为残差，这里记为 $ \vec{R}_{I,J,K} $。因此，式（4.2）可简写为：

$$ \frac{\mathrm{d}\vec{W}_{I,J,K}}{\mathrm{d}t}=-\frac{1}{\Omega_{I,J,K}}\vec{R}_{I,J,K}.\tag{4.3} $$

对全部控制体 $ \Omega_{I,J,K} $ 写出式（4.3），即可得到一组一阶常微分方程。这些方程在时间上具有双曲型性质，意味着必须从已知初始解出发进行时间推进。还必须为黏性与无黏通量提供适当的边界条件，相关内容见第 8 章。

数值求解离散控制方程组 $ (4.3) $ 时，首先需要确定如何定义控制体，以及如何相对于计算网格布置流动变量。在结构有限体积格式中，有三种基本策略：

• 单元中心格式：控制体与网格单元相同，流动变量定义在单元形心（图 4.3）。

• 采用重叠控制体的单元顶点格式：流动物理量定义在网格点（顶点、节点），控制体则为所有共享该顶点的网格单元的并集（二维中为 4 个单元，三维中为 8 个，见图 4.4）。因此，与两个相邻网格点对应的控制体彼此重叠。

• 采用对偶控制体的单元顶点格式：流动变量同样存储在网格顶点，但控制体通过连接所有共享该顶点的单元中点构造（图 4.5）。这样，各网格点被相应控制体包围，而控制体互不重叠。

第 4.2 节讨论离散格式的一般概念时，将介绍这三种方法。这里应指出，如今采用重叠控制体的单元顶点格式已很少使用，但为完整起见仍予以介绍。

需要注意，在本书的方法中，全部流动变量——守恒变量 $ (\rho, \rho u, \rho v, \rho w, \rho E) $ 及依赖变量 $ (p, T, c, \text{等}) $——均定义在同一位置，即单元中心或网格节点。这称为同位网格格式。相比之下，许多较早的压力基方法（参见第 3.1 节）采用交错网格格式，将压力与速度分量存储在不同位置，以抑制中心差分引起的解振荡。

对流通量的计算有多种选择。基本问题在于，必须知道控制体全部 $ N_{F} $ 个面上的通量，而这些位置没有直接给定流动变量。因此，必须将通量或流动变量插值到控制体面上。原则上有两种做法：

• 像中心离散格式那样采用算术平均；

• 像迎风离散格式那样采用某种有偏插值，以考虑流动方程的特征性质。

第 4.3 节将介绍最常用的对流通量离散格式，并讨论其精度、适用范围与计算量。

在控制体面上计算黏性通量时，一种常用方法是对流动物理量作算术平均。与流动变量相比，式（2.15）和（2.24）中的速度梯度及温度梯度计算更为复杂。第 4.4 节将介绍相应方法。

4.1 控制体的几何量

在讨论离散方法之前，先考察控制体 $ \Omega_{I,J,K} $ 的几何量计算，有助于理解后续内容：包括控制体积、面 m 的单位法向量 $ \vec{n}_{m} $（定义为指向外侧）及面面积 $ \Delta S_{m} $。法向量与面面积也称为控制体的度量量。下面分别讨论一般四边形控制体的二维情况和六面体控制体的三维情况。

4.1.1 二维情况

通常将平面流动视为三维问题的一种特殊情况，即解相对于某一坐标方向（例如 z 方向）对称。为体现这种对称性，并使体积、压力等量具有正确的物理单位，将全部网格单元与控制体的厚度设为常数 b。于是，二维控制体的体积等于其平面面积乘以厚度 b。四边形面积可用 Gauss 公式精确计算。因此，对于图 4.1a 所示的控制体，经一定代数运算可得：

$$ \Omega_{I,J}=\frac{b}{2}\left[(x_{1}-x_{3})(y_{2}-y_{4})+(x_{4}-x_{2})(y_{1}-y_{3})\right].\tag{4.4} $$

这里假定控制体位于 $ x-y $ 平面，z 坐标轴为对称轴。厚度 b 是任意的，为方便可取 b = 1。在二维情况下，控制体各面由直线段表示，因此单位法向量沿每个面恒定。按式（4.2）近似积分通量时，需要计算面面积 $ \Delta S $ 与相应单位法向量 $ \vec{n} $ 的乘积，即面矢量 $ \vec{S} $。

$$ \vec{S}_{m}=\begin{bmatrix}S_{x,m}\\ S_{y,m}\end{bmatrix}=\vec{n}_{m}\Delta S_{m}.\tag{4.5} $$

由于对称性，面矢量及单位法向量的 z 分量为零，因而在表达式中省略。图 4.1a 中控制体的面矢量为：

$$ \begin{array}{l}\vec{S}_{1}=b\begin{bmatrix}y_{2}-y_{1}\\ x_{1}-x_{2}\end{bmatrix},\\\vec{S}_{2}=b\begin{bmatrix}y_{3}-y_{2}\\ x_{2}-x_{3}\end{bmatrix},\\\vec{S}_{3}=b\begin{bmatrix}y_{4}-y_{3}\\ x_{3}-x_{4}\end{bmatrix},\\\vec{S}_{4}=b\begin{bmatrix}y_{1}-y_{4}\\ x_{4}-x_{1}\end{bmatrix}.\end{array}\tag{4.6} $$

随后，由式（4.5）得到面 m 的单位法向量：

$$ \vec{n}_{m}=\frac{\vec{S}_{m}}{\Delta S_{m}}\tag{4.7} $$

其中，

$$ \Delta S_{m}=\mid\vec{S}_{m}\mid=\sqrt{S_{x,m}^{2}+S_{y,m}^{2}}. $$

实际计算中，只计算并存储面矢量 $ \vec{S}_{1} $ 和 $ \vec{S}_{4} $，每个控制体 $ \Omega_{I,J} $ 均如此。面矢量 $ \vec{S}_{2} $ 及 $ \vec{S}_{3} $ 从相应相邻控制体取得，反号后指向外侧，以节省内存并减少逐点运算量。

4.1.2 三维情况

与前面的二维情况不同，三维面矢量和体积的计算存在一些困难，主要原因是控制体某个面的四个顶点通常可能不共面。这样，法向量在面上不再恒定（图 4.2）。为克服这一困难，可将控制体的六个面分别划分为两个或更多三角形，再将体积分解为四面体。若以适当方式进行这种划分，可获得在任意网格上至少具有一阶精度的离散格式 $ [3] $。当然，计算量会明显增加，因为必须在每个子三角形上分别积分通量，逐点运算量至少加倍。然而，文献 $ [3, 4] $ 表明，对于较光顺、控制体面接近平行四边形的网格，三角形划分并不会明显改善解的精度。因此，下面采用基于平均法向量的四边形面简化处理。

计算图 4.1b 所示六面体控制体的面矢量 $ \vec{S} $，最方便的方法是采用与二维四边形面积计算相同的 Gauss 公式。例如，对于面 m = 1（图 4.1b 中的点 1、5、8、4），先定义差值：

$$ \begin{aligned}\Delta x_{A}&=x_{8}-x_{1},\quad&\Delta x_{B}&=x_{5}-x_{4},\\\Delta y_{A}&=y_{8}-y_{1},\quad&\Delta y_{B}&=y_{5}-y_{4},\\\Delta z_{A}&=z_{8}-z_{1},\quad&\Delta z_{B}&=z_{5}-z_{4}.\end{aligned}\tag{4.8} $$

于是，面矢量 $ \vec{S}_{1}=\vec{n}_{1}\Delta S_{1} $ 为：

$$ \vec{S}_{1}=\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{4.9} $$

其余五个面矢量按类似方式计算。同样，只存储六个面矢量中的三个（例如 $ \vec{S}_{1} $、$ \vec{S}_{3} $ 和 $ \vec{S}_{5} $）非常方便，每个控制体 $ \Omega_{I,J,K} $ 均如此。其余面矢量 $ \vec{S}_{2} $、$ \vec{S}_{4} $ 和 $ \vec{S}_{6} $ 从相应邻接控制体取得，反号后指向外侧。式（4.8）和（4.9）给出的是平均面矢量。当面接近平行四边形，即其全部顶点共面时，近似成为精确结果。单位法向量由式（4.7）求得，其中：

$$ \Delta S_{m}=\sqrt{S_{x,m}^{2}+S_{y,m}^{2}+S_{z,m}^{2}}.\tag{4.10} $$

一般六面体的体积计算有多种精度不同的公式（例如见 [1]）。一种在多种应用中表现很好的方法，基于散度定理 [5]；该定理将某矢量场散度的体积分与其面积分联系起来。关键思想是，将控制体 $ \Omega $ 内某点的空间位置矢量（记为 $ \vec{r} = [r_{x}, r_{y}, r_{z}]^{T} $）用作这个矢量场。于是，散度定理写为：

$$ \int_{\Omega}\operatorname{div}(\vec{r})\mathrm{d}\Omega=\oint_{\partial\Omega}(\vec{r}\cdot\vec{n})\mathrm{d}S.\tag{4.11} $$

式（4.11）左侧很容易计算，它给出所求控制体 $ \Omega $ 的体积：

$$ \begin{aligned}\int_{\Omega}\operatorname{div}(\vec{r})\mathrm{d}\Omega&=\int_{\Omega}\left(\frac{\partial r_{x}}{\partial x}+\frac{\partial r_{y}}{\partial y}+\frac{\partial r_{z}}{\partial z}\right)d\Omega\\&=3\Omega.\end{aligned}\tag{4.12} $$

假定单位法向量在控制体的每个面上恒定，则式（4.11）右侧面积分可计算为：

$$ \oint_{\partial\Omega}(\vec{r}\cdot\vec{n})\mathrm{d}S\approx\sum_{m=1}^{m=6}(\vec{r}_{\mathrm{mid}}\cdot\vec{n})_{m}\Delta S_{m}.\tag{4.13} $$

其中，$ \vec{r}_{mid,m} $ 为控制体面 m 的中点。例如，

$$ \vec{r}_{\mathrm{mid},1}=\frac{1}{4}(\vec{r}_{1}+\vec{r}_{5}+\vec{r}_{8}+\vec{r}_{4}), $$

矢量 $ \vec{r}_{1} $、$ \vec{r}_{5} $、$ \vec{r}_{8} $ 和 $ \vec{r}_{4} $ 对应图 4.1b 中面 m = 1 的顶点 1、5、8、4。其他面的中点具有类似关系。式（4.13）中面 m 的面积 $ \Delta S_{m} $ 由式（4.10）计算。结合式（4.12）和（4.13），并以面矢量 $ \vec{S} $ 代替乘积 $ \vec{n}_{m} $ $ \Delta S_{m} $，最终得到：

$$ \Omega_{I,J,K}=\frac{1}{3}\sum_{m=1}^{m=6}(\vec{r}_{\mathrm{mid}}\cdot\vec{S})_{m}\tag{4.14} $$

该式给出控制体 $ \Omega_{I,J,K} $ 的体积。对于面被划分为三角形的情况，类似公式见文献 $ ^{[6]} $。

原则上，坐标系原点可移到任意位置，而不影响式（4.14）的体积计算。因此，建议将原点设在控制体的某一顶点（例如图 4.1b 中的点 1），以改善数值的尺度。这样，可以将式（4.11）至（4.14）中的 $ \vec{r} $ 替换为变换后的矢量 $ \vec{r}^{*} $，定义为：

$$ \vec{r}^{*}=\vec{r}-\vec{r}_{\mathrm{origin}}. $$

应当注意，对于各面均为平面的控制体，用式（4.14）计算的体积是精确的。