计算流体力学

用电子计算机和离散化的数值方法对流体力学问题进行数值模拟和分析的一个新分支。

简史  流体力学和其他学科一样，是通过理论分析和实验研究两种手段发展起来的。很早就已有理论流体力学和实验流体力学两大分支。理论分析是用数学方法求出问题的定量结果。但能用这种方法求出结果的问题毕竟是少数，计算流体力学正是为弥补分析方法的不足而发展起来的。

早在20世纪初，理查德就已提出用数值方法来解流体力学问题的思想。但是由于这种问题本身的复杂性和当时计算工具的落后，这一思想并未引起人们重视。自从40年代中期电子计算机问世以来，用电子计算机进行数值模拟和计算才成为现实。1963年美国的F.H.哈洛和J.E.弗罗姆用当时的IBM7090计算机，成功地解决了二维长方形柱体的绕流问题并给出尾流涡街的形成和演变过程，受到普遍重视。1965年，哈洛和弗罗姆发表“流体动力学的计算机实验”一文，对计算机在流体力学中的巨大作用作了引人注目的介绍。从此，人们把60年代中期看成是计算流体力学兴起的标志。

计算流体力学的历史虽然不长，但已广泛深入到流体力学的各个领域，相应地也形成了各种不同的数值解法。就目前情况看，主要是有限差分方法和有限元法。有限差分方法在流体力学中已得到广泛应用。而有限元法是从求解固体力学问题发展起来的。近年来在处理低速流体问题中，已有相当多的应用，而且还在迅速发展中。

流体力学基本方程  为了说明计算流体力学主要方法，需先了解流体力学运动的基本方程的性质和分类。流体力学的基本方程是在19 世纪上半叶由C.-L.-M.-H.纳维和G.G.斯托克斯等人建立的，称为纳维－斯托克斯方程，简称N-S方程。二维非定常不可压缩流体的N-S方程为：

式中u、v为沿着x、y方向上的速度分量；t为时间；p为压力；ρ为密度；ν为运动粘性系数。在不同条件下，N-S方程的数学性质也不一样。

①N-S方程描述粘性流体随时间而变的非定常运动。时间项和方程右边的高阶导数项决定方程的性质。它同二维热传导方程类似，属于抛物型方程。

②粘性流体的定常运动是将原方程中的时间项省去。此时N-S方程的性质，取决于它的高阶导数项，和拉普拉斯方程一样，为椭圆型方程。

③无粘流的欧拉方程是将 N-S方程的右边粘性项略去而得。它也适用于可压缩流体。从形式上不容易判断欧拉方程的性质。因多数无粘流动皆为无旋流动，故如将欧拉方程改用速度势表示，则二维定常可压缩气流的方程为：

(c-u)-2uv+(c-v)＝0，

式中c为声速。此式是二阶偏微分方程

A+2B+C+D+E＝0

的一般形式，其性质要看B-AC0而定。在超声速区，B-AC＞0，即u+v＞c，上式类似于波动方程，为双曲型；在亚声速区，B-AC＜0，即u+v＜c，上式便与拉普拉斯方程相同，为椭圆型。总之，流体力学的运动方程是极其复杂的非线性偏微分方程，具有各种不同的类型，而且往往还是混合型的。要全面描述流体的运动，还必须同时考虑其他方程，如连续性方程、能量方程和状态方程等。所以计算流体力学在很大程度上就是针对不同性质的偏微分方程采用和发展相应的数值解方法。

低速无粘流动数值解  在无旋条件下，低速流动的速度势满足拉普拉斯方程或泊松方程。很多平面问题利用复变函数和保角映射可以求得解析解，这是经典流体力学的重要内容。但对几何形状比较复杂的物体，必须用下述的数值解法。

①迭代解法 这是用逐步近似求解联立方程的方法，也是椭圆型微分方程的主要数值解法。此法程序简单，存储量与运算量均比较小，一般先假定一组初值，然后求每个网点上的新值。以五点格式为例，网点上的新值是邻近四点初值的平均。新值求出后，旧值还要保留，以便计算其他各点的新值。这种简单迭代收敛很慢，现已很少使用。但若稍加改进，用算出的新值冲掉旧值，并引进一个松弛因子，以加速收敛，将算出来的新值与原来的旧值加权平均，就成为50年代发展起来的逐次超松弛法。

