x

学过弹性力学的人应该都知道什么是有限元，而对学计算流体力学的来说，有限差分和有限体积法也是两种非常重要的方法。三者虽然目前形式各异，但是思想上有很多类似的地方。CFD（Computational Fluid Dynamics）中主要的三种离散方法就是他们三个。

而这篇文章主要目的是对三者进行比较，并给出三种方法计算同一个流体一维算例的过程。

一维算例：

http://s14/middle/7a1c18a84dad74419c50d&690

流体在沿着一条流线传输过程中，关于流体的某个物理量http://s2/middle/7a1c18a84dad744394191&690的输运方程为：

http://s16/middle/7a1c18a84dad744286b5f&690

这个方程也可以看成对流-扩散方程的特例，其含义就是在控制体内部：

由源产生的http://s15/middle/7a1c18a84dad7442c1dee&690的量（因为全场没有源项，故等号右端为0） =

由于流动导致的物理量http://s2/middle/7a1c18a84dad744394191&690的增加（）+由于扩散导致物理量的减少（）

另外流场还应该满足连续条件：

    http://s14/middle/7a1c18a84dad74440315d&690     

即：流入质量=流出质量

那么控制方程组就是：

http://s6/middle/7a1c18a84dad744381385&690

下面，用三种方法分别来解这个流场的物理量http://s5/middle/7a1c18a84dad7443c5a04&690值

1. 有限差分法

主要参考《计算流体力学入门》，John D.Andersion，JR .

有限差分法的做法很简单：流体的控制方程组一般都是偏微分方程，解偏微分方程组很困难， N-S方程多数都是没有解析解的，但是解线性方程组是很简单的事情。所以，差分主要目的便是将微分近似为线性运算，进而把偏微分方程组转化为线性方程组。大多数用来代替偏微分形式是基于Taylor展开得到的，利用Taylor展开，对x一阶导数差分可以近似如下三种形式：

    http://s7/middle/7a1c18a84dad744475de6&690 http://s10/middle/7a1c18a84dad7444b7039&690     

对扩散项的二阶导数，近似为：

http://s11/middle/7a1c18a84dad7444c4c3a&690

将控制方程离散成差分方程，得到：

http://s7/middle/7a1c18a84dad7444703c6&690

http://s4/middle/7a1c18a84dad7444a7473&690

这里我们把http://s2/middle/7a1c18a84dad7445c0cf1&690，然后分别一阶和二阶格式进行差分。

通过联立连续方程，可以求解出全场所有节点处流速值http://s10/middle/7a1c18a84dad74447e679&690，代入对流扩散方程中，简化为：

http://s5/middle/7a1c18a84dad744563334&690

然后联立求解全场各个点有关http://s7/middle/7a1c18a84dad74456c716&690全场节点处数值。

1. 有限体积法

主要参考《计算流体动力学分析：CFD软件原理与应用》，王福军

（其实里面也就是流体力学控制方程简介+有限体积法，整个这个名字真让人迷惑QAQ）

1.1离散

有限体积法离散的核心和有限元法一样，使用有限个离散点来代替原来整个连续的空间。把计算区域分成不重叠的计算网格，然后确定每个节点位置和节点控制体体积（也就是节点所在的网格单元）。区域几何要素主要有以下几个：

节点：需要解未知物理量的几何位置，一般在节点上定义所有的标量，下面图中的W、P、E三个点就是节点；

控制体积：应用控制方程和守恒定律的最小单位，图中灰色区域；

界面：规定了与各节点相对应的控制体积的分界面位置，图中w和e处；

网格线：连接节点的曲线簇，图中就是x轴线。

http://s8/middle/7a1c18a84dad744574997&690

1.2建立离散方程

有限体积最特殊的一步就是在控制体上积分控制方程，以便在控制体积的界面上产生离散方程。对一维算例在控制体积http://s4/middle/7a1c18a84dad7446df053&690上积分有：

http://s14/middle/7a1c18a84dad74457579d&690

然后控制体积分式等价于：

http://s12/middle/7a1c18a84dad744693eeb&690

其中http://s12/middle/7a1c18a84dad7446644eb&690 可以看成是垂直于http://s1/middle/7a1c18a84dad74465dd20&690 轴方向上断面的面积。

这里我想详细的讨论一下，在书中作者说这是完全等价的转化，并没有取任何近似，但是我有不同的看法

http://s8/middle/7a1c18a84dad74465ddb7&690

http://s1/middle/7a1c18a84dad7446abd20&690 积分式里面的被积函数http://s13/middle/7a1c18a84dad74479424c&690做积分：

http://s10/middle/7a1c18a84dad7447a9229&690

然后根据分部积分方法：

http://s15/middle/7a1c18a84dad74474e66e&690

我们可以看出，等号左端正确结果应为：

http://s13/middle/7a1c18a84dad74474bc5c&690

所以这一步其实也做了近似。否则应该加上http://s13/middle/7a1c18a84dad74475122c&690的前提。

1.3解方程

我们已经把非线性的控制方程组通过离散，转化成为了线性方程组。

http://s10/middle/7a1c18a84dad744759ff9&690

http://s5/middle/7a1c18a84dad744839b34&690

但是在我们的解里http://s11/middle/7a1c18a84dad7448c737a&690在界面处的值用节点上的值来近似，这里就是所谓的使用差分格式。

采用最煎蛋XD的方法……我们用界面两边的节点值取平均来近似，这应该是地球人都想得到的方法，方程继续化简为：

http://s4/middle/7a1c18a84dad744847b33&690

记为：

http://s1/middle/7a1c18a84dad7448535d0&690

根据连续方程，我们可以求解整个流场，那么系数http://s16/middle/7a1c18a84dad74493b4ff&690。

http://s8/middle/7a1c18a84dad744a31fc7&690值。

1. 有限单元法（暂略）