计算流体力学过渡到编程的傻瓜入门教程.docx

1、计算流体力学过渡到编程的傻瓜入门教程借宝地写几个小短文，介绍 CFD 的一些实际的入门知识。主要是因为这 里支持 Latex ，写起来比较方便。CFD ，计算流体力学，是一个挺难的学科，涉及流体力学、数值分析和计 算机算法，还有计算机图形学的一些知识。尤其是有关偏微分方程数值分 析的东西，不是那么容易入门。大多数图书，片中数学原理而不重实际动 手，因为作者都把读者当做已经掌握基础知识的科班学生了。所以数学基 础不那么好的读者往往看得很吃力，看了还不知道怎么实现。本人当年虽 说是学航天工程的，但是那时本科教育已经退步，基础的流体力学课被砍 得只剩下一维气体动力学了，因此自学 CFD 的时候也是头

2、晕眼花。不知 道怎么实现，也很难找到教学代码 那时候网络还不发达，只在教研室 的故纸堆里搜罗到一些完全没有注释，编程风格也不好的冗长代码，硬着 头皮分析。后来网上淘到一些代码研读，结合书籍论文才慢慢入门。可以 说中间没有老师教，后来赌博士为了混学分上过 CFD 专门课程，不过那 时候我已经都掌握课堂上那些了。回想自己入门艰辛，不免有一个想法 写点通俗易懂的 CFD 入门短文 给师弟师妹们。本人不打算搞得很系统，而是希望能结合实际，阐明一些 最基本的概念和手段，其中一些复杂的道理只是点到为止。目前也没有具 体的计划，想到哪里写到哪里， 因此可能会很零散。 但是我争取让初学 CFD 的人能够了解一

3、些基本的东西，看过之后，会知道一个 CFD 代码怎么炼 成的（这 “炼”字好像很流行啊）。 欢迎大家提出意见，这样我尽可能的可 以追加一些修改和解释。言归正传，第一部分，我打算介绍一个最基本的算例，一维激波管问题。说白了就是一根两端封闭的管子，中间有个隔板，隔板左边和右边的气体 状态（密度、速度、压力）不一样，突然把隔板抽去，管子内面的气体怎 么运动。这是个一维问题，被称作黎曼间断问题，好像是黎曼最初研究双曲微分方程的时候提出的一个问题，用一维无粘可压缩 Euler方程就可以描述了。这里这个方程就是描述的气体密度、动量和能量糾禹随时间的变化（-与它们各自的流量（密度流量 ，动量流量琴阿丨呼，能

4、量流量炽：.帶）随空间变化（一】）的关系。在CFD中通常把这个方程写成矢量形式这里pu pEpn puu + p pEu 十 pu进一步可以写成散度形式+ V F=0一定要熟悉这种矢量形式以上是控制方程，下面说说求解思路。可压缩流动计算中，有限体积（FVM） 是最广泛使用的算法，其他算法多多少少都和 FVM有些联系或者共通的思路。了解的FVM，学习其他高级点的算法（比如目前比较热门的间断有 限元、谱FVM、谱FDM ）就好说点了。针对一个微元控制体弧引，把Euler方程在空间积分薯必+上VFit二0用微积分知识可以得到西.I叭一阳-n3F 十 Ar - U也就是说控制体内气体状态平均值的变化是

5、控制体界面上流通量的结果。因此我们要计算的演化，关键问题是计算控制体界面上的 1 0FVM就是以这个积分关系式出发，把整个流场划分为许多小控制体，每个控制体和周围相邻的某个控制体共享一个界面，通过计算每个界面上的通 量来得到相邻控制体之间的影响，一旦每个控制体的变化得到，整个流场 的变化也就知道了。所以，再强调一次，关键问题是计算控制体界面上的 r 初学者会说，这个不难，把界面上的-插值得到，然后就可以计算| 有道理！咱们画个图，有三个小控制体 i-1到i+1 ,中间的“|表示界面，控制体i右边的界面用7表示，左边的就是I oI i-1 I i I i+1 I好下个问题：每个小控制体长度都是

