CFD网格独立性验证（Grid Independence Study）

CFD网格独立性验证（Grid Independence Study）的理论基础

概要

老师，CFD的网格独立性验证和结构FEM的格子收敛确认本质上是一样的吗？

原理上是一样的。两者都在验证"网格再精细化，解也几乎不变"这一状态。都用Richardson外推来估计网格尺寸$h \to 0$时的解，并用GCI（网格收敛指数）来量化网格误差。

那结构分析的经验能直接套用到CFD吗？

要注意几个额外条件。第一，y+（壁面第一层网格的无量纲距离）必须与湍流模型相匹配。不匹配的话，网格再细也收敛不到正确的解。第二，CFD的非线性强度大，观测到的收敛阶$p$通常低于理论值（如二阶精度的模式中$p=2$）。第三，反复迭代残差必须足够小，否则迭代误差会污染离散化误差，使GCI失效。

比结构分析复杂不少呢。有系统的标准步骤吗？

有的。Celik等人（2008）在论文《Procedure for Estimation and Reporting of Uncertainty due to Discretization in CFD Applications》中提出的五步法是事实上的行业标准。ASME期刊投稿审查也用这个标准。掌握了这五步，无论是工程实践还是论文发表都能用。

Richardson外推法

Richardson外推的基本思想是什么？

离散解$f(h)$的Taylor展开式与真实解$f_{\text{exact}}$的关系可以表示为：

其中$h$是代表单元尺寸，$p$是离散化方法的精度阶数，$C$是未知常数。比如二阶中心差分理论上$p=2$，网格减半时误差约为$1/4$。通过使用不同网格尺寸的两个或多个解，我们可以用Richardson外推法估计$f_{\text{exact}}$。

既然$p$和$C$都不知道，2个方程需要2个未知数，为什么还要3层级网格？

提得好。用3层网格（细$h_1$、中$h_2$、粗$h_3$，其中$h_1 < h_2 < h_3$）和对应的解$f_1, f_2, f_3$，我们能从方程组解出观测到的收敛阶$p$。如果只有2层，就只能假设$p$值，这样就变成了主观判断"网格细化后变化小所以可以"，而不是严格的量化。

GCI的定义和意义

GCI的计算公式是什么？和结构分析的一样吗？

公式是一样的。细网格和中等网格之间的GCI为：

各项的含义是：

• $e_a^{21} = (f_1 - f_2)/f_1$：细网格解$f_1$和中等网格解$f_2$的相对差异
• $r_{21} = h_2 / h_1$：格子细分化比（$r > 1$）
• $p$：观测收敛阶（通过下面的公式计算）
• $F_s$：安全系数。用3层网格的外推时$F_s = 1.25$，假设$p$时$F_s = 3.0$

$F_s = 1.25$和$F_s = 3.0$差别这么大，肯定要用3层网格了？

正是。$F_s = 3.0$是对"假设$p$值"这一不确定性的惩罚。只多花一层网格的计算，安全系数就能从$3.0$降到$1.25$，相当于GCI缩小到$1/2.4$，这个投资回报率很高。比如汽车外形阻力系数$C_d$的GCI小于1.5%，就可以和风洞试验数据对标。

Celik的五步法

Celik的五步法具体是什么步骤？

我来详细讲Celik等人（2008，JFE 130(7)）的五步法。

第1步：定义代表单元尺寸$h$

三维情况下：

$N$是总网格单元数，$\Delta V_i$是每个单元的体积。二维时用面积的平方根。

第2步：用3层网格求解

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

第3步：计算观测收敛阶$p$

其中$\varepsilon_{21} = f_2 - f_1$，$\varepsilon_{32} = f_3 - f_2$。如果$r_{21} = r_{32}$（等比细分化），则$q=0$，可以用简化版本：

第4步：计算Richardson外推值$f_{\text{ext}}^{21}$

第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$是壁面到第一层单元中心的距离，$u_\tau$是摩擦速度，$\nu$是运动粘度。关键是，每个湍流模型都有明确要求的y+范围：

• 壁函数模式（标准$k$-$\varepsilon$等）：$y^+ = 30$～$300$
• 壁面分解（$k$-$\omega$ SST、SA、LES）：$y^+ < 1$（理想是$y^+ \approx 1$）
• 混合壁面处理（SST自动壁面处理）：$y^+ < 5$较好

