基于有限元法的振动分析郑佳文

1、明德明德 砺志砺志 博学博学 笃行笃行基于有限元法的振动分析基于有限元法的振动分析陕西理工学院陕西理工学院 机械工程机械工程 郑佳文郑佳文1明德明德 砺志砺志 博学博学 笃行笃行内容n 有限元法简介n 一维单元 杆单元 梁单元n 二维单元与平面问题的有限元法 三角形单元 矩形单元 等参单元n 平面问题的有限元法2明德明德 砺志砺志 博学博学 笃行笃行有限元法简介 有限元法是一种可用于精确地(但近似)解决许多复杂的振动问题的数值方法。 对于基本的一维元素进行有限元分析，能得到质量矩阵与刚度矩阵和所需的力矢量，对于二维三维，元素矩阵会转换成相关的更高维的空间。使用一致的和集中质量矩阵的有限元方程并

2、结合边界条件能为复杂系统提供解释。3明德明德 砺志砺志 博学博学 笃行笃行当用有限元法来分析时，将梁分解成为梁单元，将各单元相邻的节点处的弹性位移（横向弹性位移和弹性转角）作为问题中的位置量。 列子4单元质量矩阵单元质量矩阵单元刚度矩阵单元刚度矩阵集合系统质量矩阵系统质量矩阵M系统刚度矩阵系统刚度矩阵K 梁的有限元模型明德明德 砺志砺志 博学博学 笃行笃行有限元基本思想讲一个连续弹性体看成是由若干个基本单元在节点彼此相连的组合体；由于单元较小，允许对位移分布规律做出某种假设，基于这一简化假定，回避了弹性力学中的微分方程组难以求解的困难；从而使这一无限自由的连续体问题变成一个有限自由度的离散系统

