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

借宝地写几种小短文，简介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假设的是每个小控制体内是均匀分布，也就是所谓分段常数(piecewiseconstant)，所后来来有分段线性(picewiselinear)或者分段二次分布(picewiseparabolic)，到后来ENO/WENO出来，那这个假设的多项式次数就继续往上走了。都是用多项式近似的，这是数值计算中的一种强大工具，你可以在诸多地方看到这种近似。Godunov用的四精确黎曼解，太复杂太慢，也不必要，所后来来就有多种近似黎曼解，最有名的是Roe求解器、HLL求解器和Osher求解器，都是对精确黎曼解做的简化。这个多项式的阶数是和计算精度亲密有关的，阶数越高，误差就越小。不过一般来说，分段线性就能得到不错的成果了，因此工程中都是用这个，Fluent、Fastran以及NASA的CFL3D、OverFlow都是用这个。而黎曼求解器对精度的影响不是那么大，不过对整个算法的物理合用性有影响，也就是说某种近似黎曼求解器也许对某些流动问题不合适，例如单纯的Roe对于钝头体的脱体激波会算得乱七八糟，后来加了熵修正才算搞定。上次（）说到了求解可压缩流动的一种重要算法，通用Godunov措施。其两个重要环节就是1怎么确定和2怎么计算这里我们给出第一点一种详细的实现措施，就是基于原始变量的MUSCL格式（如下简称MUSCL）。它是一种很简朴的格式，并且具有足够的精度，NASA著名的CFL3D软件就是使用了这个格式，大家可以去它的主页（）上看手册，里面空间离散那一章清晰的写着。MUSCL假设控制体内原始变量（就是）的分布是一次或者二次多项式，假如得到了这个多项式，就可以求出控制体左右两个界面的一侧的值和。我们以压力p为例来阐明怎么构造这个多项式。这里我只针对二次多项式来讲解，你看完之后肯定能自己推导出一次多项式的成果（假如你搞不定，那我对你的智商表达怀疑）。OK，开始假设，这个假设不影响最终止论，由于你总可以把一种区间线性的变换到长度为1的区间。假设压力p在控制体i内部的分布是一种二次多项式，控制体i的中心处在处，左右两个界面就是和。这里先强调一种问题，在FVM中，每个控制体内的求解出来的变量实际上是这个控制体内的平均值。因此，。我们懂得，和，等距网格状况下和处的导数可以近似表达为那么由上述三个有关a，b和c的方程，我们可以得到这样就可以得到f(x)的体现式了，由此可以算出和一般MUSCL格式写成如下形式对应我们的推导成果（二次多项式假设）。不过这不是最终形式。假如直接用这个公式，就会导致流场在激波（间断）附近的振荡。由于直接用二次多项式去迫近一种间断，会导致这样的效果。因此科学家们提出要对间断附近的斜率有所限制，因此引入了一种非常重要的修改——斜率限制器。加入斜率限制器后，上述公式就有点变化。这里是VanAlbada限制器是一种小数（），以防止分母为0。密度和速度通过同样的措施来搞定。密度、速度和压力被称作原始变量，因此上述措施是基于原始变量的MUSCL。此外尚有基于特性变量的MUSCL，要复杂一点，不过被认为适合更高精度的格式。然而一般计算中，基于原始变量的MUSCL由于具有足够的精度、简朴的形式和较低的代价而被广泛使用。OK，搞定了。下面进入第二点，怎么求。有关这一点，我不打算做详细简介了，直接使用既有的近似黎曼解就可以了，都是通过和计算得到。例如Roe由于形式简朴，而非常流行。在CFL3D软件主页（）上看手册，附录C的C.1.3。想了一下，还是把Roe求解器稍微说说吧，力争比较完整。不过不要指望我把Roe求解器解释清晰，由于这个不是很轻易三言两语说清的。Roe求解器的数学形式是这样的显然这个公式的第一项是一种中心差分形式，先前说过简朴的中心差分不可行，原因是耗散局限性导致振荡，说得通俗点就像一种弹簧，假如缺乏耗散（阻尼）它就会一直振荡。“耗散”这个术语在激波捕捉格式中是最常见的。第二项的作用就是提供足够的耗散了。这里和已经用MUSCL求得了，的定义在第一讲中已经简介了。只有是还没说过的。这个矩阵可以写成特性矩阵和特性向量矩阵的形式而，，和的详细体现式在许多书上均有，并且这里的矩阵体既有问题，因此就不写了。