CFD的网格独立性验证与结构有限元分析的网格收敛性验证是相同的吗？

基本原理（理查森外推法/GCI）是相同的，但CFD有一个附加条件，即y+值需与所使用的湍流模型相匹配。此外，由于非线性较强，其收敛阶数通常低于理论值，这是两者的主要区别。

什么是GCI（网格收敛指数）？

GCI是一种基于理查森外推法的离散化误差估计指标，由Roache提出，后由Celik等人在2008年针对CFD应用将其标准化为五步流程。通常认为，当精细网格的GCI_fine值小于百分之几时，即可判定为网格独立。

进行网格独立性验证需要多少层网格？

至少需要三层网格（粗、中、细）。如果只有两层网格，则无法估计收敛阶数p，只能进行简单的变化率比较。Celik标准流程要求使用三层网格来计算观测收敛阶数，进而得出GCI值。

y+值与网格独立性验证有何关系？

y+是壁面第一层网格的无量纲距离。如果该值不在所用湍流模型要求的范围内，即使加密网格也无法收敛到正确解。使用壁函数时，通常要求y+在30至300之间；进行壁面解析时，则要求y+小于1。

CFD网格独立性验证（网格独立性研究）

分类: 流体解析（CFD） | 综合版 2026-04-12

CFD mesh independence study convergence plot showing GCI with grid refinement

メッシュ独立性確認：格子細分化比と解の収束挙動

理论与物理

概述

🧑‍🎓

老师，CFD的网格无关性验证，和结构FEM中做的网格收敛验证基本上是一回事吗？

🎓

原理是一样的。两者都是确认“即使网格进一步细化，解也几乎不变”的工作。用Richardson外推法估计网格尺寸$h \to 0$时的解，并用GCI（网格收敛指数）量化网格误差，这些也是共通的。

🧑‍🎓

那如果结构分析中做过的话，CFD中也能直接套用吗？

🎓

有几个不同的注意事项。CFD中有三个附加条件。第一，y+是否与湍流模型匹配。如果这个不匹配，无论网格细化到什么程度，都不会收敛到正确的解。第二，由于CFD中纳维-斯托克斯方程的非线性很强，观测到的收敛阶数$p$常常低于理论值（二阶精度格式则$p=2$）。第三，如果迭代收敛的残差没有充分降低，离散化误差和迭代误差会混在一起，使得GCI失去意义。

🧑‍🎓

比结构分析更容易踩坑啊……有系统性的步骤吗？

🎓

有的。Celik等人（2008）的论文 "Procedure for Estimation and Reporting of Uncertainty due to Discretization in CFD Applications"（Journal of Fluids Engineering）是CFD中应用GCI的事实上的标准步骤。也成为了ASME期刊投稿的评审标准。掌握了这个5步流程，在实务和论文中都能使用。

Richardson外推法

🧑‍🎓

Richardson外推法，大致是什么思路呢？

🎓

将离散解$f(h)$进行泰勒展开，与真解$f_{\text{exact}}$的关系如下：

$$ f(h) = f_{\text{exact}} + C \cdot h^p + O(h^{p+1}) $$

🎓

$h$是代表单元尺寸，$p$是离散化格式的精度阶数，$C$是未知常数。例如二阶精度的中心差分理论上$p=2$，所以网格减半时误差大约变为$1/4$。利用这个关系式，使用两个或更多不同网格尺寸的解来估计$f_{\text{exact}}$，就是Richardson外推法的基本思路。

🧑‍🎓

但是$p$和$C$两个都是未知数的话，就需要两个方程对吧。所以需要三个网格等级，是这个原因吗？

🎓

没错。通过三个网格等级（细$h_1$、中$h_2$、粗$h_3$，且$h_1 < h_2 < h_3$）得到解$f_1, f_2, f_3$，解联立方程组就可以估计观测收敛阶数$p$。如果只有两个等级，就只能假设$p$，从而变成“网格细化后变化很小所以OK”这种主观判断。

GCI的定义与含义

🧑‍🎓

请告诉我GCI的公式。和结构FEM中看到的一样吗？

🎓

公式本身是一样的。细网格与中网格之间的GCI如下：

$$ \text{GCI}_{\text{fine}}^{21} = \frac{F_s \, |e_a^{21}|}{r_{21}^{p} - 1} $$

🎓

这里各变量的含义如下：

