振动理论及工程应用8-第八章-振动分析的有限元法课件

1、第8章 振动分析的有限元法 有限元法是力学模型系统上近似的数值计算方法。它先将要分析的工程结构模型，假想地分割成有限个单元，组成离散化模型。 各个单元之间在单元的外节点处互相连接起来。然后导出各单元体的运动方程，然后由各个单元的运动方程组合而形成原工程结构的有限元运动方程。 有限元法中分析的结构，是一个由有限个单元组成的与原结构非常接近的离散系统。计算所得结果的精确程度取决于单元体的划分。 有限元计算的基本过程： （1）将结构离散化，即把结构划分成离散的单元。 （2）考虑单元的性质，建立单元的质量矩阵、刚度矩阵、阻尼矩阵、载荷矩阵，推导出单元体的运动方程式。 （3）组合各单元的质量矩阵、刚度矩

2、阵、阻尼矩阵，得到整个离散系统的运动方程。 （4）解特征方程，求出频率与振型。 （5）求解动力响应问题。 有限元法中采用的单元类型很多，其形状、大小可以变化，各单元互相之间也容易连接，因此它能适应复杂的结构，也适用于各种不同的边界条件。 有限元法中的分析顺序是比较固定的，因此便于计算机计算，并有标准程序和通用程序。 取一单元体，设单元体的动能为T，应变能为U，阻尼消耗的能量为Wd，外力的势能为We。建立拉格朗日函数为 8.1 单元体的运动方程式 设q为单元体中任一点的位移矢量，qe为单元体上各节点的位移矢量，它是时间t的函数。令单元体中任一点的位移矢量q用单元体上各节点的位移矢量qe表示为 式

3、中N为形函数矩阵，它是坐标x、y、z的函数。 式中u(t)、v(t)、w(t)分别表示该点沿x、y、z方向的位移，它们都是时间t的函数。 单元体中任一点的位移矢量q又可表示为式中r为单元体积的质量单元体的动能为单元体中任一点的速度矢量可表示为 D为应力应变关系矩阵，又称为弹性矩阵，所以单元体的应变能为由弹性力学公式，应变与位移的关系为 式中B为应变位移关系矩阵，是几何矩阵，与t无关。单元体上的应力为 设单元体振动时，阻尼系数为c，则阻尼力为 。单元体上阻尼力所消耗的能量Wd为 单元体上所受的外力分为两部分，即体积力FV和表面力FS。它们的势能分别为拉格朗日函数由哈密尔顿原理，将其在时间区间(t

4、1,t2)上对L积分，并使其变分等于零，考虑到D的对称性后，有 式中qe(t1)=0，qe(t2)=0，则只剩下第二项。应用分步积分公式，上式的第二项有用同样的方法可得到令于是变分式成为 由于单元位移的变分 是任取的，所以可由式得到单元的运动方程为 由此得到 Keq、Meq、Ceq、Feq分别表示单元体的刚度矩阵、质量矩阵、阻尼矩阵和载荷矩阵 8.2 单元体的特性分析 在有限元中，形函数的作用十分重要，因为单元形状和相应的形函数确定以后，其它运算可依照标准步骤和普遍公式进行。单元上任一点的位移用节点的位移表示为8.2.1 形函数矩阵N用u、v、w表示一点在空间沿x、y、z方向的位移 式中Ni为

5、形函数，ui、vi、wi是第i个节点的位移。形函数Ni是单元内部坐标的连续函数，它所满足的条件是： (1) 在节点i处， Ni =1; (2) 在其它节点处， Ni =0; (3) 满足 Ni =1 。 用它定义的未知量u、v、w保证了相邻单元间的连续性。 为了保证收敛于精确解，形函数应包含任意线性项和使单元包含有常应变状态和刚体位移。 单元的形状越复杂，形函数的阶次就越高，单元适应能力就越强。 形函数矩阵为写成矩阵形式于是得到形函数的边界条件 在一维单元中，有二个节点，节点的位移为u1(t) u2(t) 。设单元上距1点距离为x点的位移为 式中N1(x)、N2(x)为形函数也称为插值函数，它

6、们应使点的位移满足单元的边界条件，即 满足 Ni =1 。所以对于一维单元，形函数矩阵为 设形函数Ni（x）=ai+bix，包含常数项和线性项。将边界条件代入可得形函数是用局部坐标在单元中定义的。 对于二维单元（，），正方形单元有4个节点，经推导可得其形函数为 引入新的变量0=i，0=i。其中i、i为节点i的坐标，于是上面的四个形函数可合并表示为（i = 1，2，3，4） 对于三维单元（，），正六面体，将坐标原点取在单元形心上，单元边界是六个平面 =1， =1， =1，单元有8个节点，经推导可得其形函数为（i = 1，2，3，4，5，6，7，8） 对于其它单元的形函数可参阅有关书籍。8.2.2

