计算流体力学课件完整.ppt

《计算流体力学课件完整.ppt》由会员分享，可在线阅读，更多相关《计算流体力学课件完整.ppt（221页珍藏版）》请在三一办公上搜索。

1、计 算 流 体 力 学基础,课时安排：总学时小时，小时讲课；8+8小时上机练习。主要相关前修课程计算机语言、工程流体力学、高等数学 主要内容介绍流场计算的基本概念、基本方法和简单算例,第一章 概 述 1.1 计算流体力学的发展及特点简述 流体力学研究三种方法:实验研究、理论分析和数值计算。实验研究真实可靠、是发现流动规律、检验理论和为流体机械设计提供数据的基本手段。实验要受测量技术限制，实验周期长、费用高。理论研究在研究流体流动规律的基础上，建立了流体流动基本方程。对于一些简单流动，通过简化求出研究问题的解析 解。,对于实际流动问题，通常需运用流体力学基本方程，借助于计算机求数值解（计算机数值

2、模拟）计算流体力学CFD。计算机数值模拟 数值模拟耗费小、时间短、省人力，并能对实验难以 测量的流动进行模拟，如燃烧室、转子通道内。在航空航天、核工业、天气预报、海浪和风暴潮预报 等方面有极广泛应用。在航空航天方面，可用于计算飞行器飞行过程中周围 流场（计算出升力、阻力）。计算航空发动机各部件 内部流场，以及整台发动机三维流场。目前国内有一 些使用较多的商用软件，如fluent、Star-CD、numeca等。,美国自上二十世纪八十年代后期，由于CFD方法应用，使一台发动机设计时间从10-15年降到5-8年，试验样机数 从40-50台降到10台左右。美国NASA主持建立了推进系统数值仿真系统。

3、数值模拟与实验研究、理论分析关系 三者相互依赖、相互促进 数值模拟所占份额会越来越大（计算机技术迅速发展、计算方法的不断改进）。,1.2 流场数值模拟概念 流场数值模拟概念也称为流场计算机模拟，是以计算机为手段，通过数值 计算以数据和图像显示，再现研究对象及其内在规律。数值模拟可理解为用计算机做实验。比如一个机翼绕流 问题，通过计算可得到其升力、阻力数值；绕流流线、激 波位置、流动分离、涡的生成和传播 流场数值模拟几个步骤 建立数学模型：根据流动特点建立适当的数学模型（控制 方程）确定计算方法）控制方程的离散方法：将流体力学基本方程转化成可 用计算机语言描述的形式，称为离散方程，有限元、有限差

4、 分、有限体积等。,）边界条件的处理：有无滑移、壁面等温绝热等。编制计算机程序或运用已有程序进行计算）网格生成：在流场中按一定规律分布一些点，称为网格节点。此过程通常称为前处理。无限信息空间用若干个点近似表示）流场计算：运用离散方程求出每一网格节点上气动热力参数值，如温度、压力、速度。）计算结果后处理：根据网格节点上参数值，进一步处理出需要得到的信息，如流动阻力、升力、流量、流线等。,涡轮叶片通道内三维流计算实例,压气机转静子表面压力分布,涡轮通道内速度分布,航天飞机表面网格,航天飞机表面流速矢量图,航天飞机表面温度分布,气流绕圆柱体流动压力分布,气流流过汽车,风扇流动,直升机旋翼运动,NAC

5、A0012翼型绕流流线图,翼型绕流流线图,风力机表面极限流线图,轴流叶轮计算与实验叶片表面极限流线,轴流叶轮计算与实验性能比较,轴流叶轮计算与实验流场结构比较,第二章 流体力学数值计算数学模型及定解条件本章所涉及的基本方程有两类：流体力学基本方程，基本出发点：质量守恒、动量守恒和能量守恒简化模型方程：具有流体力学基本方程的某些特性，用于对所对应的流体力学方程理论分析 21 可压缩非定常粘性流数学模型 连续方程：运动方程：能量方程：上述基本方程构成了Navier-Stokes(简称NS)方程。,在三维直角坐标系下Navier-Stokes方程为：,上述方程组不封闭，还需要补充数学关系式：）状态方

6、程：）物性系数与状态参数关系：,22 不可压缩非定常粘性流数学模型当来流数小于0.2时，为不可压流动，以下为二种不可压粘性流动控制方程。1）不可压流Navier-Stokes方程 连续方程：运动方程：能量方程：2）流函数-涡量方程:对于平面流动:,平面流动速度与流函数涡量关系:,2.3 无粘流数学模型 1）欧拉方程：2）全位势方程：上式中:为音速;3）不可压流全位势方程：,2.4 常用的模型方程 流体力学基本方程大都为复杂、非线性方程（组），从数值计算角度分析研究比较困难。并且迄今为止还没有形成成熟的理论。为了认识基本方程的数学性质，常用一些简单的线性数学方程作为替代进行研究。这些方程具有基本

