借宝地写几个小短文，介绍CFD的一些实际的入门知识。主要是因为这里支持Latex，写起来比较方便。

CFD，计算流体力学，是一个挺难的学科，涉及流体力学、数值分析和计算机算法，还有计算机图形学的一些知识。尤其是有关偏微分方程数值分析的东西，不是那么容易入门。大多数图书，片中数学原理而不重实际动手，因为作者都把读者当做已经掌握基础知识的科班学生了。所以数学基础不那么好的读者往往看得很吃力，看了还不知道怎么实现。本人当年虽说是学航天工程的，但是那时本科教育已经退步，基础的流体力学课被砍得只剩下一维气体动力学了，因此自学CFD的时候也是头晕眼花。不知道怎么实现，也很难找到教学代码——那时候网络还不发达，只在教研室的故纸堆里搜罗到一些完全没有注释，编程风格也不好的冗长代码，硬着头皮分析。后来网上淘到一些代码研读，结合书籍论文才慢慢入门。可以说中间没有老师教，后来赌博士为了混学分上过CFD专门课程，不过那时候我已经都掌握课堂上那些了。

回想自己入门艰辛，不免有一个想法——写点通俗易懂的CFD入门短文给师弟师妹们。本人不打算搞得很系统，而是希望能结合实际，阐明一些最基本的概念和手段，其中一些复杂的道理只是点到为止。目前也没有具体的计划，想到哪里写到哪里，因此可能会很零散。但是我争取让初学CFD的人能够了解一些基本的东西，看过之后，会知道一个CFD代码怎么炼成的（这“炼”字好像很流行啊）。欢迎大家提出意见，这样我尽可能的可以追加一些修改和解释。

