4.2 一般离散方法

4.2 一般离散方法

本章引言已提到控制体定义与流动变量布置的三种方法。本节将逐一详细介绍，并讨论其优缺点。

4.2.1 单元中心格式

当控制体与网格单元相同、流动变量定义在单元形心（如图 4.3 所示）时，称为单元中心格式。计算离散流动方程（4.2）时，必须提供单元各面上的对流和黏性通量 [7]。

这些通量可用以下三种方式之一近似：

1. 使用同一个面矢量，分别根据单元面左右两侧网格单元形心的变量计算通量，再对通量取平均（一般仅用于对流通量）；

2. 对单元面左右两侧网格单元形心的变量取平均，再计算通量；

3. 分别将流动物理量插值到单元面的左右两侧，再据此计算通量（仅用于对流通量）。

以图 4.3 中的单元面 $ \vec{n}_{I+1/2,J} $ 为例，第一种方法——通量平均——在二维中写为：

$$ \left(\vec{F}_{c}\Delta S\right)_{I+1/2,J}\approx\frac{1}{2}\left[\vec{F}_{c}(\vec{W}_{I,J})+\vec{F}_{c}(\vec{W}_{I+1,J})\right]\Delta S_{I+1/2,J}\tag{4.15} $$

其中，$ \Delta S_{I+1/2,J} $ 由式（4.6）和（4.7）计算。

第二种方法——变量平均——可写为：

$$ \left(\vec{F}\Delta S\right)_{I+1/2,J}\approx\vec{F}(\vec{W}_{I+1/2,J})\Delta S_{I+1/2,J},\tag{4.16} $$

控制体面 $ \vec{n}_{I+1/2,J} $ 上的守恒／依赖变量，定义为两个相邻单元内数值的算术平均，即：

$$ \vec{W}_{I+1/2,J}=\frac{1}{2}\left(\vec{W}_{I,J}+\vec{W}_{I+1,J}\right).\tag{4.17} $$

式（4.16）中的通量矢量 $ \vec{F} $，可以表示对流通量，也可以表示黏性通量。

第三种方法首先将流动物理量（通常为速度分量、压力、密度和总焓）分别插值到单元面的两侧。插值结果称为左状态和右状态（见第 4.3 节开头），两侧数值通常不同。穿过单元面的通量，随后由左、右状态之差通过某个非线性函数计算。因此，

$$ \left(\vec{F}_{\mathrm{c}}\Delta S\right)_{I+1/2,J}\approx f_{\mathrm{Flux}}\left(\vec{U}_{\mathrm{L}},\vec{U}_{\mathrm{R}},\Delta S_{I+1/2,J}\right),\tag{4.18} $$

其中，

$$ \begin{array}{l}\vec{U}_{\mathrm{L}}=f_{\mathrm{Interp}}\left(\cdots,\vec{U}_{I-1,J},\vec{U}_{I,J},\cdots\right)\\ \vec{U}_{\mathrm{R}}=f_{\mathrm{Interp}}\left(\cdots,\vec{U}_{I,J},\vec{U}_{I+1,J},\cdots\right)\end{array}\tag{4.19} $$

表示插值状态。当然，对于其他单元面，也存在类似于 $ (4.15)-(4.19) $ 的关系。

三维中采用相同近似。例如，在单元面 $ \vec{n}_{I+1/2,J,K} $（例如图 4.1b 中的 $ \vec{n}_{2} $）上，式（4.15）的通量平均变为：

$$ \left(\vec{F}_{c}\Delta S\right)_{I+1/2,J,K}\approx\frac{1}{2}\left[\vec{F}_{c}(\vec{W}_{I,J,K})+\vec{F}_{c}(\vec{W}_{I+1,J,K})\right]\Delta S_{I+1/2,J,K}\tag{4.20} $$

其中，$ \Delta S_{I+1/2,J,K} $ 按式（4.8）和（4.9）相应定义。与式（4.16）类似，变量平均写为：

$$ \begin{array}{r}{\big(\vec{F}\Delta S\big)_{I+1/2,J,K}\approx\vec{F}(\vec{W}_{I+1/2,J,K})\Delta S_{I+1/2,J,K}}\end{array}\tag{4.21} $$

