各种研究文献和学术论文中无一例外地都使用了网格独立性验证（网格无关性验证）。但要定量地评价网格独立性则并不容易。本文内容来自NPARC联盟计算流体力学（CFD）验证与确认网站，虽然年份有些久，不过仍然有较好的指导作用。

1 引言

检验模拟的空间收敛性，是一种确定CFD模拟中 离散误差阶次 （ordered discretization error） 的方法。该方法要求在两个或更多逐级加密的网格上分别开展模拟。术语 网格收敛性研究 （grid convergence study） 等同于常用术语 网格细化研究 （grid refinement study）。

随着网格不断加密（网格单元尺寸减小、流场计算域内的单元数量增加）以及时间步长不断缩小，空间离散误差与时间离散误差应分别趋于零（仅排除计算机舍入误差）。关于CFD模拟空间与时间收敛性的检验方法，详见 Roache 所著书籍，其核心基于理查德森外推法（Richardson’s extrapolation）。此处简要概述该方法。

工程人员通常希望确定由最细网格解所得工程量的误差带。然而，若CFD模拟属于一项设计研究，可能需开展数十次乃至数百次模拟，则宜选用较粗网格以节省资源。因此，也有必要评估较粗网格解的误差。

2 网格收敛性研究中的网格设置考虑

生成网格最简便的方法是：先构建一个初始细网格，其网格间距已接近用户可接受的极限——无论是受限于网格生成阶段的计算资源，还是该网格下求解收敛所需的时间成本。随后，沿各坐标方向每隔一条网格线删除一条，即可得到一级更粗的网格；重复此操作，便可生成多级逐级增粗的网格。在构建初始细网格时，可通过使各坐标方向上的网格点数满足特定关系式，预先嵌入  级更粗网格。

其中  为整数。例如，若需构建两级更粗网格（即细、中、粗三级网格），则各坐标方向上的网格点数应为  形式。不同坐标方向上的  值可取不同数值。

无需将各坐标方向的网格点数量减半，即可实现网格粗化。可采用非整数倍的网格加密或粗化策略。此类方法在某些情况下尤为必要，因网格减半可能导致数值解脱离渐近收敛区域。实施非整数倍加密或粗化时，需生成一套新网格，并保持与原始网格一致的网格生成参数。一种可行方案是选取若干关键位置的网格间距作为参考间距：其中一项应为垂直于壁面方向的网格间距，其余可选自流场边界、几何结构连接处或区域交界面等典型位置。选定加密或粗化比例后，该比例须统一应用于所有参考间距。随后，依据网格质量指标（如网格间距比限值）调整各方向上的网格点数量。表面网格与体网格均沿用原始网格的生成方法。为使离散误差能与其他误差源（如迭代收敛误差、计算机舍入误差等）能够清晰区分，网格加密比应满足 。

3 网格收敛阶数

网格收敛阶数反映了求解误差的行为特征，该误差定义为离散解与精确解之间的差值。

其中  为常数，  表征某种网格尺度，  为收敛阶数，如二阶解对应  。

CFD程序所采用的数值算法在理论上具有特定的收敛阶数；然而，边界条件、数值模型及网格质量等因素均会导致收敛阶数下降，因此实际观测到的收敛阶数通常低于理论值。

忽略式（1）中的高阶项，并对方程两端取对数，可以得到：

收敛阶数  可通过  关于  的曲线斜率确定。若存在足够多的数据点，该斜率既可通过图示直接读取，也可通过对数据点进行最小二乘拟合计算得出；但当数据点数量较少时，最小二乘拟合结果往往不够准确。

更为直接的计算  值的方法为基于三组采用恒定网格加密比  所得的解：

求解精度阶数由截断误差主项的阶数确定，其以离散尺度  为基准进行表达。局部精度阶数反映网格某一点处用于离散控制方程的差分模板所具有的精度阶数；全局精度阶数则综合考虑误差在模板作用域外的传播与累积效应。此类传播通常使全局精度阶数较局部精度阶数降低一阶。边界条件的精度阶数可比内部区域低一阶，而不影响整体全局精度。