7、方程的某些特征，称之为模型方程1）对流方程：此方程是双曲型方程，形式类同于一维欧拉方程。,2）伯格斯（Burgers）方程：是一个非线性方程，具有NS方程类似的性态，式中系数 相当于流体的粘性系数。3）对流扩散方程：这个方程和伯格斯方程同属双曲抛物型方程，但它是 线性的，比较简单。当=0时，退化成双曲型方程，当=0时，则变成抛物 型方程4）抛物型方程：,5）椭园型方程：称为泊松方程，其右端函数项f为已知；若f=0,则成为拉普拉斯方程。,25偏微分方程的数学性质及其与流体运动的关系 流体力学基本方程及模型方程属偏微分方程（组），由于方程的复杂性通常无法采用积分方法求精确解，但可将其离散进行数值求

8、解。流体力学方程（组）的数值求解需符合流动的物理规律，同时边界条件的给定也要遵循流动的物理规律，因此首先需了解方程的数学性质。251 拟线性偏微分方程组的分类拟（准）线性方程组 对于流体力学控制方程，所有最高阶偏导数项都是线性的（这些项前仅有一个系数项，系数项是变量的函数、没有最高阶偏导数与偏导数项的乘积）,拟线性方程（组）的数学性质 以下列拟线性方程组为例 式中，系数项 是x,y,u,v的函数。u,v是因变量为独立变量x,y的函数，并且u,v是x,y的连续函数。将下式：与以上四式组合在一起并写成矩阵形式可得,（2.21）,（2.22）,令矩阵A为上式的系数矩阵，即：并将A矩阵的第一列用(2.

9、23)式右侧矢量替代构成矩阵B,（2.23）,根据Gramer法则，有同理可求出du,dv,dx,dy计算：在xy平面内任一点P，过P点作一曲线ab，如果点2无限接近于P点，则：ab曲线是任意选定的，其选择不影响计算结果。但如果选择的方向使得则无法采用(2.24)计算 值 ef称为通过P点的特征线,(2.24),所谓特征线即为通过xy平面内某点P的曲线，沿此曲线方向无法确定u和v的偏导数值。因此可通过求解：确定特征线。由 展开得：进一步可得 由上式可确定xy平面内每一点的特征线斜率，从而确定特征线。如果令：,则上式可写成：即：令：，如果在xy平面内某一点有：1），则偏微分方程组(2.21)有两

10、条各不相同特征线，称方程为双曲型；2），则偏微分方程组(2.21)只有一条特征线，称方程为抛物型；3），偏微分方程组(2.21)没有特征线，称方程为椭圆型。双曲型、抛物型和椭圆型实际上是直接借用以下二次曲线性质,252 偏微分方程组分类的通用方法 以上根据Gramer法则给出了拟线性方程组类型的确定方法。下面介绍另一种方程组类型通用确定方法。为简单起见，假设方程组(2.21)右端项为0，即：定义矢量：这样式(2.29)可写成矢量形式：(2.30)或者：(2.31),（2.29）,上式可变成：上式中 矩阵的特征值决定偏微分方程组类型。如果特征值全是实数，方程组为双曲型；如果特征值全为复数，方程组

11、为椭圆型。例 二维无旋、无粘定常可压缩流，流场中有一细长体，如机翼翼型。如果在上游有一小扰动，扰动速度分量为：。根据连续方程、运动方程和能量方程可推得：为自由来流马赫数。确定以上流动的类型。,方法一：采用Gramer法则。对照式（2.21）有：而：于是：因此当流动超音时，方程组为双曲型；当流动亚音时，方程组为椭圆型,方法二：采用特征值方法 以下方程：可写成以下矢量形式：所以：由：,求出特征值 因此采用方法二计算结果与方法一相同。由两个结果比较可看出：上式中的矩阵特征值即为特征线在某一点的斜率。,253 流体力学控制方程类型及其对流场数值计算的影响 根据具体流动特点，采用的流体力学控制方程组可分

12、为:双曲型、抛物型和椭圆型。一、双曲型方程,在二维空间坐标(x,y)下有一点P，对于双曲型方程组有二条特征线通过该点，分别称为左特征和右特征。P点的影响区域仅局限于二条特征线之间下游区域，也就是说，P点产生的扰动影响在此区域可感受到，同时也只有此区域可感受到。影响P点区域仅限于二条特征线之间的上游区域，就是说此区域并且也只有此区域的扰动会影响P点。,对于控制方程为双曲型方程的流动问题可采用空间推进方法进行求解。如上图，可给定y轴上的流动参数作为初始条件，然后沿着x轴方向一步一步推进即可求得整个流场。可推得以下几种流体力学控制方程组属于双曲型控制方程组。,例1：定常无粘超音速流 超音速气流流过一

