第六章动力问题有限元法

1、第六章动力问题的有限元法6.1概述前面几章所研究的问题都属于静力问题，其特点是施加到结构上的外载荷不会使结构产生加速度,且外载荷的大小和方向不随时间变化，因而结构所产生的位移和应力也不随时间变化。本章将要研究 结构分析中另一类重要问题的有限元解法，即动力问题的有限元解法。动力学问题的特点是，载荷是 随时间变化的，因而结构所产生的位移和应力是时间的函数，结构会产生速度和加速度。由于结构本身的弹性和惯性，结构在动力载荷的作用下，往往呈现出振动的运动形态。结构振动 是工程中一个很普遍很重要的问题。有些振动对我们有利，例如，振动打桩，振动选料，有些振动对 我们有害，例如，机床的振动，仪器与仪表的振动，

2、桥梁、水坝及高层建筑在地震作用下的振动等。 因此，我们必须对振动体本身的振动特性以及它对外部激振力的响应有一个明确的认识，才能更好地 利用它有利的一面，而避免它有害的一面，设计出更好的机械和结构。振动问题主要解决两方面的问题。1. 寻求结构的固有频率和主振型，从而了解结构的固有振动特性，以便更好地利用或减少振动。2. 分析结构的动力响应特性，以计算结构振动时动应力和动位移的大小及其变化规律。6.2结构的振动方程结构的振动方程可用多种方法建立，这里我们使用达朗伯原理(动静法)，仿照前几章建立静力有限元方程的方法，来建立动力问题的有限元方程。在静力问题中用有限元法建立的平衡方程是K 二F在振动问题

3、中，对结构的各节点应用达郎伯原理所建立的振动方程仍然具有与上式相同的形式，只不 过节点位移是动位移，节点载荷是动载荷，它们都是时间的函数。上面的方程成为上式中J(t)为节点的动位移，它是时间的函数, 该时刻的节点外力 Q(t)丿构成动态平衡。在动态情况下，结构承受的载荷(集中载荷K、(t)二Q(t)(6.1)K、：(t)是t时刻的节点位移产生的弹性恢复力，它与，分布载荷)可随时间而变化，是时间的函数。按有限元方法将此种载荷移置到节点上，得到的节点载荷向量F (t)也是时间的函数。此外，结构在运动中，各点除位移:f 以外，还有速度。按照达郎佰原理，有加速度的质量应附加有惯性力载荷。如材料的密度为

4、，则结构单位体积的惯性力为说，相当于又受有另一种体积力，大小与点的加速度成比例，而方向与加速度方向相反。另外，在结构运 动过程中，还会受到周围介质和来自内部的阻力。精确地描述这种阻力的变化规律是很困难的，一般采用 阻力与速度 f ：成比例的近似线性假定，如阻力系数为卩，则单位体积的阻力为 -叫f。这对结构来说 相当于另一种体积力，大小与点的速度成比例，方向与速度方向相反。按有限元方法，用单元节点位移技进行插值表示单元内部位移。(6.2)此处形函数仍只是位置的插值函数，与时间无关，则单元内的速度和加速度分别为门二N、e(6.3)以及* 九N、e(6.4)其中GF、怯丁为单元节点的速度及加速度向量

5、。将单元惯性力_pf"与阻力_a f作为体积力，按式(3.28)移置到单元各节点，就得到相应的单元等效节点载荷向量，记为示丁和，则有F J = -jNl fdv匸= - N叫 fdv将式(6.3)、( 6.4)代入上式，有牛 J 二一 NT TNdv 二-Nt 订Ndv；二 一 mef J = JNTMN怯fdv =JNT UNdv= -c式中me = NT rNdv(6.5)称为单元质量矩阵。由于推导式(6.5)时采用了与推导单元刚度矩阵时相一致的形函数，故式(6.5 )所表示的质量矩阵也称为一致质量矩阵。而Ce = NT 叫Ndv(6.6)称为单元阻尼矩阵，由式(6.5)和(6.

