1、基于有限元法的振动分析基于有限元法的振动分析陕西理工学院陕西理工学院 机械工程机械工程 郑佳文郑佳文1内容n有限元法简介n一维单元 杆单元 梁单元n二维单元与平面问题的有限元法 三角形单元 矩形单元 等参单元n 平面问题的有限元法2有限元法简介 有限元法是一种可用于精确地(但近似)解决许多复杂的振动问题的数值方法。对于基本的一维元素进行有限元分析，能得到质量矩阵与刚度矩阵和所需的力矢量，对于二维三维，元素矩阵会转换成相关的更高维的空间。使用一致的和集中质量矩阵的有限元方程并结合边界条件能为复杂系统提供解释。3列子当用有限元法来分析时，将梁分解成为梁单元，将各单元相邻的节点处的弹性位移（横向弹性
2、位移和弹性转角）作为问题中的位置量。4单元质量矩阵单元质量矩阵单元刚度矩阵单元刚度矩阵集合系统质量矩阵系统质量矩阵M系统刚度矩阵系统刚度矩阵K 梁的有限元模型有限元基本思想讲一个连续弹性体看成是由若干个基本单元在节点彼此相连的组合体；由于单元较小，允许对位移分布规律做出某种假设，基于这一简化假定，回避了弹性力学中的微分方程组难以求解的困难；从而使这一无限自由的连续体问题变成一个有限自由度的离散系统问题。5一维单元一 杆单元 一个杆单元是从杆上划分出的一个小段，如下图所示。由于单元很小，、A均视为常量。现在就以这最简单的杆单元，推导出它的质量、刚度矩阵，等效节点力。6图1.27(1)求杆单元上任
3、意点的位移u(x,t)本来，杆单元上任意点的位移u(x,t)与节点的位移u1(t)、u2(t)之间的关系是未知的，但是，只要单元划分的足够小，那么其间的关系就无关大局。所以可以假定它们之间有简单的线性关系，即根据节点位移对单元内任意点位移进行插值：)()(),(2211tututxu（1）式中，1、2称为线性系数，与单元里点的位置有关，是x的函数。此函数与单元的形状有关，又叫形状函数。8形状函数和插值函数一样是任意的，但必须边界条件：)(),0(1tutu)(),(2tutlu只有满足此条件单元才能协调一致运动，而不致破坏系统的完整性，因此这两个条件实际上就是变形协调条件。将式（1）带入（2）
4、中，就可以得到形状函数1(x)、2(x)所满足的边界条件：（2）1 ,00 0 ,102121ll （3）9 lxxlxx21 ,1我们可以用简单的线性函数来近似，因此有：（4）代回（1）式中有：tulxtulxtxu21)1(),(（5）为此，我们已经找到了用节点位移表示单元内任意一点位移的表达式。以上边界条件确定了 xx21、由于这两个函数的任意性10（2）计算此单元的动能和势能杆单元的动能可表示成：（6）上式中，是材料的密度，A是杆单元的横截面积。11用矩阵形式表示（6）式为：（7）其中，称为单元广义速度列阵（8）所以，质量矩阵可以认为是：（9）12杆单元的势能可以写成：（10）式中，E
5、是弹性模量，（10）表示成矩阵形式为：（11）这里，所以刚度矩阵k可以表示成：（12）13(3)计算等效节点力设单元上x处作用有分布力f(x,t)，现在要把它等效成节点力tftf21,遵循等效原则，即原载荷和等效之后的节点载荷在虚位移上所做的虚功相等。其实，就是对应于广义坐标 tutu21,的广义力，为此，计算txf,所做的虚功：把上式写成矩阵形式：（13）14所以等效节点力可以写成：（14）15由杆单元组成的系统分析FR划分为三个杆单元，共有4个节点，每个节点处设一个广义坐标，写成一个阵列：u1F12（1）u2u223（2）u3u334（3）u416每个杆单元有两个广义坐标，用列阵U来改写单
6、元微分方程，单元（1）的方程改写为：表示在节点处单元i给单元j的力。17单元格（2）：单元格（3）：18系统系统的质量矩阵、刚度矩阵和广义力矩阵的质量矩阵、刚度矩阵和广义力矩阵19当由个单元的广义力列阵叠加成系统广义力列阵F时，由于引入边界条件进行处理，在右端，杆被约束，既不肯能发生位移，式中：为待求的广义坐标列阵为已知的被约束的广义坐标为已知的广义力列阵为待求约束反力2021二 梁单元 ,31tftf如下图所示，一个梁单元也是有两个节点，但是有四个自由度，每个节点处，有两种位移形式，一个是线位移，即挠度，一种是角位移。图中，tftf42,是力，是力矩。txf,是分布载荷 ,31tt tt42
7、,是对应的线位移，是对应的转角。tx,是梁单元上任意位移 x处的挠度。图12.222在静载弯曲条件下，梁单元上任意点出的挠度是x的三次方程，可写成：此方程必须满足下面的边界条件：由此可以求解处a(t)、b(t)、c(t)、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
8、)j i m x y o 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是待定常数。因三角形单元共有六个自由度，且位移函数u、v
9、在三个节点处的数值应该等于这些点处的位移分量的数值。假设节点i、j、m的坐标分别为（xi ,yi）、（xj ,yj）、（xm ,ym），代入(b)式，得：(b)uxyvxyuxyvxyuxyvxyiiijiijjjjjjmmmmmm123456123456123456 ,(c)29mmjjiimmjjiimmmjjjiiiuxuxuxyuyuyuyxuyxuyxu11121,11121,213212111 xyxyxyiijjmm由(c)式左边的三个方程可以求得 (d)其中(5-8)从解析几何可知，式中的 就是三角形i、j、m的面积。为保证求得的面积为正值，节点i、j、m的编排次序必须是逆时针
10、方向，如图5-2所示。30 图5-2 平面三角形单元 ui(Ui)um(Um)uj(Uj)vj(Vj)vi(Vi)um(Um)j i m x y o 将(d)式代入(b)式的第一式，经整理后得到 uab xc y uab xc y uab xc y uiiiijjjjmmmm12(e)31mjmjimjmjijmmjmmjjixxxxcyyyybyxyxyxyxa1111vab xc y vab xc y vab xc y viiiijjjjmmmm12ycxbaNiiii21其中同理可得若令这样，位移模式(e)和(f)就可以写为（i,j,m轮换）(5-10)（i,j,m轮换）(5-9)32
11、式中 I是二阶单位矩阵；Ni 、Nj 、Nm 是坐标的函数，它们反映了单元的位移状态，所以一般称之为形状函数，简称形函数。矩阵 N 叫做形函数矩阵。三节点三角形单元的形函数是坐标的线性函数。单元中任一条直线发生位移后仍为一条直线，即只要两单元在公共节点处保持位移相等。则公共边线变形后仍为密合。uN uN uN uvN vN vN viijjmmiijjmm eemjiNINININvuf(5-11)也可写成矩阵形式33三、应三、应 变变 xyxyuxvyuyvx 12000000bbbccccbcbcbijmijmiijjmme有了单元的位移模式，就可以利用平面问题的几何方程求得应变分量。将(
12、e)、(f)两式代入上式，即得：(g)34 Be BBBBijmBbccbiiiii1200可简写成 其中 B 矩阵叫做单元应变矩阵，可写成分块形式而子矩阵由于和bi、bj、bm、ci、cj、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 BE
13、bcbccbiiiiiiii2 112122其中 S叫做应力矩阵，若写成分块形式，有对于平面应力问题，弹性矩阵D为(5-18)(i)所以，S的子矩阵可记为（i,j,m轮换）(5-19)37 SSSiijjmm注意到（5-7）式，则有(5-21)由（5-19）、（5-20）式不难看出，S中的诸元素都是常量，所以每个单元中的应力分量也是常量。可见，对于常应变单元，由于所选取的位移模式是线性的，因而其相邻单元将具有不同的应力和应变，即在单元的公共边界上应力和应变的值将会有突变，但位移却是连续的。38单元的变形能单元的刚度矩阵，为一6X6的矩阵。单元的动能单元的质量矩阵，为一6X6的矩阵。39矩形单元
14、yu1v11u2v22u1v14u1v13x矩形单元的位移假定采用双线性插值多项式40等参单元四边形四节点等参数单元四边形四节点等参数单元1234022(,)iivNv 0i0i(,)iiuNu 001(1)(1)4iNl 四边形四节点单元位移模式四边形四节点单元位移模式:其中其中母单元母单元41 若以母单元上若以母单元上12边为例，边为例，通过映射可得在平面内任一直通过映射可得在平面内任一直线，线，12边的方程为边的方程为=-1,代入代入(,)iixNx(,)iiyNy 1211(1)(1)22xxxcd12341(1)(1)(1)(1)(1)(1)(1)(1)4iixN xxxxx 通过局
15、部座标与整体座标的映射关系把母单元变换到通过局部座标与整体座标的映射关系把母单元变换到整体座标上成为一个任意四边形用于离散结构物，它能适整体座标上成为一个任意四边形用于离散结构物，它能适合于任意曲边的形状。合于任意曲边的形状。oxy12341o11142同理可得同理可得12341(1)(1)(1)(1)(1)(1)(1)(1)4iiyN yyyyy1211(1)(1)22yyyef把以参数把以参数 代表的代表的 x方程和方程和 y方程消去方程消去 ，则得则得 x,y 所组成的直线方程所组成的直线方程 y=kx+b所以母单元上的四个边都可以通过映射在所以母单元上的四个边都可以通过映射在x,y座标
16、面座标面上得出一个任意四边形，用该四边形离散结构物。上得出一个任意四边形，用该四边形离散结构物。43四边形四节点四边形四节点 单元的应变单元的应变,00i xii yi yi xNBNNN,ii yNNy 1234eeuxvBBBBByuvyx,ii xNNxi=1,2,3,4 其中其中44其中其中iiiNNNxyxy,xx,xx,yyiiiNNNxyxy,i xii yiNNxyNNxy,yy由于形函数由于形函数 N 是是 x,y的函数，现对的函数，现对,求导，这是复合函数求导求导，这是复合函数求导把上二式写成矩阵形式把上二式写成矩阵形式45其中其中,xyJxy 1,i xii yiNNJN
17、N,Jx yy x 1,1yyJxxJ称为雅可比称为雅可比(Jacobian)矩阵矩阵 为了把为了把 Bi矩阵中的矩阵中的Ni,x和和 Ni,y化成化成 Ni,和和 Ni,，则代入下式则代入下式其中雅可比矩阵的逆阵由下式给出其中雅可比矩阵的逆阵由下式给出46也可把形函数也可把形函数Ni对的对的,偏导数写成通式偏导数写成通式:i=1,2,3,4,2,1111122i xi yiii xi yi yi xNNESDBNNNN xeeyxyDBSl 四边形四节点四边形四节点 单元的应力单元的应力对于平面应力情况对于平面应力情况47 1111eTTKBDB tdxdyBDB t J d d 11121
18、314212223243132333441424344ekkkkkkkkKkkkkkkkk 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.集中力的等
19、效节点力集中力的等效节点力:由于计算复杂，把有集中力处由于计算复杂，把有集中力处设置为节点。设置为节点。2.体枳力的等效节点力体枳力的等效节点力:两类平面问题：（1）平面应力问题对于厚度远小于其他尺寸的薄板，且弹性体仅受平行于板面沿厚度方向不变化的外力作用时，可按平面应力问题来处理49平面问题有限元法（2）平面应变问题对于长度大于其横向尺寸的弹性体，且仅平行于截面、沿其长度方向不变化的作用力时，可按平面应变问题来处理，如下图50弹性矩阵将弹性构件划分为单元时，要确定单元的种类、数目和布局。单元种类可根据需要选择。选择精度比较高的单元，单元数目就可以少一些。单元尺度小一些，但单元的数目就要增加，计算精度自然会提高，但这涉及计算的时间和成本，有时还要考虑计算机容量的限制。在应力集中的部位以及应力比较剧烈的区域，网路的划分宜稠密一些；反之，在应力较小的区域，单元的尺度可大一些，但疏密的过渡应较为平缓。一个单元的各边长度不宜过为悬殊。构件上载荷的作用点除应取为节点。各个单元的节点也同时应该是相邻单元的节点，才能保证变形协调条件。51网格划分52举例L设置材料属性，划分网格，添加约束。53求出固有频率动画动画演示演示一一阶阶模态模态5455二阶模态二阶模态56三三阶模态阶模态57四四阶模态阶模态58五五阶模态阶模态