其中，

$$ \vec{W}_{I+1/2,J,K}=\frac{1}{2}\left(\vec{W}_{I,J,K}+\vec{W}_{I+1,J,K}\right).\tag{4.22} $$

通过流动变量插值计算通量的方法，与式（4.18）类似，给出：

$$ \left(\vec{F}_{\mathrm{c}}\Delta S\right)_{I+1/2,J,K}\approx f_{\mathrm{Flux}}\left(\vec{U}_{\mathrm{L}},\vec{U}_{\mathrm{R}},\Delta S_{I+1/2,J,K}\right),\tag{4.23} $$

其中，$ \vec{U}_{L} $ 与 $ \vec{U}_{R} $ 是单元面上的插值结果。

离散流动方程（4.2）中最后一个尚待计算的项，是源项 $ \vec{Q} $。如引言所述，通常假定源项在控制体内恒定，因此使用相应单元中心的流动变量计算，即可定义为：

$$ (\vec{Q}\Omega)_{I,J,K}=\vec{Q}(\vec{W}_{I,J,K})\Omega_{I,J,K}.\tag{4.24} $$

利用上述关系，可以计算各面通量，并按式（4.2）完成对 $ \Omega_{I,J,K} $ 边界的数值积分；换言之，得到完整残差 $ \vec{R}_{I,J,K} $。第 4.3 和 4.4 节将详细介绍对流通量与黏性通量的计算。

4.2.2 单元顶点格式：重叠控制体

在单元顶点格式中，全部流动变量定义在计算网格的节点上。对于采用重叠控制体的方法，

与单元中心格式一样，网格单元仍作为控制体。不同之处是，计算所得控制体残差现在必须分配到网格点 $ [4, 8, 9] $。图 4.4 给出了这一过程。

考虑图 4.4 中的控制体 $ \Omega_{I,J} $，其由以下节点定义：

$$ \begin{array}{ccc}(i,j)&\quad(i+1,j)&\quad(i+1,j+1)\quad(i,j+1).\end{array} $$

注意，点 $ (i,j) $ 位于 $ \Omega_{I,J} $ 的左下角。例如，面 $ \Delta S_{I,J-1/2} $ 由点 $ (i,j) $ 和 $ (i+1,j) $ 定义，其对流通量近似为：

$$ \begin{array}{r}{\big(\vec{F}_{\mathrm{c}}\Delta S\big)_{I,J-1/2}\approx\vec{F}_{\mathrm{c}}(\vec{W}_{I,J-1/2})\Delta S_{I,J-1/2}.}\end{array}\tag{4.25} $$

面中点的变量，通过对定义该面的节点变量作算术平均获得，即：

$$ \vec{W}_{I,J-1/2}=\frac{1}{2}\left[\vec{W}_{i,j}+\vec{W}_{i+1,j}\right].\tag{4.26} $$

面面积 $ \Delta S_{I,J-1/2} $ 由式（4.6）和（4.7）计算。

三维情况下方法相同。例如，对于图 4.1b 中与法向量 $ \vec{n}_{3} $ 对应的面，平均变量为：

$$ \vec{W}_{I,J,K-1/2}=\frac{1}{4}\left[\vec{W}_{1}+\vec{W}_{2}+\vec{W}_{5}+\vec{W}_{6}\right].\tag{4.27} $$

假定边 1–2 沿 i 方向、边 1–5 沿 j 方向、边 1–4 沿 k 方向，并将点 $ (i,j,k) $ 对应于图 4.1b 的角点 1，则式（4.27）的平均也可写成：

$$ \vec{W}_{I,J,K-1/2}=\frac{1}{4}\left[\vec{W}_{i,j,k}+\vec{W}_{i+1,j,k}+\vec{W}_{i,j+1,k}+\vec{W}_{i+1,j+1,k}\right].\tag{4.28} $$

随后，对流通量仍由下式求得：

$$ \begin{array}{r}{\big(\vec{F}_{\mathrm{c}}\Delta S\big)_{I,J,K-1/2}\approx\vec{F}_{\mathrm{c}}(\vec{W}_{I,J,K-1/2})\Delta S_{I,J,K-1/2},}\end{array}\tag{4.29} $$