7、 应变位移关系矩阵B 将位移函数矩阵表达式，代入由弹性力学的应变与位移的关系其中子矩阵为可得到 决定了应变与形函数之间的关系对于一维单元，可求得B为即其中E为弹性模量，G为剪切弹性模量，为泊松比8.2.3 弹性矩阵D材料力学中的广义虎克定律剪应变与剪应力的关系为写成矩阵形式 对于平面应力问题，z=xz=yz=0，可得到应力与应变的弹性方程平面应力问题的弹性矩阵为平面应变问题的弹性矩阵为可得到单元体的运动方程式 将以上推导的形函数矩阵N，应变与位移的关系矩阵B，弹性矩阵D代入以下有关方程 计算的质量矩阵称为一致质量矩阵，它总是正定的。如果选择的位移函数接近真实位移，那么计算的结果比较正确，频率与

8、振型比较可靠，接近频率的上界。但是它是一个满矩阵。 8.2.4 质量矩阵质量矩阵的计算公式为 另一种质量矩阵。即将单元体的质量简单地分配于单元的节点上，每个节点上分配到质量的多少，要根据该节点所管辖的范围而定，这样得到的质量矩阵称为集中质量矩阵，它是一个对角矩阵。一般应用集中质量矩阵得到的频率比较偏低。 其中m为三角形单元体的质量，计算的频率接近实际频率的上界。 以平面问题中的三角形单元，其一致质量矩阵为 如果用集中质量矩阵，则只要将单元的质量一分为三，集中作用在三个节点上，即 对角矩阵。得到的频率比实际频率偏低。 在有限元法中，所用的坐标系是局部坐标系，由于单元在空间的局部坐标系不同。但描述

9、整个结构的运动，需要选择一个统一的坐标系，称为总体坐标系。在分析计算时，必须将单元特征的各个方程通过方向余弦矩阵转换到总体坐标系中。 设 为总体坐标系，x，y，z为局部坐标系，其方向余弦矩阵为8.3 坐标转换 式 中表示轴x与 轴之间夹角的余弦等等。则有将以上两位移分量式合并，可写成 同样位移矢量间也适用此交换关系。对于平面二维系统单元 除了某些单元在空间具有相同的方位，即其局部坐标系是平行的情况外，不同单元的矩阵L是不同的。由于l表示两个正交轴系间的一种变换，所以矩阵L是正交的，即 。 特例：在平面结构的特殊情况下，所有各单元的局部坐标系中，都有一根坐标z轴与总体坐标系的一根z轴平行。单元体

10、上任意一点M的局部坐标在总体坐标系中的关系为写为矩阵形式 称为方向余弦矩阵。于是可得到在总体坐标系中的该单元的运动方程式 有了变换矩阵之后，可以将单元的质量矩阵、刚度矩阵等换算至总体坐标系中，即由单元载荷Feq按其相应的贡献叠加，可的节点载荷矢量F。于是整个结构的运动方程式为 并将用总体坐标系的单元质量矩阵、阻尼矩阵、刚度矩阵，按其相应的贡献叠加到总质量矩阵、总阻尼矩阵、总刚度矩阵中。 为书写方便，将总体坐标系的运动方程式中的矩阵上的横线省略。 总质量矩阵、总阻尼矩阵、总刚度矩阵为例8.1 试用有限元法列出图示的简支梁的弯曲自由振动方程式。已知梁长为L，抗弯截面模量为EI，单位长度的质量为m

11、。 解: 将梁离散分为二个单元，单元长l=L/2，现取其中一个单元研究，并确定形函数N及应变与位移的关系矩阵 B。单元节点位移矢量 由于单元有四个端点位移，所以有四个待定常数a1、a2、a3、a4。写成矩阵形式为现假设单元中的点的位移函数为多项式由边界条件时得时即写为矩阵形式由此解出四个待定常数为得到即其中于是有应变与位移的关系矩阵B梁的应变应变与应力的关系为=E，即D = E，可求得单元刚度矩阵Keq。注意到单元刚度矩阵为单元质量矩阵为积分得求结构的总刚度矩阵和总质量矩阵。单元编码 单元编码如图示，取统一的总体坐标系，写出结构的位移列矩阵和单元的位移列矩阵。结构位移列矩阵结构位移列矩阵单元节