13、双圆弧机翼，在翼型前缘产生弓形激波，激波后气流仍为超音速。可以证明这种流动控制方程组为双曲型（流动可近似采用小扰动方程描述）。对于此流动可在翼型上游设初始边界ab，边界上流动参数取自由流参数，沿x方向向下游推进即可求得整个流场。,例2：非定常无粘流 对于非定常的欧拉方程组，无论流动是否超音都是双曲型（关于时间是双曲型的）。对于一维非定常流，在xt坐标系中，阴影部份为P点的影响区域；P点解由在x轴上（即初始时刻，t=0），区间ab数值确定。管道内一维波运动为曲型的一维非定常无粘流例子。,通常在流场数值计算中更多采用非定常方程时间推进求定常解。只要边界条件不随时间变化，当计算推进时间足够长时，流动

14、趋于定常、流动参数不再随时间变化，这时得到的解即为定常解。,采用非定常方法求定常解的求解过程似乎绕了弯道。实际上对于工程中的有些定常流动问题，采用定常流控制方程无法求解。比如：超音速流绕钝头体的流动，属于超音和亚音混合流动问题。超音区域流动属双曲型；亚音区域流动属椭圆型。在流场计算出来以前无法确定超音区和亚音区的分界线，同时目前还没有对于不同类型的流动都适用的求解方法。将此类定常流动控制方程加入非定常项变成非定常流控制方程，而非定常流动方程无论在亚音区还是超音区都属于双曲型方程（关于时间），因此解决了此类流动不能求解的困难。,二、抛物型方程 根据前面分析，对于抛物型方程通过任一点只有一条特征线

15、。如图2.5，假设过P点有一条垂直于x轴方向的特征线，则P点的扰动将影响特征线右边的阴影区域。抛物型方程与双曲型方程一样可采用空间推进方法求解。首先给定初始边界ac上数据，沿x方向推进即可求得边界ab和cd间的解。,例1：附面层流动 对于附面层流动，通过对NS方程进行简化处理得到适用于附面层内流动的简化方程组为抛物型。给定附面层进口边界ab和ef上数值，采用沿壁面方向空间推进即可求出整个附面层内流动。壁面采用无滑移边界条件，bc和fg两个外边界采用无粘流计算结果。采用附面层方程计算附面层内流动，需先给定附面层外边界流动参数。附面层外边界流动参数决定附面层厚度发展，附面层厚度又影响附面层外势流区

16、流动。因此附面层与势流区流动相互影响，需采用迭代方法进行流场计算，计算方法复杂目前已少有人采用。,附面层流动分析,三、椭圆型方程 对于椭圆型方程，无特征线或特征方程是虚根。流场中任意一点P的扰动会向周边任意方向传递。因此边界点的数值同样影响流场中任意一点的解。在所有边界上都要给出边界条件。,通常边界条件有以下三种类型：1）给定变量u,v数值，此类边界条件称为Dirichlet边界条件；2）给定变量u,v的导数值，此类边界条件称为Neumann边界条件；3)部份边界给定变量u,v数值，部份边界给定变量u,v的导数值，称为混合边界条件。,椭圆型方程影响区域,例：定常无粘亚音流动控制方程 该方程属椭

17、圆型方程。在此关键是流动亚音，因为对于亚音流，流场中一点的扰动理论上可向各方向传递到无限远处。如下图亚音翼型绕流，翼型上游的流线向上折转，翼型下游的流线向下折转。翼型产生的扰动引起整个流场的变化（理论上直至无穷远处）。,亚音速翼型绕流,26 流体力学问题的定解条件数学方程建立后，为确定解必须给出定解条件定解条件包括初始条件和边界条件一、初始条件初始条件就是在某一起始时刻给出流场中速度、压力、密度和温度等参数分布。对于定常问题并不需要初始条件实际计算，对于非线性方程（组）要进行迭代求解，需要初始条件作为迭代的初值。初始条件给定不影响最后结果，但初始条件的合理性会影响迭代计算收敛速度，甚至于影响收

18、敛性。,二、边界条件关于各种流动边界上要给多少个边界条件、给出哪些边界条件，目前还没有一个完善的理论。对于绝大多数工程实际中的流动问题，研究人员根据理论分析结合经验都能给出合适的边界条件。下面介绍一些常见的流动边界及边界条件。1）来流边界（进口边界）对于外流流动前方边界称为来流边界；对于内流流动，如进气道和叶轮机内流动进口截面称为进口边界。来流边界理论上应在物面上游无穷远处，在那里流动未受扰动易于给出边界条件在此边界上一般给出：总压、总温、气流角等参数,2）下游边界（出口边界）下游边界（外流流动）和出口边界（内流流动）要设定在绕流体的远下游，在那里流动通过充分掺混已比较均匀，这样有利于边界条件

19、的给定。对于亚音速流，通常给出出口边界上静压（又叫出口反压）；对超音速流，由于下游扰动对上游流动没有影响，因而不能给定出口反压。其他未确定参数如速度、密度、温度以及超音速流的静压等，则采用计算区域内部的数值外插求得。,3）壁面边界 速度的给定a)粘性流，流体在壁面边界上的速度等于壁面的运动速度，如果壁面静止，则流体速度为零，即无滑移边界条件 b)无粘流，流体在边界处的法向速度为零，而切向速度则由计算求得不再为零，即滑移边界条件 温度的给定 a)等温壁，给出壁面温度，并假设壁面处流体的温度与壁面温度相同 b)绝热壁，壁面热流量为零，即：压力的给定：壁面法向压力梯度为零，即：,第三章 有限差分近似