$e_a^{21} = (f_1 - f_2)/f_1$：细网格解$f_1$与中网格解$f_2$的相对差
$r_{21} = h_2 / h_1$：网格细化比（$r > 1$）
$p$：观测收敛阶数（由后述公式计算）
$F_s$：安全系数。使用三个等级外推时$F_s = 1.25$，两个等级且假设$p$时$F_s = 3.0$

🧑‍🎓

$F_s = 1.25$和$3.0$差别很大呢。果然还是应该做三个等级吗？

🎓

是的。$F_s = 3$是对“假设了$p$的不确定性”的惩罚。仅仅多做一个等级，安全系数就变成$1/2.4$了，所以没有理由不做三个等级。例如，如果汽车外气动中阻力系数$C_d$的GCI在1.5%以下，就可以说具备了与风洞试验比较的精度。

Celik的5步流程

🧑‍🎓

Celik的流程具体有哪些步骤？

🎓

我来解释一下Celik等人（2008, JFE 130(7)）的5步流程。

🎓

Step 1. 定义代表单元尺寸$h$

3D情况：

$$ h = \left[ \frac{1}{N} \sum_{i=1}^{N} (\Delta V_i) \right]^{1/3} $$

🎓

$N$是总单元数，$\Delta V_i$是各单元的体积。2D则使用面积的平方根。

🎓

Step 2. 用三个网格等级求解

细化比$r = h_{\text{coarse}}/h_{\text{fine}}$理想情况下$r \geq 1.3$。$r = 2$（单元数8倍）最好，但如果计算成本紧张，$r = 1.3$〜$1.5$（单元数2〜3倍）也可以。得到三个解$f_1, f_2, f_3$。

🎓

Step 3. 求观测收敛阶数$p$

$$ p = \frac{1}{\ln(r_{21})} \left| \ln\left|\frac{\varepsilon_{32}}{\varepsilon_{21}}\right| + q(p) \right| $$

$$ q(p) = \ln\left(\frac{r_{21}^{p} - s}{r_{32}^{p} - s}\right), \quad s = \text{sign}\left(\frac{\varepsilon_{32}}{\varepsilon_{21}}\right) $$

🎓

这里$\varepsilon_{21} = f_2 - f_1$，$\varepsilon_{32} = f_3 - f_2$。如果$r_{21}$和$r_{32}$相等（等比细化），则$q=0$，可以使用简化版公式：

$$ p = \frac{\ln|\varepsilon_{32}/\varepsilon_{21}|}{\ln(r)} $$

🎓

Step 4. 求Richardson外推值$f_{\text{ext}}^{21}$

$$ f_{\text{ext}}^{21} = \frac{r_{21}^{p} \, f_1 - f_2}{r_{21}^{p} - 1} $$

🎓

Step 5. 报告误差估计量

相对近似误差：$e_a^{21} = |(f_1 - f_2)/f_1|$
相对外推误差：$e_{\text{ext}}^{21} = |(f_{\text{ext}}^{21} - f_1)/f_{\text{ext}}^{21}|$
$\text{GCI}_{\text{fine}}^{21} = \frac{1.25 \, e_a^{21}}{r_{21}^{p} - 1}$

🧑‍🎓

$p$的公式是隐式的（$p$在等式两边都有），怎么求解呢？

🎓

迭代求解。将简化版的值作为$p$的初始值，更新$q(p)$再重新计算$p$，这样重复几次就会收敛。用Python的话，scipy.optimize.fsolve 可以一步搞定。如果是等比细化（$r_{21} = r_{32}$），则$q=0$无需迭代，可以直接使用简化版。

y+与湍流模型的匹配性

🧑‍🎓

y+和网格无关性验证具体有什么关系，我还不太明白……

🎓

问得好。y+是壁面第一层单元的无量纲距离，定义如下：

$$ y^+ = \frac{u_\tau \, y}{\nu}, \quad u_\tau = \sqrt{\frac{\tau_w}{\rho}} $$

🎓

$y$是壁面到第一层单元中心的距离，$u_\tau$是摩擦速度，$\nu$是运动粘度系数。关键是，每个湍流模型都有其要求的y+范围：

使用壁函数（标准$k$-$\varepsilon$等）：$y^+ = 30\text{〜}300$
壁面解析（$k$-$\omega$ SST、SA、LES）：$y^+ < 1$（理想是$y^+ \approx 1$）
混合壁面处理（SST的automatic wall treatment）：期望$y^+ < 5$

🧑‍🎓

如果是壁函数模式但y+在5左右，网格无关性验证中会发生什么？

🎓

