7.7 多重网格技术
前面章节已经说明，离散误差会随着网格间距的减小而降低。换言之，网格越细，CFD 模拟的精度通常越高。与直接法相比，迭代法所需的存储开销较小，因此更适合求解由高度细化网格产生的大规模代数方程组。
此外，第 6 章已经指出，用于连续方程和动量方程耦合的 SIMPLE 算法本身也是迭代过程。因此，只要整个迭代过程最终能够收敛到真实解，中间解并不需要非常精确。遗憾的是，随着网格不断细化，Jacobi 方法和 Gauss-Seidel 方法等迭代方法的收敛速度会迅速降低。
为了考察迭代方法的收敛速度与问题中网格单元数量之间的关系，这里考虑一个简单的二维方腔驱动流问题。图 7.5 中的小插图表明，计算区域为一个尺寸为 的方形空腔。方腔顶盖以 的速度沿 正方向运动。腔内流体为空气，并假定流动为层流。采用逐线迭代求解器，分别在 、 和 三种网格上计算该问题。
为了度量迭代序列中某个中间解与真实解的接近程度，对第 个方程采用式（7.26）定义的残差。对系统中全部 个方程取平均，即对流动问题计算区域内全部控制体取平均，可以得到平均残差 。对于给定问题，平均残差是判断迭代收敛性的一个有用指标：
如果迭代过程收敛，则平均残差 应趋于零，因为当 时，所有残差分量均满足 。对于某个给定求解变量，例如 速度分量，平均残差通常需要归一化。这样可以更容易地解释不同算例中的残差大小，也便于将其与其他求解变量的残差进行比较，例如 速度、 速度或压力的残差。这些变量本身的量级可能相差很大。
最常用的归一化方式，是将第 次迭代后的平均残差与第一次迭代后的平均残差之比作为归一化残差：
*图 7.5 使用不同网格分辨率的逐线迭代求解器所得残差下降规律。
图 7.5 给出了 动量方程归一化残差随迭代次数的变化。当所有求解变量的归一化残差都低于 时，本算例停止迭代。这里的求解变量包括速度和压力。可以看到， 网格的解在 161 次迭代后收敛，而 和 网格分别需要 331 次和 891 次迭代才能收敛。
在 CFD 程序中，可以通过调整求解参数来提高收敛速度，其中包括松弛因子。为了保持比较的一致性，这里所有求解参数均保持不变。从图中可以清楚看到残差下降的规律：残差在初期迅速降低，随后下降速率逐渐变小。图中还表明，最细网格的最终收敛速度最低。如果采用更细的网格，收敛所需时间还会进一步增加。
多重网格的基本概念
为了简化多重网格方法的说明，下面采用矩阵记号，并先回顾残差的定义。考虑由流动区域内某个守恒方程经有限体积离散后得到的线性方程组：
向量 是方程组（7.32）的真实解。
如果采用迭代方法求解该方程组，在经过若干次迭代后可得到一个中间解 。该中间解并不精确满足式（7.32），因此可按前述方式定义残差向量 ：
（7.33）
还可以将误差向量 定义为真实解与中间解之差：
将式（7.32）减去式（7.33），可得到误差向量和残差向量之间的关系：
（7.35）
在迭代过程的任意阶段，只要将中间解代入式（7.33），就可以方便地计算残差向量。可以设想用一个迭代过程求解式（7.35），从而得到误差向量。为此，将该系统写成迭代矩阵形式是有用的：
由于方程组（7.32）和（7.35）的系数矩阵 相同，迭代矩阵中的系数 与所选迭代方法的系数相同。这些方法可以是 Jacobi 方法，也可以是不带松弛或带松弛的 Gauss-Seidel 方法。不过，常向量的元素不同：
实际计算中，如果用求解原始方程组（7.32）时相同的迭代方法去求解误差方程组（7.35），并不会在收敛速度上带来差别。然而，式（7.35）十分重要，因为它揭示了误差如何从一次迭代传播到下一次迭代。其等价形式（7.36）进一步突出了迭代矩阵的关键作用。正如前面介绍松弛技术时所说明的，迭代矩阵的性质决定了误差传播速度，也决定了收敛速度。
这些性质以及误差传播随迭代技术、网格尺寸和离散格式等因素变化的数学行为，已经得到广泛研究。已有结论表明，解的误差包含一系列不同波长的分量，这些波长是网格尺寸的若干倍。迭代方法能够快速削弱短波长误差分量，即波长不超过网格尺寸若干倍的误差分量。然而，长波长误差分量随着迭代次数增加通常衰减得很慢。
这种误差行为解释了图 7.5 中观察到的趋势。对于粗网格，最长可能的误差波长，即量级与计算区域尺寸相当的误差分量，仍处在该网格的短波长范围内，因此所有误差分量都能快速降低。对于更细的网格，最长误差波长则逐渐远离能够快速衰减的短波长范围。
多重网格方法正是利用误差行为的这种内在差异，在不同尺度的网格上进行迭代。用细网格处理局部误差，用粗网格处理整体误差。粗网格计算量小，多做几步迭代的代价不高，却能明显加快收敛。
7.7.1 多重网格过程概述
下面简要说明两阶段多重网格过程的基本原理。
步骤 1：细网格迭代。 在网格间距为 的最细网格上进行迭代，得到方程组 的中间解 ，其中真实解向量为 。迭代次数应足够多，使误差中的短波长振荡分量得到有效降低，但不试图消除长波长误差分量。
该网格上解的残差向量 满足：
误差向量 可写为：
前面已经得到误差和残差之间的关系：
步骤 2：限制。 将解从网格间距为 的细网格传递到网格间距为 的粗网格上，其中 。由于粗网格的网格间距较大，原细网格上的长波长误差在新网格上表现为短波长误差，因此会快速降低。
如果采用网格间距为细网格两倍的粗网格，传递过程会更简单。这里不直接求解粗网格上的解向量 ，而是在粗网格上求解误差方程：
计算从初始猜测 开始。为了进行求解，需要残差向量和系数矩阵的值。已知细网格上的 后，必须采用适当的平均过程求得粗网格上的残差向量 。矩阵 的系数可以在粗网格上重新计算，也可以通过对细网格系数矩阵 进行某种平均或插值来得到。由于粗网格上的单次迭代计算成本较小，因此可以进行足够多的迭代，以获得误差向量 的收敛解。
步骤 3：延拓。 在得到粗网格误差向量 的收敛解后，需要将其传回细网格。需要注意的是，粗网格上的数据点少于细网格上的点数。因此，需要采用合适的插值算子，例如线性插值，在细网格中的中间点生成延拓后的误差向量 。
步骤 4：修正和最终迭代。 计算出延拓误差向量 后，可以修正细网格上的中间解：
由于长波长误差已经被消除，改进后的解更接近真实解向量 。不过，限制和延拓过程中引入了若干近似，因此仍需以改进解为初值再进行少量迭代，以消除限制和延拓过程中可能引入的误差。
上述说明针对的是两阶段过程，即一个细网格和一个粗网格。实际应用中，限制过程通常会传递到多个逐渐变粗的网格层级。随后，在返回初始网格的过程中，也要在每个层级执行延拓过程。
7.7.2 示例说明
例 7.5
考虑求解一根绝热金属杆的一维导热方程，该金属杆内部存在热生成。控制方程为：
该问题的尺寸和其他数据如下：杆长为 ，横截面积为 ，热导率为 ，热生成率为 ，两端温度分别为 和 。假设希望采用 20 个单元的网格求解该问题，此时网格间距为 ，并将其命名为网格 1。
图 7.6 用于求解该问题的 20 节点网格——网格 1。
图 7.6 给出了网格 1 以及两端标注的边界条件。这里没有必要说明该问题的离散方程如何得到，其过程类似于例 4.2。表 7.12 汇总了节点 1、2、3、……、20 处离散方程的系数。
表 7.12 各节点处离散方程的系数。
利用表 7.12 中的表达式，可得到表 7.13 中的系数数值，并构造矩阵方程 ，其中解向量 包含网格 1 各节点处的温度。
表 7.13 离散方程系数的数值。
矩阵方程为：
系数矩阵为三对角矩阵，因此可以采用 TDMA 一次扫描得到解。为了后续验证多重网格解，表 7.14 给出了该解。
表 7.14 TDMA 解：网格 1 各节点温度。
步骤 1：细网格迭代
采用 Gauss-Seidel 迭代（7.20）求解这些方程。为了启动迭代过程，简单地将所有位置的温度初始化为 。如果初始场更接近最终解，就不容易体现该方法的优势。经过 5 次 Gauss-Seidel 迭代后，解向量 为：
此时残差向量 为：
总均方根残差为 14.951。如果继续迭代，残差向量将缓慢降低，直到满足收敛判据。图 7.9 在本节末尾给出了 Gauss-Seidel 迭代的收敛模式。当采用均方根残差和小于 作为收敛判据时，最终解在 664 次迭代后获得。该收敛解当然与表 7.14 中的 TDMA 解一致。
步骤 2：限制
为了应用多重网格方法，首先必须构造粗网格。最简单的方法是构造一个单元数为原来一半的网格。图 7.7 将细网格和拟采用的粗网格从上到下画出。第一个粗网格采用 10 个单元，网格间距为 ，命名为网格 2。下一个粗网格为网格 3，由 5 个单元组成，网格间距为 。
图 7.7 用于求解该问题的网格。
如果细网格间距为 ，则采用一半单元数的网格的网格间距为 。在多重网格文献中，通常用上标表示网格间距。按这种记号，细网格上的残差向量为 。
现在需要将残差向量从细网格插值到粗网格。由于网格 2 的节点恰好位于网格 1 节点之间的中点，可以通过对 进行简单平均，得到粗网格上的残差向量 。这些值汇总于表 7.15。表中仅显示实际数值的 3 位小数。如前所述，该传递过程称为“限制”。
表 7.15 从细网格 1 传递到粗网格 2 的限制过程：细网格和粗网格残差。
限制后的网格 2 残差向量以矩阵形式写为：
注意，现在只有 10 个值。粗网格上的误差满足方程 。已经计算出向量 ，但还需要矩阵 ，才能求解该方程并得到 。在多重网格文献中，有许多方法采用精巧的插值算子来估计 。在本例中，不对系数矩阵进行插值，而是使用表 7.12 中的表达式精确计算粗网格矩阵的系数。因此，误差向量 的矩阵方程为：
现在以 为初始猜测，用 Gauss-Seidel 迭代过程求解方程组（7.39）。由于迭代现在发生在粗网格上，残差降低速度更快，而且单次迭代的计算成本低得多。经过 10 次粗网格迭代后，在第一个粗网格即网格 2 上得到误差向量 ：
由于只进行了 10 次迭代，该解只是部分收敛的，因此仍会存在残差：
其数值与限制后得到的网格 3 残差 一并列于表 7.16。
表 7.16 网格 2 上的残差和限制到网格 3 后的残差。
现在，残差 已传递到只有 5 个节点的更粗网格上，得到残差 。然后在网格 3 上求解误差 ，对应的方程组为：
其中， 的系数仍采用表 7.12 中的表达式重新计算。由于单次迭代的成本很低，可以在网格 3 上进行更多迭代，以非常有效地降低误差。经过 10 次迭代后，表 7.17 给出了误差向量 的解。
表 7.17 网格 3 上的解，即网格 3 上的误差向量 。
该粗化过程可以继续进行，但在这个说明性示例中，限制过程在 5 节点网格处停止。
步骤 3：延拓
下一步是反向进行，将误差向量从每一个粗网格传递到上一级细网格。这一过程称为延拓。可以采用线性插值或其他插值格式，根据粗网格值构造细网格值。采用线性插值时，部分示例值为：
依此类推。需要注意的是，撇号用于表示网格 2 上延拓得到的误差向量 ，以便将其与误差向量 区分开。此外，对于最靠近边界的节点，由于问题变量在边界处的值已知，因此边界上的误差为零。由式（7.41）计算得到延拓误差向量 为：
现在，用延拓误差向量 修正式（7.40）给出的网格 2 原误差向量 ：
于是得到：
在将该误差向量传递到上一级网格之前，通常要进行若干次平滑迭代。这里先进行 2 次 Gauss-Seidel 平滑迭代，得到网格 2 上经过修正和平滑后的误差向量：
随后，用这一结果通过线性插值计算网格 1 上的延拓误差向量 。具体过程与式（7.41）相同，只是将上标 替换为 ，将 替换为 ：
步骤 4：修正和最终迭代
最后，用延拓误差向量 计算网格 1 上的修正中间解 ：
因此：
将修正解（7.43b）与中间解（7.38）以及表 7.14 中的 TDMA 解比较，可以看出多重网格过程显著降低了误差。将修正解代入 ，得到均方根残差为 8.786，低于网格 1 上先前的均方根残差 14.951。
由于限制和延拓过程中包含插值误差，不能期望通过一个多重网格循环就得到真实解。为了进一步改进解，需要在细网格上再进行若干次迭代，例如 2 次，然后重复“细网格—粗网格”过程，直到收敛。这里采用三网格过程，并根据需要多次往返，以将均方根残差降至 。该多重网格循环称为三网格 V 循环。图 7.8 给出了这些过程步骤，圆圈中的数字表示各层级上的迭代次数。该图也说明了“V 循环”这一名称的来源。
图 7.8 本例采用的 V 循环多重网格过程示意图。
重复执行由 2 次细网格迭代、10 次粗网格迭代和 10 次更粗网格迭代组成的 V 循环后，得到的收敛模式如图 7.9 所示。多重网格过程快速且有效：它在 60 次细网格迭代后收敛。相比之下，普通 Gauss-Seidel 方法需要 664 次迭代才能达到相同的残差水平。
即使考虑粗网格迭代带来的额外计算量，多重网格过程所带来的数量级收敛速度提升仍然明显有利。当多重网格加速技术应用于二维和三维问题时，得到的收敛收益非常有吸引力，这也是该方法在 CFD 用户中广泛使用的原因。
图 7.9 普通 Gauss-Seidel 迭代和多重网格 Gauss-Seidel 迭代的残差下降规律。
7.7.3 多重网格循环
多重网格技术可以与任意迭代技术结合使用。在上面的简单示例中，已经说明了多重网格方法的主要概念。在实际 CFD 计算中，多重网格传递过程更加复杂，并会采用不同的粗化和细化循环，同时在不同加密层级上配合特定的限制和延拓调度方式。常见的多重网格循环包括所谓的 V 循环、W 循环和 F 循环，如图 7.10 所示。
图 7.10 不同多重网格循环策略示意图：（a）V 循环，（b）W 循环，（c）F 循环。
图 7.10a 所示的简单 V 循环由两条分支组成。计算从最细层级开始。在任意层级上的迭代都称为松弛迭代。在最细层级上进行若干次松弛扫描后，将残差限制到下一个粗层级；在该层级上松弛后，再将残差传递到下一个粗层级。该过程一直持续到最粗层级为止。在最粗层级上完成最后的松弛后，在 V 循环的上行分支执行延拓步骤，直到回到最细层级。
在 W 循环中，为了更好地降低长波长误差，会在粗层级上增加额外的限制和延拓扫描。典型模式如图 7.10b 所示。柔性循环或 F 循环与 W 循环非常相似，但其粗层级扫描模式不同，如图 7.10c 所示。
在称为完全多重网格（FMG）的方法中，计算并不从最细网格开始，而是从最粗层级开始。解依次传递到更细的网格层级。在最细层级上，延拓后的解被用作迭代过程的初始猜测。该求解过程还可以利用任意循环过程进一步加速。例如，可以在每个逐级细化层级上使用 V 循环，如图 7.11 所示。
图 7.11 完全多重网格方法中采用的循环策略。
7.7.4 多重网格方法中的网格生成
如上例所示，需要生成网格来构造粗网格。最直接的方法是合并控制体，或者在上一层级网格的节点数基础上减半并重新生成网格。对于二维结构化网格，例如图 7.12 所示的 Cartesian 网格，可以通过删除交替网格线方便地生成粗网格。这样，每 4 个细网格控制体可以构造出 1 个粗网格控制体。该方法也可以很容易地推广到三维网格，即每 8 个细网格控制体对应 1 个粗网格控制体。
图 7.12 二维 Cartesian 网格——通过删除交替网格线或合并四个控制体构造粗网格。
在上面的示例中，粗网格系统矩阵和其他所需量是利用粗网格的实际几何性质计算得到的，见表 7.12。这类多重网格过程称为几何多重网格过程。该方法的另一种变体中，为了节省计算量，不再根据网格几何重新计算系数，而是将粗网格系数近似表示为细网格方程系数的线性组合。这类多重网格方法称为代数多重网格，并广泛用于商用 CFD 求解器。Hutchinson 和 Raithby（1986）提出的加性修正多重网格（ACM）策略，也是许多 CFD 过程中常用的一种多重网格方法。
