9.4 多重网格

9.4 多重网格

多重网格是一种强有力的加速技术（参见图 3.7）。它在一系列逐渐变粗的网格上求解控制方程，然后将粗网格的解更新量合并并加到最细网格的解上。Brandt $ [19] $ 最初为椭圆型偏微分方程发展了该技术，随后 Jameson 将其用于 Euler 方程 $ [20–22] $，之后又用于 Navier–Stokes 方程 $ [11–13, 23–32] $。多重网格可以结合显式或隐式时间推进格式实现 $ [33–39] $。当前研究旨在显著提高多重网格方法处理双曲型流动问题的效率 $ [40–43] $。

多重网格的基本思路是利用粗网格，使最细网格的解更快达到定常状态。为此利用两种效应：

1. 粗网格控制体较大，允许采用更大的时间步长，同时计算量较少。确定新解的工作主要分配到粗网格，因此收敛更快，计算时间减少。

2. 多数显式、隐式时间推进及迭代格式主要能高效消除解误差的高频分量（见 10.3 节），低频分量通常很难得到阻尼。因此，初始阶段最大的误差被消除之后，定常收敛变得缓慢。多重网格可在此发挥作用：最细网格上的低频分量在粗网格上变为高频分量，从而逐级得到阻尼。结果是在整个频谱内迅速降低误差，显著加快收敛。

由此可见，多重网格的成功高度依赖于时间推进或迭代格式对高频误差分量的良好阻尼。

除几何多重网格，即基于网格的多重网格外，还可以采用代数多重网格（AMG）方法 $ [44-53] $。AMG 专为隐式格式发展，直接作用于方程组矩阵，即左端算子。其基本思路是采用粗化矩阵降低隐式算子的维数及方程数量，然后求解表示粗层的降阶方程组，得到细层解的修正量。构造粗化矩阵时，将耦合最强的方程相加，即系统矩阵非对角项最大的方程。因此，粗层的生成完全由

流动问题的物理特性决定，而不依赖于网格。AMG 的优势在于无需构造或存储粗网格拓扑，对复杂非结构网格尤其有利。

9.4.1 基本多重网格循环

应用几何多重网格之前，必须生成粗网格。标准方法是沿所有坐标方向均匀粗化，但 Mulder $ [54] $ 提出了称为半粗化的方法：每次只沿一个方向粗化，从一层粗网格到另一层时再改变方向。半粗化特别适用于控制方程沿某个空间方向具有刚性的流动，例如边界层的壁面法向。将半粗化用于 Navier–Stokes 方程的实例见文献 $ [26, 55–57] $。

根据关系式（6.1），细网格上的离散控制方程为

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

下面用下标 h 表示最细网格，h 为网格线间距，非结构网格中则为控制体的特征尺寸。从已知解 $ \vec{W}_{h}^{n} $ 出发，采用适宜迭代格式推进一个时间步后，得到新解 $ \vec{W}_{h}^{n+1} $，并用其计算新残差 $ \vec{R}_{h}^{n+1} $。为了利用粗网格改善解 $ \vec{W}_{h}^{n+1} $，执行以下三个步骤。

将解及残差传递到粗网格

通过下列插值，将解传递到粗网格：

$$ \vec{W}_{2h}^{(0)}=\hat{I}_{h}^{2h}\vec{W}_{h}^{n+1},\tag{9.17} $$

其中，下标 2h 表示粗网格¹，$ \hat{I}_{h}^{2h} $ 为插值算子。还必须将残差传递到粗网格，使其低频误差分量得到光顺。为此采用守恒传递算子，即控制体尺寸增大多少，残差也必须按相同比例增大。细网格残差还用于在粗网格上保持细网格解的精度。为此，

构造一个称为驱动函数 [19,22] 的源项，其为细网格传递来的残差与粗网格上采用初始解 $ \vec{W}_{2h}^{(0)} $（式（9.17））计算的残差之差，即

$$ (\vec{Q}_{F})_{2h}=I_{h}^{2h}\vec{R}_{h}^{n+1}-\vec{R}_{2h}^{(0)}.\tag{9.18} $$