20、及其数学性质 计算流体力学任务是将描述流体运动的偏微分方程转化成离散形式，然后在计算机上求出这些方程的解。方程的离散方法有：有限差分法、有限元法、有限体积法等 有限差分法用差商代替微商，将微分方程转化成差分方程。实现偏微分方程的离散化，以适合于计算机编程计算。3.1 差分格式基本概念对于一个二维定常问题，如果求解域如图示在直角坐标系下，变量可表示成：U(x,y)流场中任一网格节点表示为（i,j）,i=1，m；j=1，n网格点（i,j）上差分计算值表示为，它是对函数值 的近似。,空间步长时间步长,流体力学方程是偏微分方程，主要由一阶和二阶偏导数项组成 差商代替微商,3.2 常用偏导数项的差分格式

21、和精度分析 321 一阶偏导数差分格式一阶偏导数通常有。常用的有中心差分和向前、向后差分格式。1.向前差分格式 由泰勒级数：,由于,所以,截断误差,是步长的一次方，称此差分格式为一阶精度，记作：,在数值计算过程中，时间和空间步长取值都很小，因而截断误差R数值也很小，这样确保用差商替代微商有足够的精度。,忽略掉截断误差项,上式为向前差分格式，其精度为一阶。,.向后差分格式,精度为一阶,.中心差分格式,精度为二阶,中心差分比向前和向后差分离散精度高。差分格式选择：a)考虑差分格式的稳定性 b)在边界上适用性,322 二阶偏导数差分格式1.普通中心差分2.普通一侧差分格式 3.二阶混合偏导数项通常采

22、用中心差分,3.3 差分方程和相容性 差分方程：偏微分方程中的偏导数项用差商代替得到的差分形式方程 1)的差分方程 a)时间向前、空间中心的差分格式(FTCS),差分方程的截断误差：差分方程与微分方程之间存在一个误差，对于上方程：,这个差分方程具有一阶时间精度二阶空间精度,定解条件离散化,a)初始条件,离散形式,b)边界条件,离散形式,由差分方程和定解条件采用时间向前推进可求出n=2、3、4各时间层上内部节点上的函数值。,格式图：表示差分方程相邻网格节间关联性的图形。图中表示方程在该网格节点上离散，表示差分方程所涉及的网格节点。,(a)FTCS 格式(b)FTFS 格式(a)FTBS 格式,b

23、)时间向前、空间向前的差分格式(FTFS)c)时间向前、空间向后的差分格式(FTBS)差分方程相容性分析 微分方程:对应的差分方程:,截断误差为:如果有:微分方程与差分方程相容微分方程的定解条件为:对应的差分问题的定解条件:截断误差为:如果有:微分方程与差分方程定解条件相容,有限差分方法求解流体力学问题举例 在两固定平板间，流体由于平板二端压差驱动，作层状流动。这种流动称为库特（Couette）剪流。不考虑端部效应，在每个等x截面流体速度分布完全相同。,由NS方程可推得关于速度u(y)的微分方程,边界条件,求出解析解为：,有限差分方法进行求解步骤）建立基本方程和适当的定解条件）网格划分：沿y方

24、向将线段6等分，则空间步长为）偏微分方程及边界条件差分离散,离散方程,中心差分,离散边界条件,节点,节点,节点,节点,节点,）编制计算机程序进行数值计算，得结果：,这个例子虽然非常简单，但反映了流体力学问题有限差分数值计算的全过程。通常对于工程实际问题，大部份工作量是花在第步（网格生成）和第步（编制计算机程序进行数值计算）。,3.4 差分方程的收敛性差分方程收敛性定义：当步长趋于零时，差分格式的解能否趋近于微分问题的解称为差分格式的收敛性对差分网格上的任一网格节点（i,n），设差分格式在此点的解为，相应的微分问题解为二者之差为 称为离散误差。如果有：此差分格式是收敛的，即差分方程的解收敛于相对

25、应微分问题的解。否则不收敛。,条件收敛差分格式在一定的约束条件范围内是收敛的。无条件收敛指差分格式在任何条件下都收敛。关于差分方程收敛性举例。例：一维问题如下：,精确解为：,差分方程:或：因此有:,即：,进一步：,而：,因此：,即：,所以上述差分方程收敛,例2:对流模型方程 FTBS差分格式:离散误差为:精确解满足:于是有:,若:即有:于是:.将上面n-1个不等式相加得:,以上例为无条件收敛；例为在条件下收敛的差分格式。对于流体力学问题由于是多个方程组成的非线性方程组，差分格式收敛性的证明目前还是比较棘手的数学问题。目前应用比较广泛的是采用冯.纽曼（Von Neumann）方法通过差分格式稳定

