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

借宝地写几个小短文，介绍CFD的一些实际的入门知识。主要是因为这 里支持Latex，写起来比较方便。CFD，计算流体力学，是一个挺难的学科，涉及流体力学、数值分析和计 算机算法，还有计算机图形学的一些知识。尤其是有关偏微分方程数值分 析的东西，不是那么容易入门。大多数图书，片中数学原理而不重实际动 手，因为作者都把读者当做已经掌握基础知识的科班学生了。所以数学基 础不那么好的读者往往看得很吃力，看了还不知道怎么实现。本人当年虽 说是学航天工程的，但是那时本科教育已经退步，基础的流体力学课被砍 得只剩下一维气体动力学了，因此自学CFD的时候也是头晕眼花。不知 道怎么实现，也很难找到教学代码一一那时候网络还不发达，只在教研室 的故纸堆里搜罗到一些完全没有注释，编程风格也不好的冗长代码，硬着 头皮分析。后来网上淘到一些代码研读，结合书籍论文才慢慢入门。可以 说中间没有老师教，后来赌博士为了混学分上过CFD专门课程，不过那 时候我已经都掌握课堂上那些了。回想自己入门艰辛，不免有一个想法一一写点通俗易懂的CFD入门短文 给师弟师妹们。本人不打算搞得很系统，而是希望能结合实际，阐明一些 最基本的概念和手段，其中一些复杂的道理只是点到为止。目前也没有具 体的计划，想到哪里写到哪里，因此可能会很零散。但是我争取让初学CFD 的人能够了解一些基本的东西，看过之后，会知道一个CFD代码怎么炼 成的（这炼字好像很流行啊）。欢迎大家提出意见，这样我尽可能的可 以追加一些修改和解释。言归正传，第一部分，我打算介绍一个最基本的算例，一维激波管问题。说白了就是一根两端封闭的管子，中间有个隔板，隔板左边和右边的气体 状态（密度、速度、压力）不一样，突然把隔板抽去，管子内面的气体怎 么运动。这是个一维问题，被称作黎曼间断问题，好像是黎曼最初研究双 曲微分方程的时候提出的一个问题，用一维无粘可压缩Euler方程就可以 描述了。这里E = U2 + EP =（】一 1）m这个方程就是描述的气体密度,、动量*和能量随时间的变化（#）与它们各自的流量（密度流量 * 动量流量，能量流量:）随空间变化（土）的关系。在CFD中通常把这个方程写成矢量形式寄+如。这里Q = pu_ PEpuF =puu + ppEu + pu进一步可以写成散度形式器+ F = 一定要熟悉这种矢量形式以上是控制方程，下面说说求解思路。可压缩流动计算中，有限体积（FVM） 是最广泛使用的算法，其他算法多多少少都和FVM有些联系或者共通的 思路。了解的FVM，学习其他高级点的算法（比如目前比较热门的间断有 限元、谱FVM、谱FDM）就好说点了。针对一个微元控制体口可，把Euler方程在空间积分器血+ J：VF心=。用微积分知识可以得到蠹+曾=0也就是说控制体内气体状态平均值的变化是控制体界面上流通量的结果。 因此我们要计算Q的演化，关键问题是计算控制体界面上的F。FVM就是以这个积分关系式出发，把整个流场划分为许多小控制体，每个 控制体和周围相邻的某个控制体共享一个界面，通过计算每个界面上的通 量来得到相邻控制体之间的影响，一旦每个控制体的变化得到，整个流场 的变化也就知道了。所以，再强调一次，关键问题是计算控制体界面上的F。初学者会说，这个不难，把界面：上的以插值得到，然后就可以计算F：有道理！咱们画个图，有三个小控制体i-1到i+1，中间的“|表示界面，控制体i 右边的界面用；+普表示，左边的就是土I i-1 I i I i+1 I好下个问题：每个小控制体长度都是一；如何插值计算界面+普上的 以七？最自然的想法就是：取两边的平均值呗，Ql+i =捋也但是很不幸，这是不行的。那么换个方法？直接平均得到虬+;？还是很不行，这样也不行。我靠，这是为什么？这明明是符合微积分里面的知识啊？这个道理有点复杂，说开了去可以讲一本书，可以说从50年代到70年 代，CFD科学家就在琢磨这个问题。这里，初学者只需要记住这个结论： 对于流动问题，不可以这样简单取平均值来插值或者差分如果你非要想 知道这个究竟，我现在也不想给你讲清楚，因为我眼下的目的是让你快速 上手，而且该不刨根问底的时候就不要刨根问底，这也是初学阶段一种重 要的学习方法。 好了，既然目的只是为了求虬+；,我在这里，只告诉你一种计算方法，也 是非常重要、非常流行的一种方法。简单的说，就是假设流动状态在界面 +普是不连续的，先计算出界面；+普两边Q的值，w十:和二；，再由 它们用某种方法计算出F十；。上述方法是非常重要的，是由一个苏联人 Godunov在50年代首创的，后来被发展成为通用Godunov方法，著名的ENO/WENO就是其中的一种。好了，现在的问题是：1怎么确定口基和馄二；2怎么计算FL十；对于第一个问题，Godunov在他的论文中，是假设每个控制体中Q是均 匀分布的，因此W十:一。切=。中第二个问题，Godunov采用了精确黎曼解来计算F*十；。什么是精确黎 曼解，就是计算这个激波管问题的精确解。既然有精确解，那还费功夫搞 这些FVM算法干什么？因为只有这种简单一维问题有精确解，稍微复杂 一点就不行了。精确解也比较麻烦，要分析5种情况，用牛顿法迭代求解 （牛顿法是什么？看数值计算的书去，哦，算了，现在暂时可以不必看）。这是最初Godunov的方法，后来在这个思想的基础上，各种变体都出来 了。也不过是在这两个问题上做文章，怎么确定*4十1,怎么计算虬十；。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怎么计算FL十；这里我们给出第一点一个具体的实现方法，就是基于原始变量的MUSCL 格式(以下简称MUSCL)。它是一种很简单的格式，而且具有足够的精 度，NASA著名的CFL3D软件就是使用了这个格式，大家可以去它的主页(http:/cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html)上看手册， 里面空间离散那一章清楚的写着。MUSCL假设控制体内原始变量(就是7)的分布是一次或者二次多项式， 如果得到了这个多项式，就可以求出控制体；左右两个界面的一侧的值 Q*和心。我们以压力p为例来说明怎么构造这个多项式。这里我只针对二次多项式 来讲解，你看完之后肯定能自己推导出一次多项式的结果(如果你搞不定， 那我对你的智商表示怀疑)。OK，开始假设一； =1，这个假设不影响最终结论，因为你总可以把一个区间线性 的变换到长度为1的区间。假设压力p在控制体i内部的分布是一个二次多项式力=。+ 、： + .，控制体i的中心处于.；、=。，处，左右两个界面就是 Y和+ 土这里先强调一个问题，在FVM中，每个控制体内的求解出来的变量7实 际上是这个控制体内的平均值以| = 士二以心。所以，Pi = (ajc2 + 赤 + c)djc=即+ b按+ cr二=。我们知道扑i,打和执十L,等距网格情况下；一普和* + i处的导数可以近似表示为票|_广=邑邑T如汁广 2中-瓦那么户（号）=（2那+炒|至=_直= = 尸（+号）=（2愈，+刖王=十* = b +。= 十 由上述三个有关a，b和c的方程，我们可以得到A+-A-6 =-中这样就可以得到f（x）的表达式了，由此可以算出f十：和二-； 通常MUSCL格式写成如下形式口;；言=虱 + *）+（1+ *）打二直=虱一H（i 一 +（1 + 、= 1 ：；对应我们的推导结果（二次多项式假设）。但是这不是最终形式。如果直接用这个公式，就会导致流场在激波（间断） 附近的振荡。因为直接用二次多项式去逼近一个间断，会导致这样的效果。 所以科学家们提出要对间断附近的斜率有所限制，因此引入了一个非常重 要的修改一一斜率限制器。加入斜率限制器后，上述公式就有点变化。pi =瓦 + 堂（1 屈）厂 +（1 + 屈,）,=瓦堂（1屈）+ （1 +屈,）这里、是Van Albada限制器、= 心侦十亡T 。尸+。）牛是一个小数（= 11L），以防止分母为0。密度和速度通过同样的方法来搞定。密度、速度和压力被称作原始变量，所以上述方法是基于原始变量的 MUSCL。此外还有基于特征变量的MUSCL，要复杂一点，但是被认为适 合更高精度的格式。然而一般计算中，基于原始变量的MUSCL由于具有 足够的精度、简单的形式和较低的代价而被广泛使用。OK，搞定了。下面进入第二点，怎么求F十；。关于这一点，我不打算做 详细介绍了，直接使用现有的近似黎曼解就可以了，都是通过7 r ：和7二； 计算得到F+；。比如Roe因为形式简单，而非常流行。在CFL3D软件主 页（http:/cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html）上看手册， 附录C的C.1.3。想了一下，还是把Roe求解器稍微说说吧，力求比较完整。但是不要指望 我把Roe求解器解释清楚，因为这个不是很容易三言两语说清的。Roe求解器的数学形式是这样的=部（0） + F（Q鬲）-据蓦（。切-。疆）显然这个公式的第一项是一个中心差分形式，先前说过简单的中心差分不 可行，原因是耗散不足导致振荡，说得通俗点就像一个弹簧，如果缺乏耗 散（阻尼）它就会一直振荡。耗散这个术语在激波捕捉格式中是最常见 的。第二项的作用就是提供足够的耗散了。这里口二；和十：已经用MUSCL求得了，F的定义在第一讲中已经介绍 了。只有是还没说过的。A 丝牡g 这个矩阵可以写成特征矩阵和特征向量矩阵的形式A = R-AR+而|A| = R-|A|R+A，R，R一和A的具体表达式在许多书上都有，而且这里的矩阵表达有 问题，所以就不写了。十;是由厂十；、七十；和十；代入a计算得到。而厂十；、七十；和十；采 用所谓Roe平均值这才是Roe求解器关键的地方！总结一下，就是用Roe平均计算界面上的气体状态匕十，然后计算得到 俱讦，这样咯+就可以得到了。如果有时间，我后面会找一个代码逐 句分析一下。总之，计算F+；还是很不直接的。构造近似黎曼解是挺有学问的，需要对 气体动力学的物理和数学方面有较深的理解。通常，如果不是做基础研究， 你只需要知道它们的特点，会用它们就可以了，而不必深究它们怎么推导 出来的。