中心差分格式
理论:由两侧相邻点求梯度
概述
CFD中用了中心差分,温度分布出现了不可能的凹凸。是bug吗?
不是bug,而是中心差分的性质。它用两侧相邻值之差求某点梯度,二阶精度,优点是不引入人工扩散。但当流动(对流)相对扩散很强时,下游值对本点起反向作用,解逐点锯齿状振荡。判断标准是“网格佩克莱数”是否超过2。对策是加密网格,或混入迎风类格式。
离散方程
对一维稳态对流扩散方程的对流项和扩散项都用中心差分离散的形式。右边第二个系数在流动强时变为负。
网格佩克莱数条件
若系数全部非负,解就是相邻值的加权平均,不会振荡;该条件即网格佩克莱数≤2。右边是下面算例用的精确解,$Pe = uL/\Gamma$ 为全局佩克莱数。
“迎风”的思想
中心差分的振荡从20世纪50年代起就困扰着数值计算研究者。1952年Courant、Isaacson和Rees提出使用信息来源方向——上游(迎风)值的差分,人们认识到不振荡的代价是精度降到一阶。20世纪80年代Harten在数学上整理了不产生振荡的条件(TVD),只在变化剧烈处偏向迎风的带限制器格式得到普及。今天通用CFD软件的对流格式选项,就建立在这场长期讨论之上。
算例:一维对流扩散
在区间[0, 1]上左端0、右端1、全局佩克莱数50的稳态对流扩散,用中心差分(CDS)和一阶迎风(UDS)求解。精确解在右端附近急剧上升:
| 分段数 | 网格Pe | 中心差分最小值 | 中心差分最大误差 | 迎风最大误差 |
|---|---|---|---|---|
| 10 | 5.0 | −0.429 | 0.436 | 0.160 |
| 20 | 2.5 | −0.111 | 0.193 | 0.204 |
| 25 | 2.0 | 0.000 | 0.135 | 0.198 |
| 50 | 1.0 | 0.000 | 0.035 | 0.132 |
| 100 | 0.5 | 0.000 | 0.008 | 0.077 |
| 200 | 0.25 | 0.000 | 0.002 | 0.042 |
10分时,不振荡的迎风格式误差反而更小。
对。粗网格上,中心差分让本应在0到1之间的解振荡到−0.43——就像温度低于绝对零度,物理上不可能。但从网格佩克莱数降到2的25分开始振荡消失,越加密误差按网格间距平方减小。100分时中心差分0.008,迎风0.077,相差10倍。迎风格式为一阶精度,加密后误差也减得很慢。这是“稳定但粗糙”与“精确但有条件”之间的取舍。
例2:收敛阶
佩克莱数10时加密网格,中心差分的收敛阶趋近2.00,迎风趋近0.97(200分时中心差分0.0001,迎风0.0090)。
格式选择
- 由流速、网格间距和扩散系数(黏度)估算网格佩克莱数分布。
- 稳态RANS以二阶迎风类或带限制器的格式为基本。
- LES中数值扩散会抹掉湍流涡,所以用中心差分类格式(必要时混入少量迎风)。
- 对剧烈变化(激波、薄层)使用抑制振荡的限制器。
- 加密网格确认结果不变。
“漂亮的结果其实是错的”
用CFD研究管内混合时,浓度分布很漂亮也没有振荡,但混合进展比实验快得多。为了稳定用了一阶迎风,数值扩散是实际扩散的好几倍。改为二阶格式后,重现了实验中较慢的混合。不振荡的结果看起来令人放心,但可能因数值扩散而“混得过头”。越是粗糙的迎风结果,越需要改变网格细度来确认。
常见错误
错误与对策
| 错误 | 后果 | 对策 |
|---|---|---|
| 粗网格用中心差分 | 振荡、超出范围的值 | 网格佩克莱数≤2 |
| 一直用一阶迎风 | 数值扩散导致混合过度 | 二阶格式 |
| LES中用迎风 | 湍流涡被抹掉 | 中心差分类格式 |
| 把振荡误认为物理现象 | 错误结论 | 改变网格确认 |
| 没振荡就认为正确 | 漏掉数值扩散 | 确认网格收敛 |
我想了解相关内容。
可以用Richardson外推工具试算收敛阶。格式详细比较参考迎风格式,基础参考有限体积法基础,时间方向参考CFD时间积分,高精度方法参考谱方法,压力求解参考压力-速度耦合求解器。
详细
错误