面面积 $ \Delta S_{I,J,K-1/2} $ 由式（4.8）和（4.9）计算。黏性通量通常采用与对偶控制体格式相同的方法 [10, 11]，从而得到更紧致的格式（即涉及更少节点值的格式）。

分别将式（4.25）、（4.26）及（4.28）、（4.29）给出的各面贡献相加，即可得到全部网格单元的中间残差 $ \vec{R}_{I,J,K} $。为将单元残差与节点残差联系起来，采用残差分配公式作进一步近似。它本质上是一个函数，根据所有共享某节点的单元残差的加权和，计算未知的节点残差。已提出以下分配公式：

• Ni 的体积加权和 [8]；

• Hall 的无权重求和 $ [9] $；

• Rossow 的特征（迎风）加权方法 $ [12, 13] $。

截断误差的理论研究 $ [4] $ 表明，Ni 格式比 Hall 方法更精确。然而，实践中 Ni 分配公式在网格强烈扭曲和拉伸的位置会出现问题。例如，采用 O 型网格时，翼型后缘附近出现强烈压力振荡 $ [4, 14] $；只有显著增大数值黏性，才能收敛。迎风加权方法 $ [12, 13] $ 的基本思想与波动分裂格式 $ [15–18] $（参见第 3.1.5 节）很相似，但数值实现简单得多。其基本做法是，仅沿特征方向向上游传递残差。

三种方法中，Hall 分配格式被证明最稳健。Hall 格式中，某节点的残差就是所有共享该节点的单元中间残差 $ \vec{R}_{I,J,K} $ 的简单和。因此，对于图 4.4 的二维情况，有：

$$ \vec{R}_{i,j}=\vec{R}_{I,J}+\vec{R}_{I-1,J}+\vec{R}_{I-1,J-1}+\vec{R}_{I,J-1}.\tag{4.30} $$

三维情况下，同样需要将八个单元残差相加。仔细考察 $ (4.30) $ 可知，$ \vec{R}_{i,j} $ 正是穿过超单元边界的净通量，

$$ \Omega_{i,j}=\Omega_{I,J}+\Omega_{I-1,J}+\Omega_{I-1,J-1}+\Omega_{I,J-1},\tag{4.31} $$

因为内部面的通量彼此抵消。超单元也就是以点 $ (i,j) $ 为中心的“总”控制体。三维情况下，总控制体由以下单元构成：

$$ \begin{array}{l}{\Omega_{i,j,k}=\Omega_{I,J,K}+\Omega_{I-1,J,K}+\Omega_{I-1,J-1,K}+\Omega_{I,J-1,K}}\\ {\quad+\Omega_{I,J,K-1}+\Omega_{I-1,J,K-1}+\Omega_{I-1,J-1,K-1}+\Omega_{I,J-1,K-1}.}\end{array}\tag{4.32} $$

由图 4.4 可见，控制体至少重叠一个单元，这也是该格式名称的由来。

源项用相应网格节点的流动变量计算，即：

$$ (\vec{Q}\Omega)_{i,j,k}=\vec{Q}(\vec{W}_{i,j,k})\Omega_{i,j,k}.\tag{4.33} $$

采用上述定义后，式（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.34} $$

它构成一组常微分方程，必须像单元中心格式一样，在每个网格点 $ (i,j,k) $ 求解。注意，总控制体积 $ (4.31) $ 或 $ (4.32) $ 用于式 $ (4.34) $。

4.2.3 单元顶点格式：对偶控制体

在这种格式中，控制体以特定网格节点（顶点）为中心，全部流动变量存储于该节点 $ [19, 20] $。如图 4.5 所示，二维对偶控制体通过连接共享该节点的四个单元中点构造；三维中则连接八个单元的形心，形成控制体各面。另一种做法是，从一个单元形心连接到边中点（二维），再连接到相邻单元形心 $ [21] $。这样，控制体的面由法向量不同的两部分组成，这在非结构网格中很常见（见第 5 章）。然而，对于结构网格，只有在边界处采用这种定义才有必要，否则会改变表面的离散。图 4.6 展示了这一情况。

在较光顺网格的内部区域，不能期待第二种做法显著改善精度。因此，下面采用较简单的对偶控制体定义。边界处理见第 8 章。