其中，$ I_{h}^{2h} $ 为将残差从细网格传递到粗网格的限制算子。这类多重网格称为全近似存储（FAS）方法 [19]，特别适用于非线性方程，因为通过重新离散，将系统非线性传递到粗层。

1 所用术语不应按严格的几何意义理解。在非结构网格上，细网格和粗网格的控制体特征尺寸会有差异。半粗化情形也如此。

在粗网格上计算新解

粗网格上的解采用与细网格相同的方法计算，并将式（9.18）的驱动函数加到粗网格残差上，即

$$ (\vec{R}_{F})_{2h}=\vec{R}_{2h}+(\vec{Q}_{F})_{2h}.\tag{9.19} $$

因此，式（6.1）的时间推进格式可写为

$$ \frac{\mathrm{d}}{\mathrm{d}t}\vec{W}_{2h}=-\frac{1}{\Omega_{2h}}(\vec{R}_{F})_{2h}=-\frac{1}{\Omega_{2h}}\left[\vec{R}_{2h}+(\vec{Q}_{F})_{2h}\right].\tag{9.20} $$

对于显式多级格式（6.1.1 节），得到

$$ \vec{W}_{2h}^{(k)}=\vec{W}_{2h}^{(0)}-\alpha_{k}\frac{\Delta t_{2h}}{\Omega_{2h}}\left[\vec{R}_{2h}^{(k-1)}+(\vec{Q}_{F})_{2h}\right],\quad k=1,\cdots,m\tag{9.21} $$

$$ \vec{W}_{2h}^{n+1}=\vec{W}_{2h}^{(m)} $$

这与式（6.5）一致。需要注意，第一次迭代时，即式（9.21）的第一级，$ (\vec{R}_{F})_{2h} $ 与细网格传递来的残差相同，即式（9.18）的 $ (\vec{R}_{F})_{2h}=I_{h}^{2h}\vec{R}_{h}^{n+1} $。这保证粗网格解取决于细网格残差，从而保留细网格精度。

一个重要问题是粗网格空间离散的精度。由于粗网格不影响细网格解的精度，一阶格式已足够。与高阶格式相比，粗网格一阶离散具有更好的稳健性与阻尼性能，计算成本也更低。

从粗网格向细网格插值解

在粗网格上执行一个或多个时间步（迭代）后，计算相对于插值得到的初始解（式（9.17））的修正量，即所谓粗网格修正：

$$ \delta\vec{W}_{2h}=\vec{W}_{2h}^{n+1}-\vec{W}_{2h}^{(0)}.\tag{9.22} $$

将粗网格修正插值到细网格，以改善细网格解。因此，新的细网格解为

$$ \vec{W}_{h}^{+}=\vec{W}_{h}^{n+1}+I_{2h}^{h}\delta\vec{W}_{2h},\tag{9.23} $$

其中，$ I_{2h}^{h} $ 称为延拓算子。

9.4.2 多重网格策略

上述基本多重网格仅含一层粗网格。若有多层粗网格，则重复步骤 1 和 2，直到最粗网格。必须注意，粗网格驱动函数由式（9.19）的修正残差经过限制后构造。例如，粗网格 4h 的驱动函数为

$$ (\vec{Q}_{F})_{4h}=I_{2h}^{4h}(\vec{R}_{F})_{2h}^{n+1}-\vec{R}_{4h}^{(0)}=I_{2h}^{4h}\left[\vec{R}_{2h}^{n+1}+(\vec{Q}_{F})_{2h}\right]-\vec{R}_{4h}^{(0)}.\tag{9.24} $$

这样，最细网格的残差控制所有粗网格上解的精度。在最粗网格执行指定数量的时间步后，依次重复步骤 3，直到再次返回最细网格。这一过程称为锯齿形循环或 V 循环（见图 9.5a）。也可以在粗网格上执行更多循环，即图 9.5b 的 W 循环；它尤其常用于跨声速流动。对于超声速及高超声速流动，V 循环更高效。

