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

x l
)u1
t


x l
u2 t
（5）
为此，我们已经找到了用节点位移表示单元内任意一点位移的表达式。
▪ （2）计算此单元的动能和势能
杆单元的动能可表示成：
（6）
▪ 上式中，ρ是材料的密度，A是杆单元的横截面积。 用矩阵形式表示（6）式为：
其中， 所以，质量矩阵可以认为是：
（7） （8） （9）
一 单元的质量、刚度矩阵，等效节点力矢量 一 杆单元
一个杆单元是从杆上划分出的一个小段，如下图所示。 由于单元很小，ρ、A均视为常量。现在就以这最简单 的杆单元，推导出它的质量、刚度矩阵，等效节点力。
图12.1
▪ (1) 求杆单元上任意点的位移u(x,t)
本来，杆单元上任意点的位移u(x,t)与节点的位移u1(t)、u2(t)之间 的关系是未知的，但是，只要单元划分的足够小，那么其间的关系 就无关大局。所以可以假定它们之间有简单的线性关系，即根据节 点位移对单元内任意点位移进行插值：
集中质量矩阵变为：
集中质量矩阵与一致质量矩阵： 对于一般的动态问题，两者谁能得到更精确的解？ 两个质量矩阵很相似， 他们不考虑各种位移自由度的元素之间
的动态耦合。 他们的形状函数也近似，都是用静态的位移模型推导而来。 然而,由于集中质量矩阵对角，在计算时他使用更少的存储空
间。 下面的例子说明了在一个简单的振动问题中，集中和一致质量
如右图的系统，有四个 杆件， u 1(t) 、 u 2(t) 为局部坐标系的节点位移， U i 为全局坐标系下的位移
图12.3
如下图，节点位移在局部、全局坐标系中的关系：
（23）
坐标变换矩阵
其中， 因为单元的动能、势能与坐标系无关:
(24)
▪ 得到在全局坐标下的单元质量、刚度矩阵为：
类似地，根据单元在两个坐标系下的力所做的虚功相等： 得到在全局坐标系下的等效节点力：
▪ 对于基本的一维元素进行有限元分析，能得到质量矩 阵与刚度矩阵和所需的力矢量，对于二维三维，元素矩阵 会转换成相关的更高维的空间。使用一致的和集中质量矩 阵的有限元方程并结合边界条件能为复杂系统提供解释。
▪ 最后,使用MATLAB程序得到在轴向载荷下的指定节点 位移,固有振动频率和特征值分析。
本章目的
矩阵的应用。
例：杆的一致和集中质量矩阵
用一致和集中质量矩阵求如图所示两端固定杆的固有频率，用 两个杆单元建模。 解：单元的刚度和质量矩阵分别是：
质量矩阵的下标c和l分别表示一致和 集中质量矩阵。
由于该杆由两个单元建模，组合的刚度和质量矩阵如下： 方框中的部分分别与单元1和2相关。
添加边界条件U1=U3=0后，特征值问题为： 特征值w由以下方程课解： 代入已知条件得：
集中质量针对平移和旋转的元素, 假设在平均位置两侧的特定 位移表现得像个刚体而剩余的元素不参与运动。
因此这种假设不包括元素位移之间存在的动态耦合,因此产生 的元素质量矩阵是纯粹的对角矩阵。
1、杆的集中质量矩阵： 2、梁的集中质量矩阵： 旋转自由度的惯性影响被假定为0；若考虑惯性影响，有转动 惯量：
l
b： 由刚度矩阵和质量矩阵，得特征值方程：
式中w位固有频率，U1、U2分别是节点1、2的振幅，添加边 界条件：U1=0；解得：
例2：梁的自然频率 解：梁被理想化为单一单元，局部和整体的节点位移相同，如 图所示：
梁的刚度矩阵：
质量矩阵： 节点位移向量： 与端点相关的边界条件：W1=0，W3=0；解得：
求解特征值： 乘以l/2EI得： 令系数矩阵的行列式等于0得： 方程的根即梁的自然频率：
结果可以和精确解比较：
12.7——一致、集中质量矩阵
12.3节中推出的质量矩阵是一致质量矩阵，因为用于推导刚度 矩阵的位移模型也用于推导质量矩阵。
一些动态问题可以用形式简单的质量矩阵求解。最简单的质量 矩阵——集中质量矩阵，可以通过将质点指定到节点上。
10 1, 1l 0,
20 0 2l 1
（3）
▪ 以上边界条件确定了 1 x 、2 x 由于这两个函数的任意性
我们可以用简单的线性函数来近似，因此有：
1
x