6、6)可见单元质量矩阵和单元阻尼矩阵是对称的。将移置到节点上的动载荷、惯性力、阻力作为载荷，按单元叠加，得到有限元节点位移方程：mK 二F -、(me、e ce、e)或M'Cf |K' 二F(6.7)其中mM八me(6.8)e=1称为结构质量矩阵或总质量矩阵，而mC八ce(6.9)e#称结构阻尼矩阵。可见结构质量矩阵和结构阻尼矩阵分别为单元质量矩阵和单元阻尼矩阵的叠加，其叠加方法与结构刚度矩阵的形成完全一样，借助于单元定位向量，用单元集成法完成叠加过程。由于me和ce是对称的，因而叠加合成的M 1和C 1也是对称的。式（6.7）是节点位移的二阶微分方程，称为结构的动力方程式。对于

7、不同的结构，可以选用不同的 单元，有不同的形函数矩阵 N 1,但动力方程（6.7）的建立过程都是一样的。当结构不受外载荷时，；F心0，如果再忽略阻尼，则动力方程（6.7）式成为M、 KJ =0（6.10）这是系统的自由振动方程。弹性结构的振动实际上是连续体的振动，位移（f 是连续的，具有无限多个自由度。经有限元离散化处理后，单元内的位移按假定的位移形式来变化，可用节点位移插值表示。这样，连续系统的运动就离散化为有限个自由度系统的运动了。如果全部节点有N个自由度，则式（6.10）就是N阶的自由振动微分方程了。6.3 一致质量矩阵与集中质量矩阵按照式（6.5），当单元的位移插值形函数矩阵确定后，就

