1、一、有限元方法解题分析 为了说明应用有限元方法的解题步骤,以及每一步骤中的要点,下面我们以两点边值问题为例进行具体分析。 考虑两点边值问题 ,1.10, 01.2udduLpquf axbdxdxu au b 1,0,0,p xCa bpqC a b qfC a b其中我们将从Ritz法和Galerkin法两种观点出发,导出解边值问题(1.1)、(1.2)的线性有限元方法。(一)从Ritz法出发建立有限元方程1、写出Ritz形式的变分问题由变分原理可知,与边值问题(1.1)(1.2)等价的变分问题是:求1,EuH使 1*minEu HJ uJ u 1,1.32J ua u uf u(),.bb
2、aadu dva u vpquv dxf ufudxdx dx,其中积分表达式(1.3)是应用有限元法求解(1.1)、(1.2)式的出发点。 2、区域剖分 剖分原则与差分法相同,即将求解区域剖分成若干个互相连接,且不重叠的子区域,这些子区域称为单元。单元的几何形状可以人为选取,一般是规则的,但形状与大小可以不同。对于一维情形最为简单:将求解区域 剖分成若干个子区间,其节点为, a b01inaxxxxb每个单元 的长度为1,iiiexx1.iiihxx 单元在区间中分布的疏密程度或单元尺寸的大小, 可根据问题的物理性质来决定,一般来说,在物理量变 化剧烈的地方,单元尺寸要相对小一些,排列要密一
3、些。 3、确定单元基函数 有限元法与Ritz-Galerkin方法的主要区别之一,就 在于有限元方法中的基函数是在单元中选取的。由于各个单元具有规则的几何形状,而且可以不比考虑边界条件的影响,因此在单元中选取基函数可遵循一定的法则。设hV为 的有限维子空间,它的元素为1EH hux要构造hV,只需构造单元基函数 。构造单元基i函数应遵循如下原则: (1)、每个单元中的基函数的个数和单元中的节点数相同,每个节点分别对应一个基函数,本例中,单元 有两个节点,因此基函数有两个。ie(2)基函数应具有下面的性质:1,0,.jkjkjkxjk其中 是单元节点序号为k的节点。kx若取 为线性函数,则按上述
4、原则,可将 中的基函数取为jxhV显然, 中任一函数 可以表示为基函数 的线性组合,即 hVhu ix 1122,hnnuuxuxux 1111111,1,2,1,0,0,iiiiiiiiinnnnnxxxxxinhxxxxxxhxxxxxhx在别处在别处(1.4)其中, 是 在节点上的值,即12,nu uuhu1,2,hiiuxu in在单元 上, 表示为ie hux 11111 =, 1.5 ,.hiiiiiiiiiiiiuxuxuxxxxxuuhhxxx可见,单元中的近似函数由单元基函数线性组合产生,全区域的近似函数由各个单元的近似函数叠加而成。 从以上可以看出, 是满足下列条件的所有函
5、数 的集合:hVhu 2(1),;(2)10.hhhhihua buuL a bueua 、在上连续,且 ,、在 上是次数不超过 的多项式;(3)、故 是 的一个n维子空间,称为试探函数空间 称为试探函数。hV1EHhhuV4、有限元方程的形成 与Ritz法一样,以 代替 ,在 上解泛函数(1.3)的极小问题。hV1EHhV将(1.5)代入(1.3),得 ,11,21 =,.2hhhhnnijijjji jjJ ua u uf uauuuf 令0,(1.6)hjJ uu便得到确定 的线性代数方程组12,nu uu1,.1,2, . 1.7 nijijiaufjn 称(1.7)为有限元方程。显然
6、,只要我们分别算出 及,ija ,( ,1,2, ), jfi jn,就可以求解(1.7)。但在工程计算中,并不是按照上述步骤形成有限元方程的,而是首先建立单元有限元特征式(称这一过程为单元分析),然后再将单元的有限元特征式进行累加,合成为总体有限元方程(这一过程称为总体合成)。下面分步分析具体的计算方法。第一步:单元分析。注意到 1122111,21 =,(1.8)2iiiihhhhnnxxhhhxxiiJ ua u uf upuqudxfu dx 我们来计算单元 上的积分。为讨论方便,作变换ie1(1.9)iixxh并引入记号 011,NN 则在 上, 可写成iehu 011101,hii
7、iiuxNuNuuNNu或写成 ,ihuxNu(1.10)其中, 011,.TiiiNNNuuu而 可表示为hu 11,1.11ihiiiuxuuMuh式中,1/,1/.iiMhh 于是有 11221010()()()(),1.12iiiixhhxxTThhhhxTTiiiiiiTTTiiiiTpuqu dxp uuq uudxhp MMuq NuNuduhpM MqN N duuK u 这里, 101,11,1,1.13iTTiiiiiiiiii ii iKhpM MqN N daaaa 称为单元刚度矩阵,其中 1121,1110112,110111,1110(1)(1), 1.14iiii
8、iiiiiii iiiiiiiiiiii iiiiiiiah p xhhq xhdah p xhhq xhdaah p xhhq xhd 对(1.8)式右端第二项积分,同样有 1110,(1.15)iiTxihiiixTiifu dxh Nuf xhduF式中, 111101110 ,11.16TiiiiiiiiiiiiiiiFFFFhf xhdFhf xhd 称 为单元“荷载”向量。 iF根据以上分析,便有 111.1.172nnTTiiiiihiiJ uuK uuF这样,我们就得到了单元有限元特征式的一般表示形式: iiiK uF第二步:总体合成总体合成就是将单元的有限元特征式进行累加,合
9、成为总体有限元方程。这一过程实际上是将单元有限元特征式中的系数矩阵逐个累加,合成为总体系数矩阵(称为总刚度矩阵);同时将右端单元荷载向量逐个累加,合成为总荷载向量,从而得到关于 的线性代数方程组。12,nu uu为了形成总刚度矩阵,我们令 122 n ( ,) ,0010000010 1 Tniuu uuBii列 列 于是有 1,TiiiiuuuB u从而(1.17)右端第一个和式为 11121()212nTiiiiniiiTTiTuK uuBKBuu Ku其中, 111,11,1,1,()nniiiiTiiiiniiiiiiii ii iKBKBKaaaa第i-1行第i行这就是总刚度矩阵(未
10、标明的元素均为0)。对(1.17)右端第二个和式,有 11,nnTTiiiiTTiiuFuBFu b 11110,0,0,0.(1.19)nnTiiiiinTiiiiibBFFFF其中,这就是总荷载向量。这样,就可以将(1.17)式写成1.2TThJ uu Kuu b因此,有限元方程为.(1.20)Kub从总刚度矩阵和总荷载向量的形成过程可以看出, 的计算,实际上是把 中四个元素在适当的位置上“对号入座”地叠加, 的计算也是如此。我们引入 ,只是为了叙述方便,实际上,在编制程序时并不需要。K iKb iB 显然,方程组(1.20)的系数矩阵K是一个对称正定的对角矩阵,因此可采用追赶法求出u在节
11、点上的近似值 。如果我们认为这个近似解不够精确,则可以使剖分更细,即节点取得更多。这样,就产生一个收敛性与误差估计的问题。由于此问题所用的数学工具较多,本课程不做讨论。另一方面,我们以上是在单元剖分的基础上,利用Lagrange型的分段线性插值函数构造出的n维子空间 ,这样自然想到,如果不采用分段的线性插值,而采用分段的高次插值,则会得到更好的近似。12,nu uunV注1、当第一边值条件非齐次时,例如 ,则需象其它单元一样形成 上的单元刚度矩阵。但形成总刚度矩阵K时,先把当作未知量,K扩大成 矩阵,然后去掉第一行(或者一开始就不计算第一行),把第一列的第j行元素 乘以 累加到第j个方程的右端
12、后,再去掉第一列。最后仍然归结到方程(1.20),只不过右端向量因第一边值作了修改。 u a101,ex x0u 11nn0ja0,u 2212.2bnaJ upuqufu dxp b u它只是比齐次边值多了第二项。由于第二项只含有 的一次项,因此从上述泛函出发所形成的有限元方程不影响总刚度矩阵,唯一的改变量是第n个方程nu u b注2、若第二边值条件(右边值条件)非齐次,例如 ,则需从下列泛函出发:的右端要累加 .p b对于第三边值条件 u bu b则不但要修改第n个方程的右端,而且总刚度矩阵的第n行n列元素也要作适当的修改。有兴趣的同学可以自行推导。(二)、从Galerkin法出发建立有限
13、元方程从Galerkin法出发形成有限元方程的过程与前面完全一样,得到的结果也是一致的。但是从Galerkin法出发形成的有限元方程更具一般性,它不仅适用于对称正定的算子方程,而且也适用于非对称正定的算子方程。在实际问题中,主要是依据这一观点建立有限元方程。下面对这一问题作一简单陈述。由变分原理可知,与边值问题 ,1.10, 01.2udduLpquf axbdxdxu au b 等价的Galerkin形式的变分问题是: 11 , ,0, 1.21EEuHa u vf vvH 求使得我们仍用分段线性函数构成的试探函数空间代替 ,由(1.4)定义的分段线性函数是的一组基。和前面一样的方法,把hV
14、1EHhV 11,nnhiihiiiiuux vvx代入(1.21),便得到 所满足的方程组12,nu uu1,1,2, . 1.22nijijiaufjn 这和方程组(1.7)是完全一样的。 容易看出,方程组(1.22)的系数矩阵就是总刚度矩阵,在总刚度矩阵形成的过程中,注意到 11 1.23nnTTiiiiiiivK uvF而 ,iiiiuB u vB v从而有 11 nnTTiiiiiiiv BKBuv BF即0,hhv KubvV故有Kub这就是有限元方程(1.20)。由上述看出,按Galerkin法推导有限元方程更加直接方便。尤其重要的是。按这一观点推导的有限元方程,不仅适用于稳定问
15、题,而且也适用于非稳定的问题,因此它具有广泛的适用性。(三)、应用举例用有限元方法解边值问题 22,01,010,d uuxxdxuu将区间0,1等分成4个单元。 解、利用上述分析结果,我们只需构造出单元刚度矩阵和单元荷载向量,然后合成为总刚度矩阵和总荷载向量。 注意到(1.14)和(1.16): 1121,1110112,110111,1110(1)(1), 1.14iiiiiiiiiii iiiiiiiiiiii iiiiiiiah p xhhq xhdah p xhhq xhdaah p xhhq xhd 11101110 11.16iiiiiiiiiiFhf xhdFhf xhd 若将
16、 取成单元 上的中点值 则不难得到 ,p xq xf x1,iixx,iiip qf 1,11,1,11121,11126,1,1,2iiiiiiiiiiiiii ii iTTiiiiiiiaapq hKhaahFFFf其中 单元 的中点为1/4,/4,iihxi1,iixx 1/221 /8,1,1,.iixxipqf xx于是有 11219895114,111295986 42412112111212 4864iiKiiFi 如果把单元刚度矩阵 和单元荷载向量“扩大”,便得到 和 为 iK iF iK iF 1989595981,02400K 2098951,95982400K类似地,可写
17、出 34KK和。 1234100013001111,.03506464646400570007FFFF 然后进行叠加,便得到总刚度矩阵和总荷载向量: 4=19895959898951=9598989524959898959598iiKK98959519695195196952495196959598 41111 3411.3586464571277iiFF最后,考虑到约束条件 0400,10,0,uuuu即令并在K中划去首末两行和首末两列,F中划去首末两行,便得到如下线性代数方程组:精品课件精品课件!精品课件精品课件!123196951395196952 ,2951963uuu 解之得:1230.03521,0.05686,0.05052.uuu