用于振动分析的有限元方法

1、用于振动分析的有限元方法指导老师：陈益 报告人：成志斌 韩宗彪 何瑜 宁鹏内容有限元介绍单个元素的运动方程整个系统的运动方程整个系统的边界条件的加载及质量矩阵 MATLAB实例及总结单个元素的质量矩阵、刚度矩阵、力矢量及其转化有限元法简介 有限元法是一种可用于精确地(但近似)解决许多复杂的振动问题的数值方法。 对于基本的一维元素进行有限元分析，能得到质量矩阵与刚度矩阵和所需的力矢量，对于二维三维，元素矩阵会转换成相关的更高维的空间。使用一致的和集中质量矩阵的有限元方程并结合边界条件能为复杂系统提供解释。 最后,使用MATLAB程序得到在轴向载荷下的指定节点位移,固有振动频率和特征值分析。本章目

2、的 *认识用于解决不同类型振动问题的刚度和质量矩阵。 *将矩阵元素从局部坐标系变换到全球坐标系。 *装配单元矩阵和应用边界条件。 *对杆、梁元素进行静态分析。 *对杆、梁元素进行动态分析来得到固有频率和振型。 *在有限元振动分析使用一致的集中质量矩阵。 *使用MATLAB解决振动问题。有限元思想 1，实际结构被一些元素所取代,这些元素都是被假定为一个连续的结构部件即有限元，这些元素在特定点即节点上互相关联。 2，如果解决方案的各方面都选择得当,那么它可以收敛到精确的解决方案，因为组成总体结构的元素很小，在节点上的力的平衡和元素之间的位移都令人感到满意,这样整个结构(组合的元素)表现为单一实体。

3、 3，因为得到准确解很难，所以得到一个方便且逼近的近似解很有价值。元素的运动方程龙门刨铣床有限元模型三角板元素梁元素元素的运动方程 位移函数 形状函数 各点对应位移 未知节点位移数 n 动能 应变能质量矩阵刚度矩阵 位移函数 形状函数 各点对应位移 未知节点位移数 n 动能 应变能主要内容： 一，单元的质量、刚度矩阵、等效节点力 二，单元矩阵的坐标变换 三，整个系统的运动方程一 单元的质量、刚度矩阵，等效节点力矢量一 杆单元一个杆单元是从杆上划分出的一个小段，如下图所示。由于单元很小，、A均视为常量。现在就以这最简单的杆单元，推导出它的质量、刚度矩阵，等效节点力。图12.1图12.1 (1)