计算离散流动方程 $ (4.2) $ 时，需要求控制体各面上的对流和黏性通量。可以采用下列三种方法之一：

1. 使用同一个面矢量，分别根据控制体面左右节点的变量计算通量，再取平均（一般仅用于对流通量）；

2. 对控制体面左右两侧节点存储的变量取平均，再计算通量；

3. 将流动物理量分别插值到面的左右两侧，再计算通量（仅用于对流通量）。

例如，对于图 4.5 的单元面 $ \vec{n}_{i+1/2,j} $，第一种方法——通量平均——在二维中写为：

$$ \left(\vec{F}_{c}\Delta S\right)_{i+1/2,j}\approx\frac{1}{2}\left[\vec{F}_{c}(\vec{W}_{i,j})+\vec{F}_{c}(\vec{W}_{i+1,j})\right]\Delta S_{i+1/2,j}\tag{4.35} $$

其中，$ \Delta S_{i+1/2,j} $ 根据已知单元形心坐标，分别按式（4.6）和（4.7）计算。

第二种方法——变量平均——可写为：

$$ \left(\vec{F}\Delta S\right)_{i+1/2,j}\approx\vec{F}(\vec{W}_{i+1/2,j})\Delta S_{i+1/2,j},\tag{4.36} $$

控制体面 $ \vec{n}_{i+1/2,j} $ 上的守恒变量（或依赖变量），由两个相邻节点的值作算术平均获得，即：

$$ \vec{W}_{i+1/2,j}=\frac{1}{2}\left(\vec{W}_{i,j}+\vec{W}_{i+1,j}\right).\tag{4.37} $$

式（4.36）的通量矢量 $ \vec{F} $ 可以表示对流通量，也可以表示黏性通量。

第三种方法将流动物理量（通常为速度分量、压力、密度和总焓）分别插值到面的两侧。插值结果称为左状态和右状态（见第 4.3 节开头），两侧数值通常不同。随后，通过某个非线性函数，根据左、右状态之差计算控制体面通量。因此，

$$ \left(\vec{F}_{c}\Delta S\right)_{i+1/2,j}\approx f_{\mathrm{Flux}}\left(\vec{U}_{\mathrm{L}},\vec{U}_{\mathrm{R}},\Delta S_{i+1/2,j}\right),\tag{4.38} $$

其中，

$$ \begin{array}{l}\vec{U}_{\mathrm{L}}=f_{\mathrm{Interp}}\left(\cdots,\vec{U}_{i-1,j},\vec{U}_{i,j},\cdots\right)\\ \vec{U}_{\mathrm{R}}=f_{\mathrm{Interp}}\left(\cdots,\vec{U}_{i,j},\vec{U}_{i+1,j},\cdots\right)\end{array}\tag{4.39} $$

表示插值状态。当然，对于控制体其他面，同样存在类似于 $ (4.35)-(4.39) $ 的关系。

三维中采用相同近似。例如，在单元面 $ \vec{n}_{i+1/2,j,k} $（例如图 4.1b 的 $ \vec{n}_{2} $）上，式（4.35）的通量平均变为：

$$ \begin{aligned}\big(\vec{F}_{c}\Delta S\big)_{i+1/2,j,k}\approx\frac{1}{2}\big[\vec{F}_{c}(\vec{W}_{i,j,k})+\vec{F}_{c}(\vec{W}_{i+1,j,k})\big]\Delta S_{i+1/2,j,k},\end{aligned}\tag{4.40} $$

其中，$ \Delta S_{i+1/2,j,k} $ 由式（4.8）和（4.9）求得。与式（4.36）类似，流动变量平均给出：

$$ \left(\vec{F}\Delta S\right)_{i+1/2,j,k}\approx\vec{F}(\vec{W}_{i+1/2,j,k})\Delta S_{i+1/2,j,k}\tag{4.41} $$

其中，

$$ \vec{W}_{i+1/2,j,k}=\frac{1}{2}\left(\vec{W}_{i,j,k}+\vec{W}_{i+1,j,k}\right).\tag{4.42} $$

最后，流动变量插值与式（4.38）对应，给出：

