有限元理论与方法讲

1、第 13 讲(第 13 周)4.1 结构动力学问题有限元方法动力学问题在国民经济和科学技术的发展中有着广泛的应用领域。最经常遇到的是结构动力学 问题，它有两类研究对象：一类是在运动状态下工作的机械或结构，例如高速旋转的电机、汽轮机、 离心压缩机，往复运动的内燃机、冲压机床，以及高速运行的车辆、飞行器等，它们承受着本身惯 性及与周围介质或结构相互作用的动力载荷。如何保证它们运行的平稳性及结构的安全性，是极为 重要的研究课题。另一类是承受动力载荷作用的工程结构，例如建于地面的高层建筑和厂房，石化 厂的反应塔和管道，核电站的安全壳和热交换器，近海工程的海洋石油平台等，它们可能承受强风、 水流、地震以

2、及波浪等各种动力载荷的作用。这些结构的破裂、倾覆和垮塌等破坏事故的发生，将 给人民的生命财产造成巨大的损失。正确分析和设计这类结构，在理论和实际上也都是具有意义的 课题。动力学研究的另一重要领域是波在介质中的传播问题。它是研究短暂作用于介质边界或内部的 载荷所引起的位移和速度的变化，如何在介质中向周围传播，以及在界面上如何反射、折射等的规 律。它的研究在结构的抗震设计、人工地震勘探、无损探伤等领域都有广泛的应用背景，因此也是 近 20 多年一直受到工程和科技界密切关注的课题。现在应用有限单元法和高速电子计算机，已经可以比较正确地进行各种复杂结构的动力计算， 本章阐明如何应用有限单元法进行动力分

3、析。4.1.1 运动方程结构离散化以后，在运动状态中各节点的动力平衡方程如下(2-2-1)F i +Fd +P(t)=F e式中： Fi、Fd、P(t)分别为惯性力、阻尼力和动力荷载，均为向量；F e为弹性力。弹性力向量可用节点位移 和刚度矩阵 K 表示如下F e =K 式中：刚度矩阵 K 的元素 Kij 为节点 j 的单位位移在节点 i 引起的弹性力。根据达朗贝尔原理，可利用质量矩阵M 和节点加速度2 t2表示惯性力如下FiM t 2式中：质量矩阵的元素 M ij为节点 j 的单位加速度在节点 i引起的惯性力。设结构具有粘滞阻尼，可用阻尼矩阵 C 和节点速度 表示阻尼力如下tFC FdC t

4、 2式中：阻尼矩阵的元素 Cij为节点 j 的单位速度在节点 i 引起的阻尼力。 将各力代入式 (2-2-1) ，得到运动方程如下M 2 C K P(t) t 2 t(2-2-2)则运动方程可写成 C K P(t)(2-2-3)在地震时，设地面加速度为 a，结构相对于地面的加速度为 ，结构各节点的实际加速度等于 a+ ， 在计算惯性力时须用它代替式 (2-2-3)中的 。至于弹性力和阻尼力，则分别取决于结构的应变和应 变速率，即取决于位移 和速度 ，与地面加速度无关。2.2.2 质量矩阵下面用 m 表示单元质量矩阵， M 表示整体质量矩阵。求出单元质量矩阵后，进行适当的组合即 可得到整体质量矩

5、阵。组合方法与由单元刚度矩阵求整体刚度矩阵时相似。在动力计算中可采用两种质量矩阵，即协调质量矩阵和集中质量矩阵。1. 协调质量矩阵从运动的结构中取出一个微小部分，根据达朗贝尔原理，在它的单位体积上作用的惯性力为式中： 为材料的密度。在对结构进行离散化以后，取出一个单元，并采用如下形式的位移函数pi N2t2ei t 2再利用荷载移置的一般公式求得作用于单元节点上的惯性力为FieNT pidVNT NdV 22 tF im可见，单元质量矩阵为(2-2-4)m N T N dV如此计算单元质量矩阵，单元的动能和位能是互相协调的，因此叫做协调质量矩阵。2. 集中质量矩阵 假定单元的质量集中在它的节点