如果用壁函数模式但y+只有5左右，网格独立性验证时会怎样？

壁函数基于对数律（$y^+ > 30$的区间），$y^+ = 5$处于缓冲层（buffer layer），既不是对数律也不是线性粘性区，"物理上不相容"。那边的壁面剪应力会严重错误。网格越精细，y+越小，越偏离壁函数的适用范围。结果是解出现发散或非单调行为，GCI的前提（单调收敛）就破了。

所以网格独立性验证前，先要检查y+在允许范围内？

完全正确。工程现场常见的坑是"用壁函数计算，但没注意有些壁面y+只有3，后来做网格研究发现'收敛不了'就很困惑"。这等于没把体温计夹好就量温度，当然会出问题。

CFD网格独立性验证（Grid Independence Study）的数值计算方法

代表单元尺寸的定义

老师，代表单元尺寸$h$具体怎么算？非结构网格单元大小不均匀啊。

Celik方法用全局平均。三维时所有单元平均体积的立方根：

工程上很简单——只要知道总单元数$N$和总计算域体积$V_{\text{total}}$，就能快速算出$h = (V_{\text{total}}/N)^{1/3}$。求解器的输出里一般都有单元数，所以算$h$很方便。

但要是在看壁面摩擦这种局部量，用全局平均的$h$合适吗？

提得很敏锐。严格上说，只有在评估的物理量依赖的区域网格均匀细化时，全局$h$才合适。比如只在壁面附近细，远方还是粗，那细分化比$r$的实际含义就变了。最佳实践是全域均匀相似细化——保持几何比例，全体单元数均匀倍增。Fluent中Mesh > Controls > Sizing可以一次性缩放所有参数。

观测收敛阶的计算

实际CFD问题中$p$一般是多少？二阶精度方案不就是$p=2$吗？

理想情况是这样，但工程实际中$p$往往分散在0.5～3.0之间。非线性问题常见$p < 2$。粗略判断标准是：

• $p \approx 1.5$～$2.5$：正常。说明在二阶精度的渐近区内
• $p < 0.5$：网格还太粗或有bug
• $p > 3.0$：误差相消可能性大。应再加一层网格验证
• $p < 0$（$\varepsilon_{32}/\varepsilon_{21} < 0$）：非单调收敛。标准方法不适用

$p$算出来1.8，这样还能用来算GCI吗？

完全没问题。Celik方法的优势就在于不固定$p$为理论值，而是用观测值，这样对非线性问题更诚实。拿$p=2$硬套反而不如用观测的$p=1.8$靠谱。

外推值和GCI求解

具体数字例子呢？比如管道内流的压力降，3层网格得$\Delta P_1 = 245.3$Pa、$\Delta P_2 = 248.1$Pa、$\Delta P_3 = 258.7$Pa？

等比细分化$r = 2.0$的计算例子：

• $\varepsilon_{21} = f_2 - f_1 = 248.1 - 245.3 = 2.8$
• $\varepsilon_{32} = f_3 - f_2 = 258.7 - 248.1 = 10.6$
• $p = \ln(10.6/2.8)/\ln(2.0) = \ln(3.786)/\ln(2.0) = 1.331/0.693 = 1.92$
• $f_{\text{ext}}^{21} = (2.0^{1.92} \times 245.3 - 248.1)/(2.0^{1.92} - 1) = (3.784 \times 245.3 - 248.1)/2.784 = 244.3$ Pa
• $e_a^{21} = |245.3 - 248.1|/245.3 = 1.14\%$
• $\text{GCI}_{\text{fine}}^{21} = 1.25 \times 0.0114 / (2.0^{1.92} - 1) = 0.01425/2.784 = 0.51\%$

GCI 0.51%！看起来不错，这样就能说"网格独立"了？

作为全局压力降来说，GCI < 1%足够了。但要同时用多个评估量确认。比如压力降收敛了，但壁面摩擦系数或剥离点位置可能还在动。最低限度要报告"全局量（力系数、流量）"和"局部量（壁面y+、某截面速度分布）"两类的GCI。

收敛判定的陷阱

网格独立性验证容易漏掉的坑有哪些？

常见的三个坑：