3、问题。5明德明德 砺志砺志 博学博学 笃行笃行一维单元一 杆单元 一个杆单元是从杆上划分出的一个小段，如下图所示。由于单元很小，、A均视为常量。现在就以这最简单的杆单元，推导出它的质量、刚度矩阵，等效节点力。6图1.2明德明德 砺志砺志 博学博学 笃行笃行)()(),(2211tututxu7(1) 求杆单元上任意点的位移u(x,t)本来，杆单元上任意点的位移u(x,t)与节点的位移u1(t)、u2(t)之间的关系是未知的，但是，只要单元划分的足够小，那么其间的关系就无关大局。所以可以假定它们之间有简单的线性关系，即根据节点位移对单元内任意点位移进行插值：（1）式中，1、 2称为线性系数，与单

4、元里点的位置有关，是x的函数。此函数与单元的形状有关，又叫形状函数。明德明德 砺志砺志 博学博学 笃行笃行)(), 0(1tutu8)(),(2tutlu形状函数和插值函数一样是任意的，但必须边界条件：只有满足此条件单元才能协调一致运动，而不致破坏系统的完整性，因此这两个条件实际上就是变形协调条件。将式（1）带入（2）中，就可以得到形状函数1(x)、2(x)所满足的边界条件： （2） 1 , 00 0 , 102121ll （3）明德明德 砺志砺志 博学博学 笃行笃行 lxxlxx21 ,19 xx21、我们可以用简单的线性函数来近似，因此有：（4）代回（1）式中有： tulxtulxtxu2

5、1)1(),(（5）为此，我们已经找到了用节点位移表示单元内任意一点位移的表达式。以上边界条件确定了由于这两个函数的任意性明德明德 砺志砺志 博学博学 笃行笃行10（2）计算此单元的动能和势能杆单元的动能可表示成：（6）上式中，是材料的密度，A是杆单元的横截面积。明德明德 砺志砺志 博学博学 笃行笃行11用矩阵形式表示（6）式为：（7）其中，称为单元广义速度列阵（8）所以，质量矩阵可以认为是：（9）明德明德 砺志砺志 博学博学 笃行笃行12杆单元的势能可以写成：（10）式中，E是弹性模量，（10）表示成矩阵形式为：（11）这里，所以刚度矩阵k可以表示成：（12）明德明德 砺志砺志 博学博学 笃

6、行笃行 tftf21,13 tutu21,(3)计算等效节点力设单元上x处作用有分布力f ( x , t)，现在要把它等效成节点力遵循等效原则，即原载荷和等效之后的节点载荷在虚位移上所做的虚功相等。其实，就是对应于广义坐标的广义力，为此，计算txf,所做的虚功：把上式写成矩阵形式：（13）明德明德 砺志砺志 博学博学 笃行笃行14所以等效节点力可以写成：（14） 明德明德 砺志砺志 博学博学 笃行笃行15由杆单元组成的系统分析FR划分为三个杆单元，共有4个节点，每个节点处设一个广义坐标，写成一个阵列： u1F12（1）u2u223（2）u3u334（3）u4明德明德 砺志砺志 博学博学 笃行笃

7、行16每个杆单元有两个广义坐标，用列阵U来改写单元微分方程，单元（1）的方程改写为： 表示在节点处单元i给单元j的力。明德明德 砺志砺志 博学博学 笃行笃行17单元格（2）：单元格（3）：明德明德 砺志砺志 博学博学 笃行笃行18 系统的质量矩阵、刚度矩阵和广义力矩阵系统的质量矩阵、刚度矩阵和广义力矩阵明德明德 砺志砺志 博学博学 笃行笃行19当由个单元的广义力列阵叠加成系统广义力列阵F时，由于 引入边界条件进行处理，在右端，杆被约束，既不肯能发生位移， 式中： 为待求的广义坐标列阵为待求的广义坐标列阵 为已知的被约束的广义坐标为已知的被约束的广义坐标 为已知的广义力列阵为已知的广义力列阵 为

8、待求约束反力为待求约束反力明德明德 砺志砺志 博学博学 笃行笃行20 明德明德 砺志砺志 博学博学 笃行笃行 ,31tftf21tftf42, 二 梁单元如下图所示，一个梁单元也是有两个节点，但是有四个自由度，每个节点处，有两种位移形式，一个是线位移，即挠度，一种是角位移。图中，是力，是力矩。txf,是分布载荷 ,31tt tt42, 是对应的线位移，是对应的转角。tx,是梁单元上任意位移 x处的挠度。图12.2明德明德 砺志砺志 博学博学 笃行笃行22在静载弯曲条件下，梁单元上任意点出的挠度是x的三次方程，可写成：此方程必须满足下面的边界条件：由此可以求解处a (t)、b (t)、c (t)

9、、d (t)，进而挠度方程为：（16）（17）（18）明德明德 砺志砺志 博学博学 笃行笃行23上式可以写成形状函数的表示：其中，形函数分别为：梁单元的动能、势能、虚功表达式分别为：（19）明德明德 砺志砺志 博学博学 笃行笃行24式中I是横截面的惯性矩上式中：（20）（21）（22）明德明德 砺志砺志 博学博学 笃行笃行25通过上式，可以得到梁单元的质量、刚度矩阵，等效节点力：明德明德 砺志砺志 博学博学 笃行笃行26平面单元一 三角形单元 平面三角形单元 ui (Ui ) um (Um ) uj (Uj ) vj (Vj ) vi (Vi ) um (Um ) j i m x y o 明德

10、明德 砺志砺志 博学博学 笃行笃行27首先，我们来分析一下三角形单元的力学特性，即建立以单元节点位移表示单元内各点位移的关系式。设单元e的节点编号为i、j、m，如图3-2所示。由弹性力学平面问题可知，每个节点在其单元平面内的位移可以有两个分量，所以整个三角形单元将有六个节点位移分量，即六个自由度。用列阵可表示为： TmmjjiiTTmTjTievuvuvu Tiiivu其中的子矩阵（i，j，m 轮换） (a)式中 ui、vi 是节点i在x轴和y轴方向的位移。(5-7)明德明德 砺志砺志 博学博学 笃行笃行28uxyvxy123456将位移在x轴和y轴的分量设为u，v，则：式中 1、2、6是待定

11、常数。因三角形单元共有六个自由度，且位移函数u、v在三个节点处的数值应该等于这些点处的位移分量的数值。假设节点i、j、m的坐标分别为（xi , yi ）、（xj , yj ）、（xm , ym ），代入 (b) 式，得：(b)uxyvxyuxyvxyuxyvxyiiijiijjjjjjmmmmmm123456123456123456 , , , (c)明德明德 砺志砺志 博学博学 笃行笃行29mmjjiimmjjiimmmjjjiiiuxuxuxyuyuyuyxuyxuyxu11121 , 11121 , 213212111 xyxyxyiijjmm由 (c) 式左边的三个方程可以求得 (d)

12、其中(5-8) 从解析几何可知，式中的 就是三角形i、j、m的面积。为保证求得的面积为正值，节点i、j、m的编排次序必须是逆时针方向，如图5-2所示。明德明德 砺志砺志 博学博学 笃行笃行 ui (Ui ) um (Um ) uj (Uj ) vj (Vj ) vi (Vi ) um (Um ) j i m x y o 30 图5-2 平面三角形单元将 (d) 式代入 (b) 式的第一式，经整理后得到 uab xc y uab xc y uab xc y uiiiijjjjmmmm12(e)明德明德 砺志砺志 博学博学 笃行笃行31mjmjimjmjijmmjmmjjixxxxcyyyybyx

13、yxyxyxa1111vab xc y vab xc y vab xc y viiiijjjjmmmm12ycxbaNiiii21其中同理可得若令这样，位移模式 (e) 和 (f) 就可以写为（i , j , m轮换） (5-10)（i , j , m轮换） (5-9)明德明德 砺志砺志 博学博学 笃行笃行32 式中 I是二阶单位矩阵；Ni 、Nj 、Nm 是坐标的函数，它们反映了单元的位移状态，所以一般称之为形状函数，简称形函数。矩阵 N 叫做形函数矩阵。三节点三角形单元的形函数是坐标的线性函数。单元中任一条直线发生位移后仍为一条直线，即只要两单元在公共节点处保持位移相等。则公共边线变形后仍

14、为密合。 uN uN uN uvN vN vN viijjmmiijjmm eemjiNINININvuf(5-11)也可写成矩阵形式明德明德 砺志砺志 博学博学 笃行笃行33三、应三、应 变变 xyxyuxvyuyvx 12000000bbbccccbcbcbijmijmiijjmme有了单元的位移模式，就可以利用平面问题的几何方程求得应变分量。将 (e) 、(f) 两式代入上式，即得：(g)明德明德 砺志砺志 博学博学 笃行笃行34 Be BBBBijmBbccbiiiii1200可简写成 其中 B 矩阵叫做单元应变矩阵，可写成分块形式而子矩阵由于和bi 、bj 、bm 、ci 、cj 、

15、cm 等都是常量，所以矩阵B中的诸元素都是常量，因而单元中各点的应变分量也都是常量，通常称这种单元为常应变单元。 （i , j , m轮换） (3-15)(3-14)(3-13)明德明德 砺志砺志 博学博学 笃行笃行35四、应四、应 力力 eBD SD B Se D 求得应变之后，再将（3-13）式代入物理方程 ，便可推导出以节点位移表示的应力。即(5-16)(h)(5-17)令则明德明德 砺志砺志 博学博学 笃行笃行36 SD BBBSSSijmijm DE11100122对称 SD BEbcbccbiiiiiiii2 112122其中 S叫做应力矩阵，若写成分块形式，有对于平面应力问题，弹

16、性矩阵D为(5-18)(i)所以，S的子矩阵可记为（i , j , m轮换） (5-19)明德明德 砺志砺志 博学博学 笃行笃行37 SSSiijjmm注意到（5-7）式，则有(5-21) 由（5-19）、（5-20）式不难看出，S中的诸元素都是常量，所以每个单元中的应力分量也是常量。 可见，对于常应变单元，由于所选取的位移模式是线性的，因而其相邻单元将具有不同的应力和应变，即在单元的公共边界上应力和应变的值将会有突变，但位移却是连续的。明德明德 砺志砺志 博学博学 笃行笃行38单元的变形能 单元的刚度矩阵，为一6X6的矩阵。 单元的动能单元的质量矩阵，为一6X6的矩阵。明德明德 砺志砺志 博

17、学博学 笃行笃行39矩形单元v1yu11u2v22u1v14u1v13x矩形单元的位移假定采用双线性插值多项式明德明德 砺志砺志 博学博学 笃行笃行40等参单元四边形四节点等参数单元四边形四节点等参数单元1234022( , )iivNv 0i0i( , )iiuNu 001(1)(1)4iNl 四边形四节点单元位移模式四边形四节点单元位移模式: 其中其中母单元母单元明德明德 砺志砺志 博学博学 笃行笃行41 若以母单元上若以母单元上12边为例，边为例，通过映射可得在平面内任一直通过映射可得在平面内任一直线，线，12边的方程为边的方程为 = -1,代入代入( , )iixNx ( , )iiy

18、Ny 1211(1)(1)22xxxcd12341(1)(1)(1)(1)(1)(1)(1)(1)4iixN xxxxx 通过局部座标与整体座标的映射关系把母单元变换到通过局部座标与整体座标的映射关系把母单元变换到整体座标上成为一个任意四边形用于离散结构物，它能适整体座标上成为一个任意四边形用于离散结构物，它能适合于任意曲边的形状。合于任意曲边的形状。oxy12341o111明德明德 砺志砺志 博学博学 笃行笃行42同理可得同理可得12341(1)(1)(1)(1)(1)(1)(1)(1)4iiyN yyyyy1211(1)(1)22yyyef把以参数把以参数 代表的代表的 x方程和方程和 y

19、方程消去方程消去 ，则得则得 x , y 所组成的直线方程所组成的直线方程 y=kx+b所以母单元上的四个边都可以通过映射在所以母单元上的四个边都可以通过映射在x, y座标面座标面上得出一个任意四边形，用该四边形离散结构物。上得出一个任意四边形，用该四边形离散结构物。明德明德 砺志砺志 博学博学 笃行笃行43四边形四节点四边形四节点 单元的应变单元的应变,00i xii yi yi xNBNNN,ii yNNy 1234eeuxvBBBBByuvyx,ii xNNxi = 1,2,3,4 其中其中明德明德 砺志砺志 博学博学 笃行笃行44其中其中iiiNNNxyxy,xx,xx,yyiiiNN

20、Nxyxy,i xii yiNNxyNNxy,yy由于形函数由于形函数 N 是是 x, y的函数，现对的函数，现对, 求导，这是复合函数求导求导，这是复合函数求导把上二式写成矩阵形式把上二式写成矩阵形式明德明德 砺志砺志 博学博学 笃行笃行45其中其中 ,xyJxy 1,i xii yiNNJNN,Jx yy x 1,1yyJxxJ称为雅可比称为雅可比(Jacobian)矩阵矩阵 为了把为了把 Bi矩阵中的矩阵中的Ni,x和和 Ni,y化成化成 Ni,和和 Ni, ，则代入下式则代入下式其中雅可比矩阵的逆阵由下式给出其中雅可比矩阵的逆阵由下式给出明德明德 砺志砺志 博学博学 笃行笃行也可把形函

21、数也可把形函数Ni对的对的, 偏导数写成通式偏导数写成通式: i= 1,2,3,446,2,1111122i xi yiii xi yi yi xNNESDBNNNN xeeyxyDBSl 四边形四节点四边形四节点 单元的应力单元的应力对于平面应力情况对于平面应力情况明德明德 砺志砺志 博学博学 笃行笃行47 1111eTTKBDB tdxdyBDB t J d d 11121314212223243132333441424344ekkkkkkkkKkkkkkkkk 1111TTijijijkBDBtdxdyBDBt J d d l 四边形四节点单元的刚度矩阵四边形四节点单元的刚度矩阵Ke:Ke可划分成四行四列的子矩阵可划分成四行四列的子矩阵 :i, j = 1,2,3,4 明德明德 砺志砺志 博学博学 笃行笃行48对于平面应力情况对于平面应力情况: ,.,2,112211122i xj xi yj yi xj yi yj xTiji yj xi xj yi yj yi xj xNNNNNNNNEBDBNNNNNNNN TTePNp tdxdyNp t J d d l 等效节点力计算等效节点力计算: 1. 集中力的等效节点力集中力的等效节点力:由于计算复杂，把有集中力处