26、性分析来证明收敛性。,35 差分方程的稳定性及稳定性分析 差分格式的依赖区间、决定区域和影响区域 初边值问题:1）采用FTCS格式离散 2）采用FTFS格式离散 3）采用FTBS格式离散,FTCS格式 FTFS格式 FTBS格式 差分解的依赖区间和决定区域,FTCS格式 FTFS格式 FTBS格式 差分解的影响区域,同一微分问题，当采用不同的差分格式时，其依赖区间、决定区域和影响区域是不一样的。进而影响差分方程收敛性。,考察步长比 对误差传递的影响方程:其解为零,即:采用任何差分格式,若计算中不产生误差,有:假设产生计算误差:采用FTBS格式有:,上例子显示了该格式的影响区域同时显示了数值不同

27、时，计算误差对后面时间层上节点的计算值所产生的影响。），所产生的影响在数值上不会扩大。），所产生的影响在数值上会越来越大，这样随着计算向前推进，误差会将真解湮没，并最终导致数值趋于无穷，计算失败（发散）。在采用有限差分法的运算过程中，计算误差总是不可避免的（如计算机舍入误差）有些情况下，误差在传播过程中逐渐衰减；而另一些情况下，误差在传播过程中会逐渐递增、积累。若计算中产生的误差，在一定条件下逐渐衰减，那么就称这个差分格式在给定的条件下稳定，这个条件就是它的稳定准则。反之则称差分格式不稳定。,以对流扩散方程，说明稳定性概念,采用FTCS格式离散得：,假设在时刻以前的运算中不产生任何计算误差，在

28、时刻以后的运算中也不产生新的计算误差，只在时刻产生了误差。考察误差传播。这样有：,这时有：,于是：,两式相减得：,上式称为这就是误差传递方程 误差传递方程与原差分方程形式相同 上式可改写成：,此式右边第一项是由于对流项而产生的误差增长，第二项为由于扩散项而产生的误差增长。讨论这二项引起的误差增长情况。,误差沿节点分布情况可以是各种各样的。无论误差分布呈何种形态，随着计算由n 时间层向前推进，稳定格式误差应逐渐减小，而不稳定格式将逐渐增大。,假设某时刻产生的误差沿节点是振荡的，并且振幅沿节点增加,.不考虑粘性，即无扩散项，。a)结果使误差振幅随时间n单调增加,因而不稳定,b),误差随着n的增大逐

29、步减小。格式稳定如果 过大,校正会过头，称为过冲 稳定条件,2.不考虑对流项,当,当,误差随着n的增大逐步减小。格式稳定如果 过大,校正也会过头,对流项和扩散项同时存在时，它们各自所产生的误差在传递过程中将会相互影响，在一定的约束条件下，差分格式稳定,稳定性的数学定义 将差分解的误差扩展成连续函数 Z(x,t)，如果 则对应差分格式稳定。K是有限常数。考虑边界和初始值（即定解条件），稳定性定义式写成：差分问题在初始时刻或某任一时刻引入的误差为小量，此后的解与差分问题的精确解的误差也一定为小量，所以差分格式为稳定格式。,取：,则：,Von Neumann 的稳定性分析方法，又称傅氏级数（Four

30、ier）方法,考察对流方程：,格式离散：,误差传递方程：,误差是节点上的离散量，现将其扩展成空间连续量。这样节点误差的傅氏级数为:,代入离散方程得：,所以对于任意k值：,取：,误差放大因子,如果：,即差分方程稳定,考察对流方程FTBS格式稳定条件,要使：必须,即为差分方程稳定条件,对流扩散方程FTCS格式差分方程稳定性分析 误差传递方程为：,稳定性条件,误差的傅氏级数简写成,3.6 差分方程的相容性、收敛性和稳定性的关系,前面讨论了差分问题的相容性、收敛性和稳定性已经知道相容性是收敛性的必要条件发现稳定性与收敛性之间有一定的联系Lax等价定理：对一个适定的线性微分问题及一个与其相容的差分格式，

31、如果该格式稳定则必收敛，不稳定则必不收敛。适定指适当的定解条件 若线性微分问题适定，差分格式相容，则稳定性是收敛性的充分必要条件。,稳定性,线性微分问题适定,差分格式相容,收敛性,由于收敛性的证明通常比稳定性证明要难，故借助于Lax定理，可将收敛性证明转化成稳定性的证明。,第四章 模型方程的常用差分格式 41 对流方程的差分格式 1.逆风差分格式 稳定条件 截断误差格式分析,为格式,）,a相当于动量方程中的速度u 在节点i采用向后差商，因此称为逆风差分。,根据流动的物理规律，流场中某点的流动参数受上游流动影响比下游大，并且速度越大差别越大。当流动超音时，就不再受下游影响。因而从流动机理分析，采

