第十一章教案

1、第十一章 流体力学的数值研究方法计算流体力学是当今流体力学的三大分支之一。其英文名称包括Computational Fluid Mechanics，Computational Fluid Dynamics(CFD)和Numerical Simulation of Fluid Flow。l 流体力学的早期分支n 理论流体力学精确、严谨，但只能求解简单流动问题。n 实验流体力学能解决复杂几何边界的流动问题，但费时、费工、费物力，又比尺效应和速度效应。l 现在的分支n 增加了计算流体力学分支从上世纪60年代末开始，用数值方法来求解流动问题。l 计算流体力学的特点n 优点u 可解决理论上解不出，实测有

2、比较困难的流动问题u 解决前两个分支的大部分问题u 揭示新的流动现象u 优化设计n 缺点u 对精确方程的简化过程将引入误差u 流体流动的性质和趋势不能象解析解那样明显表达出来u 计算解的精确性有待于实验的验证l 计算流体力学的方法n 有限差分法n 有限元法n 有限分析法n 边界元法n 谱方法本课程主要介绍前两种方法。l 有限差分法（FD）的特点n 用差商代替为偏导数n 求解网格上的函数值n 比较成熟、简单l 有限元法（FEM）的特点n 普适方法（模块性较好）n 适用于初值问题（椭圆方程）和边值问题（近似抛物线方程），解决初值问题不如差分法n 数学物理基础是能量极值原理和分段逼近11.1 有限差

3、分法11.1.1 微分与积分的数值运算计算机的计算结果均为近似解（或数值解、离散解）。差分法的基本思想为：将微分方程按一定的差分格式离散为差分方程，然后求解其数值解。l 微分插值公式、差分格式及其精度分析首先介绍插值公式，一般有线性插值格式、二阶插值格式和中心插值格式。n 线性插值格式将函数在某一点出展开，有（11-1） 式中，为余项，记为也可以记为如果取，则。由于余项，此式的精度为二级。若在轴上划分了n个网格，网格节点从，相应的函数值为，。由式（11-1）可得（11-2） 由此得（11-3） 因此有（11-4） n 二阶插值格式将式（11-2）等式右边的第三项写出，有（11-5） 由此式可分

4、别得到1、2点的表达式。（11-6） （11-7） 将式（11-6）乘以4，减去式（11-7），可得解得（11-8） 将式（11-6）乘以2，减去式（11-7），可得解得（11-9） 将式（11-8）和式（11-9）代入式（11-5），可得（11-10） n 中心插值格式增设一个-1点，有（11-11） 将式（11-6）减去式（11-11），可得（11-12） 将式（11-6）等号右边取5项，式（11-11）等号右边取5项，两式相加，可得二阶导数的差商。（11-13） 此式是中心差分格式，有二阶精度。将式（11-12）和（11-13）带入到式（11-5）可得（11-14） 因为在处查值，用到了

5、相邻两点的函数值，故称其为中心插值公式。从上式可以看出它有三阶精度。l 积分的数值运算公式及其精度为了简单起见，采用线性插值公式来表示，可进行积分运算。（11-15） 如果采用的二阶插值公式或中心插值公式，以上定积分的数值解的精度要高一些。11.1.2 由常微分方程表示的流动问题数值解l 常微分方程的数值解常微分方程的一般形式如下式中，。将式改变形式，可得（11-16） 式中，。用取代式（11-6）中的，i取代1，式（11-6）可改写为（11-17） 将式（11-16）对求导，可得的二阶导数。略掉式（11-16）中的高阶小量，可得到计算公式（11-18） 这就是所谓的Euler公式，具有二阶精

6、度。根据这一公式，已知和值，可以算出值，最后计算出值，由此类推，可以计算出x计算域上的各y值。为了得到更高阶精度的计算公式，将式（11-11）改写为由此可得（11-19） 将式（11-19）与（11-17）相加，等式两边除以2，可得（11-20） 为了找出的表达式，仿照式（11-17）按泰勒级数展开，有将上式代入式（11-20）等式右边的第三项，可得由此可得计算公式（11-21） 与Euler公式相比，式（11-21）有更高的精度。但是不能被直接解出，公式是一个关于的隐含关系式。尽管如此，可以用迭代法解出所有的y值。式（11-21）还可以用于常微分方程组，且只有一个独立变量的场合。例如，若有则