$$ \left(\vec{F}_{c}\Delta S\right)_{i+1/2,j,k}\approx f_{\mathrm{Flux}}\left(\vec{U}_{\mathrm{L}},\vec{U}_{\mathrm{R}},\Delta S_{i+1/2,j,k}\right),\tag{4.43} $$

其中，$ \vec{U}_{L} $ 和 $ \vec{U}_{R} $ 表示面上的插值。第 4.3 节将详细介绍对流通量的几种计算方法，第 4.4 节介绍黏性通量的计算。

离散流动方程（4.2）中最后尚待计算的是源项 $ \vec{Q} $。如引言所述，通常假定源项在控制体内恒定，因此用相应网格点的流动变量计算，可定义为：

$$ (\vec{Q}\Omega)_{i,j,k}=\vec{Q}(\vec{W}_{i,j,k})\Omega_{i,j,k}.\tag{4.44} $$

利用上述关系，可以计算各面通量，并按式（4.2）完成对 $ \Omega_{i,j,k} $ 边界的数值积分，从而得到包含源项的完整残差 $ \vec{R}_{i,j,k} $。于是，每个网格点的守恒变量随时间的变化满足：

$$ \frac{\mathrm{d}\vec{W}_{i,j,k}}{\mathrm{d}t}=-\frac{1}{\Omega_{i,j,k}}\vec{R}_{i,j,k}.\tag{4.45} $$

第 6 章将介绍几种合适的求解方法。

4.2.4 单元中心格式与单元顶点格式的比较

前三节介绍了单元中心和单元顶点离散方法。这里比较三种格式，并概述围绕其相对优劣展开的、时而存在争议的讨论。

首先考察离散精度。文献 [4, 22] 的讨论表明，单元顶点格式（无论采用重叠或对偶控制体）在扭曲网格上可构造成一阶精度。在笛卡尔网格或光顺网格（即相邻单元体积变化适中、偏斜较小的网格）上，单元顶点格式可具有二阶或更高阶精度 [23]，具体取决于通量计算格式。相比之下，单元中心格式的离散误差强烈依赖网格光顺性。例如，对于图 4.7 所示的单元排列，即使函数线性变化，平均也不能给出面中点的正确值。因此，在斜率不连续的网格上，即使无限加密网格，离散误差也不会减小。文献 [4] 表明，这种零阶误差表现为等值线的振荡或折转，而单元顶点格式在相同情况下没有问题。不过，在笛卡尔网格或充分光顺的网格上，单元中心格式也可达到二阶或更高阶精度。

文献 $ [24–27] $ 对离散误差作了进一步分析。

其次比较三种方法在边界处的特性。采用对偶控制体的单元顶点格式，主要在固壁边界处遇到困难。再次观察图 4.6，可见边界处只剩约一半控制体。对各面通量积分所得残差，实际位于控制体内部，理想位置是形心；但该残差却被赋给直接位于壁面的节点。相较于单元中心格式，这种位置不匹配会增大离散误差。对偶控制体的定义也会在尖角（如后缘）处产生问题，表现为非物理的压力或密度尖峰。在坐标切口或周期边界等处（见第 8 章），还必须正确汇总控制体两部分的通量，带来进一步复杂性。全部单元顶点格式都需要额外逻辑，以确保多个网格块共用边界点的解相容。单元中心格式不存在这些问题。

采用重叠控制体的单元顶点格式，在壁面边界处理上优于对偶控制体格式，但不能与 TVD、AUSM 或 CUSP 等常用迎风离散方法结合。其离散涉及的点比单元中心和对偶控制体格式更多（三维中为 27 点而非 7 点），导致间断被抹平，并增加隐式时间离散的内存开销。

单元中心与单元顶点格式的最后一项主要差异出现在非定常问题中。如第 3.2 节所述，单元顶点格式至少需要近似处理质量矩阵 $ [28, 29] $。相反，由于残差自然对应控制体形心，单元中心格式可完全略去质量矩阵。

总体而言，在定常流场内部，对偶控制体单元顶点格式与单元中心格式的数值表现十分相近。主要差异出现在扭曲网格、边界处理和非定常流动中。后两种情况中，单元中心方法更有优势，使流动求解器的实现更直接。