CFD的时间积分法
理论：显式与隐式格式
概述
我用CFD算圆柱后的卡门涡，涡很快消失，变成了稳态流动。哪里出错了？
检查一下时间积分方法。很多软件默认用一阶隐式（后向欧拉），非常稳定，但使振荡逐渐减弱的“数值阻尼”很强。配合大时间步长，可能连本应持续的涡脱落都被削弱。显式格式时间步长受CFL条件严格限制，隐式大步长也稳定，但能否正确保持振荡取决于方法和步长。用简单振荡模型比较，差别一目了然。
每步放大率
把各方法用于振荡方程时的每步放大率 $G$。$|G| = 1$ 振幅保持，$|G| < 1$ 衰减，$|G| > 1$ 发散。
时间步长限制
显式格式要求每步流体不越过一个单元（CFL约≤1），还有扩散限制。
时间步进始于天气预报
1922年理查森尝试用手算做天气预报，6小时的预报花了几周，结果还大错特错。原因之一是当时不知道时间步长与网格的关系决定计算的稳定性。1928年库朗、弗里德里希斯和列维证明，显式格式的时间步长若不小于由网格和波速决定的值，解就会发散（CFL条件）。为避开这一限制发展了隐式方法，代价是要面对数值阻尼这个新问题。
算例：用振荡模型比较
用不同的每周期步数求解无阻尼振荡（振幅本应不变），比较10个周期后剩余振幅和每周期相位误差：
| 每周期步数 | 显式欧拉 | 隐式欧拉 | BDF2 | Crank-Nicolson |
|---|---|---|---|---|
| 10 | 发散 | 6×10⁻⁸ | 0.14 | 1.00 |
| 20 | 发散 | 8×10⁻⁵ | 0.68 | 1.00 |
| 40 | 发散（131倍） | 0.0076 | 0.944 | 1.00 |
| 80 | 发散（11.7倍） | 0.086 | 0.993 | 1.00 |
每周期相位误差：40步时后向欧拉和BDF2为2.9°，Crank-Nicolson为0.7°；80步时分别为0.7°和0.2°。
一阶隐式即使40步，10个周期后振幅也不到1%。
对。一阶隐式每步都削减一点振幅，算得越久影响越大。实际卡门涡靠流动不稳定性自己产生振荡，不一定完全消失，但强度和频率会改变。BDF2为二阶精度，阻尼小得多，是实务非稳态计算的标准。Crank-Nicolson保持振幅，但在剧烈变化处容易残留振荡。比如直径0.1 m圆柱在10 m/s流动中，涡脱落频率约20 Hz，每周期40步时时间步长1.25 ms；单元1 mm时CFL为12.5，显式无法求解，只能用隐式。
选择参考
- 估算要解析的最快变化（涡脱落、脉动等）的周期。
- 要看非稳态结果，用BDF2等二阶方法而非一阶隐式。
- 按每周期20～40步以上确定时间步长，并确认CFL。
- 把时间步长减半，确认频率和振幅不变。
- 也要确认每步内的迭代（外迭代）充分收敛。
“频率对了，幅值只有一半”
某研究桥梁断面气动振动的CFD，涡脱落频率与风洞实验吻合，但升力脉动幅值只有一半左右。时间积分用的是一阶隐式，每周期约15步。改为BDF2、每周期40步后，脉动幅值与实验相差一成以内。频率反映相位误差，幅值反映振幅衰减，频率对了不代表幅值也对。
常见错误
错误与对策
| 错误 | 后果 | 对策 |
|---|---|---|
| 非稳态用一阶隐式 | 振荡衰减或消失 | BDF2等二阶 |
| 时间步长太粗 | 相位与振幅误差 | 每周期20～40步以上 |
| 显式忽略CFL | 发散 | 确认CFL |
| 只用频率验证 | 漏掉振幅误差 | 也比较振幅 |
| 步内迭代不足 | 时间精度下降 | 确认残差 |
我想了解相关内容。
详细
错误