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

明德砺志博学笃行基于有限元法的振动分析陕西理工学院机械工程郑佳文1明德砺志博学笃行内容有限元法简介一维单元杆单元梁单元二维单元与平面问题的有限元法三角形单元矩形单元等参单元平面问题的有限元法2明德砺志博学笃行有限元法简介有限元法是一种可用于精确地(但近似)解决许多复杂的振动问题的数值方法。对于基本的一维元素进行有限元分析，能得到质量矩阵与刚度矩阵和所需的力矢量，对于二维三维，元素矩阵会转换成相关的更高维的空间。使用一致的和集中质量矩阵的有限元方程并结合边界条件能为复杂系统提供解释。3明德砺志博学笃行当用有限元法来分析时，将梁分解成为梁单元，将各单元相邻的节点处的弹性位移（横向弹性位移和弹性转角）作为问题中的位置量。列子4单元质量矩阵单元刚度矩阵集合系统质量矩阵M系统刚度矩阵K梁的有限元模型明德砺志博学笃行有限元基本思想讲一个连续弹性体看成是由若干个基本单元在节点彼此相连的组合体；由于单元较小，允许对位移分布规律做出某种假设，基于这一简化假定，回避了弹性力学中的微分方程组难以求解的困难；从而使这一无限自由的连续体问题变成一个有限自由度的离散系统问题。5明德砺志博学笃行一维单元一杆单元一个杆单元是从杆上划分出的一个小段，如下图所示。由于单元很小，ρ、A均视为常量。现在就以这最简单的杆单元，推导出它的质量、刚度矩阵，等效节点力。6图1.2明德砺志博学笃行)()(),(2211tututxu7(1)求杆单元上任意点的位移u(x,t)本来，杆单元上任意点的位移u(x,t)与节点的位移u1(t)、u2(t)之间的关系是未知的，但是，只要单元划分的足够小，那么其间的关系就无关大局。所以可以假定它们之间有简单的线性关系，即根据节点位移对单元内任意点位移进行插值：（1）式中，Φ1、Φ2称为线性系数，与单元里点的位置有关，是x的函数。此函数与单元的形状有关，又叫形状函数。明德砺志博学笃行)(),0(1tutu8)(),(2tutlu形状函数和插值函数一样是任意的，但必须边界条件：只有满足此条件单元才能协调一致运动，而不致破坏系统的完整性，因此这两个条件实际上就是变形协调条件。将式（1）带入（2）中，就可以得到形状函数Φ1(x)、Φ2(x)所满足的边界条件：（2）1,000,102121ll（3）明德砺志博学笃行lxxlxx21,19xx21、我们可以用简单的线性函数来近似，因此有：（4）代回（1）式中有：tulxtulxtxu21)1(),(（5）为此，我们已经找到了用节点位移表示单元内任意一点位移的表达式。以上边界条件确定了由于这两个函数的任意性明德砺志博学笃行10（2）计算此单元的动能和势能杆单元的动能可表示成：（6）上式中，ρ是材料的密度，A是杆单元的横截面积。明德砺志博学笃行11用矩阵形式表示（6）式为：（7）其中，称为单元广义速度列阵（8）所以，质量矩阵可以认为是：（9）明德砺志博学笃行12杆单元的势能可以写成：（10）式中，E是弹性模量，（10）表示成矩阵形式为：（11）这里，，所以刚度矩阵[k]可以表示成：（12）明德砺志博学笃行tftf21,13tutu21,(3)计算等效节点力设单元上x处作用有分布力f(x,t)，现在要把它等效成节点力遵循等效原则，即原载荷和等效之后的节点载荷在虚位移上所做的虚功相等。其实，就是对应于广义坐标的广义力，为此，计算txf,所做的虚功：把上式写成矩阵形式：（13）明德砺志博学笃行14所以等效节点力可以写成：（14）明德砺志博学笃行15由杆单元组成的系统分析FR划分为三个杆单元，共有4个节点，每个节点处设一个广义坐标，写成一个阵列：u1F12（1）u2u223（2）u3u334（3）u4明德砺志博学笃行16每个杆单元有两个广义坐标，用列阵U来改写单元微分方程，单元（1）的方程改写为：表示在节点处单元i给单元j的力。明德砺志博学笃行17单元格（2）：单元格（3）：明德砺志博学笃行18系统的质量矩阵、刚度矩阵和广义力矩阵明德砺志博学笃行19当由个单元的广义力列阵叠加成系统广义力列阵F时，由于引入边界条件进行处理，在右端，杆被约束，既不肯能发生位移，式中：为待求的广义坐标列阵为已知的被约束的广义坐标为已知的广义力列阵为待求约束反力明德砺志博学笃行20明德砺志博学笃行,31tftf21tftf42,二梁单元如下图所示，一个梁单元也是有两个节点，但是有四个自由度，每个节点处，有两种位移形式，一个是线位移，即挠度，一种是角位移。图中，是力，是力矩。txf,是分布载荷,31tttt42,是对应的线位移，是对应的转角。tx,是梁单元上任意位移x处的挠度。图12.2明德砺志博学笃行22在静载弯曲条件下，梁单元上任意点出的挠度是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)jimxyo明德砺志博学笃行27首先，我们来分析一下三角形单元的力学特性，即建立以单元节点位移表示单元内各点位移的关系式。设单元e的节点编号为i、j、m，如图3-2所示。由弹性力学平面问题可知，每个节点在其单元平面内的位移可以有两个分量，所以整个三角形单元将有六个节点位移分量，即六个自由度。用列阵可表示为：TmmjjiiTTmTjTievuvuvuTiiivu其中的子矩阵（i，j，m轮换）(a)式中ui、vi是节点i在x轴和y轴方向的位移。(5-7)明德砺志博学笃行28uxyvxy123456将位移在x轴和y轴的分量设为u，v，则：式中1、2、…6是待定常数。因三角形单元共有六个自由度，且位移函数u、v在三个节点处的数值应该等于这些点处的位移分量的数值。假设节点i、j、m的坐标分别为（xi,yi）、（xj,yj）、（xm,ym），代入(b)式，得：(b)uxyvxyuxyvxyuxyvxyiiijiijjjjjjmmmmmm123456123456123456,,,(c)明德砺志博学笃行29mmjjiimmjjiimmmjjjiiiuxuxuxyuyuyuyxuyxuyxu11121,11121,213212111xyxyxyiijjmm由(c)式左边的三个方程可以求得(d)其中(5-8)从解析几何可知，式中的就是三角形i、j、m的面积。为保证求得的面积为正值，节点i、j、m的编排次序必须是逆时针方向，如图5-2所示。明德砺志博学笃行ui(Ui)um(Um)uj(Uj)vj(Vj)vi(Vi)um(Um)jimxyo30图5-2平面三角形单元将(d)式代入(b)式的第一式，经整理后得到uabxcyuabxcyuabxcyuiiiijjjjmmmm12(e)明德砺志博学笃行31mjmjimjmjijmmjmmjjixxxxcyyyybyxyxyxyxa1111vabxcyvabxcyvabxcyviiiijjjjmmmm12ycxbaNiiii21其中同理可得若令这样，位移模式(e)和(f)就可以写为（i,j,m轮换）(5-10)（i,j,m轮换）(5-9)明德砺志博学笃行32式中I是二阶单位矩阵；Ni、Nj、Nm是坐标的函数，它们反映了单元的位移状态，所以一般称之为形状函数，简称形函数。矩阵[N]叫做形函数矩阵。三节点三角形单元的形函数是坐标的线性函数。单元中任一条直线发生位移后仍为一条直线，即只要两单元在公共节点处保持位移相等。则公共边线变形后仍为密合。uNuNuNuvNvNvNviijjmmiijjmmeemjiNINININvuf(5-11)也可写成矩阵形式明德砺志博学笃行33三、应变xyxyuxvyuyvx12000000bbbccccbcbcbijmijmiijjmme有了单元的位移模式，就可以利用平面问题的几何方程求得应变分量。将(e)、(f)两式代入上式，即得：(g)明德砺志博学笃行34BeBBBBijmBbccbiiiii1200可简写成其中[B]矩阵叫做单元应变矩阵，可写成分块形式而子矩阵由于和bi、bj、bm、ci、cj、cm等都是常量，所以矩阵[B]中的诸元素都是常量，因而单元中各点的应变分量也都是常量，通常称这种单元为常应变单元。（i,j,m轮换）(3-15)(3-14)(3-13)明德砺志博学笃行35四、应力eBDSDBSeD求得应变之后，再将（3-13）式代入物理方程，便可推导出以节点位移表示的应力。即(5-16)(h)(5-17)令则明德砺志博学笃行36SDBBBSSSijmijmDE11100122对称SDBEbcbccbiiiiiiii2112122其中[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矩形单元v1yu11u2v22u1v14u1v13x矩形单元的位移假定采用双线性插值多项式明德砺志博学笃行40等参单元四边形四节点等参数单元1234ξη022(,)iivNv0i0i(,)iiuNu001(1)(1)4iN四边形四节点单元位移模式:其中母单元明德砺志博学笃行41若以母单元上12边为例，通过映射可得在平面内任一直线，12边的方程为η=-1,代入(,)iixNx(,)iiyNy1211(1)(1)22xxxcd12341(1)(1)(1)(1)(1)(1)(1)(1)4iixNxxxxx通过局部座标与整体座标的映射关系把母单元变换到整体座标上成为一个任意四边形用于离散结构物，它能适合于任意曲边的形状。oxy12341o111明德砺志博学笃行42同理可得12341(1)(1)(1)(1)(1)(1)(1)(1)4iiyNyyyyy1211(1)(1)22yyyef