1

x l
,
2 x

x l
代回（1）式中有：
（4）
u(x, t)

(1
点处，有两种位移形式，一个是线位移，即挠度，一种是角位移。
图中， f1t, f3t 是力，
f2t, f4t 是力矩。
1t,3t 是对应的线位移，
2 t,4 t 是对应的转角。
f x,t 是分布载荷
x, t 是梁单元上任意位移 x处的挠度。
图12.2
有限元思想
▪ 1，实际结构被一些元素所取代,这些元素都是被假定为一 个连续的结构部件即有限元，这些元素在特定点即节点上 互相关联。
2，如果解决方案的各方面都选择得当,那么它可以收敛到 精确的解决方案，因为组成总体结构的元素很小，在节点 上的力的平衡和元素之间的位移都令人感到满意,这样整个 结构(组合的元素)表现为单一实体。
如图12.12所示的两端固定的阶梯梁，编写一个MATLAB程 序，命名为Program17.m，对它进行特征值分析。
Program17.m程序的输入数据如下： xl(i)=阶梯梁的长度 i xi(i)=阶梯梁的转动惯量 i a(i)=阶梯梁的横截面积 i b(i，j)=对应本地的自由度 j 和阶梯梁的自由度 i e=杨氏模量 rho=质量密度
用一致质量矩阵 用集中质量矩阵
解得： 事实上，方程的精确解析解是：
12.8 MATLAB应用举例
▪ 例12.5 阶梯轴的有限元分析 ▪ 图12.11中的阶梯轴满足一下条件：A1=16×10-4 m2 ， A2=9×10-4 m2
， A3=4×10-4 m2 ，Ei=20×1010Pa，i=1,2,3，pi=7.8×103Kg/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 阶梯梁的特征值分析的程序
▪ 杆单元的势能可以写成：
式中，E是弹性模量，（10）表示成矩阵形式为：
（10）
这里，
（11） ，所以刚度矩阵[k]可以表示成：
(3)计算等效节点力
（12）
遵循等效原则，即原载 荷和等效之后的节点载 荷在虚位移上所做的虚
功相等。
设单元上x处作用有分布力f ( x , t)，现在要把它等效成节点力 f1t , f2 t
▪
谢谢观赏
▪ 程序给出的阶梯梁的自然频率和模态的输出结果如下：
▪ 本章小结：
▪ 有限元法是一种流行的数值计算方法，准确的找到解决 复杂的实际系统的方案。我们介绍的方法适用于振动问题 。我们介绍简单的结构元素，如条，杆，梁，全局坐标系 变换的矩阵，装配的单元矩阵，和解决方案的有限元方程 推导的刚度和质量矩阵的方法。通过几个静态和动态（振 动）的例子，作为我们提出的方法的实际应用，可以知道 以有限元法为基础的MATLAB方案可以很好地解决振动问 题。
用于振动分析的有限元方法
指导老师：陈益 报告人：成志斌 韩宗彪
何瑜 宁鹏
内容
有限元介绍 单个元素的运动方程 单个元素的质量矩阵、刚度矩阵、力矢量及其转化 整个系统的运动方程 整个系统的边界条件的加载及质量矩阵 MATLAB实例及总结
有限元法简介
▪ 有限元法是一种可用于精确地(但近似)解决许多复杂的 振动问题的数值方法。
u(x, t) 1u1(t) 2u2 (t)
（1）
式中，Φ1、 Φ2称为线性系数，与单元里点的位置有关，是x的函数 。此函数与单元的形状有关，又叫形状函数。
形状函数和插值函数一样是任意的，但必须边界条件：
u(0, t) u1(t) u(l,t) u2 (t)
（2）
只有满足此条件单元才能协调一致运动，而不致破坏系统的完 整性，因此这两个条件实际上就是变形协调条件。将式（1）带 入（2）中，就可以得到形状函数Φ1(x)、Φ2(x)所满足的边界条 件：
▪ 前文中，节点没有固定，结构在节点力的作用下会发生刚体 位移。 ▪ 也就是说，矩阵[K]是奇异矩阵。通常情况下，我们希望结构 的位移为零。 ▪ 因此，我们需要添加边界条件对矩阵[M]、[K]和向量F进行约 束。
例1：杆件分析
如图：均质；长0.5m；断面截面积5e-4m^2；杨氏模量 200GPa；密度7850Kg/m^3；左端固定。
wenku.baidu.com
其实，就是对应于广义坐标 u1t , u2 t 的广义力，为此，计算 f x, t 所做的虚功：
把上式写成矩阵形式：
（13）
所以等效节点力可以写成：
mut kut f t （14）
二 梁单元
如下图所示，一个梁单元也是有两个节点，但是有四个自由度，每个节
▪ *认识用于解决不同类型振动问题的刚度和质量矩阵。 ▪ *将矩阵元素从局部坐标系变换到全球坐标系。 ▪ *装配单元矩阵和应用边界条件。 ▪ *对杆、梁元素进行静态分析。 ▪ *对杆、梁元素进行动态分析来得到固有频率和振型。 ▪ *在有限元振动分析使用一致的集中质量矩阵。 ▪ *使用MATLAB解决振动问题。
在静载弯曲条件下，梁单元上任意点出的挠度是x的三次方程，可写成：
（16）
此方程必须满足下面的边界条件：
挠度的斜率
tan
（17）
由此可以求解处a (t)、b (t)、c (t)、d (t)，进而挠度方程为：
（18）
上式可以写成形状函数的表示： 其中，形函数分别为：
（19）
梁单元的动能、势能、虚功表达式分别为：
三 全系统运动方程
经过坐标变换，各个单元的节点位移方向被统一起来，但是不同的 节点有不同的节点位移，为了便于综合出全系统的运动方程，首先要 建立全系统的节点位移向量。
每个单元的节点位移向量与全系 统的节点位移向量之间的关系：
单元节点 位移向量
长方形矩阵 由1和0组成
例如，图的12.5中的1单元，方程，变为：
式中I是 横截面的
惯性矩
上式中：
（20） （21） （22）
通过上式，可以得到梁单元的质量、刚度矩阵，等效节点力：
二 单元矩阵的坐标变换
局部坐标系：以各个单元本身的轴线为基准所设立的坐标系。便于计算节点位移。 缺点：如果整个系统里各个单元取向各异，各个节点位移方向不一致，如下图。如 何使汇交于一个节点的各个杆件的节点位移真正相等? 解决方法：进行坐标变换
a.节点2处施加1000N静态轴向外力u2，求应力 b.求系统固有频率
解： a：
平衡方程：
A=5e-4，E=2e11，l = 0.5，f2=1000，代入方程得：
u1：位移，f1：节点1处应力，添加边界条件：u1=0，解得： u2=5e-10m
由应力 与应变 的关系：
l u2 u1 表示长度的变化，l 表示应变；
▪ 3，因为得到准确解很难，所以得到一个方便且逼近的近似 解很有价值。
元素的运动方程
龙门刨铣床
三角板元 素
有限元 模型
梁 元 素
元素的运动方程
▪ 位移函数 ▪ 形状函数 ▪ ▪
各点对应位移 动能 应变能
未知节点位移数 n
质量矩阵
刚度矩阵
主要内容：
▪ 一，单元的质量、刚度矩阵、等效节点力 ▪ 二，单元矩阵的坐标变换 ▪ 三，整个系统的运动方程
把每个单元的动能相加，就得到了整个系统 的动能：
（把整个系统的动 能表示成关于节点 速度矢量的形式）
▪ 这样就得到了整个系统的质量矩阵:
类似地，考虑整个系统的势能，便可以得到整个系统的刚度矩阵：
整个系统的广义力向量： 最后得到整个系 统的运动方程：
12.6——添加边界条件
N：结构中自由节点位移的数目