时间步数

限制前及延拓后的最佳时间步数取决于时间推进格式。对于显式多级格式（6.1.1 或 6.1.2 节），通常仅在限制残差前执行一个时间步，延拓后不执行时间步。不过，将粗网格修正（式（9.22））加到细网格解 $ \widetilde{W}_{h}^{n+1} $（式（9.23））前进行光顺，可以提高多重网格的稳健性。这里采用与 9.3 节相同的中心隐式光顺，系数取常数。

另一种常用时间推进方法是隐式 LU-SGS 格式（见 6.2.4 节），其在限制前执行两次迭代可获得最佳多重网格效率 $ [35] $。延拓后的时间步数取决于空间离散。中心格式（4.3.1 节）不需要时间步 $ [35–37] $，但可以光顺解修正；迎风空间离散则应在延拓后执行一个时间步。这种 $ (2,1) $ 策略在多种流动条件下都能在稳健性与计算时间之间取得最佳效果 $ [36, 37] $。

起始网格

实际应用中，多重网格通常并不直接从最细网格开始，而是在某个粗网格上执行若干多重网格循环，然后将近似解插值到下一层较细网格（采用与延拓相同的算子），再执行若干循环，并继续向下一层较细网格插值，直到最细网格。这样，仅需适度计算成本，就能在最细网格上获得良好初始解。这一高效过程称为全多重网格（FMG）方法 $ [19] $。

传递算子的精度

限制算子（9.18）和延拓算子（9.23）必须满足一定精度要求，即 $ [58] $

$$ m_{\mathrm{R}}+m_{\mathrm{P}}>m_{\mathrm{E}},\tag{9.25} $$

其中，$ m_{R} $ 和 $ m_{P} $ 分别表示限制算子及延拓算子能够精确插值的多项式次数加 1。例如，线性插值时，$ m_{R} $ 或 $ m_{P} $ 均等于 2。此外，$ m_{E} $ 表示控制方程的阶数，因此 Euler 方程取 $ m_{E} = 1 $，Navier–Stokes 方程取 $ m_{E} = 2 $。若违反条件（9.25），限制或延拓引入的附加误差会扰动细网格解，导致多重网格收敛缓慢，甚至发散。

9.4.3 结构网格上的实现

结构网格上实现多重网格十分直接，因为沿相应坐标方向每隔一条删除一条网格线，即可生成粗网格，间距为 2h、4h 等。因此，与最细网格相比，粗网格计算成本保持较低。代表性示例见文献 $ [11–13, 20–22, 26, 33–37] $。

由图 9.6 可以看出，单元中心和节点中心（单元顶点）格式中的解插值、残差限制及修正延拓算子必须分别定义。例如，相邻两层网格每隔一个网格点相互重合，而单元中心的位置始终不同。因此，下面分别讨论单元中心及节点中心（等同于单元顶点）有限体积格式的标准传递算子。

除下面介绍的对称、仅由几何确定的限制和延拓算子外，文献 [16] 第4章还提出了偏迎风形式。迎风限制和延拓均计及流动方程的特征，可提高高超声速流动多重网格的稳健性。为节省篇幅，这里仅讨论节点中心空间离散的迎风延拓算子实现。

单元中心格式的传递算子

采用体积加权插值，将解从细网格传递到粗网格。在二维情况下，式（9.17）为（见图 9.7a）

$$ \begin{aligned}(\vec{W}_{2h}^{(0)})_{I,J}&=\frac{(\vec{W}_{h}^{n+1})_{I,J}\Omega_{I,J}+(\vec{W}_{h}^{n+1})_{I+1,J}\Omega_{I+1,J}+(\vec{W}_{h}^{n+1})_{I,J+1}\Omega_{I,J+1}}{\Omega_{I,J}+\Omega_{I+1,J}+\Omega_{I,J+1}+\Omega_{I+1,J+1}}\\&\quad+\frac{(\vec{W}_{h}^{n+1})_{I+1,J+1}\Omega_{I+1,J+1}}{\Omega_{I,J}+\Omega_{I+1,J}+\Omega_{I,J+1}+\Omega_{I+1,J+1}}.\end{aligned}\tag{9.26} $$

