← 回到首页 Skip to main content

5.2 一般离散方法

5.2 一般离散方法

本章开头已提到，定义控制体及确定流动变量位置的两种最常用方法，是单元中心格式和中位对偶格式。本节将更详细地介绍它们。

不过，首先简要讨论非结构网格流动求解器所需的基本数据结构。事实上，既灵活又能节省内存和运算次数的数据结构，是任何非结构格式的关键。可以说，网格本身缺少的结构，必须在流动求解器内部补足。至少需要以下数据：

  • 网格节点（顶点）的坐标；

  • 从单元指向网格节点的指针；

  • 从边界上的单元面指向网格节点的指针。

离散格式需要的其他数据结构，可由这些信息生成。为说明上述数据的一种存储方式，以图 5.4(a) 的四面体为例。若进一步假设面 1–2–4 位于某个边界（壁面、入口、远场等）上，则可采用以下格式：

# 节点（x、y、z）：
P1.x P1.y P1.z
P2.x P2.y P2.z
P3.x P3.y P3.z
P4.x P4.y P4.z
# 四面体：
...
P1 P2 P3 P4
...
# 边界：
...
边界类型 P1 P4 P2
...

随书源代码提供的二维非结构网格流动求解器也采用类似格式。

关于计算域边界，必须注意两点。第一，存储边界面比仅存储边界节点更方便。图 5.5 所示情形说明了原因：节点 $ P_{1} $ 为三个边界共享，节点 $ P_{2} $ 和 $ P_{3} $ 为两个边界共享，而这些边界的物理类型可能不同。因此，施加正确的边界条件可能变得十分繁琐。相比之下，一个面只能属于一个边界，例如面 $ P_{1}-P_{2}-P_{4} $ 只属于边界 1。

第二点涉及边界面节点的编号。必须按统一规则编号，例如从流动域外部观察时采用逆时针顺序，使所有面矢量（式（5.6）、（5.11）或（5.13））统一朝外或朝内。这不仅对通量积分重要，对网格生成也很重要。

Image
图 5.5 三个边界在同一角点相交，导致网格点所属边界类型不唯一。图内 Boundary 1／2／3：边界 1／2／3。

5.2.1 单元中心格式

如图 5.6 所示，若控制体与网格单元重合，且流动变量位于单元形心，则称为单元中心格式。计算离散流动方程（5.2）时，需要给出控制体面中点处的对流通量和黏性通量。对光滑网格上的二阶精度离散，这已足够；四边形面采用平均法向量。通量可按以下三种方式之一近似：

  1. 分别用单元面左右两侧网格单元形心处的数值计算通量，再取平均；两侧计算采用同一单位法向量。这一方法一般只用于对流通量。

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

  3. 根据周围单元的数值，分别在单元面两侧重构流动量，再计算通量。这一方法只用于对流通量。

以图 5.6 中单位法向量为 $ \vec{n}_{01} $ 的单元面为例，第一种方法（通量平均）在二维情况下为：

$$ \left(\vec{F}_{c}\Delta S\right)_{01}\approx\frac{1}{2}\left[\vec{F}_{c}(\vec{W}_{0},\vec{n}_{01})+\vec{F}_{c}(\vec{W}_{1},\vec{n}_{01})\right]\Delta S_{01}\tag{5.18} $$

其中，面面积 $ \Delta S_{01} $ 由式（5.6）和式（5.7）计算。

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

$$ \left(\vec{F}\Delta S\right)_{01}\approx\vec{F}(\vec{W}_{01},\vec{n}_{01})\Delta S_{01},\tag{5.19} $$

Image
图 5.6 二维单元中心格式的控制体。圆圈表示网格节点，矩形（C）表示单元中心。

其中，单位法向量为 $ \vec{n}_{01} $ 的面上的守恒变量／依赖变量，定义为两相邻单元对应数值的算术平均，即：

$$ \vec{W}_{01}=\frac{1}{2}\left(\vec{W}_{0}+\vec{W}_{1}\right).\tag{5.20} $$

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

第三种方法首先将流动量（通常为速度分量、压力、密度和总焓）分别插值到单元面的两侧。重构所得的量称为左状态和右状态（见 5.3.3 节），二者一般不同。随后，通过适当的非线性函数，根据左右状态之差计算穿过该面的通量。因此有：

$$ \left(\vec{F}_{c}\Delta S\right)_{01}\approx f_{\mathrm{Flux}}\left(\vec{U}_{\mathrm{L}},\vec{U}_{\mathrm{R}},\Delta S_{01}\right),\tag{5.21} $$