8、可用此式算出它的一致质量矩阵。例如，对 平面三节点三角形单兀，按照2.3节的叙述，形函数矩阵 N可用面积坐标表示为N二L ILj ILm其中I为二阶单位阵。将上式代入式（6.5）可得平面三节点三角形单元的一致质量矩阵etm=N r【N(ILILj IL ILmILL= PJj|ILjLiJLmLi利用积分公式(2.56)，可由上式求得ILj ILmdxdytILiLjILiLm IIL jLjIL jLm dxdytILmLjILmLm 一'1111000244111000244111000424111000424111000442111000442 一(6.11)用式（6.5）可以计

9、算出其它类型单元的质量矩阵。例如平面桁架单元的质量矩阵为而平面刚架单元的质量矩阵为140me70°式(6.11)、(6.12)及(6.13)中的M为单元质量。1 11nz 3 6：6 Y(6.12)15622L4L25413L-13L-3L2对称14001560-22L 4L(6.13)由单元一致质量矩阵叠加形成的结构质量矩阵，一般都是稀疏、带状的，但都有相当的半带宽，如果将单元质量矩阵近似作为对角型矩阵，单元叠加后的结构质量矩阵也是对角阵。这种对角型的质量矩阵称 为集中质量矩阵，而对角型的集中质量矩阵对结构的动力分析是非常有利的。可以将单元的质量以某种方 式分配在单元的节点上而得到

10、单元集中质量矩阵。例如当质量均匀分布时，将质量平均分配给各节点。按 照这种方法平面三节点三角形单元的集中质量矩阵为11其中I 1为6阶单位阵。M为单元质量。而四节点四面体单元的集中质量矩阵为m(6.14)meI(6.15)其中I为12阶的单位阵。对于矩形弯曲板单元，在不计转动质量影响的条件下，它的集中质量矩阵为(6.16)对于平面刚架单元，在不计转动质量影响的条件下，它的集中质量矩阵为1e M(6.17)m-6.4 阻尼矩阵各种工程结构的阻尼力及其产生的机理是非常复杂的。从宏观上看，阻尼有两种主要形态。一种是结 构周围粘性介质产生的阻尼，称为粘性阻尼。粘性阻尼的阻尼力一般近似与运动速度成正比。

11、另一种是结 构材料内部磨擦产生的阻尼，称为结构阻尼或材料阻尼。结构阻尼的阻尼应力一般近似地认为与弹性体的 应变速率成正比。如果假定阻尼力与运动速度成正比，那么在运动的弹性体中任意点处单位体积上作用的阻尼力为© f二"N匚ct式中:-比例常数;?材料密度;N 1单元形函数矩阵;V，单元节点速度矢量。可以将阻尼力Pd看成是一种体积力，其等效的单元节点阻尼力向量为Fd° = . NTPddvNT 订NdvJvv或写成Fde =c:因此，单元阻尼矩阵eTeC :Jill N订 Ndv 二：m(6.18)v它正比于单元质量矩阵me。如果假定阻尼应力与弹性体的应变速率成正比，

12、则阻尼应力可表示为江九 1D=扌二订DB、Je式中，一：为比例系数，D 1为弹性矩阵，B】为应变矩阵，1了为单元节点速度向量。下面推导阻尼应力的单元等效节点阻尼力向量。设单元等效节点阻尼力向量仍用Td尹表示，设节点发生虚位移f *，单元内各点产生的虚应变为 ；*，则丄B、*。e*Fd在虚位移f 上做的虚功为We*TFde单元的虚应变能为* t* TtUe = .； %dv= ,.B ；ddvvv由We =Ue得到eFde = .BTJdv 一 .BTDBdv、vv或写成Fde 二ce e因此，单元阻尼矩阵ce ： ； BtD Bd -ke(6.19)v它正比于单元刚度矩阵ke在实践中，要精确地

13、计算各种单元的阻尼矩阵是很困难的。通常在程序设计中，假定结构总体阻尼矩阵C是结构总体刚度矩阵 K与总体质量矩阵 M 1的线性组合。称为瑞利阻尼，其表达式为(6.20)Cp - M ：K其中，比例系数:及1可通过实验确定。采用瑞利阻尼近似，可以使运动方程求解大大简化，并且在程序中不必单独存贮总体阻尼矩阵。在实 际工程问题中，阻尼的作用对结构的动力响应的影响并不大，这种近似处理具有实用价值。6.5结构的自振特性和特征值问题结构的自振特性是指结构的振动频率和振型，求结构的自振频率和振型也称对结构进行模态分析，是 结构动力计算的主要内容之一。计算经验指出，结构的阻尼对结构的频率和振型的影响很小，所以求

14、频率 振型时可以不考虑阻尼的影响。此时系统的自由振动方程如式(6.10)，即K、 M、 = 0(6.21)当系统作自由振动时，各质点作简谐振动，各节点的位移可表示为 ' = cos -1(6.22)将(6.22)代入(6.21)，并消去公因子 cos t得到(KP 2M) = 0或K =2M (6.23)因此，求解(6.21)式就是寻找满足式(6.23)的，值和非零向量 心， 这种问题称为广义特征值问题。记 = ' '2， 和分别称为广义特性值和广义特征向量。式(6.23)可以写成(K- M) =0这是一个齐次的线性方程组，若要有 的非零解，系数行列式必须等于零，即K-

15、 M0(6.24)(6.25)展开此式可得K11 一 M 11K12 一M 12K1n 一 /-M 1nK 2n 一 '-M 2n=0心1 J； M n1 Kn2M n2nnK 22 一 M 21K 22 一M 22如果弹性结构的总体刚度矩阵K 1和总体质量矩阵 M 1的阶数是n,则上述行列式展开后成为的n次代数方程式，由此可以求出n个根，即n个广义特征值'i , i = 1,2,., n，从而求出结构的n个自振频率'-i = i i -1,2, ,n 。求得广义特征值 、后，就可以利用式(6.24)算得对应的广义特征向量i I它代表n个质点的振幅构成的振型。由式(6.

16、24)可知，X ?只能被确定到相差一常数因子的程度。N个特征值；和相应的n个特征向量：i 1常称为n个特征对。在弹性结构的振动分析中，结构自由度总数n往往很大，因此无法直接从上述代数方程求解广义特征值i ,i =1,2- , n。现在已研究出多种有效的广义特征值问题求解方法和相应的程序模块，供求解上述 问题使用。6.6 特征值和特征向量的性质下面讨论广义特征方程 K M 中特征值，和特征向量的性质，这些性质在结构动力 分析中是很重要的。1当K 和 M 1矩阵是实系数对称矩阵时，其特征值一定是实数。如果K1为正定矩阵时，则特征值一定是正实数，如果 K 1为半正定矩阵时，则特征值一定是非负实数。并

