有限元理论与方法讲

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

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

3、动力分析。4.1.1 运动方程结构离散化以后，在运动状态中各节点的动力平衡方程如下(2-2-1)F i + F d + P（t尸 F e式中：FFd、P分别为惯性力、阻尼力和动力荷载，均为向量；Fe为弹性力。弹性力向量可用节点位移8和刚度矩阵K表示如下Fe=K 8式中：刚度矩阵 K的元素Kj为节点j的单位位移在节点i引起的弹性力。软S根据达朗贝尔原理，可利用质量矩阵M和节点加速度 丁表不惯性力如下ft.2-M-:t2式中：质量矩阵的元素 Mj为节点j的单位加速度在节点i引起的惯性力。设结构具有粘滞阻尼，可用阻尼矩阵C和节点速度 "表示阻尼力如下日F -_C /F d 一 C 千式中：

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

5、法与由单元刚度矩阵求整体刚度矩阵时相似。在动力计算中可采用两种质量矩阵，即协调质量矩阵和集中质量矩阵。1 .协调质量矩阵从运动的结构中取出一个微小部分，根据达朗贝尔原理，在它的单位体积上作用的惯性力为p i 一3；t2式中：p为材料的密度。在对结构进行离散化以后，取出一个单元，并采用如下形式的位移函数则n铲§eP i - - N -2-ft再利用荷载移置的一般公式求得作用于单元节点上的惯性力为$*eFiein N T PidV - - N T NdV-2二 t即Fie - m 才可见，单元质量矩阵为m = N T N dV（2-2-4）如此计算单元质量矩阵，单元的动能和位能是互相协调

6、的，因此叫做协调质量矩阵。2 .集中质量矩阵假定单元的质量集中在它的节点上，质量的平移和转动可同样处理。这样得到的质量矩阵是对角线矩阵。单元集中质量矩阵定义如下：(2-2-5)式中，中为函数叼的矩阵，叼在分配Z节点i的区域内取1,在域外取0。由于分配给各节点的区域不能交错，所以由上式计算的质量矩阵是对角线的。3 .平面等应变三角形单元集中质量矩阵与协调质量矩阵设单元重量为 W,将它3等分，分配给每一节点，得到单元集中质量矩阵如下100000010000(2-2-6)W001000m =3g000100000010000001单元协调质量矩阵为Wm =3g20140141401201400140

7、1201414014012014014012(2-2-7)在单元数目相同的条件下，两种质量矩阵给出的计算精度是相差不多的。集中质量矩阵不但本 身易于计算，而且由于它是对角线矩阵，可使动力计算简化很多。对于某些问题，如梁、板、壳等。 由于可省去转动惯性项，运动方程的自由度数量可显著减少。当采用高次单元时，推导集中质量矩 阵是困难的。另外，只要离散化时保持了单元之间的连续性，由协调质量矩阵算得的频率代表结构 真实自振频率的上限。2.2.3阻尼矩阵如前所述，结构的质量矩阵M和刚度矩阵K是由单元质量矩阵m和单元刚度矩阵Me经过集 合而建立起来的。相对来说，阻尼问题比较复杂，结构的阻尼矩阵C不是由单元阻

8、尼矩阵经过集合而得到的，而是根据已有的实测资料，由振动过程中结构整体的能量消耗来决定阻尼矩阵的近似值。1.单自由度体系的阻尼单自由度体系的自由振动方程为m，5 二 c15、k、. = 0式中：m为质量；c为阻尼系数；k为刚度系数；8为变位。 上式两边除以m后得到、. 2 - . ' . J =0其中，s=Jk7m, U=c/(2mo ), t称为阻尼比，3为体系的自振频率(角频率)设初始条件为：当t=0时，8 = 8°, 6 = Vo ,符合这些初始条件的解为(2-2-8)v0 二0 .、=expi -ct % cos dt -sin dt- -d d = . 1 .2体系的

9、自振频率为 3d,其振幅随着时间而逐渐衰减。可见，阻尼对自振频率的影响是很小的，通常可取co d= co oI =0.01 0.10,般都小于0.20。根据实测资料，大多数结构的阻尼比都是很小的数，较多为5 / 82.多自由度体系的阻尼如果假定阻尼力正比于质点运动速度，从运动的结构中取出一微小部分，在它的单位体积上作 用的阻尼力为pd =-aPr = -aPN eft式中：a为比例常数；p为材料密度；N为形函数。利用荷载移置的一般公式求得作用于单元e的节点上的阻尼力如下F deT=N pddV = - ： N N dVS即F de = -C 8e而TC = N N dV =1 m(2-2-9)