三维采用相似的传递算子，求和包括构成一个粗网格单元的八个细网格控制体。

限制算子定义为一个粗网格控制体内所有单元残差之和。因此，二维时有（参见图 9.7a）

$$ (I_{h}^{2h}\vec{R}_{h}^{n+1})_{I,J}=(\vec{R}_{h}^{n+1})_{I,J}+(\vec{R}_{h}^{n+1})_{I+1,J}+(\vec{R}_{h}^{n+1})_{I,J+1}+(\vec{R}_{h}^{n+1})_{I+1,J+1}\tag{9.27} $$

三维亦同。式（9.27）中位于物理区域外部的残差 $ \vec{R}_{h}^{n+1} $ 必须置零。

式（9.23）的粗网格修正延拓有两种方法。首先，若将 $ \delta\vec{W}_{2h} $ 均匀分配给周围所有单元中心，如图 9.7b 所示，则得到零阶延拓算子。例如，

$$ (\vec{W}_{h}^{+})_{I,J+1}=(\vec{W}_{h}^{n+1})_{I,J+1}+(\delta\vec{W}_{2h})_{I,J},\quad\text{等等。}\tag{9.28} $$

第二种方法可使多重网格收敛更快，分为两步：首先像单元顶点格式那样，将 $ \delta\vec{W}_{2h} $ 插值到网格节点（见下文）；随后对节点值取平均，得到细网格单元中心的值。根据图 9.7c，可以推导出最终关系：

$$ \begin{array}{l}(I_{2h}^{h}\delta\vec{W}_{2h})_{I,J}=\frac{1}{16}\left[9(\delta\vec{W}_{2h})_{I,J}+3(\delta\vec{W}_{2h})_{I-1,J}\right.\\ \left.\quad+3(\delta\vec{W}_{2h})_{I,J-1}+(\delta\vec{W}_{2h})_{I-1,J-1}\right].\end{array}\tag{9.29} $$

对应三维表达式可用类似过程得到，即

$$ \begin{aligned}(I_{2h}^{h}\delta\vec{W}_{2h})_{I,J,K}&=\frac{1}{64}\left[27(\delta\vec{W}_{2h})_{I,J,K}+9(\delta\vec{W}_{2h})_{I-1,J,K}\right.\\&\quad+\left.9(\delta\vec{W}_{2h})_{I,J-1,K}+9(\delta\vec{W}_{2h})_{I,J,K-1}\right.\\&\quad+\left.3(\delta\vec{W}_{2h})_{I-1,J-1,K}+3(\delta\vec{W}_{2h})_{I-1,J,K-1}\right.\\&\quad+\left.3(\delta\vec{W}_{2h})_{I,J-1,K-1}+(\delta\vec{W}_{2h})_{I-1,J-1,K-1}\right].\end{aligned}\tag{9.30} $$

单元顶点格式的传递算子

由于细网格和粗网格具有共同节点，可以通过直接注入传递解，即

$$ (\vec{W}_{2h}^{(0)})_{i,j,k}=(\vec{W}_{h}^{n+1})_{i,j,k}.\tag{9.31} $$

标准中心限制算子从构成一个粗网格单元的四个细网格单元（三维时为八个）的全部节点作线性插值。根据图 9.8，二维限制残差计算为

$$ \begin{aligned}(I_{h}^{2h}\vec{R}_{h}^{n+1})_{i,j}=&(\vec{R}_{h}^{n+1})_{i,j}+\frac{1}{2}\left[(\vec{R}_{h}^{n+1})_{i-1,j}+(\vec{R}_{h}^{n+1})_{i+1,j}\right.\\&+\left.(\vec{R}_{h}^{n+1})_{i,j-1}+(\vec{R}_{h}^{n+1})_{i,j+1}\right]\\&+\frac{1}{4}\left[(\vec{R}_{h}^{n+1})_{i-1,j-1}+(\vec{R}_{h}^{n+1})_{i+1,j-1}\right.\\&+\left.(\vec{R}_{h}^{n+1})_{i-1,j+1}+(\vec{R}_{h}^{n+1})_{i+1,j+1}\right].\end{aligned}\tag{9.32} $$