4 收敛渐近区域

评估代码及计算结果的精度，要求网格足够精细，使解处于收敛的渐近区域。当网格尺度  满足不同  值与对应误差  均使常数  保持不变时，即达到收敛渐近区：

关于渐近区的另一项判据将在“网格收敛指数”一节中讨论。

5 理查森外推法

理查森外推法是一种利用一系列低阶离散解来估算连续介质解（即网格尺度为零时的解）的高阶估计值的方法。

一次数值模拟所得物理量  可用如下级数展开式一般性地表示：

其中，  为网格间距，函数  、  和  与网格间距无关。若  ，则量  被视为二阶精度。  表示网格间距趋近于零时的连续介质解。

若假定解具有二阶精度，并已在网格间距分别为  和  的两套网格上计算出  （其中  为更细密（即更小）的网格间距），则可对上述展开式写出两个方程，并求解  ，从而估计连续介质解：

其中网格加密比定义为：

理查森外推法可推广至 p-th 阶方法及任意  值的网格比（该比值不必为整数），其形式为：

传统上，理查森外推法采用网格加密比  ，此时上述公式简化为：

理论上，若  和  恰好由严格二阶方法计算得到，则上述理查森外推公式将给出  的四阶精度估计；否则，其精度为三阶。一般而言，将  视为  阶精度。理查森外推法既可用于每个网格点上的解，也可用于解的泛函量（如压力恢复系数或阻力系数）。这要求解不仅在局部具有二阶精度，且在整个计算域内亦为全局二阶精度；同时，相关解泛函也需采用一致的二阶方法计算。关于理查森外推法使用的其他注意事项（如非守恒性、舍入误差放大等）详见 Roache 所著书籍。

本研究中，设  为某一解泛函（如压力恢复系数）。此时， 表示网格间距  趋于零极限下  的估计值。 可用于替代 CFD 研究中原始估计值 ，作为更优的参考值予以报告；但须充分认识并说明其适用前提与相关限制条件。

 的另一用途是估算 CFD 计算所得  所含的离散化误差范围。以下将重点考察该用途。

 与  的差值可作为误差估计的一种参考；但该估计需结合  所附带的前述限制条件审慎评估。

研究聚焦于利用  与  构建误差估计。考察前述广义理查森外推公式，其右侧第二项可作为  的误差估计量。该公式可重写为：

其中，  为实际相对误差，定义为：

而  为  的估计相对误差，定义为：

其中相对误差  定义为：

该量不宜直接用作误差估计器，因其未考虑网格加密比  和精度阶数 ，可能导致误差被低估或高估。例如，仅通过选取接近 1.0 的网格加密比 ，即可人为地使其趋近于零。

估计相对误差  是一种误差阶次估计器；当  与  具有足够高的计算精度（即  ）时，  可较好地近似细密网格上的离散化误差。若  与  相对于  为零或极小，则  的数值可能失去意义；此时应选用其他归一化量替代  。

若需执行大量CFD计算（例如开展试验设计（DOE）研究），可考虑采用较粗网格  。此时需估算该较粗网格上的误差。理查德森外推法可表示为：

对应  的估计相对误差定义为：

6 网格收敛性指数（GCI）

Roache 提出网格收敛指数(Grid Convergence Index,GCI)  ，旨在为网格收敛性研究结果的报告提供统一方式，并可能为解的网格收敛性提供误差带。该 GCI 可利用两级网格计算得出；但为准确估算收敛阶数并验证解是否处于渐近收敛范围内，推荐采用三级网格。

严谨的数值分析应能保证：当网格分辨率趋近于零时，所得结果趋近于真实解。此时，离散化方程的解将趋近于原始控制方程的真实解。数值计算中一个关键问题是确定合适的网格分辨率水平，该水平取决于流动条件、分析类型、几何构型及其他变量。实践中常需先选定初始网格分辨率，再通过一系列网格加密操作评估网格分辨率对结果的影响，此类分析称为网格加密研究。

