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

借宝地写几个小短文，介绍 CFD 的一些实际的入门知识。主要是因为这 里支持 Latex，写起来比较方便。 CFD，计算流体力学，是一个挺难的学科，涉及流体力学、数值分析和计 算机算法，还有计算机图形学的一些知识。尤其是有关偏微分方程数值分 析的东西，不是那么容易入门。大多数图书，片中数学原理而不重实际动 手，因为作者都把读者当做已经掌握基础知识的科班学生了。所以数学基 础不那么好的读者往往看得很吃力，看了还不知道怎么实现。本人当年虽 说是学航天工程的，但是那时本科教育已经退步，基础的流体力学课被砍 得只剩下一维气体动力学了，因此自学 CFD 的时候也是头晕眼花。不知 道怎么实现，也很难找到教学代码——那时候网络还不发达，只在教研室 的故纸堆里搜罗到一些完全没有注释，编程风格也不好的冗长代码，硬着 头皮分析。后来网上淘到一些代码研读，结合书籍论文才慢慢入门。可以 说中间没有老师教，后来赌博士为了混学分上过 CFD 专门课程，不过那 时候我已经都掌握课堂上那些了。 回想自己入门艰辛，不免有一个想法——写点通俗易懂的 CFD 入门短文 给师弟师妹们。本人不打算搞得很系统，而是希望能结合实际，阐明一些 最基本的概念和手段，其中一些复杂的道理只是点到为止。目前也没有具 体的计划， 想到哪里写到哪里， 因此可能会很零散。 但是我争取让初学 CFD 的人能够了解一些基本的东西，看过之后，会知道一个 CFD 代码怎么炼 成的（这“炼”字好像很流行啊）。欢迎大家提出意见，这样我尽可能的可欢迎大家提出意见，这样我尽可能的可 以追加一些修改和解释。以追加一些修改和解释。 言归正传，第一部分，我打算介绍一个最基本的算例，一维激波管问题。 说白了就是一根两端封闭的管子，中间有个隔板，隔板左边和右边的气体 状态（密度、速度、压力）不一样，突然把隔板抽去，管子内面的气体怎 么运动。这是个一维问题，被称作黎曼间断问题，好像是黎曼最初研究双 曲微分方程的时候提出的一个问题，用一维无粘可压缩 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）上看手册， 里面空间离散那