三维细网格残差按下式汇集：

$$ (I_{h}^{2h}\vec{R}_{h}^{n+1})_{i,j,k}=(\vec{R}_{h}^{n+1})_{i,j,k}+\frac{1}{2}\mathcal{A}+\frac{1}{4}\mathcal{B}+\frac{1}{8}\mathcal{C}\tag{9.33} $$

其中，各系数为

$$ \begin{aligned}\mathcal{A}&=(\vec{R}_{h}^{n+1})_{i+1}+(\vec{R}_{h}^{n+1})_{i-1}+(\vec{R}_{h}^{n+1})_{j+1}+(\vec{R}_{h}^{n+1})_{j-1}+(\vec{R}_{h}^{n+1})_{k+1}+(\vec{R}_{h}^{n+1})_{k-1},\\\mathcal{B}&=(\vec{R}_{h}^{n+1})_{i+1,j+1}+(\vec{R}_{h}^{n+1})_{i-1,j+1}+(\vec{R}_{h}^{n+1})_{i+1,j-1}+(\vec{R}_{h}^{n+1})_{i-1,j-1}\\&\quad+(\vec{R}_{h}^{n+1})_{i+1,k+1}+(\vec{R}_{h}^{n+1})_{i-1,k+1}+(\vec{R}_{h}^{n+1})_{i+1,k-1}+(\vec{R}_{h}^{n+1})_{i-1,k-1}\\&\quad+(\vec{R}_{h}^{n+1})_{j+1,k+1}+(\vec{R}_{h}^{n+1})_{j-1,k+1}+(\vec{R}_{h}^{n+1})_{j+1,k-1}+(\vec{R}_{h}^{n+1})_{j-1,k-1},\\\mathcal{C}&=(\vec{R}_{h}^{n+1})_{i+1,j+1,k+1}+(\vec{R}_{h}^{n+1})_{i-1,j+1,k+1}+(\vec{R}_{h}^{n+1})_{i-1,j-1,k+1}+(\vec{R}_{h}^{n+1})_{i+1,j-1,k+1}\\&\quad+(\vec{R}_{h}^{n+1})_{i+1,j+1,k-1}+(\vec{R}_{h}^{n+1})_{i-1,j+1,k-1}+(\vec{R}_{h}^{n+1})_{i-1,j-1,k-1}+(\vec{R}_{h}^{n+1})_{i+1,j-1,k-1}\end{aligned}\tag{9.34} $$

上式仅列出了与 i、j、k 不同的下标。限制之前，必须将所有物理边界点和虚拟点上的残差 $ \vec{R}_{h}^{n+1} $ 置零。

粗网格修正的延拓可通过遍历粗网格点实现。在循环中，采用与限制相同的权重，将值 $ (\delta\vec{W}_{2h})_{i,j,k} $ 分配给细网格点（参见图 9.8b），再累加各部分贡献，得到细网格每个点的完整传递修正量。

迎风延拓（单元顶点格式）

原则上，迎风延拓可以用特征变量或守恒变量表示。文献 $ [16] $ 提出了一种特别高效的守恒变量实现，根据马赫数和速度方向，对式 $ (9.23) $ 的修正量作偏迎风插值。

计算成本很低，却能显著提高高马赫数流动中多重网格的稳健性 $ [16] $。

解修正量 $ \delta\tilde{W}_{2h} $ 向细网格的插值分两步。第一步，将两层网格共有的 A、B、C、D 点（见图 9.9）的修正量直接传递到细网格；第二步，将修正量插值到仅属于细网格的 e、f、$ e^{\prime}, f^{\prime} $、g 点。插值取决于相应点马赫数的正负及绝对值。此处马赫数采用逆变速度计算，例如点 e 采用公式