6、上，质量的平移和转动可同样处理。这样得到的质量矩阵是对 角线矩阵。单元集中质量矩阵定义如下：式中， 为函数 i 的矩阵， i 在分配给节点 i 的区域内取 l，在域外取 0。由于分配给各节点的区域不能交错，所以由上式计算的质量矩阵是对角线的。3. 平面等应变三角形单元集中质量矩阵与协调质量矩阵设单元重量为 W，将它 3 等分，分配给每一节点，得到单元集中质量矩阵如下1000001000W00100m3g000100000100000000001(2-2-6)单元协调质量矩阵为W m3g10020140140140(2-2-7)在单元数目相同的条件下，两种质量矩阵给出的计算精度是相差不多的。集中

7、质量矩阵不但本 身易于计算，而且由于它是对角线矩阵，可使动力计算简化很多。对于某些问题，如梁、板、壳等。 由于可省去转动惯性项，运动方程的自由度数量可显著减少。当采用高次单元时，推导集中质量矩 阵是困难的。另外，只要离散化时保持了单元之间的连续性，由协调质量矩阵算得的频率代表结构 真实自振频率的上限。2.2.3 阻尼矩阵如前所述，结构的质量矩阵 M和刚度矩阵 K是由单元质量矩阵 m和单元刚度矩阵 Me 经过集 合而建立起来的。相对来说，阻尼问题比较复杂，结构的阻尼矩阵C不是由单元阻尼矩阵经过集合而得到的， 而是根据已有的实测资料， 由振动过程中结构整体的能量消耗来决定阻尼矩阵的近似值。1.单自

8、由度体系的阻尼单自由度体系的自由振动方程为m c k 0式中： m 为质量； c 为阻尼系数； k 为刚度系数； 为变位。上式两边除以 m 后得到2 2 0其中，k /m ，c/ 2m ， 称为阻尼比， 为体系的自振频率（角频率） 。设初始条件为：当t0 时，0， v 0 ，符合这些初始条件的解为exp t 0 cos dt v00 sin dtdd1体系的自振频率为 d，其振幅随着时间而逐渐衰减。 =0.01 0.10 ，一般都小于 0.20。根据实测资料，大多数结构的阻尼比都是很小的数，较多为可见，阻尼对自振频率的影响是很小的，通常可取d 。2.多自由度体系的阻尼 如果假定阻尼力正比于质点

9、运动速度，从运动的结构中取出一微小部分，在它的单位体积上作 用的阻尼力为pdrNet式中： 为比例常数； 为材料密度； N 为形函数。利用荷载移置的一般公式求得作用于单元e 的节点上的阻尼力如下FdeTNT pddVN NdVF deC eTC N N dVm(2-2-9)可见，此时单元阻尼矩阵正比于单元质量矩阵。如果假定阻尼力正比于应变速度，则阻尼应力可表 为D DB et所以作用于单元 e 的节点上的阻尼力为其中eTFdeB d dVBT DBdVeC eBT DBdVeKe(2-2-10)可见，此时单元阻尼矩阵正比于单元刚度矩阵前面已经说过， 通常是根据实测资料，因此不是计算单元阻尼矩阵

10、，而直接计算结构的整体阻尼矩阵为瑞利 (Rayleigh) 阻尼，即K e。 由振动过程中结构整体的能量消耗来决定阻尼的近似值， 般采用如下的线性关系，并称C。CMK(2-2-11)其中的系数 和 根据实测资料决定。现在说明如何计算 和。设i和 j为两个振型。对式(2-2-11)的两边先后乘以 i ，再前乘以 Tj得到j Cij M ij Ki(2-2-12)根振型正交性再由式 (2-2-12) 得到j C i 0 j C ij mpjijij其中mpjj M j令2j2 j j(2-2-13)则Tj Cj2 j j mpj由式( 2-2-13 )得到j jj(2-2-14)j 2 j2实测两

11、个阻尼比即可求解与 。结构动力学方程主要采用振型叠加法和直接积分法。前者用到振型正交条件，但不同的振型之 间不能解耦时(在结构与地基的相互作用问题中，地基的阻尼往 往大于结构本身的阻尼，对于结 构和地基应分别给以不同的 与 值)，应采用直接积分法求解。2.2.4 结构自振频率与振型在式 (2-2-3) 中，令 P(t)0，得到自由振动方程。在实际工程中，阻尼对结构自振频率和振型的影响 不大，因此可进一步忽略阻尼力，得到无阻尼自由振动的运动方程K M 0 (2-2-15) 设结构作下述简谐运动 cos t把上式代人式 (2-2-15) ，可得到齐次方程2(K2 M ) 0 (2-2-16)在自由