17、且，特征向量V也是实向量。2 当K和M 1是对称矩阵时，广义特征方程的不同特征根、j所对应的特征向量 丄 °j /具有正交性。当对特征向量进行正则化处理后，即令 iTM i =1时，则特征向量与质量矩阵和刚度矩阵的正交性，可表达为(6.26)ij式中 为克罗内克符号。式（6.26）第一式称作正则化条件，而把满足该式的特征向量称作正则化特征向量。设是由特征值组成的对角阵，即I -1 - diag、i =1,2, ,n(6.27)又设/是由n个特征向量组成的特征向量矩阵，即(6.28)则特征向量的正交性，即式 （6.26）也可表示为(6.29)!：TKb ：.：打TM】13 .瑞雷(Ra

18、yleigh )商关于特征值'的性质，也可从瑞雷商出发去研究它。对于广义特征值问题，瑞雷商定义为(6.30)将特征值由小到大排列成如下顺序:'1 '2 -'n设一向量：；是n个特征向量的线性组合，可表示成(6.31)其中工 为式(6.28)定义的特征向量矩阵，而是由n个系数组成的向量。将式(6.31)代入(6.30)，并注意到特征向量的正交关系式(6.29),得展开后可写成a 卩 R Kaa Flant2aj迟a2i 1(6.32)可见，瑞雷商在最小特征值和最大特征值之间、_ L ? _ 'n这说明当取第i阶特征向量 时瑞雷商达到它的一个极值，该极值就是

19、与J ?对应的特征值-，特别是mi n 皿I这就是瑞雷商的极小值原理。 此原理常被工程界作为预估振型的依据去计算特征值 1，从而求得基频-11的近似值。如果a1 =a2二二am2乞m乞n-1 ,即向量、J和前m-1阶特征向量正交，则由式(6.32)看出、"心且当-m 时min戸"：“叮=爲这就是说，'m是瑞雷商在向量 U与前m-1阶特征向量正交条件下的极小值。4.移位在特征值问题的求解过程中，广泛使用所谓移位去加速计算的收敛性。对广义特征值问题，可以使 K 1作一移位，即k L K ：- M 1(6.33)然后去求解下列特征值问题总厂-M旷/为了确定原特征值问题与上

20、述特征值问题的关系，我们将式(6.33 )代入上式，得K 1 二- M C(6.34)与式(6.24)对比，显然有i -(6.35)因此，k.J - Ml :的特征向量与 K - Ml ?的特征向量是相同的，但特征值减小一个值。6.7逆迭代法求解特征值问题的方法很多，但基体上可分为两大类：变换法和迭代法。熟知的雅可比法就是一种变 换法，可用于求问题的全部特征对，具有简单和稳定的优点。迭代法有所谓正迭代和逆迭代两种。由于有 限元分析中，系统的自由度往往很高，而对工程结构进行动力分析时，通常只需了解少数几个较低阶的特 征值和相应的特征向量，逆迭代法正是适应上述要求而被广泛应用的一种方法；在求解大型