$$ M_{e}=\frac{\vec{n}\cdot\vec{ v}_{e}}{c}\tag{9.35} $$

式（9.35）的法向量 $ \vec{n} $ 可以由控制体在 A—B 方向的面矢量平均得到，也可将 A 到 B 的矢量归一化得到。随后按以下规则将修正量传递到点 e：

$$ (I_{2h}^{h}\delta\vec W_{2h})_e=\begin{cases}(\delta\vec W_{2h})_A&\text{当 }M_e>1,\\\frac12[(\delta\vec W_{2h})_A+(\delta\vec W_{2h})_B]&\text{当 }|M_e|\leq1,\\(\delta\vec W_{2h})_B&\text{当 }M_e<-1.\end{cases}\tag{9.36} $$

点 f 至 $ f' $ 也采用相同过程。迎风处理使网格间的信息交换更符合实际物理过程。向

点 g 插值更加困难。文献 $ ^{[16]} $ 直接对周围 A 至 D 点的值取平均：

$$ (I_{2h}^{h}\delta\vec{W}_{2h})_{g}=\frac{1}{4}\left[(\delta\vec{W}_{2h})_{A}+(\delta\vec{W}_{2h})_{B}+(\delta\vec{W}_{2h})_{C}+(\delta\vec{W}_{2h})_{D}\right].\tag{9.37} $$

不过，采用某种迎风加权插值更为适宜。三维中也可以用相似方式实现迎风延拓。尽管作了简化，在多个测试算例中仍获得了良好结果 $ [16] $。

9.4.4 非结构网格上的实现

与结构网格相比，非结构网格粗网格的构造复杂得多。问题在于如何从没有特定排列顺序的一组单元中构造均匀粗化的网格。此外，粗网格与细网格单元体积之比必须保持在一定范围内，二维约为 4，三维约为 8。一种解决方法是采用 9.4 节开头简述的 AMG 方法 $ [44–53] $，但几何多重网格仍应用更广，因此这里重点讨论这一方法。

粗网格生成主要有三类方法：

• 非嵌套网格方法；

• 拓扑方法；

• 控制体聚合。

上述方法的综述见文献 $ [59, 60] $。

标准限制和延拓算子仅依据几何插值。Leclercq 和 Stoufflet $ [61] $ 提出了迎风传递算子，对含强激波的流动尤其有潜力。其迎风限制和延拓先将残差或修正量变换为特征变量，作偏迎风插值后，再将限制或延拓所得值变换回物理变量。文献 $ [16] $ 给出了计算成本低得多的迎风多重网格方法，亦见 9.4.3 节及式 $ (9.36) $。

非嵌套网格

最直接的思路是生成一系列完全独立、逐渐变粗的网格 $ [62–67] $。网格之间不必有任何共同节点，因此称为非嵌套网格。但必须在所有粗网格上保留前缘、尾缘、机身头部等主要几何特征；对于复杂几何构型，这并不容易。因此，如今很少采用基于非嵌套网格的多重网格。

拓扑方法

一种方法采用基于图的算法，从细网格删除某些节点，再对剩余节点重新三角剖分 $ [68, 69] $。与非嵌套网格相比，相邻层网格具有共同节点，因此网格间插值更简单；但它仍存在非嵌套网格在几何吻合方面的缺点。

另一种拓扑方法采用网格加密 $ [29, 70, 71] $，从粗网格出发，通过单元分割生成细网格。既可对整个物理区域应用，也可仅作局部加密，例如边界层。这种方法的缺点是最细网格质量高度依赖于初始粗网格及加密过程，边交换可以部分改善这一问题 $ [72, 73] $。

生成粗网格的另一思路基于边折叠 $ [74] $，最初用于四面体网格上的无黏流动；文献 $ [75] $ 又将其扩展到混合单元网格的黏性流动。

聚合多重网格方法