7、根据式（11-21），可得（11-22） （11-23） l 非稳态管流非稳态管内流动如图11-1所示。管内的流量是随时间变化的，有阻力引起的压力降也是变化的。根据动量积分方程，有（11-24） Lp1p2图11-1 管内非稳态流动由于管内流速可近似用平均流速v代替，式（11-24）中等式右边最后一项为代入式（11-24），等式两边同除以，可得（11-25） 由于，式（11-25）可改写为（11-26） 对于过渡区流动或粗糙管流动 （11-27）对于光滑管流动（11-28） 注意，计算的公式里含有速度v。因此，计算v的步骤如下：n 根据初始值v计算出Re，若是粗糙管计算出；n 根据Re计算l，

8、若是粗糙管根据Re和计算l；n 由式（11-26）得出的表达式；n 由式（11-18）计算出新的v值；n 重复以上步骤直到时间t达到要求值为止。例题11-1 管内非稳态到稳态流动的速度变化计算一个大水池向一条长200ft、直径8in的水平圆管供水，水池的自由液面距管道进口中心线50ft高。在管道末端装有一台离心泵和一个旁通阀，泵前也有一个阀门。当泵操作时旁通阀是关闭的，流体作稳态流动，流速为19.81ft3/s。当迅速将泵前的阀门关闭，同时将旁通阀门打开时，求达到新的稳态流动状态过程中，流量随时间变化的关系。管道是市场销售的钢管，忽略热损失。解：考虑大水池液面到管道进口的Bernoulli方程

9、：得（a）由式（11-26）可得（b）根据初始，可知用式（11-27）可计算l。方程的形式为：v的初始值由下式计算。，取时间网格为50，计算过程时间期间足够大，为18s，则时间步长为s。用Euler公式计算，结果如表11.1所示。表11.1 时间步长为s时的非稳态流动计算结果No.t(s)v(ft/s)Rel1056.7520.36643.303.481060.0152130.73236.662.651060.0152941.4630.162.011060.0153952.8325.701.611060.0154865.8624.041.471060.01552711.7123.861.461

10、060.01553进一步，取时间网格数为100，则时间步长为s。用Euler公式计算，结果如表11.2所示。l 高尔夫球的运动轨迹可将高尔夫球的运动看作是球形颗粒在空气中的运动。那么，可以计算出其运动轨迹。球形颗粒在空气中运动，用表示其运动速度。颗粒所受的力包括重力和空气对它的曳力。设空气的速度表11.2 时间步长为s时的非稳态流动计算结果No.t(s)v(ft/s)Rel1056.7530.36645.123.071060.0152550.73238.482.541060.0153071.4631.462.011060.0153992.8326.261.631060.01547115.862

11、4.141.481060.015521311.7123.861.461060.01553为，球形颗粒与空气的相对速度为，球形颗粒的直径为D，质量为m。忽略浮力的影响，根据牛顿定律有（11-29） 令，则式（11-29）可改写为（11-30） 式（11-30）在迪卡尔坐标下的分量表达式如下：（11-31a） （11-31b） （11-31c） 式中，。根据以上这些公式，如果已知颗粒的初始位置和初始速度，就可以计算出颗粒的运动轨迹。例11-2 高尔夫球的运动轨迹已知高尔夫球的初始速度为120ft/s，与地平面（xy平面）的夹角为30，质量为1.5oz，直径为1.75in。假定高尔夫球没有旋转，取空

12、气粘度。解：首先计算t（a） 对于这个问题，由式（11-31）可得：（b） 式中u、v、w表示风速。对于阻力系数CD，可根据Re选取。球形颗粒的临界雷诺数Recr为9.0104，当ReRecr时，当ReRecr时，。因此，用中心差分公式计算，有（c） 令式中（d） 式（c）可改写为（e） 进一步，可以列出计算位置坐标的公式。（f） 根据题意，风速为零。因此可以认为高尔夫球在xoy平面内运动，只须计算x和z坐标就行了。如果不考虑阻力的影响，可以算出球的飞行时间为3.73s。因此可初定飞行时间为4s，划分100网格，时间步长为0.04s。由于开始时高尔夫球在地上，可令，通过计算，结果如表11.3所