6、一、如何插值计算界面 ：上的最自然的想法就是：取两边的平均值呗，但是很不幸，这是不行的。那么换个方法？直接平均得到 r - ?口 _ F卄珥+二耳+* - 3还是很不行，这样也不行。我靠，这是为什么？这明明是符合微积分里面的知识啊？这个道理有点复杂，说幵了去可以讲一本书，可以说从 50年代到70年代，CFD科学家就在琢磨这个问题。这里，初学者只需要记住这个结论： 对于流动问题，不可以这样简单取平均值来插值或者差分。如果你非要想 知道这个究竟，我现在也不想给你讲清楚，因为我眼下的目的是让你快速 上手，而且该不刨根问底的时候就不要刨根问底，这也是初学阶段一种重 要的学习方法。好了，既然目的只是为了

7、求 F ，我在这里，只告诉你一种计算方法，也 是非常重要、非常流行的一种方法。简单的说，就是假设流动状态在界面 :*是不连续的，先计算出界面=寺两边的值，1 !和7 :，再由它们用某种方法计算出一-。上述方法是非常重要的，是由一个苏联人 Godunov在50年代首创的，后来被发展成为通用 Godunov方法，著名 的ENO/WENO 就是其中的一种。好了，现在的问题是：Q 1 Q+ 11怎么确定 ：和2怎么计算P -对于第一个问题，Godunov在他的论文中，是假设每个控制体中！是均 匀分布的，因此!% = Q+i第二个问题，Godunov采用了精确黎曼解来计算尸:。什么是精确黎 曼解”就是计

8、算这个激波管问题的精确解。 既然有精确解，那还费功夫搞这些FVM算法干什么？因为只有这种简单一维问题有精确解，稍微复杂 一点就不行了。精确解也比较麻烦，要分析 5种情况，用牛顿法迭代求解（牛顿法是什么？看数值计算的书去， 哦，算了，现在暂时可以不必看）。这是最初Godunov的方法，后来在这个思想的基础上，各种变体都出来Q-峠 R了。也不过是在这两个问题上做文章，怎么确定 怎么计算 I :。Godunov假设的是每个小控制体内是均匀分布，也就是所谓分段常数(piecewise con sta nt) ，所以后来有分段线性 (picewise lin ear) 或者分段二次分布(picewise

9、 parabolic) ，到后来ENO/WENO 出来，那这个假设的多项式次数就继续往上走了。都是用多项式近似的，这是数值计算中 的一个强大工具，你可以在很多地方看到这种近似。Godu nov用的四精确黎曼解，太复杂太慢，也不必要，所以后来就有各 种近似黎曼解，最有名的是 Roe求解器、HLL求解器和Osher求解器， 都是对精确黎曼解做的简化。这个多项式的阶数是和计算精度密切相关的，阶数越高，误差就越小。不 过一般来说，分段线性就能得到不错的结果了，所以工程中都是用这个，Flue nt、Fastra n 以及 NASA 的 CFL3D、OverFlow 都是用这个。而黎曼求解器对精度的影响不

10、是那么大，但是对整个算法的物理适用性有影 响，也就是说某种近似黎曼求解器可能对某些流动问题不合适，比如单纯 的Roe对于钝头体的脱体激波会算得乱七八糟，后来加了熵修正才算搞定。上次()说到了求解可压缩流动的一个重要算法，通用 Godunov方法。其两个主要步骤就是Q 1 Q+ h1怎么确定 ：和：F2怎么计算这里我们给出第一点一个具体的实现方法，就是基于原始变量的 MUSCL格式(以下简称 MUSCL )。它是一种很简单的格式，而且具有足够的精度，NASA著名的CFL3D软件就是使用了这个格式，大家可以去它的主 页（）上看手册，里面空间离散那一章清楚的写着。MUSCL假设控制体内原始变量（就是