12、振动时，结构中各节点的振幅 不全为零，所以结构自振频率方程为2K2 M 0 (2-2-17)结构的刚度矩阵 K和质量矩阵 M都是 n 阶方阵，其中 n 是节点自由度的数目，所以上式是关于 的 n 次代数方程，由此可求出结构的自振频率1 2 3 n对于每个自振频率，由式 (2-2-16) 可确定一组各节点的振幅值i i1， i2， inT，它们互相之间应保持固定的比值，但绝对值可任意变化，它们构成一个向量，称为特征向量，在工程上通 常称为结构的振型。因为在每个振型中，各节点的振幅是相对的，其绝对值可取任意数值。在实际工作中，常用以 下两种方法之一来决定振型的具体数值：(1)规准化振型：取 i 的

13、某一项，例如取第 n 项为 1，即 in 1，于是i i1， i2， 1T(2-2-18)这样的振型称为规准化振型。(2)正则化振型：选取 ij 的数值，使这样的振型称为正则化振型。设已求得一振型 ii1, i2 , , in ，如令jiij/ in(2-2-20)则得到的 ii1, i 2, , inT 为规准化振型。如令ji ij/c(2-2-21)ciT M1/ 2i则得到的 iT 为正则化振型。mpi i M i(2-2-22)M 为集中质量矩阵时，则m1i1mpi i1i2ininm2i22ms is2s1mnini 为正则化振型时，有mpi=1T T 2k pii Kii i M

14、ii mpi(2-2-23)式中， mpi 和 kpi 分别称为第 i 阶振型相应的广义质量和广义刚度。由式（ 2-2-23 ）得/mpi(2-2-24)例 2-3 求解 K =2M 的振型，其中求解说明频率方程为求得三个自振频率为22 0.5 2114222 0.5122 ， 224 ， 32(2-2-25)(2-2-26)(2-2-27)(2-2-28)将 12 2代入式（ 2-2-16）中，得到第 1 振型必须满足的方程组如下11- 12+0=0， - 11+2 12- 13=0， 11- 12+ 13=0 联立前两个方程解出11= 13， 12= 13 取 13=1 ，得到规准化的第一

15、振型为1=1 1 1T 用同样方法得到第 2、 3振型为2=-1 0 1T3=1 -1 1T 由式（ 2-2-21 ）得到正则化振型如下1=1/ 2 1/ 2 1/ 2 T2=-1 0 1T3=1/ 2 -1/ 2 1/ 2 T2.2.5 振型叠加法求解结构的受迫振动目前，常用的求解结构受迫振动的方法有两种，即振型叠加法和直接积分法。用振型 i 的线性叠加来表示处于运动状态中的结构位移向量n 1 1 t 2 2 tn n ti i ti1用Tj M 前乘上式的两边，由于振型正交性，等式右边的n项中只剩下 ij 这一项，Tj M j t Tj M j mpj j t由此得到i tiT Mmpii

16、和 i 的初始值可表示为iT M (0) i 0 i mpiiT M (0) i 0 i mpi现在考虑下列运动方程的求解：M C K P (t)把式( 2-2-25 )代入上式，得到nnnM i i C i i K i i P(t)i 1i1i1对上式两边前乘以 Tj ，并令 C= M+K，得到n n nj M i ij M K i ij Ki i j P(t)i 1 i 1 i 1由于振型正交性，得到2 2 Tmpi ii mpi ii mpi i j P(t)由于i2 2 i i ，上式进一步化为2 1 Ti 2 i i ii2 iTj P(t)mpii 1,2,3, ,n (2-2-29)这是二阶常微分方程，这样的方程共有 n 个，它们是互相独立的。式 (2-2-29) 在形式上与单自由 度体系的运动方程相同。其解答可用数值积分方法计算，也可用Duhamel 积分计算如下：i t 1P *e i i t sin di t d(2-2-30)其中e i it i 0 cos diti 0 i i i 0 sin ditP* tiT P(t)dim pi把 i(t)代人式 (2-2-25)，即得到所需解答