13、示。表11.3 高尔夫球的运动轨迹No.t ( s )x ( ft )z ( ft )No.t ( s )x ( ft )z ( ft )10.044.1492.370112.00186.949.3420.2020.6211.26122.20202.146.9330.4040.9121.08132.40216.643.3540.6060.9129.47142.60230.638.6650.8080.6236.45152.80244.132.8861.00100.142.07163.00257.026.0671.20119.146.30173.20269.418.2481.40137.349.0

14、8183.40281.49.46791.60154.650.47193.60292.8-0.288101.80171.150.53l 压力容器的排空问题图11-2 压力容器的排空如图11-2所示，以密闭圆筒容器装有一定量的液体。容器顶部高距放液管中心高H，放液前液面高度为，压力为。对于液体，可用柏努利方程描述，由于槽内流动为准静态过程，可认为是无粘流动。对于液体上的密闭气体，可假定为等熵膨胀过程。根据连续性方程，有：（11-32） 而液面和喷嘴口之间的能量关系为：（11-33） 设初始液面高度为hi，气体压力为pi，容积为Vi，可得气体方程（11-34） 式中k为比热比。设容器中气体的质量为M

15、，则气体的密度。于是有（11-35） 由于，则式（11-34）可改写为（11-36） 将vfs和式（11-36）代入式（11-33），可得：最后可得常微分方程的表达式。（11-37） 此时可写成的形式，因此可以用式（11-18）或式（11-21）求解之。l Blasius边界层方程的数值解在第十章第二节介绍了描述边界层流动的微分方程式（10-24）的推导过程。该方程不能直接求解，现介绍数值解的方法。式（10-24）可改写为（11-38） 边界条件为（11-39） 由于式（11-38）是一个三阶微分方程，h 的取值范围是0到，不是简单的初值问题。需要进一步讨论，可以将式（11-38）改写为：（1

16、1-40a） （11-40b） （11-40c） 边界条件为，。由于是未知的，故还不能求解上式。n 解法一：可预设一个，然后求解式（11-40），算出的不一定满足。此时改变，再算，直到为止。n 解法二：对于两个不同的假定值，可以计算绘出随变化的曲线，如果，则可以按下式计算新值。（11-41） 根据计算绘出随变化的曲线，如果，则令，反之，令，然后再根据式（11-41）计算出新值。以此类推，直到后一次的与本次的之间的相对误差小于等于给定值为止。（11-42） 在实际计算中，不可能把取到无穷大。对于各，按以下的方法可以找到相应的。根据式（10-15），有（11-43） 当时，。可近似认为在边界层的外

17、边界上，这时对应的近似为。则根据式（10-12），有（11-44） 取，就可以根据式（11-41）、（11-42）找出，然后再用下面各式计算出式（11-40）各函数值。（a）（b）（c）例11-3 布拉修斯方程的求解划分200个网格，设，空气流速为u0 = 20 ft/s，温度为20，求x = 3 ft时的速度分布u(y)。解：当x = 3时，由式（11-13）可知y与h 的对应关系，即。解出不同h 对应的G值，就可以得到速度的分布u(y)，因为由式（11-43）可知，。已知，经过试算，得。取h 的步长为0.03，按照式（a）、（b）、（c）可以计算出对应的f、G和H的值。列表如下：网格序号h

18、fGH10000110.300.0149280.0995130.331549210.600.0596850.1987880.329930411.200.2378030.3936510.316862611.800.5293740.5749500.283715812.400.9224200.7297410.2291331013.001.3975690.8473810.1621651213.601.9311960.9249860.0983171414.202.5007090.9685920.0502561614.803.0889080.9892030.0214441815.403.6852800.9