11、）的分布是一次或者二次多项式， 如果得到了这个多项式，就可以求出控制体 左右两个界面的一侧的值 ：和戏！。我们以压力p为例来说明怎么构造这个多项式。这里我只针对二次多项式 来讲解，你看完之后肯定能自己推导出一次多项式的结果 （如果你搞不定, 那我对你的智商表示怀疑）。OK，幵始假设一 1 ,这个假设不影响最终结论，因为你总可以把一个区间线性 的变换到长度为1的区间。假设压力p在控制体i内部的分布是一个二次多项式+ c，控制体i的中心处于e 处，左右两个界面就是:和-0 这里先强调一个问题，在 FVM中，每个控制体内的求解出来的变量际上是这个控制体内的平均值 所以，Pi = J*?（加十吐十二

12、ax3/3 + &72+ 111/2我们知道和,等距网格情况下；=和2处的导数可以近似表示为那么尸(-1) = (2ax2 + &) |I=_j - b - a - A-尸(+壬)=(加声+圳#十斗=b -p a = 十由上述三个有关a，b和c的方程，我们可以得到a 二 2L_ _ 2- &t-A-* 決这样就可以得到f（X）的表达式了，由此可以算出和 通常MUSCL格式写成如下形式Pi二注+扑仕一上）+ +町訂 % =注一扑（1上）F + U +町訂 fe = L/3对应我们的推导结果（二次多项式假设）但是这不是最终形式。如果直接用这个公式，就会导致流场在激波（间断） 附近的振荡。因为直接用

13、二次多项式去逼近一个间断， 会导致这样的效果。所以科学家们提出要对间断附近的斜率有所限制，因此引入了一个非常重 要的修改一一斜率限制器。加入斜率限制器后，上述公式就有点变化。P； = A + 亍肛瓯）+ （14 fesJAt= +（L + fe）A7 2 弋这里：；是Van Albada 限制器耳二远予T否两:是一个小数（ 一 - 丄:），以防止分母为o。密度和速度通过同样的方法来搞定。密度、速度和压力被称作原始变量，所以上述方法是基于原始变量的MUSCL。此外还有基于特征变量的 MUSCL，要复杂一点，但是被认为适 合更高精度的格式。然而一般计算中，基于原始变量的 MUSCL由于具有足够的精

14、度、简单的形式和较低的代价而被广泛使用。OK，搞定了。下面进入第二点，怎么求 F 0关于这一点，我不打算做 详细介绍了，直接使用现有的近似黎曼解就可以了， 都是通过Z :和咗:计算得到F 。比如Roe因为形式简单，而非常流行。在CFL3D软件主 页（）上看手册，附录 C的C.1.3 o想了一下，还是把Roe求解器稍微说说吧，力求比较完整。但是不要指望 我把Roe求解器解释清楚，因为这个不是很容易三言两语说清的。Roe求解器的数学形式是这样的E出=河皿）+巩富）-枫+*|（疇* - Q希）显然这个公式的第一项是一个中心差分形式，先前说过简单的中心差分不可行，原因是耗散不足导致振荡，说得通俗点就像

15、一个弹簧，如果缺乏耗散（阻尼）它就会一直振荡。 耗散”这个术语在激波捕捉格式中是最常见的。第二项的作用就是提供足够的耗散了。这里咗和，已经用MUSCL求得了，1的定义在第一讲中已经介绍 了。只有是还没说过的。 9Q这个矩阵可以写成特征矩阵和特征向量矩阵的形式A = R- AR+而|A| = R-|A|R+A, 1：, 和丄的具体表达式在许多书上都有，而且这里的矩阵表达有问题，所以就不写了。N 是由、 -和 计*代入A计算得到。而阿+寺、叫44和皂计扌采用所谓Roe平均值这才是Roe求解器关键的地方!总结一下，就是用 Roe平均计算界面上的气体状态 门1 -，然后计算得到 计刑，这样尺叶寺就可以得到了。如果有时间，我后面会找一个代码逐句分析一下。总之，计算F还是很不直接的。构造近似黎曼解是挺有学问的，需要对气体动力学的物理和数学方面有较深的理解。 通常，如果不是做基础研究,你只需要知道它们的特点，会用它们就可以了，而不必深究它们怎么推导 出来的。