1、流流 体体 力力 学学顾伯勤 主编研 究 生 教 材 退 出中国科学文化出版社目 录流体力学基础第一篇第二篇流体动力学基本原理及流体工程退 出第三篇计算流体动力学 第一篇 流体力学基础 绪论 场论与正交曲线坐标 流体静力学 流体运动学第一章第二章第三章第四章退 出返 回 第二篇 流体动力学基本原理及流体工程 流体动力学微分形式基本方程 流体动力学积分形式基本方程 伯努利方程及其应用 量纲分析和相似原理 流动阻力与管道计算 边界层理论 流体绕过物体的流动 气体动力学基础 第五章第六章第七章第八章第九章退 出返 回第十章第十一章第十二章第三篇 计算流体动力学 计算流体动力学数学物理基础 流体动力学
2、问题的有限差分解法 流体动力学问题的有限元解法 第十三章第十四章第十五章退 出返 回第十四章流体动力学问题的有限差分解法 势流问题的数值计算 回流流动问题的数值计算 第一节第二节退 出返 回第十四章 流体动力学问题的有限差分解法 退 出返 回 第一节 势流问题的数值计算 第第1页页 求解NavierStokes方程的有限差分解法可以用速度、压力作为基本变量，也可以用涡量、流函数作为基本变量。前者称为原始变量法，后者叫涡量流函数法。原始变量法以压力和速度为基本变量，在数值计算过程中可能会出现所谓的波形压力场，从而产生虚假的模拟结果。为了解决这一问题，Patankar和Spalding于1972年
3、提出了著名的SIMPLER算法（Semi-Implicit Method for Pressure-Linked Equations），其中采用了交错网格系统，即将 、和 分别存储于三套网格（交错网格系统）的系统。原始变量法的优点是能够在求解速度的同时解出压力场，并且能推广到三维问题。但应用交错网格导致的复杂性问题至今仍未解决。涡量流函数法避免了原始变量法中压力梯度项引起的一系列问题，不必采用交错网格系统。其突出的优点是简洁方便，容易掌握，尤其适用于不可压缩流体的二维或轴对称流动问题，这也是本章将其作为重点介绍的主要原因。但要指出，此方法不能在求解流场的同时解出压力场，欲计算与压力有关的参数（
4、如可压缩流动中的流体密度）只能在流场解出之后进行。xwywp第十四章 流体动力学问题的有限差分解法 退 出返 回 第第2页页第一节 势流问题的数值计算 另外此方法不易推广至三维问题，因为对三维问题虽然可采用空间涡量来构造模型方程，但失去了求解的简洁性。本章将结合典型流动问题的模拟，对涡量流函数法进行系统阐述，并在最后对原始变量法作一简介。实际流体都是有粘性的，但流体粘性的作用只在具有比较大速度梯度的流动区域才表现出来。按照普郎特边界层理论，只有在固体表面附近速度发生剧烈变化的薄层内，流体粘性的作用才需考虑，而在这一薄层以外，流动可以当作是无粘性的。例如，离机翼较远区域的掠机翼气流流动，大楼的侧
5、墙受到均匀来流的冲击时的气流流动等都可以作为无粘性流动来处理。不可压缩无粘性流动由Euler（欧拉）方程及连续性方程所描述：xpywwxwwtwxyxxx1（14.1a）第十四章 流体动力学问题的有限差分解法 退 出返 回 第第3页页第一节 势流问题的数值计算 ypywwxwwtwyyyxy10ywxwyxxwywp0 xwywywxwywxwxwywxwxwywywxwywtyxyyxxyxxyxyyx 上述三个方程有三个因变量 ,及 ，因而方程组是封闭的。为了用较简便的方法求解速度场，将式（14.1a）对y求导并减去式（14.1b）对x的求导结果，可消去压力梯度项，得：（14.1b）（14
6、.1c）xwywyx令（14.2）0ywxwywxwtyxyx则得：（14.3a）第十四章 流体动力学问题的有限差分解法 退 出返 回 第第4页页第一节 势流问题的数值计算 0ywxwtyx据连续性方程（式（14.1c），上式中括号内的部分等于零，且可进一步写成为守恒形式：（14.3b）式（14.2）所定义的是涡矢量在xoy平面上的分量，简称涡量，因而式（14.3）为涡量守恒方程。它表明在无粘性流动中若起始时刻流场中无涡，且边界上也不产生涡，则整个流动将处处无涡。0 xwywyx0ww（14.4a）（14.4b）其中 为速度矢量。对于三维问题中涡矢量为零（简称无旋流动）的情形，按涡矢量的定义，
7、有：涡量为零的无粘性流动称为势流，其物理意义为流体微团仅有平动与变形而没有绕其自身的旋转运动，此时有：第十四章 流体动力学问题的有限差分解法 退 出返 回（14.5）第第5页页第一节 势流问题的数值计算 由数学分析可知，若 满足式（14.4b），则其三个分量 、必是某个标量函数 的偏导数：wxwywzwzwywxwzyx,02222yx对于二维问题，将式（14.5）代入式（14.1c）得：（14.6）标量 称为速度势函数。上式表明势流的速度势函数满足Laplace方程。xwywyx,02222yx在第十一章中分析二维势流时，常用流函数作为变量。流函数 的定义为：（14.7）不难证明，式（14.
8、7）的定义使连续性方程（14.1c）自动满足。将式（14.7）代入式（14.4a）得：（14.8）第十四章 流体动力学问题的有限差分解法 退 出返 回 第第6页页第一节 势流问题的数值计算 上式被称为流函数方程。在第十一章中已经证明，等势函数线与等流函数线是相互正交的曲线簇。由上面的讨论可知，对二维的有势流动，不必直接求解式（14.1）那样较复杂的控制方程，而是可以通过求解更简单的速度势函数或流函数的Laplace方程获得速度场，即计算出各处的 与 ，进而应用伯努利方程求出各节点处的流体压力。从数值求解的角度来看，求解流函数方程比求解势函数方程更简便，原因有二：（1）速度势函数仅在势流中存在，
9、对于实际粘性流体，不存在速度势函数，而流函数则在粘性流场中也存在。若将（14.7）式代入涡量定义式后得 ，即对实际粘性流体，流函数满足Poisson方程，其数值解法与Laplace方程有许多共同之处。（2）势函数方程的边界条件对无渗透的固体表面都是第二类的，即固体表面上流体的法向速度为零，而流函数方程的边界条件对其则是第一类的，即给定了边界上的流函数之值，第一类边界条件较之第二类边界条件在数值计算中容易使计算收敛。xwyw2222yx第十四章 流体动力学问题的有限差分解法 退 出返 回 第第7页页第一节 势流问题的数值计算 一、二维通道内绕圆柱体的势流如图14.1所示，设在两块相距为H的平行平
10、板之间的中心有一直径为d的圆柱。在圆柱的上游L处有一股均匀的不可压缩流体的势流进入该通道，现在利用数值计算方法来确定流场中各点的速度及流函数。首先确定计算区域与边界条件。由于对称性，只要计算四分之一区域即可，取ABCDE这一块面积作为计算区域。对于流函数的Laplace方程，这一区域的边界条件如下：ywwdyyy00)(AE边界 按 的定义AB边界 为区域的对称线，没有流体横穿此线，故此处 为常数，可以取为零；BC边界 为无渗透的固体边界，为常数，应取与AB上相同之值，故亦为零；CD边界 由于对称性，应有 ；DE边界 也是无渗透的固体边界，为常数，可取为流经该通道的总流量 。0 x2/0Hw第
11、十四章 流体动力学问题的有限差分解法 退 出返 回 第第8页页第一节 势流问题的数值计算 0AELCBODLH图14.1 绕圆柱体的流动d=0.08m，H=0.2m，L=0.12m，=1m/s w0 xyxy下一步进行区域离散化并建立离散方程。采用区域离散A方法，正方形均分网格，即 =0.01。考虑到不规则边界BC上是第一类边界条件，故不必对其边界上的节点补充离散方程，但应对与边界相邻的内节点来单独建立方程。如图14.2所示，与边界相邻的内节点有F、G、H、I和J。按照差分理论处理不规则边界的方法，对图14.3中所示的节点 ，可导得其离散方程为：ji,第十四章 流体动力学问题的有限差分解法 退
12、 出返 回 第第9页页第一节 势流问题的数值计算 IJCGFBH0ri,j+12i-1,ji,j1图14.2 与圆柱面相邻的内节点图14.3 与边界相邻点的内节点示意图babbaabajijiji1111111,121,（14.9）2222111111jxrLxxjrxLa如图14.2所示，对图中的H点（其下标变量为i与j），有：（14.10a）第十四章 流体动力学问题的有限差分解法 退 出返 回 第第10页页第一节 势流问题的数值计算 其中L1为x方向的节点总数。根据F、G、H、I和J各点上的 值，可计算出每个节点对应的a、b值，如表14.1所示。2222111111LxrjxxLrxjbj
13、i,（14.10b）节点编号i，jabF9，20.12701.0000G9，30.53591.0000H10，40.35420.3542I11，51.00000.5359J12，51.00000.1270在对称边界CD上，离散方程为：上式利用了对称性 。对其余边界可按规定的值读入。1,1,1,241jijijijijiji,1,1表表14.1 邻近边界的节点的邻近边界的节点的a、b值值（14.11）第十四章 流体动力学问题的有限差分解法 退 出返 回 第第11页页第一节 势流问题的数值计算 为编制程序的方便，可把计算区域扩大成为一个矩形，即图14.1中的AODE区域。这时可将式（14.9）作为
14、适合于整个计算区域的通用方程，不同的节点可据a、b值的不同来判断，计有三种类型的节点：（1）a1，b1，表明计算点为除FJ五个节点外的其余内节点，可取a1，b=1，此时式（14.9）即可写成为式（14.12）；（2）0a1，0b1，代表上述与曲线边界相邻的五个节点，按计算所得的a、b代入式（14.9）；（3）a0或b0，表明计算已进入被拓宽的区域（圆柱体内部），可规定取该区域中 值为零而不必解代数方程。代数方程的求解可以采用GS迭代法。在对每一节点的值作更新之前，先据节点的坐标 计算a，b之值，再据a或b的三个数值范围（大于等于1，小于1而大于0及小于0）决定其取值。在获得了满足收敛准则的流函
15、数以后，可按定义计算速度 和 ：ji,jxiw,jyiw,第十四章 流体动力学问题的有限差分解法 退 出返 回 第第12页页第一节 势流问题的数值计算 该式适用于除FJ五个节点外的其余内节点，对FJ五个内节点因其与左右或与上下两邻点的距离不相等，为了获得具有二阶截差的一阶导数表达式，可采用下列计算式（以图14.2中H点为例）：ywjijijxi21,1,xwjijijyi2,1,1,ybbbwjijijxi111,21,xaaaawjijijijyi11,12,2,1,，（14.13）（14.14a）（14.14b）0,1ji01,ji其中，。由数值计算得出的流函数线如图14.4所示。第十四章
16、 流体动力学问题的有限差分解法 退 出返 回 第第13页页第一节 势流问题的数值计算 DEFABCxywo u twinHLH2H1图14.4 通道内外掠圆柱势流等流函数图图14.5 单侧送风与回风二、受限射流的势流模型在工程中常有这样的流动：一股气流从缝口喷出，直接喷向有限空间的对壁，流动方向发生转折后从有限空间的另一开口处流出。例如空调房间中气流的运动，飞机垂直起飞时喷向地面的气流，钢化玻璃器皿在钢化过程中受到的气流的冲击等均属于这种情形。假定由喷口喷向壁面的气流速度是均匀、无旋涡的，且流体没有粘性，则有限空间内的气体流动就可以用势流方程来描述。本节中通过两个例题介绍这一类势流的计算方法。
17、第十四章 流体动力学问题的有限差分解法 退 出返 回 第第14页页第一节 势流问题的数值计算 从喷口射出的气流称为射流。当射流进入到相当大的空间以致其运动不会受到限制和影响时，称为自由射流；当射流受到四周固体壁面的限制时，称为受限射流。下面以空调工程中两种侧面送风方式为例来计算其势流的流动情况。目的在于使读者进一步熟悉势流问题的数值解法以及掌握如何建立边界条件。例例14.1 一空调房间采用上侧进风、下侧回风的单侧进、回风方式送风，其简化模型如图14.5所示。假设进风处流速 =2m/s，气流密度为常数，且出风口处的流速均匀，试用数值方法确定其势流流型，取H=3m，L3.75m，H1=H2=0.1
18、5m。inw解解：坐标的选取如图14.5所示，流函数的控制方程为Laplace方程。下面着重讨论边界条件的设置方法。由于流场中任意点处的流函数代表了过该点的流线与参考流线间区域内的流量，故对于连续的、无渗透的固体表面，其上各点流函数必处处相等，对于两个不连续的、无渗透的固体表面，流函数必相差一个常数。第十四章 流体动力学问题的有限差分解法 退 出返 回 第第15页页第一节 势流问题的数值计算 又由于参考流线位置的选择具有一定任意性，因此若取一个固体表面的流函数为零，则另一表面的流函数即为流经这两个表面间区域的总流量。本例中取EF边界上的流函数为零，则在其余固体边界上，即在ABBCCD边界上的流
19、函数即为总流量 H1=20.15=0.3。所以本例的问题就相当于第一类边界条件的二维、稳态、无内热源的导热问题。若采用均分网格和区域离散A方法，则所有的内节点方程都可以表示成为：inw1,1,1,1,41jijijijiji（14.15）该问题的计算程序与二维稳态导热问题的程序类似，此处不再列出。流函数的计算结果如图14.6所示（取x=y=0.15）。可见，在势流的假设下，从送风口进入的流体流经房间后从出风口排出，在整个计算截面上并无回流与死滞区形成。第十四章 流体动力学问题的有限差分解法 退 出返 回 第第16页页第一节 势流问题的数值计算 例例14.2 在上例中，把出风口开在与送风口相对的
20、墙上，其余条件不变，试确定势流的流场。解解:从数值计算的角度，本例与上例的区别仅在于边界条件。按上例所述，如果取任两条相邻边界上的流函数为零，则与之相对另两条相邻边界上的流函数就等于流经该区域的总流量（图14.7）。计算所得势流的流函数分布如图14.7所示。与上例相类似，在势流的假定下，计算区域内无回流区形成。实际上，由于流体粘性的作用，在上述两例的情况下，流场内均应有回流形成。第十四章 流体动力学问题的有限差分解法 退 出返 回 第第17页页第一节 势流问题的数值计算 图14.6 单侧进出风的势流模型图14.7 一侧进风，另一侧出风的势流流型三、长方形截面通道内的充分发展层流流动在采暖通风等
21、工程领域中，广泛应用长方形截面的管道输送气体。本小节中以这一流动情况为例，介绍平直通道中层流充分发展流动的数值计算方法。所谓充分发展的管内流动是指截面上的速度分布不再沿主流方向发生变化的流动阶段。xazby图14.8 长方形截面通道 o第十四章 流体动力学问题的有限差分解法 退 出返 回 第第18页页第一节 势流问题的数值计算 在平直通道内的充分发展流动区域，Navier-Stokes方程可以简化为一个导热型的方程，从而使数值计算大为简化。如图14.8所示，设相应于x、y、z三个方向的速度矢分量为 、和 ，则在流动充分发展阶段，x、y方向的速度分量 、均等于零，且z方向的速度分量 沿流动方向不
22、再发生变化。于是仅有z方向的动量方程：xwywzwxwywzw222222zwywxwzpzwwywwxwwzzzzzzyzx（14.16）由于 、均等于零，故 ，且同一横截面上的压力可以认为是常数，于是上式可简化为：xwyw022zzpywxwzzdd2222（14.17）在充分发展流动的阶段，不再沿着流动方向发生变化，而是常数，因而 不再是被求解的变量。于是式（14.17）就是一个带源项的导热型方程，而就是源项。为更清楚地说明这一点，定义下列无量纲量：zpddzpdd第十四章 流体动力学问题的有限差分解法 退 出返 回 第第19页页第一节 势流问题的数值计算 其中D为通道的某一特性尺寸（如
23、高度a或宽度b或当量直径De），于是式（14.17）可以写成为下列无量纲方程：DyYDxXzpDwWz,dd2012222YWXW（14.18）（14.19）由于对称性，可取四分之一区域作为计算对象，此时上式的边界条件为：在固体边界上W0，在对称线上速度的法向导数为零。由此，平直通道内充分发展区域流场的求解相当于求解一个广义源项为1，广义导热系数为1的导热型问题。有关求解导热问题的一些处理方法都可以用来求解这一问题。在流体力学中，确定速度场的重要目的是要计算壁面上的切应力或阻力系数。在获得了无量纲速度场W以后就能够导出W与阻力系数f之间的关系。按阻力系数及Re的定义，第十四章 流体动力学问题的
24、有限差分解法 退 出返 回 第第20页页第一节 势流问题的数值计算 可以写出：2222ddDDWDwwzpDfReememmemwfRejijijijijimFFWW,其中Wm及 分别是截面上的无量纲平均流速及有量纲的平均流速。可知，一旦得出了Wm，之值即可求出。其中平均流速可由下式计算：式中 是第 个控制容积的流动截面积。计算时取整个流动截面为计算区域，采用B方法离散，节点数为3030。由此计算所得到的部分结果及文献中的推荐值如表14.2所示，可见，两者间的符合程度是很好的。而当a/b=1时截面上的流速分布如图14.9所示。jiF,ji,ba/数值计算值文献中的推荐值156.8857.01/
25、472.4373.01/678.3480.0表14.2 长方形截面通道的f值第十四章 流体动力学问题的有限差分解法 退 出返 回 第第21页页第一节 势流问题的数值计算 采用四分之一的流动截面作为计算区域时，可以把有限个节点应用于较小的计算区域内，但在计算时间上比采用整个截面为计算区域（总节点数不变）时更长。这是因为采用全截面计算时，四条计算边界上速度均为零值，边界条件的作用最强，容易促使代数方程迭代收敛。采用四分之一区域为计算区域时，两条边界上具有的均为“绝热型”的边界条件，即速度的法向导数为零，从边界条件传人到计算区域内部的信息较少，代数方程的迭代不易收敛。对 =1/4的情形，节点数取为3
26、030，收敛条件取绝对偏差 10-5，采用逐线迭代法，当以全截面为计算区域时，只经42次迭代即收敛，而采用14区域时需经383次迭代才收敛。但若保持相同的网格疏密度，则采用全截面为计算区域时，由于节点数几乎增加四倍，迭代求解次数会大幅度地增加，总的所需计算时间就比以14区域作为计算对象时要多。ba/第十四章 流体动力学问题的有限差分解法 退 出返 回 第第1页页第二节 回流流动问题的数值计算 有回流的流动现象是用椭圆型的控制方程来描述的，这类椭圆型控制方程的求解是计算流体力学中的一个基本课题。本节讨论的重点放在不可压缩流体场的求解上。求解有回流的流场时，可以用原始变量法，也可以用涡量流函数法。
27、本节中主要介绍涡量流函数法，在末尾对原始变量法作简单的介绍。在对Navier-Stokes方程进行数值求解时，所遇到的困难主要来自对两个一阶导数项对流项 （以x方向的对流项为代表）和 项的离散。在涡量流函数法中，压力梯度项已经从控制方程中消去，建立离散方程的关键在于对流项的离散，同时如何建立合适的涡量边界条件，也是涡量流函数法计算的重要内容。而在介绍原始变量法时重点在于 的离散及所得离散方程的求解步骤上。无论是涡量流函数法还是原始变量法，控制方程都是非线性的，因而整个问题的求解必然采用迭代方法。至于在每一个迭代的层次上，代数方程的求解又有直接解法与迭代法两种，有关这方面的内容可以参阅相关文献。
28、xwwxxxpxp第十四章 流体动力学问题的有限差分解法 退 出返 回 第第2页页第二节 回流流动问题的数值计算 一、对流扩散方程的离散格式无论是原始变量法还是涡量流函数法，控制方程中都有对流项，对流项的存在是对流问题区别于导热方程的关键，因而可先讨论对流项的离散方法。由于整个对流换热问题数值解的结果还与扩散项的离散方法有关，而且有一些离散格式是把对流项与扩散项合在一起来定义的，所以在本小节中将把对流项与扩散项的离散联合起来考虑。本小节的讨论将针对最简单的模型方程，即稳态、无源项的一维对流扩散方程进行。因为这一方程既具备了对流扩散方程的基本特性，同时又具有精确解，在数值计算过程中便于把数值解与
29、精确解相比较。对流项差分格式的定义可以通过两种方式来进行：（1）规定控制容积界面上函数数值的方式。由于在进行数值计算时，所求解变量之值都是存放在节点位置上的，而采用控制容积积分法时，积分过程中要用到面上所求解变量的值，故在此方式中需要确定从节点值获得界面值的方式（也相当于选取型线），不同的插值方式代表不同的离散格式。第十四章 流体动力学问题的有限差分解法 退 出返 回 第第3页页第二节 回流流动问题的数值计算（2）规定节点上一阶导数的离散形式，由Taylor展开法来导出离散方程时采用这种定义方式。下面主要采用前一种方式来介绍五种常用的离散格式。（一）中心差分该格式规定界面上的值及界面上的一阶导
30、数均按分段线性的型线确定。将一维、稳态、无内热源方程：xxwxdddddd（14.20）对图13.6所示控制容积P作积分，对流与扩散项均采用分段线性的型线，有：wwwWeeeEwwweeePwxwxxwxw21212121（14.21）w/x PPEEWWaaa（14.22）把通过界面的流量 记为F，界面上扩散阻力的倒数（扩导）记为D，上式可写成：第十四章 流体动力学问题的有限差分解法 退 出返 回 第第4页页第二节 回流流动问题的数值计算 式中各系数为：式（14.23）表明，如果在数值计算过程的每一迭代层次中，连续性方程都能得到满足，则 仍等于各邻点系数之和。应指出，本节中一维模型方程（14
31、.20）的五种常用离散格式均类似于式（14.22），仅系数 的表达式不同。因此在讨论其它几种离散格式时仅介绍系数 的计算式11,()22EeeWwwPEWewaDFaDFaaaFFPaEWaa、EWaa、（14.23）注意到 ，这是以 为特性尺寸的Peclet（贝克莱）数，记为 。如果把 看作是运动粘性系数 ，则 也为网格Reynolds数的表达式。因而本节以下均把 称为网格Peclet数。对均分网格，式（14.22）可化为：/F Dw x xPxwReWEPPP2112112P可见，若 2时，则系数aE0，此时数值解会出现振荡，可举例说明。第十四章 流体动力学问题的有限差分解法 退 出返 回
32、 第第5页页第二节 回流流动问题的数值计算 例例14.3 在一均分网格中求解一维模型方程，设已知：=200，=100，试对 =0，l，2，3及4等五种情形按中心差分计算 值，并讨论所得之结果。解解：由上式可得EWPP2211211WEPPP图14.10 中心差分特性分析按此式计算的结果及精确解之值示于图14.10及表14.3中（一维模型方程的精确解将在下面给出）。由图可见，当大于2时，之值即小于，这是不合理的。因为这是一个无源项的对流扩散问题，之值必定介于与之间。这种在物理意义上不真实的数值解是由于的系数aE在大于2时小于零而造成的。与导热问题离散方程系数的物理意义相类似，aE、aW代表了P的
33、邻点E与W通过扩散与对流对P点所产生的影响，其值均需大于零，系数的负值会导致物理意义上不真实的解。第十四章 流体动力学问题的有限差分解法 退 出返 回 第第6页页第二节 回流流动问题的数值计算 在将对流扩散方程进行离散时，可认为对流项中的流速都是已知的，且在迭代式求解过程中取上一层次迭代所得之值，这种做法称为局部线性化，要求各个控制容积界面上的 小于2，即对各控制容积的尺度进行限制，从而导致在计算时须采用较多节点数，使所需的计算机内存与计算时间大大增加。P（二）迎风差分迎风差分的离散格式可以克服对流项采用中心差分时因 大于2而引起的上述困难。用迎风差分构造对流项的离散形式时，永远从上游（迎着来
34、流）去获得所必须的信息。采用Taylor展开法在i点上来离散通用变量的对流项时，其形式如图14.11所示，可写成：P11,0,0iiiiiiiidwwwdxxwwx（14.24）第十四章 流体动力学问题的有限差分解法 退 出返 回 第第7页页第二节 回流流动问题的数值计算 采用控制容积积分法时，迎风差分规定界面上的变量按下列方式取值：在w界面上：；式（14.24）称为第一类型迎风差分，式（14.25）称为第二类迎风差分。为表达上的简洁，并便于编制程序，采用第二类迎风差分时界面上的对流项 可写成如下紧凑格式：（14.25）Peew,0Eeew,0Wwww,0Pwww,0()w 在e界面上：；0,
35、max0,maxeEePeeeFFFw0,max0,maxwPwWwFFw（14.26a）（14.26b）类似地有：将式（14.20）对控制容积作积分时，如果对流项采用第二类迎风格式，扩散项采用中心差分（即分段线性型线），就得迎风差分的离散方程，形式仍如式（14.22）所示，其中系数计算式则为：max,0max,0EeeWwwPEWewaDFaDFaaaFF（14.27）第十四章 流体动力学问题的有限差分解法 退 出返 回 第第8页页第二节 回流流动问题的数值计算 上式表明，由于系数aE、aw永远大于等于零，因而采用迎风差分的离散方程不会出现物理上不合理的解，即无论网格Peclet数多大，数值
36、解都不会出现振荡。迎风差分虽有不会振荡的优点，也有其不合理之处，表现在：（1）迎风差分以界面上流速是大于零还是小于零来决定界面上的变量值，显得过于简单，界面上的变量值在一定的流速范围内还与流速有一定关系；（2）在迎风差分中不论 的大小，扩散项一律取中心差分，即一律以界面两侧节点的函数值之差除以两节点的间距作为界面上变量的梯度。实际上只有当 之值较小时，这样做才是合理的，这可以从后面要介绍的一维模型方程的精确解的图示中看出。为了克服这些不足，Spalding提出了混合格式。PP第十四章 流体动力学问题的有限差分解法 退 出返 回 第第9页页第二节 回流流动问题的数值计算 图14.12 系数aE，
37、aW 的关系 aE、aW具有影响系数的意义。如图14.12所示，节点 对 点的作用表现为 ，而 点对 点的作用则表现为 ，由于该两点共享一个界面，其间的扩散阻力及界面上流量的绝对值都相同（但方向相反），因而 及 间必有一定关系。所以有可能通过对规定aE的表达式而定义一种差分格式，以中心差分为例：1ii iaEi1i1iaW iaE1iaW（三）混合格式eeeeEPDFDa21121eeEPDa211wwwwWPDFDa21121wwWPDa211对i点 ，或对i+1点 ，或i1iDiDiDwe)1()(FiFiFwe)1()(PiPiPwe1FiaiaEW)()1(注意到 点的e界面即为 点的
38、w界面，故 ，于是可得：（14.28a）PiDiaiDiaeEwW11或（14.28b）第十四章 流体动力学问题的有限差分解法 退 出返 回 第第10页页第二节 回流流动问题的数值计算 可以证明式（14.28）对迎风差分及下面要介绍的其它几种差分格式也是成立的。这样，定义一个差分格式就可以通过规定 的表达式来进行。Spalding提出的混合格式的定义为：eEDa0,2eEeDaPeeEePDaP211,22eeEePDaP,20,21,maxeeeEPPDa（14.29a）（14.29b）（14.29c）（14.30）即：由这一定义式可知，当 时，混合格式采用中心差分，对流与扩散的作用同时予以
39、考虑，式（14.29b）中的“1”代表了扩散的影响，而 则为对流项的作用。当 时，扩散项的作用可以忽略，表现为在式（14.29a）及（14.29c）中略去了式（14.29b）中的“1”这一部分。当 时，E点在P的下游，因而E点对P的影响可以忽略不计，而当 时，E在P的上游，它对P的影响完全来自对流。2P12eP2P2P2P第十四章 流体动力学问题的有限差分解法 退 出返 回 第第11页页第二节 回流流动问题的数值计算（四）指数格式容易证明，若一维模型方程（14.20）具有下列边界条件：并假定、均为已知常数时具有精确解：0,0 xLLx,w wLPPLxPeeeL,1exp1exp00（14.3
40、1）；（14.32）这一精确解的图线如图14.13所示，可见，当 时，成直线分布，这是一个纯扩散问题；在 时，曲线仍与直线分布相接近；随着 的增加精确解越来越呈现出边界层类型问题的特性，即在 至 的大部分范围内上游的值 占优势，仅在靠近区间端点的一个薄层内 值才由来流值迅速增大到边界值。由图还可以看出，在 的绝对值较大时，界面上的扩散作用（即曲线在界面上的斜率）已比较小，当 时，区间中截面上的扩散作用已接近于零。0eP2ePeP0 xLx eP10eP第十四章 流体动力学问题的有限差分解法 退 出返 回 第第12页页第二节 回流流动问题的数值计算 为便于将其它格式与精确解进行比较，我们把这一精
41、确解也表示成式（14.22）的形式，这样格式的比较可通过比较 来进行。为此把式（14.20）改写成为：上式括号中 的代表了对流项的作用，而 则是扩散作用的度量。式（14.33）表明，对于一维、稳态、无内部源项的对流扩散问题，通过任一截面由对流、扩散所传递的总量是相等的。将式（14.33）对控制容积P作积分（参见图13.6），得：eEDa0ddddxwxw xddwexwxwdddd（14.33）（14.34）将精确解的表达式（14.32）代入式（14.34），并注意到：对e界面对w界面 eeeELPPxwP,0wwePLPPxwP,0第十四章 流体动力学问题的有限差分解法 退 出返 回 第第1
42、3页页第二节 回流流动问题的数值计算 经整理可得：1expexp1exp1exp11expexpwwwWeeEwweeePPPFPFPFPPF（14.35）为了将此式写成式（14.22）的形式，在上式左端方括号内加入 项，其中 原有的两项结合各自构成 的系数，于是有：式（14.36）称为指数格式，它是相应于一维模型方程精确解的格式，对于其它情形，它仅是一种近似的离散格式。)(wweeFFFFweFF与wE及weWEPeewWeeEFFaaaPPFaPFa,1expexp,1exp（14.36）（五）乘方格式由于指数格式在计算上很费时间，Patankar提出了与指数格式十分接近的乘方格式，其计算
43、工作量比指数格式要小得多。乘方格式的定义为：0,101.01,1001.01,010,1055eEeeeEeeeeEeeeEeDaPPDaPPPDaPPDaP（14.37）第十四章 流体动力学问题的有限差分解法 退 出返 回 第第14页页第二节 回流流动问题的数值计算 还可以写成：至此，已介绍了五种关于对流扩散方程的差分格式。其中混合格式、指数格式与乘方格式在给出定义时，对流与扩散作用是联合考虑的，而中心差分与迎风差分中，扩散项与对流项是分别定义的，且扩散项均取为中心差分。这五种格式的定义都可以用 的计算式来表示，在表14.4中作了归纳。进一步的理论分析表明，在这五种格式中，实际上只需规定 时
44、的 ，的部分可由格式的定义及其特性推演而得。在图14.14中画出了 （）的曲线。如图所示，乘方格式与指数格式十分接近，而对所讨论的一维模型方程，混合格式可以看成是精确解的一种很好的逼近。0,max0,1.01max5eeeEPPDa/EeaD0P/EeaD0P/EeaD0P（14.38）最后要指出，所介绍的五种格式中，对流项的中心差分不具有迁移特性，其所构成的离散方程只能用于 的情形，当 时数值解就要振荡，这称为条件稳定。2P2P第十四章 流体动力学问题的有限差分解法 退 出返 回 第第15页页第二节 回流流动问题的数值计算 而迎风差分、指数格式及乘方格式都具有迁移特性，由它们所构成的离散方程
45、无论 多大都不会引起数值解的振荡，这种特性称为绝对稳定。混合格式在 时虽然不具有迁移特性，但在此区域数值解不会发生振荡，因而总体上也是绝对稳定的。进一步的分析可以证明，凡是对流项的差分格式不具有迁移特性的，这种格式一定是条件稳定的；凡对流项具有迁移特性的，其对流扩散问题的差分方程一定是绝对稳定的。P2P采用涡量流函数法时，通过引入涡量与流函数，把压力梯度项从动量方程中消去，从而避免了求解压力场的一些困难。本节中将对直角坐标中的对流换热问题导出涡量流函数法的控制方程。为使讨论具有一定的通用性，考虑自然对流换热的情形。图14.15 导出涡量流函数的直角坐标wxxywyg二、涡量流函数法的控制方程选
46、取坐标如图14.15所示（重力与x轴成 角）。此时稳态的对流换热的控制方程可以表示成为：第十四章 流体动力学问题的有限差分解法 退 出返 回 第第16页页第二节 回流流动问题的数值计算 0yxwwxy2222cosywxwgxpywwxwwxxxyxx2222sinywxwgypywwxwwyyyyyxyTyxTxyTwxTwcyxp（14.39a）（14.39b）（14.39c）（14.39d）上述方程组中，流体的密度是温度的函数。但是在自然对流情形下，密度的变化对流动影响最大的是体积力项中的部分，由于正是密度的变化引起了浮升力，从而驱使流体流动，因此如果把密度随温度的变化仅限在体积力项，既
47、抓住了物理现象的本质又可以使计算工作得以简化。从而引入了Boussinesq假设：在动量控制方程中，除体积力项中的密度外，其它的物性（包括其它项中的密度）均为常数，并且密度与温度的关系满足：第十四章 流体动力学问题的有限差分解法 退 出返 回 第第17页页第二节 回流流动问题的数值计算 rrTT 1rTrrTcoscos()cos()cossinsin()sinrrrrrrrrrggTTggTTgggTT (cos)(sin)effrrppgxgy（14.40）（14.41）其中，为参考温度，为与 相应的密度，为体积膨胀系数。于是式（14.39）中的体积力可以写成：定义：则有：称为有效压力或折
48、算压力。将式（14.41）代人式（14.39b）、（14.39c），令非体积力项中的 均为 ，可得：cos,sineffeffrrppppggxxyyeffpr2222cos1ywxwTTgxpywwxwwxxrreffrxyxx（14.42a）第十四章 流体动力学问题的有限差分解法 退 出返 回 第第18页页第二节 回流流动问题的数值计算 2222sin1ywxwTTgypywwxwwyyrreffryyyxeffpprrxwywyx（14.42b）以下为书写简便，仍记为 ，记为 ，记为 ，并定义：（14.43）将式（14.42a）对y求导，减去（14.42b）对x求导后的结果，有：ywxw
49、ywxwxwywxwywxwywywxwywxwyxwwywxwxwwxwxwywwywywyxwwxwywyxyxyyxxyxyyxxyyyyyxyxxyxyxxxx222222左端第十四章 流体动力学问题的有限差分解法 退 出返 回 故得：第第19页页第二节 回流流动问题的数值计算 222222222222222222sincossincossin1cos1yxxTyTgxwywyxwywxxTyTgxwyxwxxTgyxpywyywxyTgyxpyxyxyyxx右端xTyTgyxywxwyxsincos2222（14.44a）其守恒形式为：xTyTgyxywxwyxsincos2222（
50、14.44b）式（14.44a）或（14.44b）称为涡量方程。对于不计及体积力的强制对流，只要令上式右端的体积力项等于零即可。将流函数的定义式代入涡量定义式，可得流函数的微分方程：2222yx（14.45）第十四章 流体动力学问题的有限差分解法 退 出返 回 第第20页页第二节 回流流动问题的数值计算 而能量方程（14.39d）的守恒形式为：式（14.44b）、（14.45）及（14.46）就构成了用涡量流函数法来计算流动与换热问题的控制方程。可定义下列无量纲形式：yTcyxTcxyTwxTwppyx2/,/,/,/(/),/(/),/(/),xxyyrrxx Lyy LwwLwwLTTLT