12、点位移转换矩阵同理，单元节点位移转换矩阵其中0和I都是22阶子矩阵两单元的刚度矩阵分别为总刚度矩阵其中代入上式得同样质量矩阵也可按类似方法叠加得到。简支梁整个结构系统的刚度矩阵和质量矩阵。 因简支梁在x = 0和x = l处，点的位移为零，即v1=v3=0，因此,可以删去第一行、第五行及第一列、第五列 于是得简支梁的振动方程式为将l=L/2代入，得到矩阵形式8.4 固有频率及主振型 用有限元法求得结构的自由振动方程式为设结构做简谐振动，其解为或 式中q0是位移qSe的振幅矢量，p是系统的固有频率，将其代入振动方程，消去sinpt得到 用有限元法计算时，质量矩阵总是对称的正定矩阵，而刚度矩阵则因

13、为在考虑在考虑单元体的位移模式时，考虑了刚体位移和常应变状态，即各单元体都被看成是不受约束的，所以刚度矩阵是半正定的，系统含有刚体运动模态。为了以后计算方便，可将系统分为约束系统和无约束系统，对于无约束系统要从其中消去刚体运动模态。 对于无约束系统，由于刚度矩阵是奇异的，所以可以给出一些位移矢量，并满足 它们对应于零频率，一个完全自由的系统，具有六个刚体运动模态，所以有六个零频率，对应这些零频率的特征矢量，相互之间应满足对M的正交关系，则 利用正交关系，可以找到各位移之间的变换关系。经过坐标变换，可以从运动方程式中消去刚体运动模态。并使得总刚度矩阵变为对称正定矩阵，可以求解广义特征值和特征矢量

14、。 由于总刚度矩阵是大型稀疏矩阵，常采用子空间迭代法，行列式搜索法、雅可比方法、QL方法和迭代法等进行计算。例8-2 求解例8-1所示简支梁的自由振动频率。解：由例8-1，已解出简支梁的振动方程式为 由上述计算可见，单元划分的太粗影响计算结果的精确度。 本题的精确解为经计算可得例8-3 某人行中承式钢管混凝土拱桥全长126m,主桥跨中心线间距110m，桥面全宽7m，全桥钢结构净重约248t。桥高8m。本桥按超静定结构设计，承载主体为钢管混凝土拱肋，拱肋两端锚固于拱座基础内，全桥采用两根平行、对称的主拱肋，主拱肋采用等截面单圆钢管结构。平面框架由横梁、纵梁及桥面板组成，其中横梁为主要承载结构，纵

15、梁为次要承载结构。荷载直接作用于桥面，经横梁传递至吊杆，吊杆两端分别锚于主拱肋和横梁上，垂直与主拱肋及桥面，所有吊杆均采用不锈钢拉杆。桥台采用实体式钢筋混凝土桥台,基础采用钢筋混凝土扩大基础。求该桥的前阶固有频率和振型。 解：根据桥梁的特点，采用多种单元形成混合力学模型，其中、横梁、纵梁、横撑、立柱、主拱肋、立柱横联钢管采用空间梁单元，吊杆采用三维杆单元，桥面采用空间壳单元。共计6920个节点，8763个单元，有限元模型如图8-8(a)所示。计算结果如图8-8(b)(f)所示。有限元模型 桥面竖向反对称振型(固有频率0.707347Hz) 桥面竖向对称振型(固有频率1.295 Hz) 两拱及桥

16、面微扭振型(固有频率1.345 Hz) 桥面横向对称振型(固有频率1.741 Hz) 桥面竖向反对称振型(固有频率2.372 Hz) 通过建立典型的中承式钢管混凝土拱桥空间有限元模型，利用软件计算了其前阶固有频率和相应的振型，从而对其动力学性能加深了了解，对中承式钢管混凝土拱桥的设计和施工提供了一定的帮助。 8.5 系统的响应 用有限元法解系统的响应问题就是解动力方程组 目前普遍使用的有两种方法，一种是振型叠加法(最低几阶振型)，另一种是逐步积分法(高频振型)。 在振型叠加法的计算中，假定结构的响应能用前s个较低的振型（sn来描述,这时可以将所求的个主振型依此排列，构成一个ns阶的截断振型矩阵。 当然也可以用正则振型构成截断振型矩阵。 一、振型叠加法(最低几阶振型)， 用截断振型矩阵进行坐标变换