使用者需明确区分两类情形：一类是数值结果趋近于某个渐近数值解；另一类是数值结果趋近于控制方程的真实物理解。理想情况下，随着网格不断加密、分辨率持续提高，计算所得解的变化幅度应逐渐减小，并趋近于某一渐近值（即真实的数值解）。然而，该渐近值与控制方程的真实物理解之间仍可能存在误差。

Roache 提出了一套用于统一报告网格加密研究结果的方法论。其基本思想为：将任意网格加密试验的结果近似关联至采用二阶方法进行网格加倍时所预期的结果。该 GCI 基于广义理查德森外推理论导出的网格加密误差估计器。无论是否实际采用理查德森外推法提升精度，均推荐使用该指标；在某些情况下，即使理论前提未被严格满足，亦可采用。该指标旨在量化网格收敛性的不确定性。

GCI 表征计算值偏离渐近数值解的百分比，反映解与渐近值之间的误差带范围，同时指示若进一步加密网格，解可能发生的变化量。较小的 GCI 值表明当前计算已处于渐近收敛范围内。

细网格上的 GCI 定义为：

式中  为安全系数。网格加密可在空间维度或时间维度上实施。对于两级网格对比，推荐取  ；对于三级及以上网格对比，推荐取  。在正式报告中建议采用较高的安全系数，该取值对实际误差具有较强保守性。

当某项设计或分析工作需开展大量CFD模拟（例如DOE研究）时，可考虑采用较粗网格  。此时需量化该较粗网格上的误差。较粗网格的GCI定义为：

每个网格层级所得解均须位于计算解的收敛渐近范围内，此为关键前提。该条件可通过分析三套网格计算所得的两个网格收敛指数（GCI）值予以验证：

6.1 所需网格分辨率

若已知目标精度水平，并已有网格分辨率相关研究结果，即可据此估算实现目标精度所需的网格分辨率。

6.2 独立坐标方向加密与混合阶方法

网格加密比  对定常解在所有空间坐标方向  上具有同等适用性；对非定常解，则在时间方向  上同样适用。若该假设不成立，可分别计算各坐标方向的网格收敛指数，再将其叠加，获得整体网格收敛指数。

6.3 有效网格加密比

若生成了更细或更粗的网格，但尚不确定应采用何种网格加密比，可计算一个有效网格加密比：

其中  表示该网格所用总节点数，  表示流场域的维数。该有效网格加密比同样适用于非结构化网格。

6.4 网格收敛性研究示例

以下示例展示了上述流程在网格收敛性研究中的具体应用。该CFD分析旨在确定某进气道的压升恢复系数（pressure recovery）。流场计算在三套网格上完成：后一套网格在  和  坐标方向上的节点数均为前一套的两倍，而  方向节点数保持不变。由于流动在  方向上呈轴对称，较细网格可视为其相邻较粗网格在该方向上的两倍加密。下表列出了各套网格的基本参数及对应解所得的压升恢复系数。所有计算结果均已达到充分迭代收敛。间距栏表示经最细网格间距归一化后的网格间距。

下图展示了压升恢复系数随网格间距变化的趋势。随着网格间距减小，压升恢复系数逐渐趋近于零网格间距对应的渐近值。

由上述结果可确定实际观测到的收敛阶数：

理论收敛阶数为  。二者差异主要源于网格拉伸、网格质量、解的非线性特性、激波、湍流模拟以及其他可能因素。

现可利用两套最细网格实施Richardson外推法，估算零网格间距下的压升恢复系数：

随后可计算细网格解对应的网格收敛指数。由于采用三套网格估计收敛阶数  ，故取安全系数  。网格1与网格2之间的GCI为：

网格2与网格3之间的GCI为：

此时可验证，所得解处于渐近收敛范围内：

该值约等于1，表明所得解完全处于渐近收敛范围内。

依据本项研究，压力恢复系数估计为  ，误差带为  或  。

1. Roache, P.J., Verification and Validation in Computational Science and Engineering , Hermosa Publishers, Albuquerque, New Mexico, 1998.
2. https://www.grc.nasa.gov/www/wind/valid/tutorial/spatconv.html

”