19、973340.0075942016.004.2846391.0000000.002224根据以上数据，可得x = 3 ft时的速度分布曲线。现在确定边界层的厚度。因为，当或=0.99，根据上表的计算值可知，带入式（11-13）可得：故有11.1.3 由偏微分方程表示的稳态流动问题的数值解l 概述采用有限差分法摸拟流体流动问题的基本过程包括：n 建立模型方程，确定相应的初始条件和边界条件；n 将求解区域划分成网格区域并确定计算节点，即区域离散化；n 利用差分公式代替模型方程的各微分项，使微分方程转化为由节点流动参数所表示的代数方程组，即方程离散化；n 采用适宜的数学方法和计算程序求解该代数方程组

20、，即算法设计；从而获得各节点上流体速度及相关参数的近似值。l 势流的数值解n 建立模型方程正如第七章讨论的那样，对于平面势流，流函数满足拉普拉斯方程，用式（7-14）表示。即图11-3 计算网格系统1，1 2，1 3，1 4，11，2 2，2 3，2 4，21，3 2，3 3，3 4，31，4 2，4 3，4 4，4n 计算域的划分将计算域划分成网格，对于矩形计算域，网格如图11-3所示。其中，网格线的交点称为节点，计算节点为（2，2）、（3，2）、（2，3）和（3，3），其余的是边界点。一般边界点的参数值或者以参数值为变量的函数关系是已知的，需要求解的是计算节点的参数值。n 方程离散化对于式

21、（7-14）表示的势流流动，可用中心差分格式即式（11-13）代替其中的二阶偏导数，最后得：（11-45） 式中， 中心节点处的流函数值； 中心节点相邻右边节点处的流函数值； 中心节点相邻左边节点处的流函数值； 中心节点相邻上边节点处的流函数值； 中心节点相邻下边节点处的流函数值。上式可以化简为：（11-46） 式中。对于如图11-3所示的计算域的每个计算节点，可按式（11-46）写出具体表达式。（11-47a） （11-47b） （11-47c） （11-47d） 如果边界条件是第一类边界条件，则式（11-47）是一个只有4个未知数的代数方程组，有4个方程，故此代数方程组有解。这样，可以根据

22、节点（1，1）、（1，2）、（1，3）、（1，4）、（2，1）、（2，4）、（3，1）、（3，4）、（4，1）、（4，2）、（4，3）和（4，4）上的已知值，求解出计算节点上的值。n 求解代数方程组代数方程组的求解方法可分为直接解法和迭代法两大类。对于大型代数方程组，直接解法诸如高斯（Gauss）消元法及其变形等应用较少，因其所需的计算量和内存都较大，而且还存在由于计算机舍入误差带来的累积误差等问题。迭代法在所需计算内存和时间方面都可望比直接法省，程序设计也简单，而且更重要的是，实际流动问题的模型方程大多数是非线性的，离散后得到的是非线性代数方程组，方程中的系数可能都是未知变量的函数，因此，整

23、个问题的求解本身就带有迭代的性质：即首先假定流场内各节点上的未知变量值，并根据假定值计算方程系数，然后求解方程获得未知变量在各节点上的改进值，再由改进值重新计算方程系数，并获得未知变量进一步的改进值。如此反复，直到所有节点上未知变量前后两轮的解不再变化或仅在允许的误差范围内变化时，才算获得代数方程的收敛解。因此，求解实际流动问题（尤其是非线性问题）的离散方程主要采用的是迭代法。l 粘性不可压缩层流的数值解在第5章和第9章讨论过圆形直管内流体的层流流动，其速度分布是轴对称的，并且只与半径有关。因此，可用精确的解析方法求解该问题。然而，对于矩形管中的层流流动问题，由于速度分布于横截面上的两个坐标有

24、关，虽然可用N-S方程来描述，但解析法已无能为力。然而，用数值研究方法可以解决这个问题，这体现了计算流体力学的优越性。下面从简单的流动着手，讨论矩形管内层流流动的数值解法。n 模型方程的建立取管道的轴向坐标为x，其它坐标如图11-4所示。由于是定粘度、不可压缩、稳态层流流动，、为常量，由式（6-26）可将x方向的动量方程可简化为：图11-4 矩形管及计算坐标zyAB（11-48） 式中，由于是直管流动，为常量，式（11-48）是泊松方程。解这个方程可得到轴向速度在横截面上的分布值。为了使计算值具有代表性，可将式（11-48）无量纲化，这样得到的结果对具有相似横截面的管道（例如，）有普遍的指导意