②时间相关法 这是用非定常方程求解定常问题的方法，常用于求解N-S方程和欧拉方程等。虽然用的是非定常方程，但所解的并不是非定常问题。根据给定的初始条件以及随时间改变的约束条件，非定常问题是研究流动随时间的演变过程。这种非定常行为和给出的初值很有关系。然而时间相关法的初值，原则上是随意选取的，只是须满足定常问题所规定的边界条件。在求解过程中，流动随时间的变化并不代表真实的物理过程。当时间足够长后，未知函数值逐步与时间无关，便渐近趋于定常解。所以时间相关法实际上也是一种迭代法，时间变量只不过是用来记录迭代的次数而已。

③交替方向隐式法 流体力学的应用问题，往往是二维和三维的空间问题。由于稳定性的要求，时间步长受维数的限制，维数愈高，要求时间步长愈小，计算工作量也愈大。50年代中期D.W.裴斯曼和J.道格拉斯等人提出所谓交替方向隐式法，以加快计算速度。如在二维非定常方程中，第一步先对x的导数用隐式差分，而y方向的导数则用前一个的数值。第二步对 y的导数用隐式差分，x方向的导数则用第一步算出来的数值。这一方法的优点是稳定性好，有足够的二阶精度，所产生的差分方程是三对角矩阵方程，便于求解。

④有限基本解法 解位势流动的一种数值方法。航空工业中的低速飞机设计采用位势理论计算各种气动力参数，就是求解二维或三维拉普拉斯方程。在经典流体力学中，用基本解的叠加来解拉普拉斯方程的做法是很成功的。这种方法的要点是，用源、汇、偶极子的分布代替机翼和机身对流场的影响。它们的强度由边界条件确定，结果需要求解积分方程。对一些简单情况可以求解，对一般情况则比较困难。高速电子计算机的出现使这种积分方程的数值解法也有了突破。其主要思想是把积分方程离散化，积分方程代表源、汇等奇点在空间连续分布的总和。例如，若把机翼和机身表面，分割成若干个小单元，每个单元上的奇点强度取平均值。把这些奇点的总和叠加起来，就得出流场总的效应。因此，它用有限项的求和来代替积分，而最后要解的是一组代数方程。由于基本解都是具有奇点的函数，所以这种方法又称为有限奇点法或鳞片法。(见有限基本解法)

跨声速流动数值解  跨声速流动的流场是既有亚声速区又有超声速区的一种混合流场。在不考虑粘性影响和小扰动的情况下，定常二维速度势方程是混合型的，即V+＝0，式中V是来流马赫数Ma与的复杂函数。V＞0是亚声速区(椭圆型)，而V＜0为超声速区(双曲型)。美国的E.M.穆曼和J.D.科尔在1971年首先采用混合差分格式，并运用松弛法成功地解出定常小扰动速度势方程。混合差分格式就是在亚声速区用中心差分格式，所有邻近网点上的条件都会影响计算点，而在超声速区，则用迎风格式，因为上游迎风网点正好是双曲型波动方程的依赖区。(见跨声速流数值计算)

超声速流动数值解  在超声速流动中，主要问题是如何处理激波。用数值方法处理超声速流场中的激波现有两种方法。一是激波捕捉法，另一是激波装配法。激波捕捉法对激波本身并不需作任何特殊处理，只是在计算公式中，直接或间接地引进“粘性”项，自动算出激波的位置和强度，以“捕捉”激波。其中又有所谓人工粘性和格式粘性两种方法。人工粘性方法是J.von诺伊曼和R.D.里希特迈尔于1950年首先提出的，它是以真实粘性流体的物理理论为基础的一种自动处理激波近似方法。该法是在激波层内，人为地加入粘性项，使激波间断变成光滑的过渡区。近年来，在超声速流动中得到广泛的应用。格式粘性是通过某种差分格式间接地引入粘性项拉克斯格式。拉克斯-文德霍夫格式和麦克马克格式都具有类似的效果。激波装配法是把激波仍当作间断面来处理，激波前后要满足激波跳跃条件。但是在普通坐标中，它的实现很困难。一般采用坐标变换，使激波位置(此时是未知的)和一个坐标轴重合，然后把激波看作内边界。这种处理是比较精确的，但也是很麻烦和不方便的。最好的办法是把激波捕捉法和激波装配法结合起来。例如在流场外围的离体激波用激波装配法，在流场内的激波用激波捕捉法。(见超声速无粘绕流数值解)

粘性流动数值解法  可参见纳维－斯托克斯方程数值解、边界层方程数值解法和湍流数值计算等。

参考书目 朱幼兰等著：《初边值问题差分方法及绕流》，科学出版社，北京，1980。 P.J.罗奇著，钟锡昌、刘学宗译：《计算流体力学》，科学出版社，北京，1983年。 P.J.Roache，Computation Fluid Dynamics，Hermosa Pub.，Albuquerque，1972.