4、求杆单元上任意点的位移u(x,t)本来，杆单元上任意点的位移u(x,t)与节点的位移u1(t)、u2(t)之间的关系是未知的，但是，只要单元划分的足够小，那么其间的关系就无关大局。所以可以假定它们之间有简单的线性关系，即根据节点位移对单元内任意点位移进行插值：)()(),(2211tututxu（1）式中，1、 2称为线性系数，与单元里点的位置有关，是x的函数。此函数与单元的形状有关，又叫形状函数。形状函数和插值函数一样是任意的，但必须边界条件：)(), 0(1tutu)(),(2tutlu只有满足此条件单元才能协调一致运动，而不致破坏系统的完整性，因此这两个条件实际上就是变形协调条件。将式（

5、1）带入（2）中，就可以得到形状函数1(x)、2(x)所满足的边界条件： （2） 1 , 00 0 , 102121ll （3） 以上边界条件确定了 xx21、 lxxlxx21 ,1由于这两个函数的任意性我们可以用简单的线性函数来近似，因此有：（4）代回（1）式中有： tulxtulxtxu21)1 (),(（5）为此，我们已经找到了用节点位移表示单元内任意一点位移的表达式。 （2）计算此单元的动能和势能杆单元的动能可表示成：（6） 上式中，是材料的密度，A是杆单元的横截面积。用矩阵形式表示（6）式为：（7）其中，（8）所以，质量矩阵可以认为是：（9） 杆单元的势能可以写成：（10）式中，E

6、是弹性模量，（10）表示成矩阵形式为：这里，所以刚度矩阵k可以表示成：（11）（12）(3)计算等效节点力设单元上x处作用有分布力f ( x , t)，现在要把它等效成节点力 tftf21,遵循等效原则，即原载荷和等效之后的节点载荷在虚位移上所做的虚功相等。其实，就是对应于广义坐标 tutu21,的广义力，为此，计算txf,所做的虚功：把上式写成矩阵形式：（13）所以等效节点力可以写成： tftuktum （14）二 梁单元 ,31tftf如下图所示，一个梁单元也是有两个节点，但是有四个自由度，每个节点处，有两种位移形式，一个是线位移，即挠度，一种是角位移。图中，tftf42, 是力，是力矩。

7、txf,是分布载荷 ,31tt tt42, 是对应的线位移，是对应的转角。tx,是梁单元上任意位移 x处的挠度。图12.2在静载弯曲条件下，梁单元上任意点出的挠度是x的三次方程，可写成：此方程必须满足下面的边界条件：由此可以求解处a (t)、b (t)、c (t)、d (t)，进而挠度方程为：（16）（17）（18）tan挠度的斜率上式可以写成形状函数的表示：其中，形函数分别为：梁单元的动能、势能、虚功表达式分别为：（19）式中I是横截面的惯性矩上式中：（20）（21）（22）通过上式，可以得到梁单元的质量、刚度矩阵，等效节点力：二 单元矩阵的坐标变换局部坐标系：以各个单元本身的轴线为基准所设

8、立的坐标系。便于计算节点位移。缺点：如果整个系统里各个单元取向各异，各个节点位移方向不一致，如下图。如何使汇交于一个节点的各个杆件的节点位移真正相等?解决方法：进行坐标变换 如右图的系统，有四个杆件， u 1(t) 、 u 2(t) 为局部坐标系的节点位移，U i 为全局坐标系下的位移图12.3如下图，节点位移在局部、全局坐标系中的关系：坐标变换矩阵（23）其中，因为单元的动能、势能与坐标系无关:(24) 得到在全局坐标下的单元质量、刚度矩阵为：类似地，根据单元在两个坐标系下的力所做的虚功相等：得到在全局坐标系下的等效节点力：三 全系统运动方程 经过坐标变换，各个单元的节点位移方向被统一起来，

9、但是不同的节点有不同的节点位移，为了便于综合出全系统的运动方程，首先要建立全系统的节点位移向量。每个单元的节点位移向量与全系统的节点位移向量之间的关系：长方形矩阵由1和0组成单元节点位移向量单元节点位移向量 例如，图的12.5中的1单元，方程，变为：把每个单元的动能相加，就得到了整个系统的动能：（把整个系统的动能表示成关于节点速度矢量的形式） 这样就得到了整个系统的质量矩阵:类似地，考虑整个系统的势能，便可以得到整个系统的刚度矩阵：整个系统的广义力向量：最后得到整个系统的运动方程：12.6添加边界条件 前文中，节点没有固定，结构在节点力的作用下会发生刚体位移。 也就是说，矩阵K是奇异矩阵。通常

10、情况下，我们希望结构 的位移为零。 因此，我们需要添加边界条件对矩阵M、K和向量F进行约束。 N：结构中自由节点位移的数目例1：杆件分析 如图：均质；长0.5m；断面截面积5e-4m2；杨氏模量200GPa；密度7850Kg/m3；左端固定。 a.节点2处施加1000N静态轴向外力u2，求应力 b.求系统固有频率 解： a： 平衡方程： A=5e-4，E=2e11，l = 0.5，f2=1000，代入方程得： u1：位移，f1：节点1处应力，添加边界条件：u1=0，解得：u2=5e-10m 由应力 与应变 的关系： 表示长度的变化， 表示应变； b： 由刚度矩阵和质量矩阵，得特征值方程： 21

11、luu ll 式中w位固有频率，U1、U2分别是节点1、2的振幅，添加边界条件：U1=0；解得： 例2：梁的自然频率解：梁被理想化为单一单元，局部和整体的节点位移相同，如 图所示： 梁的刚度矩阵： 质量矩阵： 节点位移向量： 与端点相关的边界条件：W1=0，W3=0；解得： 求解特征值： 乘以l/2EI得： 令系数矩阵的行列式等于0得： 方程的根即梁的自然频率： 结果可以和精确解比较：12.7一致、集中质量矩阵 12.3节中推出的质量矩阵是一致质量矩阵，因为用于推导刚度矩阵的位移模型也用于推导质量矩阵。 一些动态问题可以用形式简单的质量矩阵求解。最简单的质量矩阵集中质量矩阵，可以通过将质点指定

12、到节点上。 集中质量针对平移和旋转的元素, 假设在平均位置两侧的特定位移表现得像个刚体而剩余的元素不参与运动。 因此这种假设不包括元素位移之间存在的动态耦合,因此产生的元素质量矩阵是纯粹的对角矩阵。 1、杆的集中质量矩阵： 2、梁的集中质量矩阵： 旋转自由度的惯性影响被假定为0；若考虑惯性影响，有转动惯量： 集中质量矩阵变为： 对于一般的动态问题，两者谁能得到更精确的解？ 两个质量矩阵很相似， 他们不考虑各种位移自由度的元素之间的动态耦合。 他们的形状函数也近似，都是用静态的位移模型推导而来。 然而,由于集中质量矩阵对角，在计算时他使用更少的存储空间。 下面的例子说明了在一个简单的振动问题中，

13、集中和一致质量矩阵的应用。 集中质量矩阵与一致质量矩阵： 例：杆的一致和集中质量矩阵 用一致和集中质量矩阵求如图所示两端固定杆的固有频率，用两个杆单元建模。解：单元的刚度和质量矩阵分别是： 质量矩阵的下标c和l分别表示一致和 集中质量矩阵。 由于该杆由两个单元建模，组合的刚度和质量矩阵如下： 方框中的部分分别与单元1和2相关。 添加边界条件U1=U3=0后，特征值问题为： 特征值w由以下方程课解： 代入已知条件得：用一致质量矩阵用集中质量矩阵 事实上，方程的精确解析解是： 解得：12.8 MATLAB应用举例 例12.5 阶梯轴的有限元分析 图12.11中的阶梯轴满足一下条件：A1=1610-

14、4 m2 ， A2=910-4 m2 ， A3=410-4 m2 ，Ei=201010Pa，i=1,2,3，pi=7.8103Kg/m3， i=1,2,3，l1=1m，l2=0.5m，l3=0.25m。编写一个MATLAB程序解决以下问题。 a，在在载荷p3=1000N下u1，u2，u3的位移 b，阶梯轴的固有频率和模态12.8 MATLAB应用举例 解决方案：阶梯轴的刚度矩阵和质量矩阵如下所示： 在载荷p3的作用下系统的平衡方程如下所示：解得：b.由下面方程可求得特征值K已在（E.4）中得到，M可由下式得到：由MATLAB求解（E.3）和（E.5）的程序如下所示例12.6 阶梯梁的特征值分析的程序如图12.12所示的两端固定的阶梯梁，编写一个MATLAB程序，命名为Program17.m，对它进行特征值分析。Program17.m程序的输入数据如下：xl(i)=阶梯梁的长度 ixi(i)=阶梯梁的转动惯量 ia(i)=阶梯梁的横截面积 ib(i，j)=对应本地的自由度 j 和阶梯梁的自由度 ie