有限元方法专题知识专业知识讲座.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,*,本文档所提供的信息仅供参考之用，不能作为科学依据，请勿模仿。文档如有不当之处，请联系本人或网站删除。,*,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,本文档所提供的信息仅供参考之用，不能作为科学依据，请勿模仿。文档如有不当之处，请联系本人或网站删除。,*,*,*,单击此处编辑母版文本样式,第二级,单击此处编辑母版标题样式,本文档所提供的信息仅供参考之用，不能作为科学依据，请勿模仿。文档如有不当之处，请联系本人或网站删除。,有限元法是求解偏微分方程问题的一种重要数值方法，它的基础分两个方面：一是变分原理，二是剖分插值,从第一方面看，,有限元法是Ritz-Galerkin方法的一种变形,它提供了一种选取“局部基函数”的新技巧，从而克服了Ritz-Galerkin方法选取基函数的固有困难,从第二方面看，,它是差分方法的一种变形,差分法是点近似,，它只考虑在有限个离散点上函数值，而不考虑在点的邻域函数值如何变化；,有限元方法考虑的是分段（块）的近似,因此有限元方法是这两类方法相结合，取长补短而进一步发展了的结果在几何和物理条件比较复杂的问题中，有限元方法比差分方法有更广泛的适应性,1,7.,两点边值问题的有限元方法,本节以两点边值问题为例，并从Ritz法和Galerkin法两种观点出发来叙述有限元法的基本思想及解题过程.,7.1 基于Ritz法的有限元方程,考虑两点边值问题,其中，,2,1.写出Ritz形式的变分问题,与边,值问题,(7,.1,)、(7,.2,),等价的变分问题是：,求,，使,其中，,（7.3）,式(7.3)是应用有限元法求解边值问题(7,.1,)、(7,.2,)的出发点.,3,2,.,区域剖分,剖分原则与差分法相同，即将求解区域剖分成若干个互相连接，且不重叠的子区域，这些子区域称为,单元,单元的几何形状可以人为选取，一般是规则的，但形状与大小可以不同对于一维情形最为简单,.,将求解区间,分成若干个子区间，其节点为,每个单元 的长度为 .,单元在区间中分布的疏密程度或单元尺寸的大小，可根据问题的物理性质来决定，一般来说，在物理量变化剧烈的地方，单元尺寸要相对小一些，排列要密一些.,4,设 为 的有限维子空间，它的元素为 .,要构造 ，只需构造单元基函数 .构造单元基函数所遵循的原则是：,其中，是单元节点序号为 的节点.,(1)每个单元中的基函数的个数和单元中的节点数相同，每个节点对应一个基函数，本例中，单元 有两个节点，因此基函数有两个.,(2)基函数应具有性质,5,3.确定单元基函数,有限元法与Ritz-Galerkin方法的主要区别之一，就在于,有限元方法中的基函数是在单元中选取的,由于各个单元具有规则的几何形状，而且可以不必考虑边界条件的影响，因此,在单元中选取基函数可遵循一定的法则,6,7,8,9,(7.6),令,4,.,形成有限元方程,便得到确定 的线性代数方程组,称式（7.5）为,有限元方程,.,10,(7.8),(7.7),值得注意的是，在实际计算中，并不是按照上述步骤形成有限元方程的，而是,先,进行,单元分析,，即在单元上建立有限元特征式，然后,再,进行,总体合成,，即将各单元的有限元特征式进行累加，合成为有限元方程具体过程如下：,第一步：,单元分析,.注意到,作变换,11,(7.10),并引入记号,其中，.于是,或写成,(7.9),其中，.从而有,12,(7.11),这里,(7.12),称为,单元刚度矩阵,，其中,(7.13),13,(7.16),(7.14),式中,(7.15),将式（7.11）、（7.14）代入式（7.7），便有,对式（7.7）右端第二项积分，有,这样，我们就得到了单元有限元特征式的一般表示形式：,14,于是有,第二步：,总体合成,.总体合成就是将单元上的有限元特征式进行累加，合成为总体有限元方程.,这一过程实际上是将单元有限元特征式中的系数矩阵（称为单元刚度矩阵）逐个累加，合成为总体系数矩阵（称为总刚度矩阵）；同时将右端单元荷载向量逐个累加，合成为总荷载向量，从而得到关于的线性代数方程组为此，记,从而式(7.16)右端第一个和式为,15,(7.17),其中,（未标明的元素均为0）这就是,总刚度矩阵,对式(7.16)右端第二个和式，有,其中,这就是,总荷载向量.,16,从总刚度矩阵和总荷载向量的形成过程可以看出，的计算，实际上是把 中四个元素在适当的位置上“对号入座”地叠加，的计算也是如此我们引入 ，只是为了叙述方便，实际上，在编制程序时并不需要显然，方程组(7.18)的系数矩阵 是对称正定的三对角矩阵，因此可采用追赶法求出 在节点上的近似值 .,(7.18),其这样，就可将式(7.16)写成,因此，有限元方程为,17,7.,两点边值问题的有限元方法,本节以两点边值问题为例，并从Ritz法和Galerkin法两种观点出发来叙述有限元法的基本思想及解题过程.,7.1 基于Ritz法的有限元方程,7.2 基于Galerkin法的有限元方程,从Galerkin法出发形成有限元方程的过程与前面完全一样，针对边值问题(7.1)、(7.2)所得到的结果也是一致的但是从Galerkin法出发形成的有限元方程,更具一般性,，它不仅适用于对称正定的算子方程，而且也适用于非对称正定的算子方程，所以我们今后主要是依据这一观点建立有限元方程,18,与边,值问题,(7,.1,)、(7,.2,),等价的,Galerkin,变分问题是：,求 ，使得,(7.19),其中,仍用分段线性函数构成的试探函数空间 替代 ，将,代入(7.19)，则得到 所满足的线性代数方程组,(7.20),这和方程组(7.6)是完全一样的.,19,与,容易看出，方程组,(,7.,20),的系数矩阵就是总刚度矩阵在总刚度矩阵形成的过程中，注意到,(7.21),而,从而有,即,故有,这就是有限元方程(7.18).,20,由上述看出，按Galerkin法推导有限元方程更加直接方便尤其重要的是按这一观点推导的有限元方程，不仅适用于,定常的,微分方程定解问题,，而且也适用于,不定常的,微分方程定解问题,，因此具有广泛的适应性,例7.1 用有限元方法解边值问题,将区间0,1等分成4个单元,.,解,利用上述分析结果，我们只需构造出单元刚度矩阵和单元荷载向量，然后合成为总刚度矩阵和总荷载向量,21,注意到(7.13)和(7.15)，并将 形成单元,上的中点值 则不难得到,其中，单元 的中点为,于是有,22,如果把单元刚度矩阵 和单元荷载向量 “扩大”，便得到 和 为,类似地，可写出 和 .,23,然后进行叠加，便得到总刚度矩阵和总荷载向量：,24,依边界条件 即 在 中划去首末两行和首末两列，在 中划去首末两行，便得到如下线性代数方程组：,解之，得,25,8.二维椭圆,边值问题的有限元方法,用有限元方法求解二维椭圆边值问题的过程与两点边值问题的有限元方法大体相同，只是在具体处理时比一维情形更复杂些考虑如下椭圆型方程的第一边值问题,：,(8.1),(8.2),其中，是 平面的一个有界域，其边值 是分段光滑的简单闭曲线.,以下我们从Galerkin法出发，叙述有限元求解问题(8.1)、(8.2)的全过程,26,与边,值问题,(8,.1,)、(8,.2,),等价的,Galerkin,变分问题是：,求 ，使得,(8.3),其中,8.1 区域剖分,正如,前章,所言，对高维区域的剖分与对一维区域的剖分有很大不同.对,一维区域,无论作哪一种剖分，其,单元仍然是一个区间,，,对不同的剖分只是区间长度不同而已,.,对高维区域而言，不同的剖分其,单元的形状各异,，如对二维区域，剖分后的子区域可以是,三角形、矩形或四边形,.限于篇幅，本书只讨论剖分后所得的子区域是三角形的情况，这种剖分称为,三角形剖分,.,27,将区域 划分成有限个三角形单元，剖分方法见前章，那里曾假定剖分的单元应是,锐角三角形,现在我们去掉这一限制，,只假定不同的单元是无重叠的内部，且单元的顶点不是其它单元边的内点,当然还要,尽量避免出现大钝角的三角形,在物理量变化剧烈的地方，单元要划分得细密一些，变化缓和的地方，划分得稀一些,划分好单元之后，要对单元和节点进行编号设 是区域中的单元总数，将全区域中的单元统一编号，单元号记为 全区域中的节点也要按一定的顺序统一编号，记全区域中共有 个节点，节点号记为,节点编号的一般原则是尽可能使同一单元内的节点号比较接近以后可以看到，单元内节点序号的差值决定了总体系数矩阵的带宽,28,8.2 确定单元基函数,与一维情形一样，为了构造试探函数空间 我们只需在每个单元上构造插值基函数这里，我们仅考虑三角形单元上的线性插值函数为了便于后面积分的计算，我们先将直角坐标转换为面积坐标.,1.,面积坐标及有关公式,(1)面积坐标的定义,设,是以,为顶点的,任意三角形单元，面积为,，我们,规定,的次序按逆时针方向排列,在 中(图8.1)任意一点,的位置，可用它在直角坐标系,中的两个坐标值 来确定.,图(8.1),29,如果我们过点,作与三个顶点的连线，形成三个小三角形,那么一旦 的值确定后，这三个三角形的面积也就有了确定的值；反之，这三个小三角形的面积确定之后，点 也就有了确定的位置，由此可见，三角形单元中任意一点 的位置，除了可用直角坐标 来确定外，还可以用连接点,与,与三个节点所形成的小三角形的面积来确定,用 分别表示这三个小三角形的面积，显然,30,令,(8.4),则称这三个比值 为点 的,面积坐标,.由定义可知，,所以，并不是互相独立的，其中任意一个面积坐标都可以用另外两个面积坐标来表示，而且它们与直角坐标系的选取方法是无关的，这也是采用面积坐标的一个优点,显然，三个节点 的面积坐标是,节点,节点,节点,31,(,2,)面积坐标,与直角坐标的关系,我们知道,于是，有,其中,32,由此可得到面积坐标与直角坐标之间的如下转换关系：,(8.6),(8.7),从上述关系中可以看出，面积坐标 和直角坐标,之间是线性变换的关系，它实际上是将 平面上的任意形状的三角形变换到 平面上的直角三角形单元经过这种变换，使得在任意三角形区域上的积分问题转化为在直角边为1的直角三角形区域上的积分问题，所以在计算上会带来很大的方便,33,(,3,)面积坐标函数对直角坐标的偏导数,设面积坐标函数为 是 的函数，由复合函数的求导法则，有,注意到式(8.6)，可得,34,所以面积坐标函数对直角坐标的偏导数为,(8.8),(4)面积坐标的积分,单元分析中的积分，由于基函数几乎无例外地均采用多项式函数，被积函数一般都是以幂函数形式出现的，因此在单元分析中经常考虑的是如下典型形式的积分,其中 是任意非负整数,35,以面积坐标替换直角坐标，并利用重积分变量替换公式，不难算出,(8.9),为证明此式，只需注意到变换的Jacobi行列式,以及积分关系式,(8.10),36,容易看出，在一维有限元分析中，由式(7.8)给出的变换 正与上面讨论的面积坐标相当如果说变换(8.6)将 平面上的任意形状的三角形单元变换到 平面上直角边为1的直角三角形单元，那么变换(7.8)则把 轴上的线段单元 变换到 轴上的参考单元0,1由此可见，在一维的情况下，与公式(8.9)相应的积分公式为,(8.11),其中,为线段单元 的长度,37,2.构造单元基函数,任一个三角形单元 上，可唯一确定一个线性插值函数,其中，为三角形单元顶点 处给定的函数值，,为 处相应的单元节点基函数，它们都是线性,函数，且满足条件,根据前面的讨论容易验证：,即面积坐标正好是三角形单元上线性插值函数的基函数，于,是在任一三角形单元 上,38,(8.12),称式(8.12)为单元形状函数将每一个单元上构造的函数,合并恰来就得到 在整个区域 上的分块近似函数,.由于每个节点对应一个基函数，所以整个区域 上共,有 个基函数.容易证明，由 所生成的试,探函数空间 是 的 维子空间，中的任一函数 均,可表为,其中，是 在节点 处的值.,39,8.3 单元分析,令,(8.13),下面利用式,(8.6),、,(8.8),对式,(8.13),右端的积分逐项进行分析.,设单元,为,在三个顶点上，,的值分别为,和,记,为 的面积，,40,则在单元 上,于是式,(8.13),右端的第一个积分可化为,41,其中,(8.14),称 为,单元刚度矩阵,.,同理，可将式(8.13)右端的第二个积分化为：,42,其中,(8.15),称 为单元载荷向量.,综合以上分析，我们得到,8.4 总体合成,令,(8.16),我们知道，单元刚度矩阵或单元载荷向量中的元素下标值对应于单元节点序号所谓,总体合成,，就是将这些序号转,43,换为总体节点序号，然后把这个元素加到总体系数矩阵某个位置上去这个位置的行列序号，正是相应的总体节点序号必须指出，对于每个单元中的节点号统一按逆（或顺）时针方向排列，但是它们在总的节点编号中就不一点按原来的次序排列了,设已经给出了单元节点编号与总的节点编号的对应关系，令,对单元 有,其中，是一个 的矩阵.若单元 的节点 在总节点编号中的序号为 ，则 的第一行第 个元素为1，其余元素均为0；的第二行、第三行分别也只有一个非零元素，其值为1，其位置由单元 的节点 在总编号中的序号来决定,44,同理，可以写出,于是式(8.16)右端的第一项成为,其中,(8.17),就是,总刚度矩阵,，表示对 个单元求和.,显然，只是单元刚度矩 的九个元素在,总节点编号下重新排列和“扩展”的结果，而总刚度矩阵,则是将各个单元的“贡献”叠加起来.它是一个 的对称、正定矩阵.,45,对于式(8.16)右端的第二项，同样得,其中,(8.18),就是,总载荷向量,，它是每个单元上单元载荷向量贡献的叠加,这样由式(8.3)可知，对 ,有,46,故得到 所满足的线性代数方程组,(8.19),方程组(8.19)的系数矩阵 对称、正定，故方程组(8.19)有,唯一解 求得它们后，就有,8.5 边界条件的处理,1.第一边值条件,若第一边值条件为非齐次的,(8.20),则应像内点一样，在界点也引进基函数,引进的方法及计算公式同内点完全一样,47,假定内节点和边界点的总个数为 ，则与上面的推导完全一样，可得到 所满足的线性代数方程组,或写成,(8.21),要得到问题(8.1)、(8.20)的解 在内部节点上的近似值，则可按如下方法对方程组(8.21)进行处理,假设有 个边界节点，个内节点为叙述方便，假定在总体节点编号时，,把边界 上的节点排在最前面,，于是为了得到内点所满足的有限元方程，可先从方程组(8.21)中去掉前 个方程，即：,48,(8.22),然后用边值 代入左端相应项，并移至,右端，便得到,有限元方程,，它是一个具 个未知数、个方程的线性代数方程组,若记,其中，是 的方阵，是 的方阵；都是 的，,都是 的则方程(8.23)可写成如下形式,49,(8.24),其中 是从 中划去头 行 列元素而得到的.显然，若边界条件是齐次的，则此方程组就是式(8.19),若边界 上的节点不是排在前 个,，则在 中划去相应的 行 列，在 中划去相应的 行，然后用边值,代入左端相应项，并移至右端，同样得到内点所满足的线性代数方程组,以上的约束处理方法，将 改 为时，要重新存储矩阵的元素，从而给编制程序带来麻烦因此在实际计算时，应把方程组（8.24）改写为,(8.25),50,其中,是 阶单位矩阵，是由边界 上节点 的值所组成当边界条件为齐次时，这一方程组的系数矩阵保留了 的阶数它与方程组(8.24)等价,2.第二和第三边值条件,假定问题中给定的是下列边值条件之一：,(8.26),(8.27),其中，都是给定的连续函数.,因为它们是,自然边值条件,，可像内点一样在界点形成有限元方程,只不过双线性形式 因边值条件不同需要修改.,51,例如，对Poisson方程(8.1)与边值条件(8.26)，相应的双线性形式为,变分式的右端修改为,于是有,(8.28),其中，是第 个单元如果它的边界 至少有一条边是在,的边界 上，则记这部分边界为 显然，若 是“内部的元”，即 的三条边都不在 上，则在(8.28)中将不出现,线积分项,52,为了书写方便，我们记这种情况下的 为 ，即空集，,在空集上的积分为零,在这种情况下，只需对单元载荷向量作如下修改：,(8.29),其中,53,如果 是非空集，可设 是 边（如果是其它边，只需将下列计算公式稍作修改，如果是两条或三条边，则应计算类似的两个或三个积分），令 ，引入参数 为 上的弧长，对应 有 ，对应 有 ，这样可得,这样就得到了单元载荷向量,如果对Poisson方程(8.1)与边值条件(8.27)，相应的双线性形式为,(8.30),变分式的右端项为,54,于是，有,(8.31),这时单元载荷向量也具有(8.29)的形式单元刚度矩阵修改为,其中,与(8.14)的形式完全一样，而 是由(8.31)左端的,线积分,产生的三阶矩阵，其表达式为,(8.32),55,其中,如果问题给出的是,混合边值条件,：在 的一部分 上给出第一边值，在其余部分 上给出第二边值这时可用前述方法分别处理 和 ，,和 的,交界点应取作剖分节点,，把它当做第一边值或第二边值点处理均可,56,