10、可见，此时单元阻尼矩阵正比于单元质量矩阵。如果假定阻尼力正比于应变速度，则阻尼应力可表 为od = -Pd = Pdb §eft所以作用于单元e的节点上的阻尼力为一TFde = 1fB 8d dV = -P (BT DBdV Se = -C 8e其中C =BT DBdV 6e = K e(2-2-10)可见，此时单元阻尼矩阵正比于单元刚度矩阵Keo前面已经说过，通常是根据实测资料，由振动过程中结构整体的能量消耗来决定阻尼的近似值， 因此不是计算单元阻尼矩阵，而直接计算结构的整体阻尼矩阵Co 一般采用如下的线性关系，并称为瑞利(Rayleigh)阻尼,即C = ：M K(2-2-11)

11、其中的系数”和3根据实测资料决定。现在说明如何计算 a和3。设小和扪为两个振型。对式(2-2-11)的两边先后乘以 小，再前乘以 . 得到C C G =a(|)T MG +P(|TKG(2-2-12)根振型正交性再由式(2-2-12)得到6： CG=0(i#j)c g =© +%j2 m加(i = j)其中mpj = 6：M G令" Pu 2 = 2 j , j(2-2-13)则由C由=21j0 j mpj由式(2-2-13)得到二；-： ji = - (2-2-14)j 2 , j 2实测两个阻尼比即可求解a与印结构动力学方程主要采用振型叠加法和直接积分法。前者用到振型正

12、交条件，但不同的振型之 间不能解耦时(在结构与地基的相互作用问题中，地基的阻尼往往大于结构本身的阻尼，对于结构和地基应分别给以不同的a与3值)，应采用直接积分法求解。2.2.4结构自振频率与振型在式(2-2-3)中，令P(t)=0,得到自由振动方程。在实际工程中，阻尼对结构自振频率和振型的影响 不大，因此可进一步忽略阻尼力，得到无阻尼自由振动的运动方程K S+M S=0(2-2-15)设结构作下述简谐运动S= Gcosat把上式代人式(2-2-15),可得到齐次方程2(K M )小=0(2-2-16)在自由振动时，结构中各节点的振幅中不全为零，所以结构自振频率方程为 一 2 一 一K -O M

13、 =0(2-2-17)结构的刚度矩阵K和质量矩阵M都是n阶方阵，其中n是节点自由度的数目，所以上式是关于«的n次代数方程，由此可求出结构的自振频率31 W 32W CO3W W cdn对于每个自振频率，由式(2-2-16)可确定一组各节点的振幅值%=41,由2,，* in：,它们互相之间应保持固定的比值，但绝对值可任意变化，它们构成一个向量，称为特征向量，在工程上通 常称为结构的振型。因为在每个振型中，各节点的振幅是相对的，其绝对值可取任意数值。在实际工作中，常用以 下两种方法之一来决定振型的具体数值：(1)规准化振型：取 叼的某一项，例如取第 n项为1,即6 in=1,于是中 i=

14、%1, %2,，1：(2-2-18)这样的振型称为规准化振型。(2)正则化振型：选取* ij的数值，使这样的振型称为正则化振型。设已求得一振型d=%1,礼，Wn T,如令则得到的伞=加的2,的inT为规准化振型。如令c 二 hT MGTM G =1(2-2-19)ji=ij1/2* 1/"in(2-2-20)(2-2-21)则得到的G=%1,©i2，的inT为正则化振型。mpi6T M(2-2-22)当M为集中质量矩阵时，则而1mpi =1；1m2i2s=1当6为正则化振型时，有kpi = GTKGmpi=1,T 2_加以M ()imn2=% mpi(2-2-23)式中，m

15、pi和kpi分别称为第i阶振型相应的广义质量和广义刚度。 由式（2-2-23）得例2-3求解Ka32中的振型，其中2-10K = -14-1 ，“-1 2 一,i k k pi / mpi2-10M 14 1JT2 一(2-2-24)频率方程为2 0.5助2K - 0 M =-1010一 2._40 1 0_ _ 212 0.5 8求得三个自振频率为22.2。1 =2,缶2=4,83=6将012 =2代入式(2-2-16)中，得到第1振型必须满足的方程组如下11 12 + 0=0 , 41+2力2弧=0,可12 + ©13=0联立前两个方程解出如=4,13,12= 413取丸=1 ,

16、得到规准化的第一振型为1=111T用同样方法得到第 2、3振型为2=-101T3=1-11T由式(2-2-21)得到正则化振型如下彳=1/ 21/ .21/ 2T2=-101T3=1/ 2-1/ 21/ 2T2.2.5振型叠加法求解结构的受迫振动目前，常用的求解结构受迫振动的方法有两种，即振型叠加法和直接积分法。用振型9的线性叠加来表示处于运动状态中的结构位移向量 nB=旷1。)十也"2(t )+十 6n”n(t)=£ 拈i(t)(2-2-25)i =1用M前乘上式的两边，由于振型正交性，等式右边的n项中只剩下i=j这一项，即6M S=,(t N：M 脑=mjj(t)由此得

17、到(2-2-26)(2-2-27)(2-2-28)一二2mpi)和可的初始值可表示为6 M S(0)i 0 =mpi“，、GT M 9(0) i 0 =mpi现在考虑下列运动方程的求解：M B+C B+K B=P (t)把式(2-2-25)代入上式，得到 nnnM Z /； +c£/+k£ a、= P(t)i 1i 衽i =1对上式两边前乘以端,并令C =a M + 3 K,得到nnn一 T. . TTTZ (J)j M 6Q +£ (J)j (o(M + pK 州i +£ ()j K Ni =(J)j P(t) i ii 1i i由于振型正交性，得到2

18、2Tmpi +(cc + %i 加村 +mi mpi =加 P(t)由于ot +P2 =2，iO ,上式进一步化为21 T q +2，©q= 心 P(t)(i =1,2,3,，n )(2-2-29)mpi这是二阶常微分方程，这样的方程共有n个，它们是互相独立的。式(2-2-29)在形式上与单自由度体系的运动方程相同。其解答可用数值积分方法计算，也可用Duhamel积分计算如下：i t二一1P * . e _ i i't_ sin di t i d .(2-2-30)+ e-5 4Tli(0 JcosSdit；0 i -0sin.it, 'di其中P (t )=叱 P (t)di m pi把