所谓聚合多重网格是一种高效的非结构网格方法，最先由 Lallemand $ [76] $、Lallemand 等人 $ [77] $ 及 Koobus 等人 $ [78] $ 提出，随后被多位研究者采用 $ [31, 79–87] $。该方法通过将细网格控制体与邻居合并生成粗网格，得到由尺寸逐渐增大、形状不规则的多面体单元组成的网格，如图 9.10 所示。聚合技术保留边界表面的完整离散，相比前述非结构多重网格方法，这是一个显著优势。不过，目前聚合多重网格的实现主要基于中位对偶单元顶点格式（5.2.2 节）；其在单元中心格式（5.2.1 节）中的应用见文献 $ [79] $。

通过体积聚合生成粗网格

节点中心格式的体积聚合按以下步骤进行：

1. 构造所谓种子点列表。种子点是选定用来聚合周围控制体的网格点。列表可以包含构成近似极大独立集的点 $ [81] $，也可以直接包含当前网格层的所有点。

2. 遍历全部种子点。

3. 若种子点尚未聚合，将与其以边相连、且尚未聚合的全部最近邻点聚合到一起。

4. 检查粗化比，即一个粗网格控制体中包含的细网格控制体数。若小于 4（三维时为 8），则加入已经聚合的最近邻点的邻居，前提是这些点未归属于其他种子点，直到达到最佳粗化比。文献 [27] 建议优先聚合与至少两个（三维时为三个）已聚合最近邻点相连的二阶邻点。

5. 若列表中仍有种子点，返回步骤 2。

6. 消除孤立单元。这些控制体由于没有未聚合邻居而未能聚合。可以将孤立单元并入粗化比最小的相邻控制体，从而使粗网格控制体尺寸分布更加规则 [31]。

重复上述过程，直到生成全部粗网格。

为保持网格各向同性，体积聚合必须从边界开始。三维时可能需要用户根据边界形状指定聚合方向。为克服这一困难，Okamoto 等人 $ [83] $ 提出了称为全局粗化的另一算法，采用基于边着色的全局划分方法，生成边的独立集。第二步，将共享独立集中某条边的全部控制体聚合，反复执行，直到达到指定粗化比。该方法无需指定初始种子点或聚合方向，也可处理任意单元类型。Moulitsas 和 Karypis $ [88, 89] $ 还提出另一方法，采用多层次思想优化给定的粗网格质量指标。

聚合多重网格的问题

对无黏流动，聚合多重网格的实现较简单。Euler 方程通常离散为各控制体面上的通量，因此控制体形状再复杂也无妨，实际使用的是平均面矢量。此外，一阶空间格式仅需相邻控制体中的流动量，这些信息在粗网格上易于获得。

对于黏性流动，在任意形状控制体上离散扩散通量不再简单，问题在于面的中点处梯度的计算（见文献 $ [31] $ 的讨论）。另一个更加

严重的困难是延拓算子的精度要求。式（9.25）提示需要 $ m_{P} = 2 $，因为标准限制算子，即残差求和，仅能达到 $ m_{R} = 1 $。但在粗网格上构造线性插值并不容易。

为避开一阶精度延拓算子的问题，Mavriplis $ [59] $ 提出采用常值延拓，即同一聚合粗网格控制体内的全部细网格点得到相同解修正，并对黏性通量额外缩放，系数为 $ 0.5^{\text{(level}-1)} $，梯度也同样缩放。不过，这种方法不能获得最佳多重网格效率。

Haselbacher $ [31] $ 建议在粗网格上保留细网格的黏性通量离散，并像最细网格一样施加边界条件，还提出了分片线性延拓算子。对于标量 U，可写为

$$ I_{2h}^{h}\delta U_{2h}=(\delta U_{2h})_{i}+\Psi_{i}(\nabla\delta U_{2h})_{i}\cdot\vec{r}_{i,j}.\tag{9.38} $$

式（9.38）的梯度 $ (\nabla\delta U_{2h})_{i} $ 采用 5.3.4 节式（5.55）的线性最小二乘重构计算。限制器函数 $ \Psi_{i} $ 的值按 5.3.5 节介绍的 Barth–Jespersen 限制器计算。