其中，

$$ \begin{aligned}\vec{U}_{\mathrm{L}}&=f_{\mathrm{Rec}}\left(\cdots,\vec{U}_{2},\vec{U}_{0},\cdots\right)\\\vec{U}_{\mathrm{R}}&=f_{\mathrm{Rec}}\left(\cdots,\vec{U}_{1},\vec{U}_{0},\cdots\right)\end{aligned}\tag{5.22} $$

表示重构后的状态。

当然，控制体其他面也有类似关系。上述近似同样可用于三维情况；此时面矢量 $ \vec{S} $ 分别由式（5.11）或式（5.13）计算。

本节引言已指出，描述单元的基本数据结构必须作适当扩充，才能支持离散方法。前文可见，数值运算主要利用单元（控制体）的面，以及相邻单元中心处的数值。因此，空间离散自然适合采用基于面的数据结构。对网格中的每个面，该结构存储以下信息（见图 5.6）：

  • 指向共享该面的两个单元的指针，用于访问这两个单元 $ (C_{0}, C_{1}) $ 的流动变量；

  • 面矢量 $ (\vec{S}_{01} = \vec{n}_{01} \Delta S_{01}) $，其方向必须统一朝外或朝内；

  • 分别从两个单元的形心指向面中点 $ M_{01} $ 的矢量，即 $ (C_{0}-M_{01}) $、$ (C_{1}-M_{01}) $，用于将流动变量准确地插值到面上。纯四面体网格不需要这些矢量，因为可以使用简单的外推公式 [2, 26]（参见式（5.44））。

因此，通量积分（例如式（5.19））可通过遍历网格中所有内部面和边界面的循环实现：

$$ \begin{aligned}&\text{遍历面 }\mathrm{face}=1,\ldots,\mathrm{nfaces}:\\&\quad I=\texttt{pointer\_to\_left\_cell}(\mathrm{face}),\\&\quad J=\texttt{pointer\_to\_right\_cell}(\mathrm{face}),\\&\quad(\vec F\Delta S)_{IJ}\approx\vec F(\vec W_{IJ},\vec n_{IJ})\Delta S_{IJ},\\&\quad\vec R_I=\vec R_I+(\vec F\Delta S)_{IJ},\\&\quad\vec R_J=\vec R_J-(\vec F\Delta S)_{IJ}.\\&\text{循环结束。}\end{aligned} $$

循环结束并加上源项 $ \vec{Q}_{I}\Omega_{I} $ 后，即得到所有单元的最终残差 $ (\vec{R}) $。若改为遍历单元，效率会较低，因为面矢量必须存储两次，通量也要计算两次，边界处除外。此外，采用完全相同的面矢量 $ \vec{S}_{IJ} $ 计算流入控制体 $ \Omega_{I} $ 和 $ \Omega_{J} $ 的部分通量，控制方程的守恒性质便会自动得到保证。

5.2.2 中位对偶单元顶点格式

在单元顶点格式中，流动变量位于网格节点（顶点）。将共享某一节点的全部单元的形心、面中点和边中点相连，即形成中位对偶控制体。四面体的构造见图 5.7(a)，六面体见图 5.7(b)。这样，每个网格节点周围形成一个多面体包络；图 5.8 以二维混合单元网格示意这一构造。该多面体可视为对偶网格，格式因此得名。值得注意的是，中位对偶有限体积离散与采用线性单元的 Galerkin 有限元格式等价，例如见文献 [37]。

计算离散流动方程 $ (5.2) $ 时，需要在控制体表面上对对流通量与黏性通量积分。因此，

Image
图 5.7 四面体 (a) 和六面体 (b) 的中位对偶格式：部分控制体及其面（阴影）。P 为网格节点，C 为单元形心（式（5.17）），F 为面形心（式（5.8）或（5.9）），M 为边中点。阴影区域分别表示对应于边 $ P_{1}-P_{3} $ 或 $ P_{1}-P_{5} $ 的控制体面的一部分。
Image
图 5.8 二维中位对偶格式的控制体。$ C_{1} $、$ C_{2} $ 等为单元中心；$ P_{1} $、$ P_{2} $ 等为网格节点。对应于边 $ P_{0}-P_{1} $ 的面用粗线表示。

