流体力学课件（第十五章）.ppt

流体力学,退出,中国科学文化出版社,目录,流体力学基础,第一篇,第二篇,流体动力学基本原理及流体工程,退出,第三篇,计算流体动力学,第一篇流体力学基础,绪论场论与正交曲线坐标流体静力学流体运动学,第一章,第二章,第三章,第四章,退出,返回,第二篇流体动力学基本原理及流体工程,流体动力学微分形式基本方程流体动力学积分形式基本方程伯努利方程及其应用量纲分析和相似原理流动阻力与管道计算边界层理论流体绕过物体的流动气体动力学基础,第五章,第六章,第七章,第八章,第九章,退出,返回,第十章,第十一章,第十二章,第三篇计算流体动力学,计算流体动力学数学物理基础流体动力学问题的有限差分解法流体动力学问题的有限元解法,第十三章,第十四章,第十五章,退出,返回,第十五章流体动力学问题的有限元解法,有限元法的基本思想与区域离散化有限元法中代数方程的建立二维边值问题有限元法求解举例有限分析法介绍,第一节,第二节,退出,返回,第三节,第四节,第十五章流体动力学问题的有限元解法,退出,返回,第一节有限元法的基本思想与区域离散化,第1页,在上一章中，对求解流体动力学问题的有限差分方法进行了比较仔细的讨论。有限差分法的优点是原理简单，便于实施，对非线性比较强的对流换热问题有比较好的适应性。其弱点是对复杂几何形状的区域适应性比较差，采用近二十年发展起来的网格生成技术后可以克服这一弱点，但增加了计算工作量。而有限元法对复杂几何形状的区域具有较强的适应性。因而这两种方法已在扩散方程的求解中得到广泛的应用，同时研究工作正在向流场求解方面深入展开。近二十年内还发展起一种称为有限分析法的数值方法。本节中将对有限元法和有限分析法的基本思想作简要介绍，并以扩散方程为例说明其实施过程的主要步骤。,用有限元法求解物理问题时，总的解题步骤仍如图13.3所示，它与有限差分法的区别主要在于区域离散化的方式不同，建立代数方程所依据的原则或方法不同，以及由此而引起的代数方程求解方法的不同。本节中主要讨论有限元法的基本思想及其区域离散化（包括插值函数）等问题，下一节中再讨论代数方程的生成。,退出,返回,第一节有限元法的基本思想与区域离散化,第2页,假设我们要在某一个区域内求解一个偏微分方程（例如二维稳态的导热方程）：,图15.1三角形单元网格,一、有限元法的基本思想,其中表示一个微分算子，例如对直角坐标中的二维稳态导热问题，就代表。首先把求解区域划分成许多子区域，称为单元。对二维问题，单元的形状可以是矩形，三角形或四边形，单元形状的这种多样性使有限元法对求解区域的几何形状有很好的适应性。,（15.1）,在每个单元中取定几个点作为节点，例如对三角形单元一般取其三个顶点作为节点。然后对于单元内的被求函数的局部变化特征作出假设，也就是选定型线或插值函数，一般选用多项式作为插值函数。例如对三角形单元（图15.1），可设：,第十五章流体动力学问题的有限元解法,第十五章流体动力学问题的有限元解法,退出,返回,第一节有限元法的基本思想与区域离散化,第3页,式中，系数，及可以用三个节点的值来表示，符号表示被求函数的近似表达式。显然这一多项式不可能恰好是式（15.1）真正的解，即若把它代入到该式中，则其右侧不会等于零。将不等于零的部分作为余量，记为R，利用加权余量法，要求余量在某种意义上为最小，即要求：式中，A为区域的面积，W为权函数。上式要求余量R与权函数在A区域上的内积为零。由于近似解是用末知节点上的函数值来表示的，因而式（15.2）给出了这些未知节点上函数值之间的代数关系式，即有限元法的离散方程式。求解这些方程，可得到有限元法的数值解。在有限元法中导出离散方程的方法较多，有变分法，最小二乘方法，加权余量法等，其中加权余量法应用范围较广，本书中仅介绍这一种方法。,（15.2）,第十五章流体动力学问题的有限元解法,退出,返回,第一节有限元法的基本思想与区域离散化,第4页,二、有限元法区域的离散化在二维问题中最常用的单元为三角形单元，其三个顶点可作为节点，如图15.2(b)所示，是二维问题的线性单元。它对不规则区域的适应性较好（参见图15.1）。,图15.2一维和二维线性单元,对于一维问题，单元都是直线段，每个单元上的节点数取决于所选定的型线。若选用线性函数作单元上的函数逼近，则在该单元上只需两个节点，线性方程中的两个未知量（截距与斜率）可由这两个节点上的未知函数值确定，具有这种特性的单元称为线性元，如图15.2(a)所示。如果在直线段的单元上选取多于两个以上的节点，就为非线性单元。本章中仅介绍线性单元。,第十五章流体动力学问题的有限元解法,退出,返回,第一节有限元法的基本思想与区域离散化,第5页,在二维问题中最常用的单元为三角形单元，其三个顶点可作为节点，如图15.2(b)所示，是二维问题的线性单元。它对不规则区域的适应性较好（参见图15.1）。,图15.3稀疏矩阵及带宽,将一个求解区域划分成许多相连接又不重叠的子区域的过程就是区域离散化，子区域就是单元。一般先将子区域分成四边形与三角形的组合，然后再将四边形细分成三角形。要注意不能把一个三角形的顶点取在任一相邻三角形一条边的中间位置上。,单元的尺寸及疏密程度据物理问题的性质及对计算精度的要求而定。一般而言，物理量变化剧烈的地方单元的尺寸要小一些、排列要密集一些。有限元法的计算精度受到单元内最长边与最短边长度比的影响，尽量不要把三角形划分成钝角三角形，因为那样会使长短边之比增加而使计算精度下降。,第十五章流体动力学问题的有限元解法,退出,返回,第一节有限元法的基本思想与区域离散化,第6页,区域离散化过程中需给单元与节点编号。整个计算区域中的单元是统一编号的，单元号的字符用e表示，其值从1开始，顺序增加。与有限差分法不同，有限元法中的节点有两个编号，即单元节点号（局部的）和总体节点号。在单元中节点一般用i，j和k（或1，2，3）按逆时针方向编号（图15.2(b)），整个计算区域内的节点则按一定的顺序统一编号。总体节点编号的原则是尽可能使同一单元内各节点的编号相近，因为单元节点号的差值，决定了所形成的代数方程系数矩阵的特性。有限元法所生成的代数方程的系数矩阵是一个稀疏矩阵，即系数矩阵中有相当多的元素为零。如果总体节点编号合适，同一单元中各节点的编号相差较小，可以使非零元素相对集中地分布在系数矩阵的对角线附近。如图15.3所示，从对角线到非零元素所在区边界之间的距离称为带宽。在用直接解法求解代数方程组时（有限元法所生成的代数方程组多用直接解法求解），在采用一定的处理方法后，只需把非零元素输入计算机进行计算。在一定的总节点数下，带宽越窄，需送入计算机的非零元素越少，所占用的计算机内存就越小。,第十五章流体动力学问题的有限元解法,退出,返回,第一节有限元法的基本思想与区域离散化,第7页,可以证明，如果计算区域中每个单元三个节点编号数的最大值为R，则带宽与成正比。图15.4所示为同一计算区域节点的两种编号方式，方式(a)可比方式(b)节省一半以上的计算机内存。,第十五章流体动力学问题的有限元解法,退出,返回,第一节有限元法的基本思想与区域离散化,第8页,三、单元的插值函数对于一维问题的线性单元（图15.2(a)），设插值函数为由1、2两节点上的函数值、可得的表达式为：式中为书写方便略去了表示近似值的符号。将该式代入式（15.3）得,（15.3）,（15.4）,（15.5a）,这里h是线性单元的长度。显然这一插值函数对计算区域中的各个单元都是适用的。取其中任一单元为e，则上式可写成为：,（15.5b）,（15.5c）,其中：,和称为单元e的形状函数。,第十五章流体动力学问题的有限元解法,退出,返回,第一节有限元法的基本思想与区域离散化,第9页,对于三角形线性元（图15.2(b)），若讨论导热问题，可假设单元中的温度为x，y的线性函数，则有：其中待定常数、及可由节点上的温度值来表示。为此将三个节点i，j和k的坐标代入上式（区域离散化后，各个节点的位置坐标及单元的面积均为已知），将所得的三元一次代数方程组写成矩阵形式，有：,（15.6）,（15.7）,第十五章流体动力学问题的有限元解法,退出,返回,第一节有限元法的基本思想与区域离散化,为书写方便，令：,第10页,利用矩阵求逆的方法，可得a1，a2及a3：,第十五章流体动力学问题的有限元解法,退出,返回,第一节有限元法的基本思想与区域离散化,第11页,再将分母中的行列式展开：可证明之值等于该三角形单元面积的两倍。把三角形单元的面积记为，则有：,（15.8）,将a1、a2及a3的表达式代入（15.6），可得单元函数的插值计算式：,（15.9a）,其中、为单元形状函数，其计算式为：,（15.9b）,第十五章流体动力学问题的有限元解法,退出,返回,第一节有限元法的基本思想与区域离散化,第12页,在获得了单元形状函数后就可以构造总体的形状函数，即构造整个求解区域上被求函数的一种近似形式。我们以一维问题为例来说明。如图15.5(a)所示，该计算区域有三个单元，图中在括号内的数字表示单元编号，水平线以上的数字表示单元节点，水平线以下的数字表示总体节点号。注意到同一个节点的相邻的两个单元中的局部编号是不同的。例如总体编号为2的节点在单元中为编号2而在单元中则为1。在每个单元中两个插值函数是线性函数，如图15.5(b)所示。把各个单元内的插值函数叠加起来，有：,式中上角标表示单元编号，下角标为该单元中节点的编号。考虑到：于是得：,（15.10）,第十五章流体动力学问题的有限元解法,退出,返回,第一节有限元法的基本思想与区域离散化,第13页,式中N1，N2，N3及N4就是相对于节点1，2，3及4的整体形状函数，而上式也就是在整个求解区域中函数的近似表达式。注意每个整体形状函数在相应的节点上都取得“1”的值，而在其余节点上为零，这使上述的整体函数近似表达式能够满足在不同节点上取得该节点函数值的要求。在图15.5(c)中画出了各个单元中的局部形状函数，它们的叠加所形成四个整体形状函数如图15.5(d)15.5(g)所示。,第十五章流体动力学问题的有限元解法,退出,返回,第二节有限元法中代数方程的建立,第1页,本节主要介绍有限元法中离散方程建立的原理与过程，包括Galerkin（伽辽金）加权余量法的原理，单元矩阵的生成，总体矩阵合成及边界条件的处理等内容。,一、加权余量法加权余量法是获得微分方程近似解的一种有效方法。假设微分方程（15.1）中的未知函数可以近似地表示成为：,（15.11）,式中，是一些所选定的线性独立的函数（即其中任意一个函数都不能由其它函数经过线性运算而得出），Ci为未知的变量，m为未知变量的个数。加权余量法要求将（15.11）代入微分方程后所得到的余量在整个计算区域上与所选定的权函数的内积为零，即满足：,（15.12）,也就是在某种平均意义上要求余量为零。,第十五章流体动力学问题的有限元解法,退出,返回,第2页,第二节有限元法中代数方程的建立,加权函数的不同选择导致多种不同的加权余量法，其中应用较广的Galerkin余量法选定权函数。因而Galerkin余量法要求：即要求余量（误差）的加权平均值在整个计算区域上应等于零。这样m个积分式就产生m个代数方程，从而可以解出m个未知量。由于式（15.11）中的形状函数都是坐标的线性函数，而导热问题控制方程的最高阶导数为二阶，线性函数的二阶导数为零，因此不能直接将式（15.11）代入上式进行积分计算。为了克服这一困难可对式（15.13）作分部积分，对于一维问题，在区域上有：,（15.13）,（15.14a）,式中，u相当于式（15.13）中的权函数（即Wi），而dv则相当于微分算子L。分部积分的结果就可以把包括在算子符号内的函数的导数降低一阶。对于二维及三维的问题，分部积分相当于应用Gauss降维定理。设有一空间区域，其体积为，表面积为A，则Gauss降维定理为：,第十五章流体动力学问题的有限元解法,退出,返回,第3页,第二节有限元法中代数方程的建立,式中，u相当于权函数；相当于微分算子，为矢量运算；n为边界外法线上的单位矢量。,（15.14b）,图15.6分段线性近似解,分部积分引起了两点变化：（1）对函数的近似表达式的要求降低了，只要求一阶导数存在即可；（2）引入了边界条件，式（15.14a）中积分的上下限就是引入了边界条件。由于式（15.14a）右端第二项的积分在计算时是按单元分段进行的，因而只需在每个单元上近似解的一阶导数存在且在单元边界上函数值连续即可。例如对一维问题的近似解式（15.10），可用分段线性函数近似表示（图15.6）,第十五章流体动力学问题的有限元解法,退出,返回,第4页,第二节有限元法中代数方程的建立,在得出了分部积分的表达式后，可将所假定的近似解（如式（15.11）所示）代入进行积分计算，对每一个整体形状函数都要在整个计算区域内作积分，每完成这样一个积分就得出一个代数方程式。但是注意到每个整体形状函数实际上只在一个单元或相似的几个单元之内值不为零，在其它单元上其值均为零（如图15.5(d)(g)所示），且在整体形状函数不等于零的单元内，整体形状函数之值就等于该单元的形状函数之值。因而整体形状函数在整个求解区域内的积分形成代数方程的过程为：对每个单元按单元形状函数作积分，然后把共享一个节点的各单元的积分结果按一定方式相加。由于每个单元的单元形状函数都是一样的，因而只要对一个代表性单元作积分即可。对每个单元积分所形成的该单元内的代数方程的矩阵称为单元矩阵，其实施过程称为单元分析。把共享一个节点的各单元的结果叠加以形成总体代数方程的过程叫作总体合成。单元分析与总体合成是有限元法中建立离散方程的两个重要步骤。,第十五章流体动力学问题的有限元解法,退出,返回,第5页,第二节有限元法中代数方程的建立,二、单元矩阵的生成下面我们以一维扩散问题为例来说明单元分析与总体合成的过程。设式（15.1）中的微分算子为一二阶导数，即，任一单元e的形状函数为，则对该单元内的积分有：即对每个单元，余量的加权平均值也等于零。式中x1，x2为单元e的两个端点。对该式做分部积分，可得：,（15.15a）,（15.15b）,利用的表达式（15.5b），有：,（15.15c）,第十五章流体动力学问题的有限元解法,退出,返回,第6页,第二节有限元法中代数方程的建立,式中，符号表示行矢量，表示列矢量。将式（15.15c）代入式（15.15b），得：对上式右端的第一项计算如下：,（15.15d）,这里已把图15.5(b)中单元形状函数的两个端点值代入。若采用下列符号：,（15.15e）,第十五章流体动力学问题的有限元解法,退出,返回,第7页,第二节有限元法中代数方程的建立,则式（15.15d）可写成：该式即单元e的代数方程，矩阵及称为单元矩阵。在获得了单元代数方程后就应将共享一个节点的各个单元的矩阵元素按一定规则叠加，以形成整个计算区域的代数方程。在叠加过程中，每个单元代数方程（式（15.15f）右端第一项的叠加结果，除了整个计算区两端单元外，其余都将相互抵消为零。而这涉及到该问题的边界条件。,（15.15f）,三、总体矩阵的合成与边界条件的处理关于边界条件引入到总体代数方程中去的方法在介绍了总体合成规则以后再讨论。为了便于叙述，以具有三个单元的一维扩散问题为例说明。总体合成的步骤如下：,第十五章流体动力学问题的有限元解法,退出,返回,第8页,第二节有限元法中代数方程的建立,（1）列出各单元节点的局部编号与总体编号间的对照表。对如图15.5a所示的三个单元的区域，这一对应关系如表15.1所示；,（2）将每个单元方程中矩阵元素的下标由局部节点号换成为总体节点号；（3）将每个单元的矩阵元素按其经过转换的下标号插入在到相应的总体矩阵的相应位置上。本例中有四个节点，因而矩阵K就为44方阵，而有关边界条件及源项的两列矢量均有四个分量；,表15.1局部节点编号与总体节点编号的对照,第十五章流体动力学问题的有限元解法,退出,返回,第9页,第二节有限元法中代数方程的建立,（4）将位于总体矩阵同一位置上的各个单元矩阵的元素相加，形成总体矩阵的元素。这一叠加过程表明总体矩阵中的元素代表了共享一个节点的各单元所作出的贡献之和。,对所举的例子，总体合成所得结果如下列矩阵方程所示：,式中L为计算区域的长度。为书写简便起见，把近似函数直接写为。等号右端的两个列矢量都是已知值，是代数方程组中的常数项，可合并记为列矢量R。于是此总体代数方程组便可表示为下列矩阵形式：,第十五章流体动力学问题的有限元解法,退出,返回,第10页,第二节有限元法中代数方程的建立,式中，n为计算区域中的总节点数；矩阵K可称为热导矩阵，相当于有限差分法中系数等所形成的矩阵，反映了热传导的机理对温度场的影响；列矢量R称为热载荷矢量，代表了热边界条件、内热源项等的影响。以下讨论如何把边界条件应用于矩阵方程（式（15.16）中。对第一类边界条件，可以采用有限差分法中对某点赋给定值的方法。例如设区域的左端点为给定值，则可以用一个很大的数，如1020，乘K中的第一个对角元，然后将R中的第一个元素以来代替，结果形成以下矩阵方程：,（15.16）,第十五章流体动力学问题的有限元解法,退出,返回,第11页,第二节有限元法中代数方程的建立,显然求解可得。对于给定热密度的边界条件，只要把给定的之值代入R即可。对第三类边界条件，若设区域的右端点为第三类边界条件，则有：式中对流换热系数、周围流体温度均为已知。考虑到导出总体矩阵方程时所依据的控制方程式相当于导热系数为1的情形，为协调起见，这里也取（时分析的方法完全一样），因而有。显然只要把R第四个元素中的用来代替，并把中的用来代替即可把对流边界条件考虑在内。,至此，可得采用Galerkin加权余量法导出有限元法离散方程的步骤：1写出所求解问题的控制方程及其边界条件的表达式；2将Galerkin加权余量法的表达式作一次分部积分，得出降维后的积分表达式；,第十五章流体动力学问题的有限元解法,退出,返回,第6页,第二节有限元法中代数方程的建立,3区域离散化，剖分单元，确定节点及其编号；4选定单元插值函数，导出单元形状函数；,5单元分析，将所选定的近似函数代入降维后的积分表达式，导出单元矩阵；6总体合成，将单元矩阵的元素按一定的方式叠加形成总体有限元方程；7处理边界条件。至于所形成的代数方程的求解与有限差分法中的代数方程求解一样，也有直接解法与迭代法两大类。在有限元法中，只要节点的总体编号合适，可以把总体代数方程组系数矩阵中的非零元素集中在对角线附近区域，同时采用压缩存储的技术，可以使总体热导矩阵中需要存储的元素个数远远小于该矩阵的元素总数，因而多采用直接解法。有关压缩存储技术的内容可参见相关文献。,第十五章流体动力学问题的有限元解法,退出,返回,第1页,第三节二维边值问题有限元法求解举例,例15.1设有两块形状为等边三角形的平行平板，相距为h，其间充满了一种粘性不可压缩的流体，粘度为。该两平板以速度w向前运动（图15.7(a)），同时粘性流体则从板中被挤出（轴承润滑是这一问题的工程背景）。由流体力学可知，此时两平板中压力场的分布由Poisson方程描述：,为了使读者进一步熟悉有限元法的求解步骤，本节中以一个二维边值问题为例给出有限元法求解的全过程。为便于手工计算，单元数取得很少，主要目的是使读者进一步熟悉单元矩阵的生成及总体矩阵的合成过程。,第十五章流体动力学问题的有限元解法,退出,返回,第2页,第三节二维边值问题有限元法求解举例,同时，在平板的边界上流体压力为零。三角形的边长为。试用有限元法确定三角形区域中任一点的流体压力值。解：为能进行手工计算在流场（平面流场，沿高度方向压力是均匀的）中只取四个节点，即三角形区域的三个顶点和一个内点4，其坐标为。并且为使计算结果能应用于内部各点，这里x，y暂不取定具体数值。节点的总体编号及单元编号如图15.7(b)所示。单元矩阵的元素为：,式中表示下标i作ijk的轮换。计算形函数的导数项，并代入上式，得：,第十五章流体动力学问题的有限元解法,退出,返回,对单元3有：,第3页,第三节二维边值问题有限元法求解举例,式中为单元面积。将三个单元矩阵中的节点号均转换为总体节点号，如表15.2所示。,于是，对单元1有：对单元2有：,第十五章流体动力学问题的有限元解法,退出,返回,第4页,第三节二维边值问题有限元法求解举例,将上述三个单元方程的元素装入到总体方程矩阵中：在上式中为表明元素来自哪一个单元，在其右上角用小括号标明了单元号。注意式中每个元素均已为整体坐标。,下面来确定边界条件。由于上式中1、2、3节点处的流体压力，和均为已知（等于零），因而采用消行修正法后只得出一个关于的方程：或,第十五章流体动力学问题的有限元解法,退出,返回,第5页,第三节二维边值问题有限元法求解举例,其中，、和的值计算如下：,代入P4的计算式得：而的关系可由三角形面积的计算公式求得：,第十五章流体动力学问题的有限元解法,退出,返回,第6页,第三节二维边值问题有限元法求解举例,代入P4计算式并整理有：可知在给定了后即可求得P4。有限元法是50年代从固体力学中发展起来的，60年代逐渐推广到传热学与流体力学的计算中。虽然有限元法在区域的离散（常称剖分）代数方程的形成等方面工作量比较大，但由于它对不规则区域的适应性强而得到重视，目前已是计算流体力学中应用较广的一种数值方法，并且随着计算机技术的高速发展，将在计算流体动力学中发挥更大的作用。,第十五章流体动力学问题的有限元解法,退出,返回,第1页,第四节有限分析法介绍,有限分析法是80年代初发展起来的一种数值计算方法。在计算流场方面它可以克服在高Re数下的有限差分法容易振荡或发散的缺点。本节中以二维稳态导热问题为例来介绍它的基本思想。用有限分析法求解边值问题的主要步骤如下：（1）把求解区域用一系列与区域边界平行的网格线进行离散，两根网格线的交点为计算节点，每一个节点与其相邻的四个网格组成一个单元，即每一个单元由一个内点和八个邻点组成（图15.8），图中有阴影线的部分即为P单元。应指出，这样定义的单元在每个相邻单元间有一部分区域是相互重叠的。（2）在单元内将控制方程的非线性项局部线性化，并对单元边界上未知函数的变化型线作出选择，把所选定型线表达式中的常数或系数项用单元边界节点上的值来表示。（3）找出上述条件下单元内未知函数的分析解（相当于边界上未知函数为给定的第一类边界条件）。,第十五章流体动力学问题的有限元解法,退出,返回,第2页,第四节有限分析法介绍,（4）选用所假定的分析解找出其内节点及八个邻点的未知函数之间的关系式。即建立该内点的离散方程。（5）逐一地对求解区域内每一节点建立起离散方程，对非第一类边界条件的边界节点，利用与该节点相邻单元中的分析解建立一个补充方程。（6）求解所建立的代数方程组，在获得了求解区域每一个点上的函数值后再进行其它量（如边界上法向导数）的计算。一个单元上九个节点的命名方式如图15.9所示。下面研究在这一单元上如何建立相应于Laplace方程=0的离散方程。,首先确定单元边界上用以逼近被求变量的近似函数的形状。一般有三种选择。以右边界为例来写出：（1）分段线性型线，以EC为原点分为上、下两支，即：,（15.17a）,第十五章流体动力学问题的有限元解法,退出,返回,第3页,第四节有限分析法介绍,（2）二项式分布：（3）多项式指数函数：其中，系数可采用边界上三点的未知值、表示，而系数A则在求解过程中确定，上述函数都是单调变化的。如果在求解前能预期到相邻两边界点间被求函数不是作单调变化时，应在该两节点间设置更多节点，以使每两邻节点间被求函数呈单调变化。,（15.17b）,（15.17c）,对于分段线性及二项式分布这两种型线采用多项式拟合法来导出系数、与边界节点函数值间的关系。把坐标原点设在EC上，则对分段线性的型线有：由此得：,（15.18a）,第十五章流体动力学问题的有限元解法,退出,返回,第4页,第四节有限分析法介绍,对二项式分布有：由此得：可计算当、和均已知时，这一单元上的分析解。这是第一类边界条件下矩形域上的Laplace方程，利用分离变量法可求得其分析解为：,（15.18b）,（15.19）,其中：,（15.20a）,（15.20b）,第十五章流体动力学问题的有限元解法,退出,返回,第5页,第四节有限