言归正传，第一部分，我打算介绍一个最基本的算例，一维激波管问题。说白了就是一根两端封闭的管子，中间有个隔板，隔板左边和右边的气体状态（密度、速度、压力）不一样，突然把隔板抽去，管子内面的气体怎么运动。这是个一维问题，被称作黎曼间断问题，好像是黎曼最初研究双曲微分方程的时候提出的一个问题，用一维无粘可压缩Euler方程就可以描述了。
[tex]\begin{equation*}\left\{\begin{array}{c}\frac{\partial \rho}{\partial t} + \frac{\partial (\rho u)}{\partial x}=0\\\frac{\partial \rho u}{\partial t} + \frac{\partial (\rho uu +p)}{\partial x}=0\\\frac{\partial \rho E}{\partial t} +\frac{\partial (\rho E u +pu)}{\partial x}=0\end{array}\right.\end{equation}[/tex]
这里
[tex]E=\frac{1}{2}u^2+e[/tex]
[tex]p=(\gamma-1)\rho e[/tex]
这个方程就是描述的气体密度[tex]\rho[/tex]、动量[tex]\rho u[/tex]和能量[tex]\rho E[/tex]随时间的变化（[tex]\frac{\partial}{\partial t}[/tex]）与它们各自的流量（密度流量[tex]\rho u[/tex]，动量流量[tex]\rho uu+p[/tex]，能量流量[tex]\rho Eu+pu[/tex]）随空间变化（[tex]\frac{\partial}{\partial x}[/tex]）的关系。

在CFD中通常把这个方程写成矢量形式
[tex]\begin{equation*}\frac{\partial \mathbf{Q}}{\partial t} + \frac{\partial \mathbf{F}}{\partial x}=0 \end{equation}[/tex]
这里
[tex]\mathbf{Q}=\left[\begin{array}{c} \rho \\ \rho u \\ \rho E \\\end{array}\right][/tex]
[tex]\mathbf{F}=\left[\begin{array}{c} \rho u \\ \rho uu+p \\ \rho E u +pu \\\end{array}\right][/tex]
进一步可以写成散度形式
[tex]\begin{equation*}\frac{\partial \mathbf{Q}}{\partial t} + \nabla\cdot{ \mathbf{F}}=0 \end{equation}[/tex]
一定要熟悉这种矢量形式

以上是控制方程，下面说说求解思路。可压缩流动计算中，有限体积（FVM）是最广泛使用的算法，其他算法多多少少都和FVM有些联系或者共通的思路。了解的FVM，学习其他高级点的算法（比如目前比较热门的间断有限元、谱FVM、谱FDM）就好说点了。

针对一个微元控制体[tex][a,b][/tex]，把Euler方程在空间积分
[tex]\int_a^b \frac{\partial \mathbf{Q}}{\partial t}dx+\int_a^b \nabla\cdot\mathbf{F}dx =0[/tex]
用微积分知识可以得到
[tex]\frac{\partial \mathbf{\bar{Q}}}{\partial t}+ \frac{\mathbf{F}_b – \mathbf{F}_a}{\triangle x} =0[/tex]
也就是说控制体内气体状态平均值的变化是控制体界面上流通量的结果。因此我们要计算[tex]\mathbf{Q}[/tex]的演化，关键问题是计算控制体界面上的[tex]\mathbf{F}[/tex]。

FVM就是以这个积分关系式出发，把整个流场划分为许多小控制体，每个控制体和周围相邻的某个控制体共享一个界面，通过计算每个界面上的通量来得到相邻控制体之间的影响，一旦每个控制体的变化得到，整个流场的变化也就知道了。

所以，再强调一次，关键问题是计算控制体界面上的[tex]\mathbf{F}[/tex]。

初学者会说，这个不难，把界面[tex]f[/tex]上的[tex]\mathbf{Q}_f[/tex]插值得到，然后就可以计算[tex]\mathbf{F}_f[/tex]。有道理！

咱们画个图，有三个小控制体 i-1到i+1，中间的“|”表示界面，控制体i右边的界面用[tex]i+\frac{1}{2}[/tex]表示，左边的就是[tex]i-\frac{1}{2}[/tex]。

| i-1 | i | i+1 |


好下个问题：每个小控制体长度都是[tex]\triangle x[/tex]如何插值计算界面[tex]i+\frac{1}{2}[/tex]上的[tex]\mathbf{Q}_{i+\frac{1}{2}}[/tex]？

最自然的想法就是：取两边的平均值呗，
[tex]\mathbf{Q}_{i+\frac{1}{2}}=\frac{\mathbf{Q}_i+\mathbf{Q}_{i+1}}{2}[/tex]

但是很不幸，这是不行的。

那么换个方法？直接平均得到[tex]\mathbf{F}_{i+\frac{1}{2}}[/tex]？
[tex]\mathbf{F}_{i+\frac{1}{2}}=\frac{\mathbf{F}_i+\mathbf{F}_{i+1}}{2}[/tex]

还是很不行，这样也不行。

我靠，这是为什么？这明明是符合微积分里面的知识啊？

这个道理有点复杂，说开了去可以讲一本书，可以说从50年代到70年代，CFD科学家就在琢磨这个问题。这里，初学者只需要记住这个结论：对于流动问题，不可以这样简单取平均值来插值或者差分。如果你非要想知道这个究竟，我现在也不想给你讲清楚，因为我眼下的目的是让你快速上手，而且该不刨根问底的时候就不要刨根问底，这也是初学阶段一种重要的学习方法。

好了，既然目的只是为了求[tex]\mathbf{F}_{i+\frac{1}{2}}[/tex]，我在这里，只告诉你一种计算方法，也是非常重要、非常流行的一种方法。简单的说，就是假设流动状态在界面[tex]i+\frac{1}{2}[/tex]是不连续的，先计算出界面[tex]i+\frac{1}{2}[/tex]两边[tex]\mathbf{Q}[/tex]的值，[tex]\mathbf{Q}^{-}_{i+\frac{1}{2}}[/tex]和[tex]\mathbf{Q}^{+}_{i+\frac{1}{2}}[/tex]，再由它们用某种方法计算出[tex]\mathbf{F}_{i+\frac{1}{2}}[/tex]。上述方法是非常重要的，是由一个苏联人Godunov在50年代首创的，后来被发展成为通用Godunov方法，著名的ENO/WENO就是其中的一种。

好了，现在的问题是：
1 怎么确定[tex]\mathbf{Q}^{-}_{i+\frac{1}{2}}[/tex]和[tex]\mathbf{Q}^{+}_{i+\frac{1}{2}}[/tex]
2 怎么计算[tex]\mathbf{F}_{i+\frac{1}{2}}[/tex]

对于第一个问题，Godunov在他的论文中，是假设每个控制体中[tex]\mathbf{Q}[/tex]是均匀分布的，因此[tex]\mathbf{Q}^{-}_{i+\frac{1}{2}}=\mathbf{Q}_i[/tex]
[tex]\mathbf{Q}^{+}_{i+\frac{1}{2}}=\mathbf{Q}_{i+1}[/tex]

第二个问题，Godunov采用了精确黎曼解来计算[tex]\mathbf{F}_{i+\frac{1}{2}}[/tex]。什么是“精确黎曼解”，就是计算这个激波管问题的精确解。既然有精确解，那还费功夫搞这些FVM算法干什么？因为只有这种简单一维问题有精确解，稍微复杂一点就不行了。精确解也比较麻烦，要分析5种情况，用牛顿法迭代求解（牛顿法是什么？看数值计算的书去，哦，算了，现在暂时可以不必看）。

这是最初Godunov的方法，后来在这个思想的基础上，各种变体都出来了。也不过是在这两个问题上做文章，怎么确定[tex]\mathbf{Q}^{-,+}_{i+\frac{1}{2}}[/tex]，怎么计算[tex]\mathbf{F}_{i+\frac{1}{2}}[/tex]。

Godunov假设的是每个小控制体内是均匀分布，也就是所谓分段常数(piecewise constant)，所以后来有分段线性(picewise linear)或者分段二次分布(picewise parabolic)，到后来ENO/WENO出来，那这个假设的多项式次数就继续往上走了。都是用多项式近似的，这是数值计算中的一个强大工具，你可以在很多地方看到这种近似。

Godunov用的四精确黎曼解，太复杂太慢，也不必要，所以后来就有各种近似黎曼解，最有名的是Roe求解器、HLL求解器和Osher求解器，都是对精确黎曼解做的简化。

这个多项式的阶数是和计算精度密切相关的，阶数越高，误差就越小。不过一般来说，分段线性就能得到不错的结果了，所以工程中都是用这个，Fluent、Fastran以及NASA的CFL3D、OverFlow都是用这个。而黎曼求解器对精度的影响不是那么大，但是对整个算法的物理适用性有影响，也就是说某种近似黎曼求解器可能对某些流动问题不合适，比如单纯的Roe对于钝头体的脱体激波会算得乱七八糟，后来加了熵修正才算搞定。

[未完，待续]