原则上应分别计算每个局部面（例如图 5.7(a) 中的 $ F_{1}-M_{13}-F_{2}-C $）的通量。不过，只有三阶及更高阶精度的离散才必须如此 [38, 39]。对最常采用的二阶格式，可以假设围绕某一边的全部面上的流动变量相同，随后利用边两端节点的变量与梯度，在边中点计算通量。这样，便可为每条边定义平均单位法向量和总面面积。以图 5.8 中的边 $ P_{0}-P_{1} $ 为例，其平均法向量为：

$$ \vec{n}_{01}=\vec{n}_{L}+\vec{n}_{R},\tag{5.23} $$

校注：原书直接相加两个局部法向量，却称结果为平均单位法向量，存在归一化疑点。计算时应由合成面矢量 $\vec S=\vec n_L\Delta S_L+\vec n_R\Delta S_R$ 得到方向，并检查 $\vec n=\vec S/|\vec S|$ 与面面积的配套定义；原文给出的两局部面积直接相加在非共面情况下也应核对。

总面面积为 $ \Delta S_{01} = \Delta S_{L} + \Delta S_{R} $。三维情况也一样，平均法向量由共享该边中点的所有局部面求和得到，见图 5.9。二维面矢量 $ (\vec{S} = \vec{n} \Delta S) $ 由式（5.6）计算。三维局部面始终为四边形，可以将其分成三角形并使用式（5.11），也可以采用式（5.13）的简化处理；对光滑网格，后一种处理已足够。单元形心及面中心分别由式（5.17）、（5.8）或（5.9）求得。

随后可采用以下三种方法之一计算通量：

  1. 用边两端节点的数值分别计算通量，再取平均；两侧采用同一平均单位法向量。这通常只用于对流通量。

  2. 对边两端节点的变量取平均，再用平均变量计算通量。

  3. 根据周围节点的数值，分别在控制体面的两侧重构流动量，再计算通量。这只用于对流通量。

Image
图 5.9 三维中位对偶单元顶点格式中，对应于边 ij 的总面面积和平均单位法向量。

通量计算在形式上与单元中心格式相同。因此，式 $ (5.18)-(5.22) $ 也适用于中位对偶格式。若采用上述每条边对应一个平均单位法向量的方法，则空间离散最高效的实现方式是采用基于边的数据结构。该结构为网格中的每条边存储以下信息（参见图 5.9）：

  • 指向定义该边的两个节点的指针，用于访问控制体 $ \Omega_{i} $ 和 $ \Omega_{j} $ 的流动变量；

  • 面矢量 $ (\vec{S}_{ij} = \vec{n}_{ij} \Delta S_{ij}) $，其方向必须统一朝外或朝内；

  • 从节点 i 指向节点 j 的边矢量，用于将流动变量插值到面上，即解的重构。也可在计算过程中直接由节点坐标求出边矢量。采用人工耗散的标准中心格式不需要这一矢量（5.3.1 节）。

于是，通量积分（例如式（5.19））可通过遍历网格全部边的循环实现：

$$ \begin{aligned}&\text{遍历边 }\mathrm{edge}=1,\ldots,\mathrm{nedges}:\\&\quad i=\texttt{pointer\_to\_left\_node}(\mathrm{edge}),\\&\quad j=\texttt{pointer\_to\_right\_node}(\mathrm{edge}),\\&\quad(\vec F\Delta S)_{ij}\approx\vec F(\vec W_{ij},\vec n_{ij})\Delta S_{ij},\\&\quad\vec R_i=\vec R_i+(\vec F\Delta S)_{ij},\\&\quad\vec R_j=\vec R_j-(\vec F\Delta S)_{ij}.\\&\text{循环结束。}\end{aligned} $$

循环结束并加上源项 $ \vec{Q}_{i}\Omega_{i} $ 后，便得到所有节点的最终残差 $ (\vec{R}) $。相比于逐个控制体分别对通量求和，这种方法高效得多，因为每个平均面矢量只存储一次，每条边也只访问一次，而非两次。此外，由于采用完全相同的平均面矢量 $ \vec{S}_{ij} $ 计算流入控制体 $ \Omega_{i} $ 和 $ \Omega_{j} $ 的部分通量，质量、动量和能量保持严格守恒。

5.2.3 单元中心格式与中位对偶格式的比较

单元中心格式和中位对偶格式各自的优缺点一直存在争论，主要原因是：针对实际构型，缺乏对两种方法在精度、计算时间及内存方面的公平比较。这里从以下方面归纳支持或反对各方法的主要论据：

  • 精度；

  • 计算量；

  • 内存需求；

  • 灵活性。