• 迭代残差收敛不足：粗网格残差$10^{-4}$就停了，但细网格要更多迭代才能达到同级别。要所有网格都收敛到至少$10^{-5}$，否则迭代误差会污染离散化误差。
• 忽视非定常效应：定常分析残差降不下来，可能是物理本身非定常（如Karman涡街）。这时要改非定常，用时间平均值计算GCI。
• 边界层网格不一致：整体细化但充气层（inflation layer）的第一层高$\Delta y_1$和生长率没改，导致y+变化不一致。壁面应按第一层高度指定，用生长率控制，确保与全体细化保持相容。

CFD网格独立性验证（Grid Independence Study）的实务应用

五步实践流程

具体工作中应该怎么操作？

工程实用的五步清单：

1. 确定评估量：最开始就决定用什么指标判断"网格独立"（阻力系数、出口温度、压力降等）。至少要全局量+局部量各一个。
2. 制作基础网格：先用中等网格算一遍，检查y+分布。如果不符合湍流模型要求，先改充气层，再做研究。
3. 生成3层网格：以基础网格为参考，用$r \approx 1.3$～$2.0$生成粗、细两层。全域相似细化很关键。
4. 充分收敛求解3层网格：残差降至$10^{-5}$以下。评估量的监测值要趋于平坦。
5. 用Celik法算GCI并报告：汇总$p$、$f_{\text{ext}}$、$e_a$、GCI的表格。

评估量的选择

就看全局量比如阻力系数不行吗？

不够。圆柱周围流，阻力系数$C_d$收敛了，后流涡结构或剥离点位置可能还在变。推荐的评估量组合：

网格细分化比的设计

$r = 1.3$和$r = 2.0$计算成本差距大吗？

三维网格单元数是$r^3$倍，差距很大：

• $r = 1.3$：单元数$1.3^3 = 2.2$倍。3层累计粗网格的$2.2^2 = 4.8$倍
• $r = 1.5$：单元数$1.5^3 = 3.4$倍。3层累计$3.4^2 = 11.4$倍
• $r = 2.0$：单元数$2.0^3 = 8$倍。3层累计$8^2 = 64$倍

Celik推荐$r \geq 1.3$。$r < 1.3$网格差异太小，容易被舍入误差淹没。$r = 1.5$是"精度和成本平衡"的人气选择。

Fluent、STAR-CCM+、OpenFOAM的具体实现

各软件怎么具体做网格独立性验证？

各软件的实践诀窍：

Ansys Fluent：Mesh > Adapt > Uniform可做均匀细化，但这是六面体网格各方向二分，$r=2$固定。非结构网格最好在Meshing端用Size参数全局缩放。Fluent 2024R2起加了Adjoint Solver的网格适应，但GCI计算仍需手工均一细化。

Simcenter STAR-CCM+：Parts > Mesh Operations只需改"Base Size"，网格自动相似重建。这是STAR-CCM+最大的优势，3层网格生成轻而易举。还能写Java宏"改Base Size → 重网格 → 求解 → 读结果"全自动，甚至有人直接写宏一口气算到GCI。

OpenFOAM：blockMesh改单元数很直接。snappyHexMesh的maxLocalCells、maxGlobalCells、refinementLevel调参控制细分化。GCI一般用Python脚本，从postProcessing读评估量，自动算$p$和GCI，比较流畅。

品质保证检查清单

报告或评审时要证明"确实做过网格独立性验证"，需要什么清单？

• 3层及以上网格的求解结果
• 全部网格用同样离散格式、同样湍流模型
• 全部网格迭代残差达$10^{-5}$以下
• 全部网格y+符合湍流模型要求范围
• 细分化比$r \geq 1.3$
• 观测收敛阶$p$已计算，与理论值对比
• GCI用全局和局部评估量都报告了
• 如非单调收敛，有明确处理说明
• 渐近域检查：$\text{GCI}_{\text{coarse}}/(\text{GCI}_{\text{fine}} \times r^p) \approx 1$

CFD网格独立性验证（Grid Independence Study）的软件对比

自动网格独立性验证功能对比

有软件能自动算GCI吗？

一体化自动GCI的工具还很少。但各软件的网格细化和自动化能力可比较：

脚本化的易用性

如果要全自动从3层网格算出GCI，哪个软件最顺手？

我觉得OpenFOAM最容易。三个理由：

1. 字典文件是纯文本，sed和Python一键替换参数
2. postProcessing自动输出评估量，取数据无痛
3. shell脚本"网格生成→求解→结果取→GCI计算"一体贯通

