第一部分,我打算介绍一个最基本的算例,一维激波管问题。说白了就是
一根两端封闭的管子,中间有个隔板,隔板左边和右边的气体状态(密度、
速度、压力)不一样,突然把隔板抽去,管子内面的气体怎么运动。这是
个一维问题,被称作黎曼间断问题,好像是黎曼最初研究双曲微分方程的
时候提出的一个问题,用一维无粘可压缩 Euler 方程就可以描述了。
这里
这个方程就是描述的气体密度 、动量 和能量 随时间的变化( )
与它们各自的流量(密度流量 ,动量流量
,能量流量
)随空间变化( )的关系。
在 CFD 中通常把这个方程写成矢量形式
这里
进一步可以写成散度形式
一定要熟悉这种矢量形式
以上是控制方程,下面说说求解思路。可压缩流动计算中,有限体积(FVM)
是最广泛使用的算法,其他算法多多少少都和 FVM 有些联系或者共通的
思路。了解的 FVM,学习其他高级点的算法(比如目前比较热门的间断有
限元、谱 FVM、谱 FDM)就好说点了。
针对一个微元控制体
,把 Euler 方程在空间积分
用微积分知识可以得到
也就是说控制体内气体状态平均值的变化是控制体界面上流通量的结果。
因此我们要计算 的演化,关键问题是计算控制体界面上的 。
FVM 就是以这个积分关系式出发,把整个流场划分为许多小控制体,每个
控制体和周围相邻的某个控制体共享一个界面,通过计算每个界面上的通
量来得到相邻控制体之间的影响,一旦每个控制体的变化得到,整个流场
的变化也就知道了。
所以,再强调一次,关键问题是计算控制体界面上的 。
初学者会说,这个不难,把界面 上的 插值得到,然后就可以计算 。
有道理!
咱们画个图,有三个小控制体 i-1 到 i+1,中间的“|”表示界面,控制体 i
右边的界面用
表示,左边的就是
。
| i-1 | i | i+1 |
好下个问题:每个小控制体长度都是 如何插值计算界面
上的
?
最自然的想法就是:取两边的平均值呗,
但是很不幸,这是不行的。
那么换个方法?直接平均得到
?
还是很不行,这样也不行。
我靠,这是为什么?这明明是符合微积分里面的知识啊?
这个道理有点复杂,说开了去可以讲一本书,可以说从 50 年代到 70 年
代,CFD 科学家就在琢磨这个问题。这里,初学者只需要记住这个结论:
对于流动问题,不可以这样简单取平均值来插值或者差分。如果你非要想
知道这个究竟,我现在也不想给你讲清楚,因为我眼下的目的是让你快速
上手,而且该不刨根问底的时候就不要刨根问底,这也是初学阶段一种重
要的学习方法。
好了,既然目的只是为了求
,我在这里,只告诉你一种计算方法,也
是非常重要、非常流行的一种方法。简单的说,就是假设流动状态在界面
是不连续的,先计算出界面
两边 的值,
和
,再由
它们用某种方法计算出
。上述方法是非常重要的,是由一个苏联人
Godunov 在 50 年代首创的,后来被发展成为通用 Godunov 方法,著名
的 ENO/WENO 就是其中的一种。
好了,现在的问题是:
1 怎么确定
和
2 怎么计算
对于第一个问题,Godunov 在他的论文中,是假设每个控制体中 是均
匀分布的,因此
第二个问题,Godunov 采用了精确黎曼解来计算
。什么是“精确黎
曼解”,就是计算这个激波管问题的精确解。既然有精确解,那还费功夫搞
这些 FVM 算法干什么?因为只有这种简单一维问题有精确解,稍微复杂
一点就不行了。精确解也比较麻烦,要分析 5 种情况,用牛顿法迭代求解
(牛顿法是什么?看数值计算的书去,哦,算了,现在暂时可以不必看)。
这是最初 Godunov 的方法,后来在这个思想的基础上,各种变体都出来
了。也不过是在这两个问题上做文章,怎么确定
,怎么计算
。
Godunov 假设的是每个小控制体内是均匀分布,也就是所谓分段常数
(piecewise constant),所以后来有分段线性(picewise linear)或者分
段二次分布(picewise parabolic),到后来 ENO/WENO 出来,那这个假
设的多项式次数就继续往上走了。都是用多项式近似的,这是数值计算中
的一个强大工具,你可以在很多地方看到这种近似。
Godunov 用的四精确黎曼解,太复杂太慢,也不必要,所以后来就有各
种近似黎曼解,最有名的是 Roe 求解器、HLL 求解器和 Osher 求解器,
都是对精确黎曼解做的简化。
这个多项式的阶数是和计算精度密切相关的,阶数越高,误差就越小。不
过一般来说,分段线性就能得到不错的结果了,所以工程中都是用这个,
Fluent、Fastran 以及 NASA 的 CFL3D、OverFlow 都是用这个。而黎
曼求解器对精度的影响不是那么大,但是对整个算法的物理适用性有影
响,也就是说某种近似黎曼求解器可能对某些流动问题不合适,比如单纯
的 Roe 对于钝头体的脱体激波会算得乱七八糟,后来加了熵修正才算搞
定。
上次(http://gezhi.org/node/399)说到了求解可压缩流动的一个重要
算法,通用 Godunov 方法。其两个主要步骤就是
1 怎么确定
和
2 怎么计算
这里我们给出第一点一个具体的实现方法,就是基于原始变量的 MUSCL
格式(以下简称 MUSCL)。它是一种很简单的格式,而且具有足够的精
度,NASA 著名的 CFL3D 软件就是使用了这个格式,大家可以去它的主
页(http://cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html)上看手册,
里面空间离散那一章清楚的写着。
MUSCL 假设控制体内原始变量(就是 )的分布是一次或者二次多项式,
如果得到了这个多项式,就可以求出控制体 左右两个界面的一侧的值
和
。
我们以压力 p 为例来说明怎么构造这个多项式。这里我只针对二次多项式
来讲解,你看完之后肯定能自己推导出一次多项式的结果(如果你搞不定,
那我对你的智商表示怀疑)。
OK,开始
假设
,这个假设不影响最终结论,因为你总可以把一个区间线性
的变换到长度为 1 的区间。
假设压力 p 在控制体 i 内部的分布是一个二次多项式
,控制体 i 的中心处于
处,左右两个界面就是
和 。
这里先强调一个问题,在 FVM 中,每个控制体内的求解出来的变量 实
际上是这个控制体内的平均值
。
所以,
。
我们知道
, 和
,等距网格情况下
和
处的导数可以近
似表示为
那么
由上述三个有关 a,b 和 c 的方程,我们可以得到
这样就可以得到 f(x)的表达式了,由此可以算出
和
通常 MUSCL 格式写成如下形式
对应我们的推导结果(二次多项式假设)。
但是这不是最终形式。如果直接用这个公式,就会导致流场在激波(间断)
附近的振荡。因为直接用二次多项式去逼近一个间断,会导致这样的效果。
所以科学家们提出要对间断附近的斜率有所限制,因此引入了一个非常重
要的修改——斜率限制器。加入斜率限制器后,上述公式就有点变化。
这里 是 Van Albada 限制器
是一个小数(
),以防止分母为 0。
密度和速度通过同样的方法来搞定。
密度、速度和压力被称作原始变量,所以上述方法是基于原始变量的
MUSCL。此外还有基于特征变量的 MUSCL,要复杂一点,但是被认为适
合更高精度的格式。然而一般计算中,基于原始变量的 MUSCL 由于具有
足够的精度、简单的形式和较低的代价而被广泛使用。
OK,搞定了。下面进入第二点,怎么求
。关于这一点,我不打算做
详细介绍了,直接使用现有的近似黎曼解就可以了,都是通过
和
计算得到
。比如 Roe 因为形式简单,而非常流行。在 CFL3D 软件主
页(http://cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html)上看手册,
附录 C 的 C.1.3。
想了一下,还是把 Roe 求解器稍微说说吧,力求比较完整。但是不要指望
我把 Roe 求解器解释清楚,因为这个不是很容易三言两语说清的。
Roe 求解器的数学形式是这样的
显然这个公式的第一项是一个中心差分形式,先前说过简单的中心差分不
可行,原因是耗散不足导致振荡,说得通俗点就像一个弹簧,如果缺乏耗