这有助于更深入地理解各格式固有的问题，并为预期应用范围选择合适的格式。

精度

在三角形／四面体网格上，单元中心格式的控制体数及自由度数约为中位对偶格式的两倍／六倍 $ [37] $。在典型的四面体与棱柱混合网格上，单元中心格式的未知量约为中位对偶格式的三倍。这似乎表明，在同一网格上，单元中心格式比单元顶点离散更精确。不过，单元中心格式每个残差所包含的通量数明显少于中位对偶格式（四面体网格上约为 3 对 7），这可能损害精度。因此，尚无明确证据表明哪种格式更优。

中位对偶格式在拉伸的三角形和四面体网格上存在一个特殊问题。以图 5.10 为例，图中网格由直角三角形铺成，黏性流动的固壁附近常采用这种网格。

Image
(a)
Image
(b)
图 5.10 拉伸直角三角形网格上的中位对偶控制体 (a) 与包容对偶控制体 (b) 的比较。
Image
Image
Image
图 5.11 锐角三角形 (a) 与钝角三角形 (c) 的部分包容对偶网格（虚线）$ [41] $。包容圆是能够包含该三角形的最小圆。对钝角三角形，其圆心位于最长边上。图内 Containment circle：包容圆。

图 5.10(a) 中可见，面 $ \Delta S_{ij} $ 相对于边 ij 严重倾斜。但空间离散格式大多假设通量方向与面正交，Riemann 求解器尤其如此。因此会引入误差，对一阶格式，这种误差尤为显著 [39]。采用所谓包容对偶控制体（containment dual）可改善这一状况 [40]。如图 5.11 所示，包容对偶方法用最小包容圆／球的中心代替单元形心来定义面，从而得到与四边形网格相同的控制体（图 5.10(b)）。注意，像 $ i j^{\prime} $ 这样的对角边没有对应的面面积。这需要增加预处理工作，但可明显改善解的精度 [41]。当然，也可直接在边界层内使用四边形或六面体单元。网格引起误差的更多讨论见文献 [22, 23]。

中位对偶格式的另一个固有问题，是物理域边界处的离散。边界处只剩约一半控制体（参见图 4.6）。沿各面计算通量积分，所得残差应位于控制体内部，理想情况下位于其中心；实际却将残差赋给直接位于边界上的节点。与单元中心格式相比，这种位置不一致会增大离散误差，对固壁尤为不利。对偶控制体的定义还会在尖角处（例如尾缘）造成问题，表现为压力或密度出现非物理峰值。周期边界还存在进一步的困难（见 8.8 节），必须正确累加控制体两部分的通量。

控制体形心与残差存储节点的位置不一致，对中位对偶格式还会产生另一后果：非定常流动中出现质量矩阵。3.2 节开头已讨论这一点。单元中心格式的优点，是可从方程中消去质量矩阵而不损害解的精度；中位对偶格式则必须对质量矩阵进行特殊处理 $ [42, 43] $。

计算量

评价两种格式的计算量，首先应考虑通量积分。前文已知，单元中心格式遍历单元面，中位对偶格式则遍历边。由于两者在界面上计算通量的过程很相似，单元面数与边数之比便给出了计算量之比。在四面体网格上，每个共享面只计一次时，单元面数约为边数的两倍，因此同一网格上单元中心格式的计算量约为中位对偶格式的两倍 $ [37] $。不过，在包含棱柱单元的混合网格上，单元中心方法更有竞争力。在六面体网格上，面数等于边数，除边界处理外，两者计算量相当。

内存需求

从内存需求看，与中位对偶格式相比，单元中心格式在四面体网格上需存储约六倍的流动变量，在标准混合网格上约为三倍。此外，两种格式对每个单元面或每条边都要存储两个整数和三个实数，分别对应指针和面矢量。单元中心格式还需为每个面存储两个指向面中点的矢量，即六个实数；中位对偶格式则可仅用节点坐标完成相关计算，所需数值明显更少。因此，平均而言，单元中心格式所需计算机内存超过中位对偶方法的两倍。

网格生成与自适应

当单元界面不相容时，例如图 3.4 中字母 F 所在的位置，单元中心格式具有明显优势。与中位对偶方法不同，它无需采用复杂且昂贵的特殊过程计算界面通量。因此，在网格生成、网格自适应和网格运动中均具有更大的灵活性。采用嵌套 笛卡尔 网格的流动求解器也利用了这一优点。

Back to top