第六章 动力问题的有限元法

第六章动力问题的有限元法6.1概述前面几章所研究的问题都属于静力问题，其特点是施加到结构上的外载荷不会使结构产生加速度，且外载荷的大小和方向不随时间变化，因而结构所产生的位移和应力也不随时间变化。本章将要研究结构分析中另一类重要问题的有限元解法，即动力问题的有限元解法。动力学问题的特点是，载荷是随时间变化的，因而结构所产生的位移和应力是时间的函数，结构会产生速度和加速度。由于结构本身的弹性和惯性，结构在动力载荷的作用下，往往呈现出振动的运动形态。结构振动是工程中一个很普遍很重要的问题。有些振动对我们有利，例如，振动打桩，振动选料，有些振动对我们有害，例如，机床的振动，仪器与仪表的振动，桥梁、水坝及高层建筑在地震作用下的振动等。因此，我们必须对振动体本身的振动特性以及它对外部激振力的响应有一个明确的认识，才能更好地利用它有利的一面，而避免它有害的一面，设计出更好的机械和结构。振动问题主要解决两方面的问题。1.寻求结构的固有频率和主振型，从而了解结构的固有振动特性，以便更好地利用或减少振动。2.分析结构的动力响应特性，以计算结构振动时动应力和动位移的大小及其变化规律。6.2结构的振动方程结构的振动方程可用多种方法建立，这里我们使用达朗伯原理(动静法)，仿照前几章建立静力有限元方程的方法，来建立动力问题的有限元方程。在静力问题中用有限元法建立的平衡方程是K6二F在振动问题中，对结构的各节点应用达郎伯原理所建立的振动方程仍然具有与上式相同的形式，只不过节点位移是动位移，节点载荷是动载荷，它们都是时间的函数。上面的方程成为K6(t)二Q(t)(6.1)上式中暫(t)为节点的动位移，它是时间的函数，K5(t)是t时刻的节点位移产生的弹性恢复力，它与该时刻的节点外力(t)构成动态平衡。在动态情况下，结构承受的载荷(集中载荷，分布载荷)可随时间而变化，是时间的函数。按有限元方法将此种载荷移置到节点上，得到的节点载荷向量F(t)也是时间的函数。此外，结构在运动中，各点除位移f以外，还有速度及加速度。按照达郎佰原理，有加速度的质量应附加有惯性力载荷。如材料的密度为P,则结构单位体积的惯性力为-P(fj。这对结构来说，相当于又受有另一种体积力，大小与点的加速度成比例，而方向与加速度方向相反。另外，在结构运动过程中，还会受到周围介质和来自内部的阻力。精确地描述这种阻力的变化规律是很困难的，一般采用阻力与速度L比例的近似线性假定，如阻力系数为“，则单位体积的阻力为-卩f。这对结构来说相当于另一种体积力，大小与点的速度成比例，方向与速度方向相反。按有限元方法，用单元节点位移e进行插值表示单元内部位移。(6.2)f=N5e此处形函数仍只是位置的插值函数，与时间无关，则单元内的速度和加速度分别为(6.3)门=N5eP=N以及（6.4）其中、幺为单元节点的速度及加速度向量。将单元惯性力-Pfk阻力-af乍为体积力，按式（3.28）移置到单元各节点，就得到相应的单元等效节点载荷向量，记为2和则有P卩b=-NTpfdvp将式（6.3）、（6.4）代入上式，有pTpN=-JNtpNdv-meNt山Ni)dv=-JNt山Ndv备=-ce2式中me=NTpNdv(6.5)称为单元质量矩阵。由于推导式（6.5）时采用了与推导单元刚度矩阵时相一致的形函数，故式（6.5）所表示的质量矩阵也称为一致质量矩阵。而(6.6)Ce=JNt卩Ndv称为单元阻尼矩阵，由式（6.5）和（6.6）可见单元质量矩阵和单元阻尼矩阵是对称的。将移置到节点上的动载荷、惯性力、阻力乍为载荷，按单元叠加，得到有限元节点位移方程：K8=F瓦（me&e+ce5e）e=1或MS+C5+K5=F（6.7）其中M=Xme（6.8）e=1称为结构质量矩阵或总质量矩阵，而C=区ce（6.9）称结构阻尼矩阵。可见结构质量矩阵和结构阻尼矩阵分别为单元质量矩阵和单元阻尼矩阵的叠加，其叠加方法与结构刚度矩阵的形成完全一样，借助于单元定位向量,用单元集成法完成叠加过程。由于me和ce是对称的，因而叠加合成的Im和IC也是对称的。式（6.7）是节点位移的二阶微分方程，称为结构的动力方程式。对于不同的结构，可以选用不同的单元，有不同的形函数矩阵In，但动力方程（6.7）的建立过程都是一样的。当结构不受外载荷时，仏=0，如果再忽略阻尼，则动力方程（6.7）式成为M8+K8=0（6.10）这是系统的自由振动方程。弹性结构的振动实际上是连续体的振动，位移f是连续的，具有无限多个自由度。经有限元离散化处理后，单元内的位移按假定的位移形式来变化，可用节点位移插值表示。这样，连续系统的运动就离散化为有限个自由度系统的运动了。如果全部节点有N个自由度，则式（6.10）就是N阶的自由振动微分方程了。6.3一致质量矩阵与集中质量矩阵按照式（6.5）,当单元的位移插值形函数矩阵确定后就可用此式算出它的一致质量矩阵。例如，对平面三节点三角形单元，按照2.3节的叙述，形函数矩阵bv可用面积坐标表示为N二ILILILijm其中I为二阶单位阵。将上式代入式（6.5）可得平面三节点三角形单元的一致质量矩阵me=NtpNdxdytIL=pHILiILjILdxdytmIL=pJJILLILLILLiiijimILLILLILLjijjjmILLILLILLmmmmmidxdyt利用积分公式（2.56），可由上式求得me=-3012014011401201400140120114014012001401401(6.11)用式（6.5）可以计算出其它类型单元的质量矩阵。例如平面桁架单元的质量矩阵为1-61-21-31-6(6.12)而平面刚架单元的质量矩阵为1400156M022L4L2me420700005413L0-13L-3L2对称14001560-22L4L(6.13)式(6.11)、(6.12)及(6.13)中的M为单元质量。由单元一致质量矩阵叠加形成的结构质量矩阵一般都是稀疏、带状的，但都有相当的半带宽，如果将单元质量矩阵近似作为对角型矩阵，单元叠加后的结构质量矩阵也是对角阵。这种对角型的质量矩阵称为集中质量矩阵，而对角型的集中质量矩阵对结构的动力分析是非常有利的。可以将单元的质量以某种方式分配在单元的节点上而得到单元集中质量矩阵。例如当质量均匀分布时，将质量平均分配给各节点。按照这种方法平面三节点三角形单元的集中质量矩阵为me=-3其中H为6阶单位阵。M为单元质量。而四节点四面体单元的集中质量矩阵为me=刃(6.14)(6.15)其中I为12阶的单位阵。对于矩形弯曲板单元，在不计转动质量影响的条件下，它的集中质量矩阵为0100me=-4010001(6.16)00对于平面刚架单元，在不计转动质量影响的条件下，它的集中质量矩阵为me=-26.4阻尼矩阵(6.17)各种工程结构的阻尼力及其产生的机理是非常复杂的。从宏观上看构周围粘性介质产生的阻尼，称为粘性阻尼。粘性阻尼的阻尼力一般近似与运动速度成正比。另一种是结构材料内部磨擦产生的阻尼，称为结构阻尼或材料阻尼。结构阻尼的阻尼应力一般近似地认为与弹性体的应变速率成正比。如果假定阻尼力与运动速度成正比，那么在运动的弹性体中任意点处单位体积上作用的阻尼力为(P=apff=apN*ddt阻尼有两种主要形态。一种是结式中a比例常数；p材料密度；Lv单元形函数矩阵;单元节点速度矢量。dt可以将阻尼力Pd看成是一种体积力其等效的单元节点阻尼力向量为Fd=BJNtPdv二aJJJNtpNdv6dFe=c5d因此，单元阻尼矩阵或写成Ce二aJJJNTpNdv=ame(6.18)它正比于单元质量矩阵me。如果假定阻尼应力与弹性体的应变速率成正比，则阻尼应力可表示为。D?8二卩DB6eddt式中，卩为比例系数，d为弹性矩阵，b为应变矩阵，0为单元节点速度向量。下面推导阻尼应力的单元等效节点阻尼力向量。设单元等效节点阻尼力向量仍用仗表示，设节点发生虚位移6*，单元内各点产生的虚应变为d*，则*=B6*。Fe在虚位移6d上做的虚功为W二6*TFeed单元的虚应变能为=*tdv=6*tBJBtdvdd由W二U得到eeFd=B!Btdv=bB!BTDBdv6ed或写成Fd因此，单元阻尼矩阵=ce6ece二训BtDBdv二Bke(6.19)它正比于单元刚度矩阵ke在实践中，要精确地计算各种单元的阻尼矩阵是很困难的。通常在程序设计中，假定结构总体阻尼矩阵c是结构总体刚度矩阵k与总体质量矩阵Im的线性组合。称为瑞利阻尼，其表达式为(6.20)C=aM+卩K其中，比例系数及P可通过实验确定。采用瑞利阻尼近似，可以使运动方程求解大大简化，并且在程序中不必单独存贮总体阻尼矩阵。在实际工程问题中，阻尼的作用对结构的动力响应的影响并不大，这种近似处理具有实用价值。6.5结构的自振特性和特征值问题结构的自振特性是指结构的振动频率和振型，求结构的自振频率和振型也称对结构进行模态分析，是结构动力计算的主要内容之一。计算经验指出，结构的阻尼对结构的频率和振型的影响很小，所以求频率振型时可以不考虑阻尼的影响。此时系统的自由振动方程如式(6.10)，即K8+MS=0(6.21)当系统作自由振动时，各质点作简谐振动，各节点的位移可表示为8二cos3t(6.22)将(6.22)代入(6.21),并消去公因子COSt得到(K2M)(M=0或K2M0(6.23)因此，求解(6.21)式就是寻找满足式(6.23)的2值和非零向量V，这种问题称为广义特征值问题。记九二2，九和分别称为广义特性值和广义特征向量。式(6.23)可以写成(K-XM)二0(6.24)这是一个齐次的线性方程组，若要有0的非零解，系数行列式必须等于零，即(6.25)K-XM|=0展开此式可得K-XMK-XMK-XM111112121n1nK-XMK-XMK-XM222122222n2nK-XMK-XMK-XMn1n1n2n2nnnnn，如果弹性结构的总体刚度矩阵Ik和总体质量矩阵Im的阶数是n则上述行列式展开后成为X的n次代数方程式，由此可以求出n个根，即n个广义特征值X,i=1,2,.,n，从而求出结构的n个自振频率i卩=jXC=1,2,n)。求得广义特征值X后，就可以利用式(6.24)算得对应的广义特征向量务,它代表n个质点的振幅构ii成的振型。由式(6.24)可知，务只能被确定到相差一常数因子的程度。N个特征值X和相应的n个特征ii向量务常称为n个特征对。在弹性结构的振动分析中，结构自由度总数n往往很大，因此无法直接从上述代数方程求解广义特征值九,i=1,2,，n。现在已研究出多种有效的广义特征值问题求解方法和相应的程序模块，供求解上述i问题使用。6.6特征值和特征向量的性质下面讨论广义特征方程K二灯M中特征值九和特征向量务的性质，这些性质在结构动力分析中是很重要的。1当Ik和Im矩阵是实系数对称矩阵时，其特征值一定是实数。如果Ik为正定矩阵时，则特征值九一定是正实数，如果h为半正定矩阵时，则特征值九一定是非负实数。并且，特征向量也是实向量。2.当Ek和m是对称矩阵时，广义特征方程的不同特征根九、所对应的特征向量1*ijij有正交性。当对特征向量务进行正则化处理后，即令TM(H=1ii时，则特征向量与质量矩阵和刚度矩阵的正交性，可表达为式中ij(6.26)为克罗内克符号。式(6.26)第一式称作正则化条件，而把满足该式的特征向量称作正则化特征向量。设1a是由特征值组成的对角阵，即(i=1,2,n)(6.27)又设Bl是由n个特征向量组成的特征向量矩阵，即=)e12n则特征向量的正交性，即式(6.26)也可表示为IkbLa!mMJl(6.28)(6.29)3.瑞雷(Rayleigh)商关于特征值九的性质，也可从瑞雷商出发去研究它。对于广义特征值问题，瑞雷商定义为P(V)=(6.30)将特征值由小到大排列成如下顺序：九XX12n设一向量U是n个特征向量的线性组合，可表示成(6.31)(6.32)|)=Ya务=loKaiii=1其中b为式(6.28)定义的特征向量矩阵，而是由n个系数组成的向量。将式(6.31)代入(6.30)，并注意到特征向量的正交关系式(6.29),得P)出展开后可写成工a2Xp()=-4=Ya2ii=1可见，瑞雷商在最小特征值和最大特征值之间Xp()X1n并且，当务=务时ip()=p()=Xii这说明当务取第i阶特征向量务时瑞雷商达到它的一个极值，该极值就是与务对应的特征值X，特iii别是minp()=p()=X11这就是瑞雷商的极小值原理。此原理常被工程界作为预估振型的依据去计算特征值X，从而求得基频11=込的近似值。i、i女口果a-a-=a(2mn-1),即向量务和前m-1阶特征向量正交，则由式(6.32)看出12m-1Xp()Xmn且当务=务时mminp()=p()=Xmm这就是说，X是瑞雷商在向量6与前m-1阶特征向量正交条件下的极小值。m4.移位在特征值问题的求解过程中，广泛使用所谓移位去加速计算的收敛性对广义特征值问题，可以使Ik作一移位卩，即(6.33)kLk-rm然后去求解下列特征值问题为了确定原特征值问题与上述特征值问题的关系，我们将式(6.33)代入上式，得(6.34)Ik!|/=(卩+耳為K|/与式(6.24)对比，显然有九二卩+耳UUk.(6.35)iii因此，Lb的特征向量与k血=九Dmh的特征向量是相同的，但特征值减小一个卩值。6.7逆迭代法求解特征值问题的方法很多，但基体上可分为两大类：变换法和迭代法。熟知的雅可比法就是一种变换法，可用于求问题的全部特征对，具有简单和稳定的优点。迭代法有所谓正迭代和逆迭代两种。由于有限元分析中，系统的自由度往往很高，而对工程结构进行动力分析时，通常只需了解少数几个较低阶的特征值和相应的特征向量，逆迭代法正是适应上述要求而被广泛应用的一种方法；在求解大型特征值问题的子空间迭代法中也要用到逆迭代法。现对此法作一讨论。设(0)是某特征向量的一个初始近似向量，并假定九二1。于是可求出式(6.23)右端的惯性力(0)(0)h由于(0)是任意假定的，一般不能满足平衡条件，即K(0)Lb(0)而按照静力平衡的关系k(i)Lb(o)可解得另一近似向量通常是一个比(o)h更好的近似特征向量。通过反复迭代就能得到满足精度要求的特征向量。在逆迭代法中，假定k是正定的。迭代公式为K(k+1)=%仏(k=0,1,2,)(6.36)由于特征向量各分量的模是任意的，只是它们之间保持一定的比例。为确定起见，在每次迭代中，对求出的(k+i)h进行正则化处理。利用正则化条件(6.26),可得正则化了的(k+ij,即(k+1)(k+1)a(k+1)HMX(6.37)如果所选初始迭代向量不与第一阶特向量关于矩阵Im正交，即1(0)hNh0(k+1)11关系式(6.36)和(6.37)是逆迭代的基本公式，然而在计算机中执行时，按如下方式进行更为有效。假设Y(0)=Im(0)按k=0,1,2,，进行下列迭代计算：Ik1X、(k+1)a=,Y(k)、Y(k+1)=p(k+1)=X/+1)鼻6*丁丿)Hx(k+1)Y(k+1)a2(k+1)TY(k)(k+1)(6.38)(k+1)式中，若Y(0)h0，则当kT8时1Y(k+1)p6(k+1)九1在迭代过程中，由(6.38)的第三式可求得由瑞雷商卩6(k+i)a给出的特征值九的近似值。该近似值可用来确定迭代法所达到的精度。当条件=10-2sp(k)-p(k-1)满足时，(6.39)(6.40)为了保证特征向量准确到r位精度，则应要求特征值准确到2s位精度，正如式(6.39)所表示的。这个要求可用瑞雷商来证明。现在来证明逆迭代的收敛性。设初始迭代向量是n个特征向量的线性组合，可表示为&(0血十a2iii=1式中系数a不能为零。应用式(6.36)和(6.37)进行迭代，并注意到关系式(6.26)和(6.24),可将第k次迭代后的向量(k+1)(k+1)(k+1)工aiu九i(i7(=12丿工匚丫2、九2k丿i(6.41)如果假设九1冬，将上式分子和分母同乘以以k，则得(k+1)、.丿i(九)1丿/i。由上式也可看出，(6.42)当当1为特征值时，b(k+1)趋于它们对应的特征向量的线性组合。由上式可以看出，迭代收敛的速度与比值九J有关，比值愈小收敛愈快。当九汎=0.99999时，收敛可能很慢；当九汎=0.01收敛就可能很快。如果我们采用移位，将九移至12121原点附近但仍保持正值，则收敛速度会很快。移位不但可用来加快收敛速度，而且可用来求中间特征对。如果作一移位卩，依据式(6.34),使n-卩为负值，而n-卩为接近于零的正值，则逆迭代的结果就可求出特征对匕,%。同112222样可以用移位来求出其它特征对。当用移位法依次自小至大求特征对时，有时很难收敛到指定的特征对，由于计算中的舍入误差很可能收敛到已求出的特征对，尤其当存在重特征值时更是如此。为了不收敛到已求出的特征对，可以利用特征向量的正交性，在开始选择初始迭代向量时就进行正交化处理。具体说，在求出m-l阶特征向量后，可以构造一新的迭代向量(0)(oJ-0aiii=1(6.43)其中b爲是任意向量，&是待定常数，它们由b爲与如,如,L关于质量矩阵m的正交条件确定。由此，应取ai(0)1,2,，m-1)(6.44)由于!X(0)為前m-1阶特征向量关于Im正交，则用它作为初始向量进行逆迭代可指望求得第m阶特征对(X,G。为了防止计算中的舍入误差而使(0)偏离正交方向，还应在每一步开始都对mm行正交化处理。以上迭代过程称为克莱姆一史密特(Gram-Schmidt)正交化。例6.1图6.1示一三层刚架结构，各层的楼面质量分别为mi=180t,m2=270t，m3二360t；各层的侧移刚度分别为k=98MN；m,k?=196MNm,k?=294MNm。求刚架的固有频率和型。解该刚架的质量矩阵和刚度矩阵为100_m=180X10301.50002r1-10k=98X106-13-20-25我们先用行列式法直接求这一问题的特征值。由式(6.25)0-2(a)1一申一1K一2M=98x1061-13-1.5申18098x103令式(a)矩阵行列式展开、化简并令其等于零、得到三次方程，解之得9=0.3514,9=1.607,9=3.542123从而得到=13.831s,=29.61s,=43.91s而由式(6.23)可得三个主振型向量，用振型矩阵表示为1.0000_b=0.6486-0.6070-2.5420.3018-现在用式(6.38)的逆迭代法来解同一问题。为以后迭代计算方便先求出刚度矩阵的逆矩阵588x1061152552222采用初始迭代向量/門k丿、1.0_10、180x1031.52.0一当k=1时x(1)Lk11Y(0)L180588x10322.59卜2.516.59】(1)Y(0)/p=卜2.516.59】1.5,588x1032.022.524.75,180=197.9797918.6857250.7542980.548580】T以下各次循环做法相同，现将结果载入下表6.1内。从上表看出，p閃比Y（k）嫩敛得快些经过5次迭代就已经达到了要求，p閃唱、因此由得到u=1（i+1）(.742622、1.000000、=0.648587,0.2241910.3018912(1+1)（1+1）Y1x0.742622k卜(k)Y（k）P(k)P(k)-p(k-1)Y（k）P(k)122.522.5197.97980.685728可见第一振型的结果与前面用行列式法求得的结果相符合。表6.116.5924.75180.7542980.548580212.411612.4116191.63930.03308530.7301558.2972812.44590.7321733.977217.954420.467946312.628512.6285191.36640.00142630.7398908.2475412.37130.7248213.860557.721100.452317412.667612.6676191.35370.00006660.7420438.228312.34250.7229993.834177.668340.449196512.675912.6759191.35300.00000320.742528.2236112.33540.7225753.828447.656960.448524612.677612.6776191.35300.00000020.7426228.2225312.33380.7224833.827247.654480.4483816.8用逐步积分法求结构动力响应概述求结构的动力响应，在数学上就是要求出运动方程(6.7)的解答。式(6.7)是一个二阶常系数微分方程组，可以用数值积分的方法对方程直接求解，即按时间增量At逐步求解运动微分方程，直至反应终了，这一方法称作逐步积分法，由于在数值积分前不对运动方程进行变换，所以又称作直接积分法。逐步积分法既可用于求解线性结构体系问题一一在整个动力反应过程中Ik,m,c矩阵保持不变的问题；也可用于求解非线性结构体系的问题一一k,m,c矩阵随动力反应的过程而变化的问题。此处，我们只讨论线性结构体系的问题。逐步积分法求解运动微分方程的基本思路是：把连续的时间过程离散为ti,”一有限个点，对于运动微分方程mC胡LK8=F只要求它们在上述每个时间离散点上得到满足，也就是说，最终求解得到的只是位移、速度和加速度在有限个时间离散点上的值，而不是连续函数幺!&t)!匕。2.在每个时间间隔At内，假定位移、速度和加速度符合某一简单的关系。而At的选择，要求保证计算的稳定性与精确性。从这样一个基本思路出发，形成了多种逐步积分的方法，如中心差分法；线性加速度法；Wilson-0法；Newmark法等等。逐步积分法的优点在于，计算的精确度可随At的减少而提高，当大型结构体系受到历时较短的脉冲荷载作用时，由于要考虑高阶振型的反应使计算振型的数目较多，而逐步积分法求解动力反应时，可省去特征方程的求解工作，在经济上、效率上可能是有利的。但是，逐步积分法的缺点，在于它不能同时给出结构的振型和自振频率，而是直接给出结构的反应过程，并且当反应的时间历程较长，At取的较短时，计算工作量也是相当可观的。法Wilson-0法实质上是线性加速度法的推广，线性加速度法假定从时刻t到t+At的时间内加速度是线性变化的。而Wilson-0法则假定在时间区间t到t+At内加速度是线性变化的，如图6.1,其中01.0。当0=1.0时，Wilson-0法法就退化为线性加速度法。Wilson证明，要达到无条件稳定，必须选用01.37,通常采用0二1.4。图6.2Wilson-0法按Wilson-0法的假定，在时间区间t到t+At内加速度的数学表达式为i(t+tT0At(6.45)式中，T为时间的增量，0T8(t)+(1-召8t)(6.60)02At202At20由式(6.57)得10(8(t+0At)-8-(t)二8-(t+At)-8-(t)将上式代入(6.58)和(6.59)式得8(t+At)二8(t)+At(8-(t+At)+8-(t)(6.61)At28(t+At)=8(t)+At8(t)+(8(t+At)+28(t)(6.62)6Wilson-0法的计算步骤归纳如下：(1) 初始计算形成总刚度矩阵Ik和总质量矩阵m，而总阻尼矩阵C=aM+卩K不存贮。(2) 形成初始值8(0),8(0)，并计算-(0)。选取时间步长At和0(一般取0=1.4)，计算积分常数。得63宀0Ataa43At1Aa=,a=,a=2a,a=,a=0,a=2,a=1,a=,a=At20(OAt)21OAf21324e5e6e72*6形成有效刚度矩阵K=K+aM+aC01对K进行三角分解K=LtDL1. 对每一时间步长，循环计算计算t+0At时刻的有效载荷向量F(t+OAt)=F(t)+O(F(t+OAt)F(t)+M(a6(t)+a8(t)+2S(t)02+C(a6(t)+26(t)+aS(t)13求解t+OAt时刻的位移向量6(t+OAt)LtDL6(t+OAt)=F(t+OAt)(3)计算t+At时刻的加速度、速度和位移向量S(t+At)=a(6(t+OAt)6(t)+aS(t)+a8(t)456S(t+At)=S(t)+a(S(t+At)+S(t)76(t+At)=6(t)+AtS(t)+a(S(t+At)+2S(t)8例6.2有一简单的结构系统，不考虑阻尼影响，它的运动方程是20_卩(t)62_6(t)0、+016-(t)J2丿246(t)J2丿、10一试用Wilson-O法计算此结构系统从0到3.5的位移响应。已知结构初始状态是静止的，受突加不变载荷作用。解：选择时间步长At。为保证计算精度，通常取0.1，其中T为结构系统自由振动的最小周期。对于复杂结构系统，一般高频部分对结构动力响应的贡献不大，可取低频部分的最小周期。本例的最小周期T=2.8s，取At=0.28s，计算12个时间步长。将“(0)=0,5(0)=0初始条件代入运动方程，求得加速度初始条件6(0)=010T。取O=1.4，计算Wilson-O法的常数，得a=39.0a=7.65a=15.3012a=0.196a=27.9a=10.9345a=1.14a=0.14a=0.0131678形成有效刚度矩阵_622084.02_+39.0=2410243.0K=K+aM+aC=01对于每一个时间步长，需先求解K5(t+OAt)二(t+OAt)式中，卡(t+OAt丄00V+、1002001(39.0(t)+15.3S(t)+2&(t)然后计算加速度、速度和位移i(t+Aty=27.9(t+OAt)-6(t)-10.9S(t)-1.14S(t)i(t+At)L6(t)+0.14(6(t+At)+6(t)匕(t+At)=6(t)+0.286(t)+0.0131(6(t+At)+26(t)12个时间步长按以上诸式计算的结果如下表所示：tAt2At3At4At5At6At7At8At9At10At11At12At610.006050.05250.1960.4900.9521.542.162.672.922.822.331.54620.3661.342.643.924.885.315.184.613.823.062.522.29法Newmark法也是一种线性加速度法的推广。它采用如下两个基本假定6(t+At)-6(t)+(1-Y)6(t)+Y6(t+At)At(6.63)At26(t+At)=6(t)+6(t)At+(12P)S(t)+2卩6(t+At)(6.64)其中0和Y是参数，根据积分的精度和稳定性要求确定。1111当Y=和卩=时，关系式(6.63)、(6.64)相当于线性加速度法。当丫=和卩=4时，相当于常平均加速度法。已经证明，当Y工0.5和00.25(0.5+丫)2时，Newmark法是无条件稳定的。根据上述两个假定，以及t+At时刻的运动方程，就可以求得t+At时刻的位移、速度和加速度的解。将6(t+At)作为基本未知量，由式(6.64)求得(6.65)(6.66)(6.67)6(t+At)=丄(6(t+At)-6(t)-丄6(t)-注-1)6-(t)0At20At20将式(6.65)代入式(6.63)得6(t+At)=0(6(t+At)一6(t)+(1-0)6(t)+(1-?0)At6(t)0At020将式(6.65)和(6.66)代入(t+At)时刻的运动方程M6(t+At)+C6(t+At)+K6(t+At)=F(t+At)得到只有一个基本未知量6(t+At)的方程式中(6.68)K5(t+At)二F(t+At)(6.69)1卩At2F(t+At)=F(t+At)+M(丄(t)+吕6(t)+(亠-l)(t)卩At2卩At2卩(6.70)+C二5(t)+(I-1)8(t)+(丄-1)AtS-(t)卩At卩2卩解方程(6.68)得到5(t+At)后，再代入式(6.65)和(6.63)求出5(t+At)和5(+At)。为了便于编写计算机程序，将Newmark法的计算步骤综述如下：1.初始计算(2)(1)形成刚度矩阵Ik,质量矩阵Im,阻尼矩阵C二M+卩K不存贮。形成初始向量5(0)、5(0)，计算(0)。(3)选择时间步长At和参数卩,Y。计算积分常数，得a=21a卩At3(4)=P-形成有效刚度矩阵At(y=2(卩-2),a=At(1-y),a=yAt67K=K+aM+aC01(5)对K进行三角分解K=LtDL。2.对每一个时间步长循环计算在时刻t+At的有效载荷向量F(t+At)=F(t+At)+M(a5(t)+a5(t)+a5(t)+C(a5(t)+a5(t)+a5(t)023145求解t+At时刻的位移向量LtDL5(t+At)=斤(t+At)计算t+At时刻的加速度和速度5(t+At)=a(5(t+At)-5(t)-a5(t)-a5(t)0235(t+At)=5(t)+a5(t)+a5(t+At)67可以看出，Wilson-0法和Neumark法的计算步骤是相同的，计算公式的形式也是一样的，只是系数有所不同，故可以在一个有限元程序中同时存在这两种直接积分法供用户选择。例6.3用Newmark计算例题6.2的结构系统的响应。取0=0.25,丫=0.5。解：选择At二0.28s。将=0代入运动方程,得*(0)=10,10卜计算纽马克法的常数a二51.00a=7.141a=14.32a二1.003a=1.004a=0.005a=0.146a=0.147形成有效刚度矩阵为-220+51.0401108-2-255对于每一个时间步长,先求解下列方程“山(t+At)=$(t+Aty式中01(51.05(t)+14.3然后计算时间At2At3At4At5At6At7At8At9At10At11At12At510.006730.05040.1890.4850.9611.582.232.763.002.852.281.40520.3641.352.684.004.955.345.134.483.642.902.442.31用以上诸式,按12个时间步长依次计算的结果如下6.9用振型叠加法求动力响应在直接积分法中,将整个动力响应历程划分成若干个时间步长。为保证计算精度,必须把每个时间步长的时间间隔At取得比较小，每个时间步长都要求解一次线性方程组(t+At)=51.0(5(t+At)-(t)J)-14.35(t丄1.05(t)5(t+At)=5(t)+0.145(t)+0.145(t+At)“4(t+At)=$(t+Aty或“鸟(t+0At)=(t+BAt)j尽管对于线性弹性结构，有效刚度矩阵可以在初始计算时就三角分解好，但每个时间步长进行前消回代仍然要花费不少机时。一般来讲，计算短时间的动力响应，使用直接积分法是很有效的，但是，要进行长时间的动力响应计算，直接积分法就变得不经济了。如果在逐步求解开始之前，先将运动方程进行变换，把有效刚度矩阵化成对角矩阵，那么每一时间步长的计算就非常简单，使逐步求解变得比较经济。振型叠加法在逐步求解之前；先将运动方程进行变换化简。考虑到弹性结构无阻尼自由振动k怡=九Imh时的振型向量6具有对质量矩阵和刚度矩阵的正交性，即时Im=5VTk血川ijiij式中，5ij0,i丰j1,i=j因此，采用振型矩阵“=“,e，,e作为变换矩阵，总体位移向量匕（t）可表示为12n匕(t)=eh(t)(6.72)或5(t)=x(t)G+x(t)G+xG1122nn其物理意义是结构的总体位移向量可以用结构振型i=1,2,n的线性组合来表示。由i于振型矩阵1在求解结构自由振动，即求解特征值问题时已得到，因此计算位移向量暫（t）的问题变成了计算未知向量h（t）的问题。需要特别指出的是，振型矩阵与时间无关。将（6.72）式代入运动方程，并左乘k得mm儿恥）+側卜風（t）+側k儿）=側曲）（6.73）考虑到振型的正交性，则斜m=b式中九W211九0W20q=2-20.0九W2nnWi固有频率，等于匹i=1,2,.,n;(6.74)Iinxn阶单位矩阵。如果采用瑞利阻尼，即c=aM+pK那么振型对阻尼矩阵也是正交的，其表达式为%TC%=0i丰jij命2eg112eg220(6.75)2egnn/(t)=Mb(t)式中，g称为振型阻尼系数i=l,2,n。i因此式(7.73)可写成(6.76)x(t)+Mb儿施(t)+hKx(t)=务(t)由于卜儿1匕都是对角阵，(6.76)式是一组相互独立的二阶微分方程，可写成:x(t)+2wgx(t)+w2x(t)二丫(t)i二1,2,n(6.77)iiiiiii式中Y(t)=GtR)(6.78)ii(6.77)式是二阶常微分方程，可以用直接积分法来求解，也可以用杜哈梅(Duhamel)积分来求得，如x.(t)=.Y.(t)e_gi-(t_T)sin0(t-t)dT+e-g巴t(asinQt+PcosQt)(6.79)iQ0iiiiiii式中，Q=qV1-g2。a和P由初始条件确定。杜哈梅积分一般采用数值积分来计算。iiiii向量“的初始值，即t=0时的x(0)值，可由下式求出(6.80)i=1,2,nx(0)=Gtm脸(0)iix(0)=Gtm(0)ii上式是由(6.72)式匕(0)=Dh(0)和振型正交性m血=l得出。知道各个时间步的向量&(t)后，就可以求出该时间步的位移向量。匕(t)=)J!x(t)从上面的讨论看出，振型叠加法与直接积分法之间的区别在于振型叠加法在开始逐步积分之前进行了基的变换，即从有限元的坐标基变换成广义特征问题1K怡=九1M山的特征向量基。从数学上讲，由于n个特征向量所张成的空间和有限元n个节点位移所张成的空间是相同的，所以两种解法的结果应该是相同的。然而实际使用的振型叠加法是近似计算这是因为要计算阶数n很高的广义特征值问题的全部特征对(九，)i=1,2,n是很ii困难和很费机时的，并且求解广义特征值问题对于低频率和相应的振型有较好的精度，而对于高频率及其相应的振型计算精度就很差。通常在振型叠加法中仅考虑部分低频率及其相应振型，而略去其余高频的影响。实际采用振型叠加法计算弹性结构的动力响应时，做了如下改动。结构的总体位移向量用p个低频振型的组合来表示，即匕(t)=Dh(t)(6.81)式中1将(6.81)式代入运动方程，再左乘6)，并考虑振型的正交性，最后得到一组p个独立的二阶常微分方程x(t)+x(t)+w2x(t)二丫(t)i二1,2,p(6.82)iiiii初始条件是(6.80)式。各个时间步的位移向量匕(t)=l*仁(t)。nxppxl由于在振型叠加法中略去了高频响应，因此适用于如地震响应，低频载荷引起结构振动等激发振型较少，高频响应不重要而需计算较长响应时间的问题。而对于象冲击等那样激发振型较多，只要计算较短响应时间的问题，通常采用直接积分法为宜。这里补充一点。如何确定瑞利阻尼b=a(M+plK的比例系数和0。对弹性结构有阻尼自由振动测定振幅的衰减，可以求得不同频率的振型阻尼系数。例如，知道频率时的振型阻尼系数2和频率时的振型阻尼系数g。由(6.75)式1122=GtC)=Gt(am+pEk|)=a+02iiiiiii得联解2g二a+p21112g=a+022222(g_g)a=U-212(6.83)(_)(+)p2(g1_g)1p=-2-1(_)(+)2121例题6.4试用振型叠加法计算例题6.2的结构系统的位移响应。解：首先计算结构自由振动的振型矩阵“。由于结构系统的运动方程为20r、-6_2rs01+=01824s111011212丿它的广义特征方程为本例比较简单，现采用代数方法计算。根据广义特征多项式det-九二2X2-14九+20二0丿可以求出两个根，即特征值九=2和九=5。12分别将X和X代回广义特征方程，可以求得对应的特征向量非零解12再进行规范化处理于是，振型矩阵山31：22片312：3按(6.76)式求得以特征向量为基向量的运动方程上式是二个不耦合的运动方程，即X+2x=-2r113迂X+5x=-10.1223根据结构系统的初始条件乞(0)=(0)L,可按(6.80)式求得x(0)=01X(0)=01x(0)=0x(0)=022二个不耦合的运动方程的解析解是x=2(1-cos、：2t)13因此,代入(6.81)式得122；3-23-cos迈t)/2IJ1+cosV5t3取At=0.28，用上式算出各时间步长的位移如下时间At2At3At4At5At6At7At8At9At10At11At12At510.0060.0450.1700.5201.051.722.3382.8613.0522.8012.1301.157520.3791.422.794.125.045.334.985