32、用向后差商比中心差商和向前差商稳定性好（实际上这时采用后两种差分格式是不稳定的）。,）,为格式,2.Lax-Wendroff格式 由泰勒展开式,采用Von Neumann方法分析可得，这种格式的稳定条件为。由泰勒级数分析可得截断误差为,利用微分方程,因此有：,稳定条件 截断误差3全隐格式上述两种差分格式都有较严格的稳定性条件。为了扩大差分格式的稳定范围，还可以构造隐式格式。隐式格式就是在差分方程中，n+1时间层上有多个节点函数值出现。显式格式就是在差分方程中，n+1时间层上有只有一个节点函数值出现。,采用中心差商逼近式中一阶和二阶空间导数，并略去高阶小项，得到：,格式恒稳，即无条件稳定,如果：

33、,则得差分方程：,显式格式时间推进求各时间层的节点函数值过程直截了当。隐式格式求n+1时间层节点i的函数值，涉及到相邻节点i-1和i+1的函数值。,对流方程：,42 扩散方程差分格式 1.古典格式稳定条件为：2.三层全隐式格式,43 对流扩散方程差分格式 1.中心显式差分格式 稳定条件：2.逆风差分格式,稳定条件：3.全隐差分格式,4.4 计算实例,假设在两相距1m的无限大平板间充满水，平板原来都处于静止状态，在某一时刻t=0，上平板突然以恒定速度平动，求在任意时刻t两板间水的速度分布。由于两板无限大，可忽略端部效应，这样每一个等x截面速度分布相同。因此只需求x=0截面的速度分布。,1）建立控

34、制方程和定解条件 运用非定常二维不可压粘性流的基本方程，根据此流动具体情况对方程进行简化，最后可得到流动的控制方程：,边界条件,初始条件,2）网格划分 采用均匀网格，网格点数为m，边界节点分别为1和m，那么网格空间步长为,3）控制方程和定解条件的离散化采用古典格式进行差分离散（时间向前、空间中心），得,4）编制和调试计算机程序,稳定条件为,边界条件的离散形式,初始条件的离散形式,编制程序要注意程序的条理性、可读性、以及通用性等。要有好的条理性，程序则要基于结构化思想进行设计，即一个程序要分成若干块，每一块赋予各自功能，块与块间逻辑关系清晰。可读性与条理性是紧密相关的，为了增加程序可读性，还要在

35、程序中适当加入说明语句；以及程序的外观布局、变量所采用的符号都要有所考虑。通用性表现在两方面，一是避免程序的局限性，比如对于现在的例子，程序中网格节点数m不要是一个确定的数,而作为变量，这样通过赋值语句对m赋值，可方便地给出不同网格节点数的 差分解。其二是在保证程序对确定流场计算模拟功能的前提下，最好也能兼顾计算其他相近流场。,START,网格划分,初始给件,n=1,边界条件,内点值计算,结果输出,否,n=n+1,程序流程图,程序调试要求上机前阅读参考源代码完成程序调试测试稳定性条件修改程序，加入图形显示功能,45 多维问题的几种常用差分格式 实际问题多为空间二维和三维的问题，上面所述的针对空

36、间一维差分方法向二维和三维的推广并没有太大的理论难度。本节以扩散方程为例，将一维差分方法向二维推广，其他情况处理方法类似。1.加权平均差分格式,二维扩散方程为,相应的定解条件为,采用均匀网格：1.加权平均差分格式,其中,当时格式是显式的；当时格式是隐式的。对于二维问题，隐式格式直接求解需要解一个大型稀疏代数矩阵，比较麻烦且耗费机时，后面要进一步介绍对此类问题的处理。,稳定性分析,差分方程,放大因子：,令：,则：,）：格式无条件稳定,）：,稳定条件：,）：,或,2.交替方向隐式格式（ADI方法）对于多维问题，采用隐式格式要求解大型稀疏矩阵，而采用显式格式稳定性限制又较严格。交替方向隐式格式综合了

37、显式和隐式格式特点，它的基本思想是将差分计算分成两步。第一步，在一个方向上（比如x方向）是隐式的，而另一个方向是显式的；第二步则两个方向交换，即在第一个方向是显式的，而另一方向为隐式。在二步计算中，由于只有一个方向是隐式，这样每一步求解的方程组都是三对角方程组，所以求解过程大为简化。同时格式的稳定性条件比之于显式格式也会大为放宽。因为计算在二个方向上交替进行，所以叫做交替方向隐式格式（Alternating Direction Implicit method,ADI）,对扩散方程：,采用交替方向隐式格式，得到差分表达式为：,在进行n+1/2时间层计算时，涉及到x方向三个节点函数值要同时求得，因