STAR-CCM+的Java宏特别强，Base Size一改全网格自动重建+重新求解，工作流最整洁。Fluent的Journal文件自动化能力也行，但网格重建依赖Meshing工具，工作流有点分散。

CFD网格独立性验证（Grid Independence Study）的前沿研究

自适应网格细分化（AMR）的关系

用了AMR（自适应网格）就不用验证网格独立性了吧？

不行。AMR自动加细梯度大的地方，提高计算效率，但"解不依赖网格"是另一回事。AMR是效率工具，不是精度保证。用了AMR也要最后验证结果的网格无关性。

可AMR网格不均匀啊，Celik法的前提和它矛盾吗？

确实矛盾。实际做法有(1)改AMR的refinement level上限生成3层、(2)用加权平均算代表$h$、(3)只在评估量密集区算GCI，等等。都不是完美方案，但总比没验证好。

LES/DNS的特殊问题

LES或DNS也能用GCI验证吗？

这是CFD网格独立性验证最难的部分。LES中滤波宽度（≈网格尺寸）是物理模型的一部分。改网格就是改方程本身。RANS的事情完全不适用。

LES的网格敏感性评估通常用：

• Pope准则：分解的湍流动能占比$M = k_{\text{resolved}}/(k_{\text{resolved}} + k_{\text{sgs}}) > 0.8$
• 能谱：能看到$-5/3$的惯性区远离网格截断频率
• 二点相关：相关长度远大于单元尺寸

DNS理论上用Kolmogorov长度$\eta$以下网格就独立了，但成本$\propto Re^{9/4}$，高Re数基本不可行。

机器学习的格子敏感性预测

最近有AI预测"这个网格够不够"的研究吗？

有两个方向。一是用粗网格解+局部网格特征（单元尺寸、长宽比、y+等）训练Surrogate Model，预测细网格解变化，省去细网格计算。二是从历史网格收敛数据库学习，对新问题推断需要的网格密度。但这些还停留在研究阶段，实务中还是Celik标准为准。

CFD网格独立性验证（Grid Independence Study）的故障处理

收敛阶$p$异常低（$p < 0.5$）

算出来$p = 0.3$，什么问题啊？

$p$太小的原因和对策：

• 还没进渐近区：粗网格太粗。换更精细的基础网格重来。
• 迭代不充分：细网格残差可能只降到$10^{-3}$。要压到$10^{-6}$。
• y+不匹配：3层网格的y+分布有的进壁函数范围，有的出范围。检查所有网格的y+分布。
• 格式不一致：是不是用了风上差分（一阶精度）？全用二阶以上精度的同一格式。
• 评估点在奇点：角点、再附着点等奇特性位置的值GCI会低。改到光滑区。

GCI过大

GCI算出15%，这是网格不够吧？

对，当前细网格离散化误差还太大。对策：

• 再细化：计算能力允许就最直接。
• 升高精度格式：一阶升二阶，二阶用MUSCL/QUICK等。$p$提升后GCI会显著下降。
• 聚焦投资：把网格密度集中在评估量相关的区域（壁面附近、尾流等），非关键区保持粗。
• 改评估量：最大局部值（最大壁面剪应力等）GCI常很大。改面积平均值看行不行。

非单调收敛的处理

3层数据$f_1 = 10.2$、$f_2 = 10.5$、$f_3 = 10.3$振动了，怎么办？

$\varepsilon_{32}/\varepsilon_{21} < 0$，非单调收敛，Celik法不适用（$p$无定义）。处理方案：

1. 检查迭代：残差真的全收敛了吗？定常求解残差有没有振荡？
2. 检查流物理：是不是物理本身非定常（涡街、弛豫等）？改非定常求解，用时间平均值。
3. 加第4层网格：可能某层有异常，添加网格找单调对子。
4. 保守报告：$p$求不出就用$e_a^{21}$（相对近似误差）作"离散化误差保守估计"报告。AIAA-G-077允许这个。

网格独立性验证的全流程我明白了。关键是：先确认y+，再3层网格，用Celik法算GCI，多个评估量齐报告。

完全正确。还有最后一点——网格独立性验证只是Verification（能否正确求解方程）的一部分。要得到可信的CFD，还要加Validation（用的物理模型对不对）。两条腿都要走。