21、特征值问题的 子空间迭代法中也要用到逆迭代法。现对此法作一讨论。设：X (0)是某特征向量的一个初始近似向量，并假定 =1。于是可求出式(6.23)右端的惯性力乍丄(i) M汶(叭由于：x(0) ?是任意假定的，一般不能满足平衡条件，即K X(0) F (0)'而按照静力平衡的关系K Ex -匸(0) 1可解得另一近似向量»通常"X是一个比'X(0) 1更好的近似特征向量。通过反复迭代就能得到满 足精度要求的特征向量。在逆迭代法中，假定 K 是正定的。迭代公式为KkX(k1)'M：Xk* k=0,1,2,(6.36)由于特征向量各分量的模是任意的，只

22、是它们之间保持一定的比例。为确定起见，在每次迭代中，对求出的：X(k 1) '须进行正则化处理。利用正则化条件(6.26),可得正则化了的(k ° ，即(k 1):(6.37)汶(k 1) /1肢(5M欧(5”2如果所选初始迭代向量不与第一阶特向量关于矩阵M 1正交，即则当k 一. ，时x(0)匚鼻 0；x(k1)Cj关系式(6.36)和(6.37)是逆迭代的基本公式，然而在计算机中执行时，按如下方式进行更为有效。假设Y-M :X(0)按k =0,1,2;，进行下列迭代计算：k ：X(5=y(k)'代(小飞(2和,(Y(k -1) /(6.38)Y(k*)= M 3t

23、X(k*)X(k*)F V(k)Y(k 1)=(X(ry Y(k*)2 J式中，若Y(0) y=0，则当k时Y(k卩M爲/"X(k1)；j、在迭代过程中，由(6.38)的第三式可求得由瑞雷商t :X (k T"给出的特征值'1的近似值。该近似值可用来确定迭代法所达到的精度。当条件满足时，迭代就终止，令r(k) _ J(k4)(k)2s-；=10(6.39)i是最后 一次迭代，则有.1X(l d)】,.X(f1厂1；x(i p ：T Y(l 1<? 2(6.40)为了保证特征向量准确到 可用瑞雷商来证明。现在来证明逆迭代的收敛性。设初始迭代向量是n个特征向量的线

24、性组合，可表示为r位精度，则应要求特征值准确到2s位精度，正如式(6.39)所表示的。这个要求X(°)?八了岸i =1a：i =1式中系数ai不能为零。应用式(6.36)和(6.37)进行迭代，并注意到关系式(6.26)和(6.24)，可将第k次迭代后的向量汶(k 1):和”X(k 1)写成x(k1)一ak h I i2ai22 ai2k 'i(6.41)(k =1,2,)n 2k人丿如果假设岂J岂乞n，将上式分子和分母同乘以k i，则得Z aiX(k 1) : 心4、丁2Z ai11.卜' 丿J(6.42)显然，当 k ',X(k °'叫

25、I由上式也可看出，当'1为特征值时，"X(k ° f趋于它们对应的特征仁2有关，比值愈小收敛愈快。当向量的线性组合。由上式可以看出，迭代收敛的速度与比值仁匕=0.99999时，收敛可能很慢；当J / =0.01收敛就可能很快。如果我们采用移位，将-1移至原点附近但仍保持正值，则收敛速度会很快。移位不但可用来加快收敛速度，而且可用来求中间特征对。如果作一移位J，依据式(6.34)，使=人-4为负值，而 J =丸2 -卩为接近于零的正值，则逆迭代的结果就可求出特征对(扎2,饰2)。同样可以用移位来求出其它特征对。当用移位法依次自小至大求特征对时，有时很难收敛到指定的特征

26、对，由于计算中的舍入误差很可能 收敛到已求出的特征对，尤其当存在重特征值时更是如此。为了不收敛到已求出的特征对，可以利用特征向量的正交性，在开始选择初始迭代向量时就进行正交 化处理。具体说，在求出m-1阶特征向量后，可以构造一新的迭代向量m 4X(0)fx(0)m n(6.43)i经其中(0) 是任意向量，ai是待定常数，它们由、X(0) 与汀2，2 1关于质量矩阵M丨的正交条 件确定。由此，应取aj =i T M iX(0) p i = 1,2, ,m -1(6.44)由于"X(0) 与前m-1阶特征向量关于 M】正交，则用它作为初始向量进行逆迭代可指望求得第m阶特征对m,m?。为

27、了防止计算中的舍入误差而使1X(0)1偏离正交方向，还应在每一步开始都对X叫进行正交化处理。以上迭代过程称为克莱姆史密特(Gram-Schmidt)正交化。例6.1 图6.1示一三层刚架结构，各层的楼面质量分别为m = 180t, m? = 270t, m3 = 360t ；各层的侧移刚度分别为 匕=98MN.m,k2 = 196MNm, k3 = 294MNm。求刚架的固有频率和型。解该刚架的质量矩阵和刚度矩阵为M 】=180汉103-100101.50-002_一1-101k =98心06-13-2-0-25我们先用行列式法直接求这一问题的特征值。由式(6.25)1110 1K 】-o2

28、M 1=98X106-13-1.5®-20-25 - 2® 一图6.1甲=3式中18098 10令式(a)矩阵行列式展开、化简并令其等于零、得到三次方程，解之得0.3514,% -1.607,加3.542从而得到、=13.831 s, 2 = 29.61 s,3 =43.91 s而由式(6.23)可得三个主振型向量，用振型矩阵表示为,.0000 1.0000 1.0000 b】=0.6486-0.6070-2.542'0.3018-0.67972.440 J现在用式(6.38)的逆迭代法来解同一问题。为以后迭代计算方便先求出刚度矩阵采用初始迭代向量115'2

29、21221.0、= d.0，1.01.0则 &(0) = M Xx(0) =180 X103* 1.52.0当k=1时22.5；x匚 K JV(0180588 "03*16.5 »9Y九M汶（叭=型2 22.524.75 卜 588I 18 1.0122.5 16.5 9.5 冷8810= 197.9797920225122.5 16.5 9! 24.75 180丨18 JY二 0.685725 0.754298 0.548580T以下各次循环做法相同，现将结果载入下表6.1内。得到从上表看出，"k）比Y（k）收敛得快些，经过5次迭代就已经达到了要求，&q

30、uot;k）； -12。因此由0.742622'p j1.000000'0.4816550.6485870.7426220.224191 ,0301891,% = 宀1 1x（i 1）:T Y（l 1） ?2可见第一振型的结果与前面用行列式法求得的结果相符合。表6.1kX（k）Y（k）p(k)|p(k)_ p(2)V(k)p(k)122.522.5197.97980.68572816.5924.75180.7542980.548580212.411612.4116191.63930.03308530.7301558.2972812.44590.7321733.977217.95

31、4420.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.74262

32、28.2225312.33380.7224833.827247.654480.4483816.8用逐步积分法求结构动力响应6.8.1 概述求结构的动力响应，在数学上就是要求出运动方程(6.7)的解答。式(6.7)是一个二阶常系数微分方程组，可以用数值积分的方法对方程直接求解，即按时间增量.：t逐步求解运动微分方程，直至反应终了，这一方法称作逐步积分法，由于在数值积分前不对运动方程进行变换，所以又称作直接积分法。逐步积分法既 可用于求解线性结构体系问题一一在整个动力反应过程中K 1，M 1，C 1矩阵保持不变的问题；也可用于求解非线性结构体系的问题一一K 1，M 1, C 1矩阵随动力反应的过程

33、而变化的问题。此处，我们只讨论线性结构体系的问题。逐步积分法求解运动微分方程的基本思路是：1 把连续的时间过程离散为12,tn有限个点，对于运动微分方程MJ C；J K、 =F只要求它们在上述每个时间离散点上得到满足，也就是说，最终求解得到的只是位移、速度和加速度在有限个时间离散点上的值，而不是连续函数<(t)?'(t):z (t):?°2 在每个时间间隔 氏内，假定位移、速度和加速度符合某一简单的关系。而氏的选择，要求保证计算的稳定性与精确性。从这样一个基本思路出发，形成了多种逐步积分的方法，如中心差分法；线性加速度法；Wils on法;Newmark法等等。逐步积分

34、法的优点在于，计算的精确度可随氏的减少而提高，当大型结构体系受到历时较短的脉冲荷载作用时，由于要考虑高阶振型的反应使计算振型的数目较多，而逐步积分法求解动力反 应时，可省去特征方程的求解工作，在经济上、效率上可能是有利的。但是，逐步积分法的缺点，在于它 不能同时给出结构的振型和自振频率，而是直接给出结构的反应过程，并且当反应的时间历程较长，厶t取的较短时，计算工作量也是相当可观的。6.8.2 Wilson-二法Wils on-二法实质上是线性加速度法的推广，线性加速度法假定从时刻 t到r =t的时间内加速度是线 性变化的。而 Wilson-二法则假定在时间区间t到i ：A内加速度是线性变化的，

35、如图 6.1，其中r _ 1.0 o 当二-1.0时，Wils on -寸法法就退化为线性加速度法。 Wils on证明，要达到无条件稳定，必须选用二_ 1.37 , 通常采用，=1.4 o图 6.2W ilson-r 法按Wilson-二法的假定，在时间区间t到t . ：t内加速度的数学表达式为1(t .)*=(t)('(t”：t) 一、( t)J(6.45)t式中，为时间的增量，0 _ _将式(6.45)积分，得速度表达式2(t )彳(t) (t)，(小)(6.46)和位移表达式3(6.47)L(t ) L(t)k *L(t).2 . .t；现在的问题是如何根据t时刻的位移、速度和

36、加速度和上述位移、速度和加速度关系式以及运动方程来求tt时刻的位移、速度和加速度。将 -t代入式(6.46)及(6.47)，得到t . ：t时刻的速度、位移计算式、(t * n ：t) =、(t)丄 n ：t(、(t u ：t)、(t)(6.48)2k彳a aa n22 ' (t n ：t) =、(t) rt、(t) 寸 2，t2(、.(t r .：t)2 ' (t)(6.49)6选择、；(t * "t)为基本未知量。由式(6.49)和(6.48)，可以得到用、：(t"t)表示的、：(t ut)及、(t珂)的算式，即(t F)*(t ”t) _、(t) -走

37、、(t) -2、(t)(6.50)l (t F)几三(、(t n ：t) -f (t) -2f (t) t (t)(6.51)t2t 时刻的运动方程是M (t u ：t)C、(t u ：t) .K、.(t n ：t)二F(t rt)(6.52)F(t Tt)二F(t)班F(t ：t)-F(t)将(6.50), (6.51)代入(6.52)式，得到只有一个未知量 、(t”t)的方程组K?、(t n ：t) =F?(tn t)GQ式中Mzkm贏©(6.55) R(trt)=F(t)T(F(t：t)-F(t)M(、(t)、(t) 2、.(t)(Mt)uAt3-6At叫(t)"2&

38、#174;求解(6.54)式，可以得到、(tn：t)。再将它代入式(6.51)和(6.50)得到、(t nt)和、(t n ：t)将.=-：t 代入式(6.45)，(6.46)和(6.47)。得到、(t ：t)、(t t)和、(t t)为' (t ：t) =、(t)丄(、(t F) -f (t)"" A + (t) =、(t) 、(t).：t 亍(、(t rt) -、(t)12 At2、(t ：t)彳(t) 违、(t)(t) 芦(、(t rt) -、)为编程方便，可以进一步简化运算过程。将式(6.50)代入(6.57)得-63"、(t 切 2 N (t r

39、t) -f.(t)r- (t)(t)6l2At2B2At2日由式(6.57)得1 (、(t"t) -、(t) * (t t) -、(t)将上式代入(6.58)和(6.59)式得''上t、(t "5 列沁P (t)2At、(t M 二、(t)：t、(t)(、(t ：t)2、(t)6(6.53)(6.54)(6.56)(6.57)(6.58)(6.59)(6.60)(6.61)(6.62)Wils on -二法的计算步骤归纳如下:1初始计算(1)形成总刚度矩阵 K 和总质量矩阵 M 1,而总阻尼矩阵C二M 订K不存贮。(2)形成初始值 f (0), f (0)，

40、并计算(0);。(3)选取时间步长 比和B (般取日=1.4 )，计算积分常数。得“冷,卄芥彳七1,*n ：t:- o:- 2,3. ：t，一：打，J5,氏 =1 一 一，一:辽：2丁6J72(4)形成有效刚度矩阵($)对K进行三角分解K?二K： 0M： dCK? =LtDL2 对每一时间步长，循环计算(1) 计算t t时刻的有效载荷向量F?(t rt)二F(t) d(F(t rt) -F(t) M(： °、(t) ： 2、(t) 2、.(t)C(：i、(t) 2、(t)：3)(2)求解tc.-t时刻的位移向量'；(t ”t)L TDL、(t 八切二F?(t r ：t)(3)

41、计算t At时刻的加速度、速度和位移向量' (t：t)4(、(t rt) -、(t)： 5( (t) ： 6卜(t)、(t ：t)二 ' ： 7(、(t 讥)、(t)p(t 十At) =p(t)十甜卩 +叫(死 +甜) +2,) 例6.2有一简单的结构系统，不考虑阻尼影响，它的运动方程是2 0"&(说6-26(t)010试用Wilson- v法计算此结构系统从 用。1$2(t)i二24鸟住)|0到3.5的位移响应。已知结构初始状态是静止的，受突加不变载荷作J0J解：选择时间步长 氏。为保证计算精度，通常取 4 <0.1，其中T为结构系统自由振动的最小周期

42、。T对于复杂结构系统，一般高频部分对结构动力响应的贡献不大，可取低频部分的最小周期。本例的最小周 期T =2.8s，取 4 -0.28s,计算12个时间步长。(0)'0 10t。将1(0) ； = 0, L(0)；= 0初始条件代入运动方程，求得加速度初始条件取v -1.4，计算Wilson -31法的常数，得:0 =39.0-0.196= 7.65：3：'6= 1.14-27.9= 0.14:2 =15.3:-5 二-10.9：8 =0.0131形成有效刚度矩阵斤二K： °M239.0'.10j -243.0 一对于每一个时间步长，需先求解!?、(t n t

43、)二诚(tn“i 05 T 。1式中，卡?(t +日盘)，+(39.06(t)+15.36(t)+2§(t)然后计算加速度、速度|WJ0J卫仃和位移T(t：t)27.9(、(t u ：t) _、(t) -10.9、. (t) 1.14、(t)(t：t) ；=、(t)0.14(、(tt) 、(t)(t：t)匚、(t)0.2&、(t)0.0131(、(t：t)2(、(t)12个时间步长按以上诸式计算的结果如下表所示:tAt2也t3 At4也t5也t6 AtMt8也t9也t10也t11也t12也t0.006050.05250.1960.4900.9521.542.162.672.9

44、22.822.331.54620.3661.342.643.924.885.315.184.613.823.062.522.29法Newmark法也是一种线性加速度法的推广。它采用如下两个基本假定:、(t t)二、(t)(1 - )、(t)、(t ：t) ：t(6.63)、(t ：t) =、(t) (t) ：t (1-2 J、(t)2 (t ：t)：t2(6.64)其中一：和是参数，根据积分的精度和稳定性要求确定。1 1 1 1当韶=2和时，关系式(6.63)、(6.64)相当于线性加速度法。当吋二鼻 和？ =4时，相当于常平 均加速度法。已经证明，当_0.5和戸二0.25(0.5 )2时，N

45、ewmark法是无条件稳定的。根据上述两个假定，以及t 氏时刻的运动方程，就可以求得时刻的位移、速度和加速度的解。将 ( ：t)作为基本未知量，由式(6.64)求得(t 切二(、(t ：t) -、(t)(t) -(尹-1)b (t)(6.65)将式(6.65)代入式(6.63)得y丫 丫(t M 二- t(、(t：t)-、(t)(1 -、(t)(1 - 2W (6.66)将式(6.65)和(6.66)代入(t：t )时刻的运动方程M、(t ：t) C、(t ：t) K、(t：t) = F(rt)(6.67)得到只有一个基本未知量 '")的方程式中? 1eKrM fc(6.69

46、)(6.70)F?(t：t)二F(t ：t) M( +、(t) (t)(右-1)b (t)yYYC:，tf (t) G -1)b (t)(2： -1)讥、(t)解方程(6.68)得到、(t . ：t)后，再代入式(6.65)和(6.63)求出、(t t)和、(t t)。为了便于编写计算机程序，将Newmark法的计算步骤综述如下:1初始计算形成刚度矩阵K ,质量矩阵 M 1,阻尼矩阵C=:MK不存贮。(2)形成初始向量r (0)、(o)，计算:(0)。(3)选择时间步长:t和参数-,。计算积分常数，得：01厂1：t2-2),>6 二：t(1 一 )，： 7 二:t形成有效刚度矩阵K =K

47、 ： °M ： JC(5)对 K 进行三角分解2 对每一个时间步长循环(1)计算在时刻t ： =t的有效载荷向量F(t：t)二F(t：t) M(：。、.(切：2( (t) ： 3、(t)C(： £ (t) ： 4(t) ： 5' (t)(2) 求解t t时刻的位移向量LTDL、(t P二?t 沏(3) 计算t 氏时刻的加速度和速度' (t：t)二：0(、(t：t) - ' (t) -： 2、(t) -： 3、(t)、(t ：t)壬(t)： 6(t) ： 7' (t：t)可以看出，Wilson-二法和Neumark法的计算步骤是相同的，计算公式

48、的形式也是一样的，只是系数有所不同，故可以在一个有限元程序中同时存在这两种直接积分法供用户选择。例6.3用Newmark计算例题6.2的结构系统的响应。取0.25，= 0.5。解：选择.Vt =0.28s。将1(0)二0,1(0)丄0代入运动方程，得 1(0)1 0,1OT。计算纽马克法的常数：0= 51.0:7.14 工2= 14.3：'3= 1.00： 4 = 1.00 :- 5=0.00：'6= 0.14：7 =0.14形成有效刚度矩阵为从 62 4251.00 _ 108 -21 _|1-255对于每一个时间步长，先求解下列方程X'2 O-K?1 (t ：t)；R(t ：t)?0 (51.0】(t* 14.3 ($ 1.°'：)1：(t ：t) 151.0(1 (t ：t)T(t)：) 一14.3： (t)1.0l (t)1(t =t) 、(tV 0.14 (tV 0.14 (t =t)'用以上诸式，按12个时间步长依次计算的结果如下时间2