38、而与一维隐式格式相同；同样对于n+1时间层也涉及到y方向三个相邻节点函数值。利用Von Neumann 方法进行稳定性分析，格式无条件稳定。交替方向法另一形式,上两式合并可得,因此相当于全隐式格式中加入一高阶小量，此格式是无条件稳定,向三维推广,上式表明：对于三维问题，先在x方向采用隐式格式进行第一步计算；再在y方向采用隐式格式进行第二步计算；最后在z方向采用隐式格式进行第三步计算。这种格式是无条件稳定的。,3.时间分裂格式,改写成,忽去高阶小量,分解成,由泰勒级数,这种格式构造的基本思想是将多维问题化为几个一维问题。,相当于求解：,稳定性条件：,时间分裂差分格式由于其构造和计算方法简单而有较

39、广泛应用,*46 数值效应,第五章 不可压流场的数值计算 流动分类：定常和非定常流；有粘和无粘流；可压和不可压流。从数值计算角度考虑，将流场计算分为不可压流计算和可压缩流计算。对于可压缩流，无论其为有粘、无粘、定常或非定常都可以采用同一类方法进行计算。但可压缩流的计算方法一般不能用于求解不可压流。,51 不可压无粘流场计算的流函数涡量法 511 基本方程推导二维不可压粘性流,对于无粘流，且忽略体积力,流场无旋,上式在二维直角坐标系下展开形式,泊松方程,512 泊松方程的差分求解,中心差分,若假设,若,512 泊松方程的差分求解,差分方程求解方法 要求解(m-2)*(n-2)个代数方程组成的方程

40、组 1)采用常规的代数方程组求解方法来求解，比如:高斯消去法 这种直接方法求解，如果节点数很多，也就是代数方程个数很多时，就会很耗费计算机内存和机时。2)迭代计算方法 方法简单，求解速度快,且节省计算机内存。,一、黎曼方法 迭代方法相当于利用非定常方程求定常解。泊松方程是一个定常方程,要求解这个方程，可以建立一非定常方程,构造方程的FTCS差分格式,若:,黎曼迭代公式,黎曼迭代的计算步骤 1）将连续函数表示的边界条件转化为离散形式,左边界：右边界：,下边界：上边界：,2）给定初场，即给定初始时刻(n=1)假想的内部节点上的 函数值（可采用边界值线性插值）,3）采用黎曼迭代公式由左到右、由下而上

41、逐点，或者由 下而上、由左到右逐点进行迭代。,4）收敛程度判断和输出结果。,二、点松驰法,松驰因子当1时称为超松驰；而=1即为黎曼迭代。如果迭代计算过程不收敛，则需减小。但越小收敛速度越慢。超松驰的收敛速度快。,三、对于非正方形网格（）,1)黎曼迭代公式为,2)点松驰迭代公式,四、线松驰法,黎曼迭代和点松驰迭代法都是点迭代法，即计算是逐点进行的。在计算(i,j)节点函数值时，周围四个节点(i-1,j)、(i+1,j)、(i,j-1)、(i,j+1)上函数值都作为已知。线松驰法把与(i,j)节点的同一行或同一列上的相邻两个节点函数值不再当成已知，从而需要与(i,j)节点函数值同时求出。即计算是逐

42、行或逐列进行的。,四、线松驰法,例如在同一列中进行线松驰迭代，黎曼迭代计算变成 1)列迭代,2)行迭代,线松驰法达到规定的误差所需要的迭代次数大约是黎曼法的一半。但每一次的迭代计算，线松驰法要解一个三对角矩阵。,513 追赶法求解三对角方程组,求解方程组,i=1:,i=2:,i=3:,i=N-1:,i=N:,.,设,写成矩阵形式为,方程组求解,1)消去过程（追的过程）,.,2）回代过程。,SUBROUTINE BANFAC(N)DIMENSION A（100），B(100)，C（100），F（100），XM（100）C 求解三对角矩阵C 追赶过程 NP=N-1 DO 1 J=1,NP JP=J

43、+1 A（JP)=A（JP)/B（J)B(JP)=B(JP)-A（JP)*C（J)F(JP)=F(JP)-A(JP)*F(J)CONTINUEC 回代过程 XM(N)=F(J)/B(N)DO 2 J=1,NP JA=N-J XM(JA)=(F(JA)-C(JA)*XM(JA+1)/B(JA)2 CONTINUE RETURN END,附加:将此程序改成C语言程序并进行测试,5.1.4 计算举例（内置方形体的突然扩张通道流）,如下图，流体由一小孔流入突然扩张通道，通道内放置一方形体，假设流动无粘、无旋、不可压且定常，求流场中速度分布。YL1=0.1,YL2=2.0,XL1=2.0,XL2=2.0

44、,XL3=3.0。,根据流动特点可得流动控制方程为：,边界条件：）沿ABCDEF是一条流线，给定；）沿IHG也是一条流线，给定。）由于出口截面离绕流体较远，因而可设,网格划分如图中所示采用矩形网格，在此要注意划分的网格要尽可能使流动边界与某一条网格线重合，以便于离散化边界条件的给定。本算例采用x方向70等分，y方向20等分即可以满足这个要求，这时 方程离散 方程采用中心差分离散，并构造点松驰迭代计算公式,边界条件离散,AB段：,BC段：,段：,段：,段：,段：,段：,段：,因为,要求:1)阅读源程序，并对各模块加说明文字；2)完成程序调试，运用Tecplot 画出速度矢量图和流函数等值线图；3