25、义。定义以下的无量纲变量：（11-49a） （11-49b） （11-49c） 将式（11-49）带入式（11-48）可得无量纲方程。（11-50） 式（11-50）是一个泊松方程，它描述了矩形管横截面上无量纲速度之间的关系。无量纲坐标的变化范围为（a）（b）（c）对于固体边界，；对于对称边界，同理可得，。根据计算出的无量纲速度分布，可以计算出矩形管道流的平均无量纲速度。（d）根据式（11-49c）可得平均速度与平均无量纲速度之间的关系。（e）由此可知，只要已知无量纲速度分布，就可以根据式（d）计算出平均无量纲速度，若改变流量qv，可以根据式（e）计算出压力降dp/dx，将其带入式（11-49

26、c）可以计算出实际流速分布v(y, z)。需要指出的是，以上结果只能用于层流的计算。n 计算域的划分由于流动相对x轴是对称的，所以，可选图11-4中四分之一的部分，即第一象限的这部分过流面积作为计算域。边界条件如图11-5所示。n 方程的离散化n 求解代数方程组例11-4 计算矩形截面管的流速分布计算矩形管道的流速分布，已知A = B = 1 ft，采用每条边10个网格，则。用式（c）迭代计算223次达到计算精度。计算结果如下图所示。0.000 0.000 0.0000.000 0.000 0.0000.1130.1110.1010.0830.0520.1930.1880.1700.1730.

27、0830.2450.2390.2150.1700.1010.2740.2660.2390.1880.1110.281 0.274 0.2450.193 0.1130.00011.2 有限元法11.2.1 有限元法的基本原理有限元法的理论基础包括泛函分析，变分法，函数逼近理论和矩阵论等。本节只讨论有限元的基本原理及其应用等有关内容。l 有限元的基本思想n 基本原理变分原理和伽辽金法n 具体方法u 采用差分离散处理网格的思想：把要分析的区域分割成很多很小的子域单元（Element），单元的顶点称为节点。这些单元的集合代表计算区域。u 在单元内用一个简单的函数插值逼近函数来逼进单元的真解。避免了伽辽

28、金法和利兹法在选取满足全域的函数时的困难。u 选取插值函数后，用变分法代入泛函表达式，再利用泛函在子域内的极值条件，建立单元方程组，如果定解问题找不到或不存在变分泛函，用加权余量法找出方程组。u 将每个单元的方程联合起来，组成整个计算域的总方程组，得到总“刚度”矩阵。u 求解方程组，得到各节点的函数值。l 有限元法和利兹法（伽辽金法）的比较n 有限元法与利兹法（伽辽金法）的关系如图11-6所示。利兹法的基础是变分法，伽辽金法的基础是加权剩余法。有限元方程组是按变分原理或者加权余量法求得的。其分析过程与利兹法或者伽辽金法很相似，但是存在着某些差异。n 有限元法中，变量的总体泛函表达式由局部反函表

29、达式汇集而成，即作用在全域上的总体泛函是由作用在分散子域上的泛函汇集而成。n 采用伽辽金法时，可将有限元插值函数作为权函数处理。n 总体微分方程由局部微分方程汇集而成。n 在有限元法中，近似函数的选取与利兹法或伽辽金法类似。基本差异在于：近似函数建立在子域上，避免了满足总体边界条件的全域的近似函数的选取。图11-6 有限元法与利兹法（伽辽金法）的关系积分方程（近似）边界元法有限基本解法物理定律偏微分方程（定解问题）（近似）有限差分法配置法矩量法最小二乘法伽辽金法加权余量法精确解析解泛函表达（变分法）欧拉方程利兹法（近似）（离散）（自对半算子）有限元法11.2.2 一维有限元法l 有限单元的划分

30、n 一维：直线n 二维：三角形，矩形，任意四边形n 三维：四面体，矩形立方体，六面体，环形体l 插值函数n 三角形式的多项式n 指数形式的多项式n 拉格朗日多项式n Hermite多项式l 单元的坐标n 不同单元可以取不同的坐标。n 直角坐标，自然坐标l 一维有限元模式n 子域划分把内部区域分为有限个子域（Element），连接各单元的点称为节点。如图11-7所示的节点为1、2、3、7等。n 编号图11-7 一维有限元的节点与编号(e)1 2ba1 23 4 5(1) (2) (3)(4)有限元的编号分节点编号，单元编号和总体编号三种。图11-7a中的1、2、3、4和5就是总体节点编号，用表示

