盖驱动腔体流动 — 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$ 的等间距网格上使用了多重网格法进行计算。这种方法因为不需要压力方程，所以适合2D问题，但难以扩展到3D。现在的商用CFD求解器主流是速度-压力形式（原始变量法）。
涡结构与雷诺数依赖性
改变Re数，流场会怎么变化？
Re数决定性地支配流场的性质。我们分阶段来看：
- Re=100：一次涡位于方腔中心稍偏上的位置。粘性占主导，涡的形状接近对称的圆形。二次涡几乎看不见。
- Re=400：一次涡中心接近几何中心。左下和右下开始出现小的二次涡（角涡）。
- Re=1000：一次涡中心移动到 $(0.5313, 0.5625)$ 附近。左下二次涡（Bottom-Left, BL）和右下二次涡（Bottom-Right, BR）可以明确确认。
- Re=5000：二次涡进一步成长。左上也开始出现小的三次涡（Top-Left, TL）。一次涡形状向上偏移，接近椭圆形。
- Re=10000：所有角涡显著成长。定常解的存在处于临界状态，微小的扰动就可能使其过渡到非定常。
诶，Re=10000也存在定常解吗？不会变成湍流吗？
在2D情况下，直到Re=10000，Ghia都得到了定常解。不过收敛非常慢。根据Erturk等人(2005)的研究，Re=10000的定常解需要 $512 \times 512$ 网格和数十万步的迭代。一般认为从Re=20000左右开始会发生Hopf分岔，向周期性流动过渡。在3D情况下，可能在更低的Re数（大约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$。因为可以用一个标量值验证，所以很方便，但速度剖面的整体比较是更严格的验证。
各项的物理意义
- 对流项 $u \partial u / \partial x$：流体粒子以其自身速度输运动量的项。Re数大时对流胜过粘性，产生急剧的速度梯度。方腔上部被顶盖拖拽的流体向下潜入的现象与此对应。
- 粘性扩散项 $\nu \nabla^2 u$：相邻流体层间剪切引起的动量扩散。Re数小时粘性占主导，速度场变得平滑。Re=100时涡呈圆形就是因为这个。
- 压力梯度项 $-\frac{1}{\rho}\partial p / \partial x$：为满足不可压条件的压力场。在方腔流中没有入口出口，因此压力与涡旋转产生的离心力相平衡。
假设条件与适用范围
- 不可压流动（马赫数 $\ll$ 0.3）
- 牛顿流体（恒定粘度）
- 2D流动假设（展向无限长）
- 顶壁与侧壁角部存在速度不连续（数学奇点）
验证数据可视化
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（Semi-Implicit Method for Pressure-Linked Equations）是标准选择。每次迭代求解一次压力修正方程，更新速度场。如果适当设置欠松弛因子，可以稳定收敛。对于Re=1000左右，速度欠松弛0.7、压力欠松弛0.3左右可以作为起点。
PISO在什么场景下使用？
PISO（Pressure-Implicit with Splitting of Operators）适用于非定常分析。每个时间步进行两次压力修正，以确保时间精度。在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数。Re=100〜400的话等间距就足够了，但Re=1000以上时壁面附近会产生急剧的速度梯度，因此向壁面方向按几何级数加密的不等间距网格更有效。典型做法是，膨胀比（expansion ratio）取1.05〜1.1，壁面第一个网格厚度约为整体的1/200。
网格收敛性验证怎么做？
なった
詳しく
報告