壁函数以对数律（$y^+ > 30$的区域）为前提，$y^+ = 5$正处于缓冲层（buffer layer）中间，是对数律和粘性底层的线性律都不成立的“物理上不匹配”的区域。在那里壁面剪切应力会出现很大误差。网格越细化，y+越低，越偏离壁函数的适用范围，解会表现出发散行为，或者出现非单调收敛。结果导致GCI的前提（单调收敛）崩溃。

🧑‍🎓

也就是说，作为网格无关性验证的“前提条件”，必须首先确认y+在湍流模型的预期范围内，对吧。

🎓

没错。现场常见的情况是“用壁函数计算了，但部分壁面的y+在3左右，没注意到这点就去做网格研究，然后苦恼‘不收敛！’”。这就像没正确夹好体温计就测量，然后说“不对劲”一样。

GCI公式各项的物理含义

相对近似误差 $e_a^{21}$：两个网格间解的变化程度。这个值越小，解对网格的依赖性就越小。但仅$e_a$小还不够，如果收敛阶数$p$不正常，则可能是“碰巧得到了相同值”。
网格细化比 $r_{21}$：粗网格与细网格的代表尺寸之比。$r$越大可靠性越高，但计算成本会增加$r^d$倍（$d$：维数）。3D分析中$r=2$则成本增加8倍。
观测收敛阶数 $p$：如果得到的值接近离散化格式的理论精度阶数，则说明处于“渐近区”（asymptotic range），外推结果可信。如果$p$与理论值偏差很大，则可能怀疑存在编程错误、迭代收敛不充分、奇点影响等。
安全系数 $F_s$：设定为使GCI具有接近“95%置信区间”的含义。Roache推荐$F_s=1.25$，但这是使用三个等级外推求得$p$时的值。

Richardson外推法的适用限制

处于渐近区：网格足够细，高阶误差项$O(h^{p+1})$可以忽略。网格太粗时$p$会不稳定。
单调收敛：$\varepsilon_{32}/\varepsilon_{21} > 0$（向同一方向变化）。非单调收敛或振荡收敛时无法应用标准流程。
相同的离散化格式：三个等级都使用相同的数值格式（对流项离散化、梯度计算、压力-速度耦合方法）。
充分的迭代收敛：每个网格的解已迭代充分收敛（残差在$10^{-5}$〜$10^{-6}$以下为参考）。迭代误差大于离散化误差时，外推就失去意义。
无奇点：如果评估量包含再附着点或角部的压力奇点，收敛将不符合理论预期。

数值解法与实现

代表单元尺寸的定义

🧑‍🎓

老师，代表单元尺寸$h$具体是怎么计算的呢？非结构网格的话，单元尺寸不是参差不齐吗？

🎓

Celik的流程中使用全局平均值。3D情况下使用全部单元的平均体积的立方根：

$$ h = \left[ \frac{1}{N} \sum_{i=1}^{N} \Delta V_i \right]^{1/3} $$

🎓

实务中，“只要知道总单元数$N$，就可以从计算域总体积$V_{\text{total}}$近似求得$h = (V_{\text{total}}/N)^{1/3}$”，所以只要求解器输出中有单元数，马上就能计算。

🧑‍🎓

观察的是局部评估量（如壁面摩擦等），却用全局平均值定义$h$，这样可以吗？

🎓

很敏锐的指正。严格来说，只有在评估量所依赖的区域被均匀细化的前提下，全局$h$才是合适的。例如，如果只细化壁面附近而远处保持粗网格，那么细化比$r$实际上会发生变化。最佳实践是对整个区域进行相似细化（保持几何比例不变，均匀增加所有单元数量）。在Fluent中，可以通过 Mesh > Controls > Sizing 将所有参数一并缩放。

观测收敛阶数的计算

🧑‍🎓

在实际的CFD问题中，$p$大概会是多少呢？二阶精度格式的话就是$p=2$吗？

🎓

理想情况下是这样，但在实务中，$p$的值常常分散在$p = 0.5$〜$3.0$左右的范围内。非线性问题中$p < 2$的情况很多。大致参考如下：

$p \approx 1.5$〜$2.5$：正常。可以判断处于二阶精度的渐近区。
$p < 0.5$：网格仍然太粗，尚未进入渐近区，或者存在错误。
$p > 3.0$：可能发生了偶然抵消。应再增加一个等级进行确认。
$p < 0$（$\varepsilon_{32}/\varepsilon_{21} < 0$）：非单调收敛。无法应用标准流程。