31、；而（1）、（2）、（3）、（4）就是单元编号，用符号（e）表示，因此，单元节点用符号表示。一维有限元的单元如图11-7b所示，其中1和2表示单元的节点。n 总体节点编号与单元节点编号的关系总体节点编号与单元节点编号满足布尔矩阵关系。如图11-17所示，可得到它们之间的关系表达式。用布尔矩阵表示，有总表达式为：（11-57） 式中，是布尔矩阵（逻辑矩阵），表达式为：反之，也有（11-58） 式中，。对于一维有限单元，总体节点编号与单元节点编号的关系见表11.4。表11.4 总体节点编号与单元节点编号的关系单元号 e总体节点编号m1m2112223334445n 局部和总体插值函数u 单元上的线

32、性插值对于局部单元，如果则函数u可近似用线性函数表示。将边界条件代入上式，可得单元上的线性插值函数。（11-59） 上式可以简化为（11-60） 式中，。称为插值函数，是节点的函数值，正是需要求解的值。u 插值函数的特性（11-61） u 单元上的二次插值设，仍用式（11-60）的形式表示，则插值函数为：（11-62） 上式的插值函数满足式（11-61）给出的插值函数的特性。u 总体插值函数所有局部插值函数的集合就是总体插值函数。l 总刚合成与有限元方程n 用利兹法导出有限元方程下面以一维微分方程为例说明有限元方程的导出过程。(1) 问题(2) 单元划分单元划分如图11-7所示。表11.4给出

33、了总体节点编号与单元节点编号之间的关系。(3) 对每个单元用变分原理，有泛函(4) 选取插值函数(5) 推导有限元方程（11-63） 令 ，式（11-63）可以改写为：（11-64） 对式（11-64）变分，并令其等于0，有可得单元的有限元方程。（11-65） 式中，是节点上的值，即待求解的值。(6) 总体有限元方程的合成有两种方法可以合成总体的有限元方程，一种是布尔矩阵法，另一种是对号入座法。(7) 总体有限元方程合成步骤总结u 按单元节点的总体编码将单元的刚度矩阵扩大成与总体刚度矩阵同阶的矩阵；u 将扩大后的单元刚度矩阵相加形成总体刚度矩阵；u 叠加是按照总体节点编号在总体刚度矩阵的相应位

34、置进行的。注意：扩大的单元刚度矩阵可以看作是该单元对总体刚度矩阵的贡献，叠加起来就是所有单元刚度矩阵对总体刚度矩阵的总贡献。n 用伽辽金法导出有限元方程(1) 问题为了便于与精确解比较，现选取一个可解析解的常微分方程作为定解问题。(2) 单元划分如图11-7所示，取4个单元，共5个节点。(3) 写出伽辽金弱表达式（11-67） （11-68） （11-69） (4) 选出插值函数故， ， (5) 推求局部有限元方程将插值函数及其微分式代入式（11-69），可得局部有限元方程。故式中若取插值函数为线性插值函数，将计算域等步长划分，得，将上式积分可得单元刚度矩阵，由此得单元有限元方程。有(6) 总

35、体矩阵合成将以上单元有限元中的未知数对号入座，可得总体有限元方程的总刚度矩阵。（11-70） (7) 计算结果解代数方程组式（11-70），可得速度在各x处的值，将其与精确解比较可知误差在1%以内（见表11.5）。表11.4 有限元法的数值计算结果与精确解的比较总体节点编号的计算值的精确解100.99861.00020.51.63481.64831.02.69612.71841.54.46074.48152.07.36907.38911.2.3 场域的划分及插值函数l 场域划分的基本要求n 把求解域分解成有限个单元。插值函数为线性的时，有3节点三角形单元，4节点四边形单元。插值函数为高次方时，有6节点三角形单元，8节点四边形单元。n 划分单元的基本要求单元划分的好坏