盖驱动腔流 — CFD基准验证指南
盖驱动腔流的理论基础
概要
老师，盖驱动腔流是什么基准？我只知道"这是CFD验证中要做的事情"...
简单来说，盖驱动腔流是CFD代码验证中最广泛使用的基准问题。一个正方形的盒子（腔）的上壁以恒定速度水平滑动，导致内部产生循环流。计算这个流场，然后与Ghia、Ghia和Shin（1982）使用高精度差分法计算的参考解进行比较，以检查求解器是否正确实现。
为什么这个问题这么有名？也有更简单的流动，对吧？
很好的问题。有三个原因。首先，形状极其简单（正方形、仅有墙壁），所以网格生成中没有歧义。其次，能得到非平凡的流场——循环涡、角部的二次涡等，包含对流与粘性的相互作用。与管内流等平行流不同，离散化方案的质量在结果中清晰可见。第三，Ghia的论文公开了Re=100到10000的7个条件下的速度分布数值表，使定量比较变得容易。
支配方程式
支配腔流的方程是什么？
二维、不可压、定常Navier-Stokes方程。成对求解连续方程和动量方程：
边界条件如下：
- 上壁（盖）：$u = U$、$v = 0$（以恒定速度滑动）
- 下壁、左壁、右壁：$u = 0$、$v = 0$（无滑移条件）
雷诺数定义为 $\text{Re} = UL/\nu$。其中$L$是腔的边长，$U$是盖的速度，$\nu$是动力粘度系数。
我听说也可以用流函数-涡度法求解。
Ghia的原始论文就是这样做的。使用流函数$\psi$和涡度$\omega$来消除压力项：
Ghia在$129 \times 129$的等间距网格上使用多重网格法进行了计算。这种方法不需要压力方程，适合二维问题，但难以扩展到三维。当前的商用CFD求解器主要采用原始变量法（速度-压力形式）。
涡结构与雷诺数依赖性
改变雷诺数时，流场如何变化？
雷诺数决定了流场的性质。让我们逐步看：
- Re=100：主要涡位于腔中心略上方。粘性占主导，涡形状对称且圆润。二次涡几乎看不到。
- Re=400：主要涡中心接近几何中心。左下角和右下角开始出现小的二次涡（角涡）。
- Re=1000：主要涡中心移至$(0.5313, 0.5625)$附近。左下二次涡（BL）和右下二次涡（BR）清晰可见。
- Re=5000：二次涡进一步增长。左上出现小的三次涡（TL）。主要涡形状向上偏移，接近椭圆。
- Re=10000：所有角涡都显著增长。定常解的存在处于临界点，微小扰动可能导致非定常转变。
等等，Re=10000也有定常解？不会变成湍流吗？
在二维情况下，Ghia给出了Re=10000的定常解。但收敛非常缓慢。Erturk等人（2005）的研究表明，Re=10000的定常解在$512 \times 512$网格上需要数十万步迭代。据报告，Re≈20000时会发生Hopf分岔，流动转为周期性。在三维情况下，非定常性在低得多的雷诺数（约Re=1000附近）就会出现。
Ghia参考数据
具体用什么数据来验证？
Ghia等人在通过腔中心的竖直线上的水平速度$u$和通过中心的水平线上的竖直速度$v$中各给出17个点的值。Re=1000时竖直中心线上的$u$速度分布是最常用的数据：
| $y$位置 | $u/U$ (Ghia, Re=1000) |
|---|---|
| 1.0000 | 1.00000 |
| 0.9766 | 0.65928 |
| 0.9688 | 0.57492 |
| 0.9609 | 0.51117 |
| 0.9531 | 0.46604 |
| 0.8516 | 0.33304 |
| 0.7344 | 0.18719 |
| 0.6172 | 0.05702 |
| 0.5000 | -0.06080 |
| 0.4531 | -0.10648 |
| 0.2813 | -0.27805 |
| 0.1719 | -0.38289 |
| 0.1016 | -0.29730 |
| 0.0703 | -0.22220 |
| 0.0625 | -0.20196 |
| 0.0547 | -0.18109 |
| 0.0000 | 0.00000 |
从这个表可以看出，$y \approx 0.17$处$u/U \approx -0.38$取最小值。这是主要涡下侧反向流最强的位置。将自己的CFD结果与此分布重叠绘制，所有点的相对误差都在1%以内，就表明求解器的实现值得信赖。
主要涡中心的流函数值也可以用于验证吗？
当然可以。Re=1000时主要涡中心的流函数值为$\psi_{\min} = -0.1189$。这是单个标量值，易于验证，但全速度分布的比较是更严格的验证。
验证数据的可视化
Re=1000时竖直中心线上的$u$速度分布比较。网格密度和离散化方案的影响的定量显示。
| 评估项 | Ghia参考值 | 计算值 (128×128, 2次CD) | 相对误差 [%] | 判定 |
|---|---|---|---|---|
| $\psi_{\min}$（主要涡中心） | -0.1189 | -0.1188 | 0.08 | PASS |
| $u_{\min}$（最大反向速度） | -0.3829 | -0.3821 | 0.21 | PASS |
| $v_{\max}$（水平中心线） | 0.3769 | 0.3758 | 0.29 | PASS |
| 主要涡中心$x$坐标 | 0.5313 | 0.5313 | 0.00 | PASS |
| 主要涡中心$y$坐标 | 0.5625 | 0.5625 | 0.00 | PASS |
判定标准：相对误差 < 1%：■ 优秀，1～5%：■ 可接受，> 5%：■ 需要检查
盖驱动腔流的数值计算方法
压力-速度耦合算法
求解腔流时，SIMPLE和PISO怎么区分使用？
求定常解就用SIMPLE（半隐式压力链接方程法）。每次迭代都求解一次压力修正方程，更新速度场。适当设置下松弛系数可以稳定收敛。Re=1000左右的情况下，速度下松弛系数0.7、压力下松弛系数0.3是不错的起点。
PISO什么时候用？
PISO（含分裂算子的隐式压力）是非定常分析用的。每个时间步进行两次压力修正，确保时间精度。当Re=10000以上需要非定常解时使用。仅求Re=1000的定常解的话，SIMPLE的计算效率更好。
| 算法 | 用途 | 压力修正次数/迭代 | 腔流中的推荐场景 |
|---|---|---|---|
| SIMPLE | 定常 | 1 | Re ≤ 5000的定常解 |
| SIMPLEC | 定常 | 1（强化版） | 收敛缓慢时的替代 |
| PISO | 非定常 | 2 | Re ≥ 10000的非定常解 |
| COUPLED | 定常/非定常 | 耦合求解 | 高Re数的高速收敛 |
对流项离散化方案
一阶迎风法和二阶中心差分结果差别很大吗？
在腔流中，方案选择的影响很明显。来看具体例子。Re=1000、使用相同网格（64×64）的情况下：
| 离散化方案 | 网格 | $\psi_{\min}$ | $u(x=0.5, y=0.1719)$ | $\psi$误差 [%] |
|---|---|---|---|---|
| 一阶迎风（UDS） | 32×32 | -0.1102 | -0.348 | 7.31 |
| 一阶迎风（UDS） | 64×64 | -0.1158 | -0.365 | 2.61 |
| 二阶中心差分（CDS） | 32×32 | -0.1170 | -0.372 | 1.60 |
| 二阶中心差分（CDS） | 64×64 | -0.1185 | -0.380 | 0.34 |
| 二阶中心差分（CDS） | 128×128 | -0.1188 | -0.382 | 0.08 |
| QUICK | 64×64 | -0.1186 | -0.381 | 0.25 |
一阶迎风法由于人工扩散（人工粘性）很大，32×32时误差超过7%。实际工作中，Re=1000的情况下，至少要用二阶精度方案加64×64或更细的网格。QUICK方案（三阶精度上游内插）的话，64×64就能达到足够精度。
不应该用一阶迎风法吧？
对于验证目的不太合适。但是收敛性非常好，可以用于初始解的推算，或者二阶方案发散时的调试。"先用一阶迎风法计算出流场的大致形状，再把结果作为初值换用二阶方案"这种分阶段方法在实务中也常用。
网格设计与收敛验证
网格是用等间距还是非等间距好？
取决于雷诺数。Re=100～400时等间距就够了，但Re=1000以上会在壁面附近产生陡峭的速度梯度，所以采用向壁面方向按几何级数密化的非等间距网格更有效。通常膨胀比（expansion ratio）为1.05～1.1，壁面最内层单元厚度约为整体的1/200。
怎么验证网格收敛性？
至少用三个不同尺度的网格（比如32×32、64×64、128×128）求解相同问题，检查关键物理量（$\psi_{\min}$、中心线速度等）的变化。用Richardson外推法评估收敛次数，计算GCI（网格收敛指数）是标准方法。
这里$f_1, f_2, f_3$分别是密、中、粗网格的结果，$r$是网格密度比（通常为2），$F_s$是安全系数（三个网格级别时为1.25）。使用二阶精度方案的话，收敛次数$p$应约为2。如果$p$与2相差很大，说明有问题。
验证数据的可视化
网格密度和离散化精度组合的精度比较（Re=1000）。
| 评估项 | Ghia参考值 | 计算值 (64×64, QUICK) | 相对误差 [%] | 判定 |
|---|---|---|---|---|
| $\psi_{\min}$ | -0.1189 | -0.1186 | 0.25 | PASS |
| $u_{\min}$ at $y$=0.1719 | -0.3829 | -0.3815 | 0.37 | PASS |
| 收敛次数 $p$ | 2.00（理论） | 1.95 | 2.50 | PASS |
| GCI (64→128) | < 1% | 0.42% | — | PASS |
判定标准：相对误差 < 1%：■ 优秀，1～5%：■ 可接受，> 5%：■ 需要检查
盖驱动腔流的实际应用
OpenFOAM中的实现
用OpenFOAM求解腔流怎么做？从最初的步骤开始教我。
其实OpenFOAM随附了tutorial案例cavity。在$FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity目录下有现成的，先从那里开始比较好。但默认是20×20网格，Re=10非常低，用于验证需要定制。
Re=1000情况下的设置关键点：
blockMeshDict：把单元数改成(128 128 1)。Z方向是1个单元（二维计算）。transportProperties：nu 0.001（$U=1$ m/s、$L=1$ m时）U边界条件：movingWall处设fixedValue uniform (1 0 0)，其他3面为fixedValue uniform (0 0 0)，frontAndBack为emptyp边界条件：全壁面zeroGradient，用1个单元设定基准压力fixedValue uniform 0- 求解器：定常时
simpleFoam，非定常时icoFoam/pisoFoam
fvSchemes怎么设置？
divSchemes中使用Gauss linearUpwind grad(U)就能得到二阶精度上游内插。Gauss linear（中心差分）也可以，但粗网格容易振荡。Gauss upwind（一阶迎风）不推荐验证用。laplacianSchemes用Gauss linear corrected实现二阶精度扩散项。
Ansys Fluent中的实现
Ansys Fluent怎么设置？
Fluent的步骤是这样的：
- 几何：用SpaceClaim或DesignModeler创建1m×1m的正方形
- 网格：在Meshing中应用Mapped Face Meshing（结构网格），128×128分割。用Edge Sizing给壁面添加偏置可提高精度。
- General设置：Pressure-Based、Steady、2D
- Viscous Model：Laminar（Re=1000是层流）
- Materials：密度$\rho = 1$ kg/m³、粘度$\mu = 0.001$ Pa·s
- Boundary Conditions：上壁设为Moving Wall（1 m/s），其他为Stationary Wall
- Solution Methods：SIMPLE、Momentum用Second Order Upwind
- Residuals：收敛到$10^{-6}$
怎么取出结果？想和Ghia数据比较。
在Fluent中Solution > Results > Plots > XY Plot，绘制竖直中心线（$x=0.5$）上的$u$速度和水平中心线（$y=0.5$）上的$v$速度。"Write to File"可以输出成CSV格式的数值数据，再用Excel或Python与Ghia数据重叠绘制。
结果验证步骤
计算完成后要检查什么？想要检查清单。
按以下5个步骤验证是标准程序：
- 残差确认：连续方程和动量残差是否收敛到$10^{-6}$以下
- 质量守恒确认：全壁面质量流量合计是否为零（闭合系统）
- 竖直中心线的$u$分布：与Ghia的17个点比较，全点相对误差目标1%以内
- 水平中心线的$v$分布：同样比较
- 主要涡中心位置和$\psi_{\min}$：确认与Ghia值的一致
能在多个Re数下验证就更好了。Re=100正确也不表示Re=1000一定没问题——对流方案的缺陷在高Re数时才会显现。至少要用Re=100和Re=1000的两个条件进行验证。
验证数据的可视化
OpenFOAM simpleFoam的验证结果（Re=1000、128×128网格、linearUpwind方案）。
| 评估项 | Ghia参考值 | 计算值 | 相对误差 [%] | 判定 |
|---|---|---|---|---|
| $u$ at $y$=0.9766 | 0.6593 | 0.6580 | 0.20 | PASS |
| $u$ at $y$=0.5000 | -0.0608 | -0.0602 | 0.99 | PASS |
| $u$ at $y$=0.1719 | -0.3829 | -0.3818 | 0.29 | PASS |
| $\psi_{\min}$ | -0.1189 | -0.1187 | 0.17 | PASS |
| BL涡 $\psi$ | -1.75e-5 | -1.72e-5 | 1.71 | PASS |
判定标准：相对误差 < 1%：■ 优秀，1～5%：■ 可接受，> 5%：■ 需要检查
盖驱动腔流的软件比较
CFD求解器间的结果比较
不同CFD求解器计算腔流结果会不一样吗？
正规的求解器在适当设置下，误差应该控制在Ghia值的1%以内。出现差异通常是"使用方法"的问题，方案选择、网格品质、收敛判定设置都会影响结果。看一下不同求解器的交叉验证例子：
| 求解器 | 网格 | 方案 | $\psi_{\min}$ (Re=1000) | 与Ghia相对误差 [%] |
|---|---|---|---|---|
| OpenFOAM (simpleFoam) | 128×128 | linearUpwind | -0.1187 | 0.17 |
| Ansys Fluent | 128×128 | 2nd Order Upwind | -0.1188 | 0.08 |
| STAR-CCM+ | 128×128 | 2nd Order Upwind | -0.1187 | 0.17 |
| COMSOL (FEM) | P2+P1, 64×64 | — | -0.1188 | 0.08 |
| FEniCS (FEM) | Taylor-Hood, 64×64 | — | -0.1189 | 0.00 |
基于有限元法（FEM）的COMSOL和FEniCS用Taylor-Hood单元（P2速度+P1压力）就能在比较粗的网格上达到高精度。这是因为速度场的多项式次数较高。基于有限体积法（FVM）的求解器要达到同等精度需要更细的网格，但每个单位计算成本的效率更高。
选择指南
最后，腔流验证用哪个求解器最好？
应根据目的区分使用：
- 教育学习：OpenFOAM的cavity教程最适合。免费，而且直接编辑配置文件能加深理解。
- 自有代码验证：直接和Ghia参考数据比较。评估求解器性能时在多个Re数上验证。
- 商用求解器初期验证：新导入的求解器"先计算腔流"是业界惯例。推荐Re=100、1000、5000三个条件。
- 研究目的（高精度）：用FEM求解器（FEniCS、deal.II）配合高阶单元，粗网格也能达到高精度。
验证数据的可视化
各求解器的$\psi_{\min}$比较（Re=1000、同等网格密度）。
| 评估项 | Ghia参考值 | 最佳值 | 相对误差 [%] | 判定 |
|---|---|---|---|---|
| OpenFOAM $\psi_{\min}$ | -0.1189 | -0.1187 | 0.17 | PASS |
| Fluent $\psi_{\min}$ | -0.1189 | -0.1188 | 0.08 | PASS |
| FEniCS $\psi_{\min}$ | -0.1189 | -0.1189 | 0.00 | PASS |
判定标准：相对误差 < 1%：■ 优秀，1～5%：■ 可接受，> 5%：■ 需要检查
盖驱动腔流的前沿研究
三维腔流
二维腔流已理解了。扩展到三维后怎么样？
三维腔（立方体，上面全部作为盖滑动）中，沿跨向的端壁影响加入。二维中看不到的Taylor-Görtler式纵向涡在端壁附近出现，流场变得复杂得多。Ku、Hirsh和Taylor（1987）或Koseff和Street（1984）的风洞实验数据用作三维验证的参考。
计算成本增加多少？
128×128的二维约16,000个单元，而128×128×128的三维有约200万个单元。内存和计算时间都是数量级的增加。三维验证作为最初步骤太重了，应该先通过二维验证，再进入三维。
高雷诺数与湍流转变
Re=10000以上怎么样？需要湍流模型吗？
二维情况下，Re=10000～20000附近发生Hopf分岔，定常解变得不稳定。产生周期性涡振荡，进一步升高Re数转为混沌行为。三维中远低的Re数（约Re=1000～3000）就会出现三维不稳定性。
求解高Re数腔流的方法有：
- DNS（直接数值模拟）：解所有湍流尺度。Re=10000左右的三维需要$512^3$以上网格。
- LES（大涡模拟）：直接解大涡，小涡用模型。DNS所需网格的1/10～1/100。
- RANS：仅求定常湍流统计量。标准$k$-$\varepsilon$模型对腔的再循环流精度低。RNG $k$-$\varepsilon$或SST $k$-$\omega$效果更好。
格子Boltzmann方法方法
听说格子Boltzmann法（LBM）也能计算腔流。
LBM不直接求解Navier-Stokes方程，而是在网格上模拟粒子的碰撞与移流。腔流也是LBM的标准基准，用D2Q9模型（二维9速度模型）能很好地再现Ghia值。
LBM的优势是易于并行化。每个网格点的计算是局部的，与GPU的兼容性非常好。近年NVIDIA推进的GPU版LBM求解器（如WaLBerla、Palabos）的高速腔流计算已有报告。不过高Re数的稳定性需要采用MRT（多松弛时间）模型等优化。
验证数据的可视化
各方法·Re数的验证结果总结。
| 方法/Re数 | Ghia参考值 $\psi_{\min}$ | 计算值 | 相对误差 [%] | 判定 |
|---|---|---|---|---|
| FVM, Re=100 | -0.10342 | -0.10340 | 0.02 | PASS |
| FVM, Re=400 | -0.11393 | -0.11380 | 0.11 | PASS |
| FVM, Re=1000 | -0.11890 | -0.11880 | 0.08 | PASS |
| LBM-BGK, Re=1000 | -0.11890 | -0.11850 | 0.34 | PASS |
| FVM, Re=5000 | -0.12200 | -0.12170 | 0.25 | PASS |
判定标准：相对误差 < 1%：■ 优秀，1～5%：■ 可接受，> 5%：■ 需要检查
盖驱动腔流的故障排查
中心线速度与Ghia值不符
计算收敛了，但与Ghia参考值相差5%以上。什么原因？
常见情况。逐步检查：
- 是否用了一阶迎风法：数值扩散大，涡变钝。改用二阶以上方案。
- 网格太粗：Re=1000至少要64×64，理想128×128以上。
- 收敛判定太松：残差降到$10^{-6}$。$10^{-3}$时速度分布会定量改变。
- Re数计算错误：OpenFOAM设kinematic viscosity$\nu$，Fluent设dynamic viscosity$\mu$。$\nu = \mu/\rho$换算错常见。
- 数据提取位置不精确：Ghia数据严格在$x=0.5$竖直线。偶数分割网格$x=0.5$可能在单元边界，需要插值。
用二阶方案、64×64网格，还是偏差2～3%...
64×64偏差2～3%是网格不足引起的正常范围。试试升到128×128，误差应该降到1%以下。计算GCI来定量评价网格依赖性是最确实的方法。
高雷诺数不收敛
Re=5000时，残差振荡不收敛。怎么办？
高Re数对流强，收敛困难。试这些对策：
- 逐步提高Re数：先在Re=100收敛→用该结果作Re=400的初值→→Re=1000→...逐步升高。直接跳到Re=5000，初值太远容易发散。
- 降低下松弛系数：SIMPLE速度下松弛0.7→0.5，压力下松弛0.3→0.2。收敛慢但稳定。
- 细化网格：粗网格高Re数单元雷诺数过大，方案不稳定。
- 改用非定常法：SIMPLE定常不收敛时，用PISO时间推进到定常。或者物理上高Re解本身是非定常的（Re>8000）。
处理角部奇异性
盖与侧壁的角上压力异常大。正常吗？
这是数学上的奇点，正常现象。上壁角处速度从$(U,0)$瞬间变为$(0,0)$。严格解中压力和涡度在此无穷发散。网格越细，该处压力值越大——这不是"未收敛"，而是奇点附近的正确行为。
实务应对：
- 验证不涉及角部：Ghia参考数据本身避开了角区。
- 平滑盖速度：某些求解器（如FreeFEM++）可对盖速度在角附近平滑处理（如$u(x) = 16x^2(1-x)^2 U$）。奇点消失但无法直接比较Ghia值。
- 别过度细化角部网格：集中网格于角也不能改进解。浪费计算资源。
所以奇点"不用管"，中心线速度分布的验证就够了？
正是。腔流验证看的是全局流场再现能力，不是局部奇点。中心线速度分布、一次涡位置和强度、二次涡出现——这些都对上参考值，求解器实现就信得过。
验证数据的可视化
故障排查前后精度对比例。
| 情况 | Ghia $\psi_{\min}$ | 计算值 | 相对误差 [%] | 判定 |
|---|---|---|---|---|
| 一阶迎风, 32×32 | -0.1189 | -0.1102 | 7.31 | FAIL |
| 一阶迎风, 128×128 | -0.1189 | -0.1172 | 1.43 | MARGINAL |
| 二阶CDS, 64×64 | -0.1189 | -0.1185 | 0.34 | PASS |
| 二阶CDS, 128×128 | -0.1189 | -0.1188 | 0.08 | PASS |
判定标准：相对误差 < 1%：■ 优秀，1～5%：■ 可接受，> 5%：■ 需要检查
更详细
错误