45、）改变松驰因子，测试其对收敛性影响。4）改成小孔流出,52 不可压粘性流场计算,不可压粘性流的基本方程为,在二维直角坐标系下展开得,5.2.1 不可压粘性流求解的流函数涡量法,1基本方程,引入流函数和涡量消去上方程中的压力项，可得：,如果需要求得压力场,采用时间向前、空间中心差分得：,2涡量边界条件,设壁面以速度运动，由于粘性无滑移条件，壁面上沿x方向速度分量为,又根据无渗透条件：沿壁面v恒为0，因而有,由此得,上面的涡量边界值计算公式是一阶精度。下面构造更高精度的涡量边界值计算公式。,如果设,为待定系数,在边界上：,由泰勒展开式得：,比较得：,开始,n=1,输入内部节点和边界节点 初始值,n

46、=n+1,由(5.32)式计算内部节点 初始值,由(5.35)或(5.37)式计算边界节点 值,由(5.31)式计算内部节点 值,由(5.35)或(5.37)式计算边界节点 初始值,由(5.32)式计算内部节点 值,输出结果,收敛,是,否,由(5.33)式计算,流函数涡量法计算流程,平板驱动方腔内流场计算,如图所示的水槽，水槽横截面为方形，边长1m。底部为一活动平板，以恒定速度U0沿x方向平移。槽内水在平板的带动下一起运动。设流动雷诺数Re极小，这时动量方程(5.31)中惯性力项比粘性力项小得多可忽略，则涡量控制方程变为,这种流动称为Stokes流，控制方程为二维抛物型方程。采用FTCS格式差

47、分离散得：,如果：并令,流函数控制方程：,流场速度分布计算,流函数和涡量边界条件,）底部边界：,）顶部边界：,）左边界：,）右边界：,要求 1）考察稳定性条件;2）不忽略对流项进行计算并进行结果比较;3)将二阶精度涡量边 界条件改成一阶精度.,5.2.2 不可压粘性流求解的原始变量法,1基本方程 对于不可压流，控制方程为,流函数涡量法推广到三维流动较困难。原始变量法直接以速度和压力为变量采用方程进行数值计算方法。二维和三维流动均可采用 计算过程较流函数涡量法复杂 原始变量法分为：人工压缩方法、投影法和法等。下面介绍人工压缩方法。,不可压非定常流连续方程中没有关于时间的偏导数项，与动量方程性质不

48、同，方程组不易求解。对于求定常解，可考虑在连续方程中加一非定常项。连续方程中加入：最合适。待解的变量：动量方程中有；方程组中尚缺 项，连续方程变为,c为音速，伪压缩方程,在二维直角坐标系下展开得,2基本方程离散,虽然压力分布不均匀，但若根据上式计算压力梯度，则都等于零。为了解决这一问题，选用交错网格对上述方程组进行差分离散。,压力梯度项离散 上方程中有x方向和y方向的压力梯度，如采用中心差分，即,将y方向运动方程在所示一类点上展开:,将连续方程在用表示节点上展开,将x方向运动方程在所示一类点上展开：,边界条件 对于粘性流，在固壁边界上速度等于壁面的运动速度。1)计算时，涉及到点 的函数值，此点

49、在求解域外不存在。为此通常采用线性外插法确定一虚假值，即,2)计算时，涉及到点的函数值，此点也在求解域外不存在。采用类似方法可得：,第六章 可压缩流场的数值计算,61 可压缩无粘流的差分计算 611 一维欧拉方程的显示格式,一维欧拉方程为,欧拉方程的差分格式可直接用于NS方程。虽然粘性流的控制方程（NS方程）比无粘流控制方程（欧拉方程）要复杂得多，其实质是比无粘流多了粘性项，通常不影响数值格式的构造，并且由于粘性项存在会增加格式的稳定性。着重介绍无粘流的差分计算方法，然后再进一步推广到粘性流的计算。,1.MacCormack格式,二步显式差分格式，第一步称为预报步，第二步称为修正步，第一步采用

50、的是一阶精度向后差分，第二步采用一阶精度向前差分，差分格式总精度为空间二阶精度。该格式广泛用于求解欧拉方程和NS方程，并在航空航天工业的流场计算中得到大量应用。MacCormack格式的不足之处是在激波附近数值解中有数值振荡，因而在实际使用时要加入人工粘性项。关于人工粘性在下面的多步龙格库塔格式中将作介绍。,2.多步龙格库塔格式,加人工粘性项,欧拉方程可采用三步、四步、五步这样一些多步龙格库塔格式，在此介绍四步龙格库塔格式。方程可离散成,其中和分别为二阶和四阶人工粘性项。,则取决于节点i的压力二阶导数值，因为：,在激波和滞止点附近压力变化较骤烈，也相应较大，引入的二阶人工粘性就增大。,612