wilson法和newmark法的理论过程要点

1、第三章离散化结构动力方程的解法（2013.4.24）立1绪言对于一个实际结构，由有限元法离散化处理后， 应用瞬时最 小势能原理可导出动力方程mku + clu + kku = f（t）（3.1）这里，；u：，、iu：，、iu：，及【f（t”分别表示加速度、速度、位移及所 作用的外力矢量，他们都是与时间有关的。从数学的角度来看，式（3.1）是一个常系数的二阶线性常 微分方程组，对于它的求解原则上弁无困难。 但是，由于m、c 和仅的阶数非常高，使得式（3.1）的求解必须花费很大的代价， 便促使人们去寻求一些效率高的近似计算方法。目前，用于求解式（3.1）的方法，大致可分为两大类。一是坐标变换法，它

2、是对结构动力方程式（3.1）,在求解之前，进行模态坐标变换，实际上就是一种ritz变换，即把原物理 空间的动力方程变换到模态空间中去求解。现在，普遍使用的方法是模态（振型）迭加法，即用结构的前q阶实际主模态集（主振型阵）构成坐标变换阵进行变换。通过这一变换，实现降阶， 求较好的近似解，而且，还用解除耦合的办法，简化方程的计算。 还有一种所谓假设模态法，即是用一组假设模态，构成模态坐标变换阵进行变换，获得一组降阶的而不解耦的模态基坐标方程。 显然，这种方法的计算精度，取决于所假设的模态。用ritz矢量 法求解的近似模态作为假设模态，可得到满足要求的精度。二是直接积分法，它是对式（3.1）在求解之

3、前，不进行坐 标变换，直接进行数值积分计算。这种方法的特点是对时域进行离散，将式（3.1）分为各离散时刻的方程，然后，将该时刻的加速度和速度用相邻时刻的各位移线性组合而成， 于是，式（3.1） 就化为一个由位移组成的该离散时刻上的响应值， 通常又称为逐 步积分法。线性代数方程组的解法与静力时刻的位移来线性组合，就导致了各种不同的方法。主要有中央差分法，houbolt方法，wilson-8法和newmark方法等。3.2模态（振型）迭加法设有n个自由度的系统，在外力 作（毋的作用下，常常被激 起较低阶的一部分模态（即振型），而绝大部分高阶模态被激起 的分量很小，一般可忽略不计。例如，在地震载荷作

4、用下，通常, 只有最低的二阶，三阶模态起主要作用。所以，对于这样的一些 问题，采用模态迭加法是有效的。设有式（3.1）的n阶动力方程，起主要作用的是其前 q阶 模态，通常取q= no按ritz变换，则可将式（3.1）中的用前 q个模态的线性组合来表示，即qu二 yl 1 丫2 2yq q- jyj j=y（3.2）其中，附为结构的已知的保留主模态矩阵， 而yqx是维的模 态基坐标矢量，它形成了一个q维的模态空间。它表示在y中, 各阶主模态所占有的成分的多少。假定户已用第二章所述的某一方法解出，再将式（3.2）代 入（3.1）,弁左乘以，可得*一”1 * *my cy ky:f（3.3）式中_*

5、_ _t _ _m = tm*_ _t _ k = tk_*.thc = tc f = tf显然，式（3.3）是一个q阶的微分方程组。由于q n ,所以， 它比式（3.1）的n阶就小的多了，实现了降阶，因而也就容易 求解多了。一,、 . . * 一- , ，,、一 一 、 、 一 .若展开上述的m的表达式，根据主模态（主振型）关于m 的表达式，根据主模态的（主振型）关于m的正交性质，可知*mj = 0（ = j）所以，m是一个对角阵。同理可知k也是一个对角阵。然 而，在一般的情况下，c*是一个非对角阵，即在模态空间中， 系统的的阻尼一般是耦合的。因此，式（3.3）是一个完全解耦的动力学方程。但

6、是，它是一个已降阶的q阶的动力方程，可使用后面即将介绍的直接积分法求解。当系统的阻尼为比例阻尼时，即c可以表示为 _* n*c = m k（3.4）则c为对角阵。此外，若系统的阻尼是一般的的线性阻尼，弁 非比例阻尼，但是只要结构的固有频率不相等， 而且不十分接近， 则可用舍去c阵中的非对角元来实现c的对角阵，也不会引 起太大的误差。在上述两种情况下，可以获得对于模态坐标的完全解耦的动 力学方程。即式（3.3）是q个独立的方程，每个方程只包含一 个未知量，相互之间不耦合。因而式（3.3）可按单自由度的动力学方程写为* -* ,mii yi cii y kii yi = fi (t)(i = 1,

7、2,.q)(3.5)或2 yi+2 i i yi+ i y = l (t)(i = i,2,.q)(3.6)其中23。i = cii / mii , fi (t) = fi (t) / mii 。式(3.6)可用直接积分法计算，或用duhamel积分求得其解为1 t二yi(t)fi ( )e i i sin i di 0(3.7)-e- i yai sin丁i t bi cos-i t( i = 1,2,. q)式中，ei =切i ji 2，而ai , b由初始条件y0 = (m*)1 tmu0(3.8)y0=(m*)t： tmu。得出的yi0v yi0决定。由于有阻尼的存在，由初始条件所激发

8、的振动，随时间的增长而衰减以致消失。因此，常可不计式(3.7)中的第二项，即是由初始条件激发的自由衰减振动。计算出yi(t)后，便可利用式(3.2),计算出物理坐标的响应u(t)。数学计算步骤可归纳如下：第一步：根据结构的离散化模型，建立系统的 m,k以及f(t),弁进行结构的固有特性分析，即求解特征值问题(k- 2m) =。求出前q阶特征对(包,电) , (i q)第二步：形成模态阵n=电曲.%,弁建立模态基坐 标下的动力方程.v、2；-% jy =力仕)(i = 1,2,., q)其中fi(t)=mtf。,而n =曲丁叫说。根据实验结果或经验 m数据确定各阶主振动中的比例阻尼1。第三步：求

9、解主模态基坐标的动力方程，有tvi(t)=i。fid)e-sin7i(i)dt ,其中,io77i = i j- 2i。第四步：进行坐标变换后，求得动力响应u ；= ：t：y)3.3模态假设法上节所述的模态迭加法，是用系统的真实主模态组成的模态 矩阵，再对系统的物理坐标进行模态坐标变换，从而在主模态空间中得到降阶弁解耦的动力学方程，这样来实现简化计算。而这里提出的假设模态法，则是用一组假设模态矩阵， 对系统的物理 坐标进行模态坐标转换，从而在模态空间中得到一组只降阶的动 力学方程。若令假设模态矩阵为1n，而mn,进行坐标变换，即x t =门匕2 t :m1(3.10)把它代入式(3.1),弁左

10、乘阵，则可得到降阶的动力学方程为(3.11)i m ： q： ic i q ik i q - q t )其中im j = m m,lk* = tk悔,c* = 6tc ,q（t）=6（t）。它们分别对应于假设模态坐标iq）的质量矩阵、阻尼矩阵、刚度矩阵与广义力列阵。因为矩阵中的各列都是假设模态，它们 一般不具有正交性，所以，m 、ic 和k 都不是对角阵。于是，方程（3.11）是不能解耦的方程组，但它却是比式（3.1）的阶数要低得多了。显然，对式（3.11）采用直接积分法求解，将比对式 （3.1）求解要简便得多。这是假设模态法的优点。假设模态法的计算精度，很显然地是取决于假设模态阵中模 态假设

11、的好坏与质量。因此，应用假设模态法能否成功的关键在 于确定出一个适宜的假设模态矩阵。在第五章中，我们介绍了几种构造假设模态的方法。实际上，在2.9中介绍的rayleigh-ritz分析，可认为是一种 假设模态法。它的作用，在于降低方程的阶数，简化计算。它的 基本思想是，事先假定出若干近似的特征矢量，然后按照这些特征矢量的最佳线性组合，而算得前若干阶特征值的近似值。显然,运用这种方法时，其计算精度与事先假定的特征矢量的近似程度 和数量有关。按照ritz变换的思想，找到了近似的特征矢量xj后，（i=0,q）,即有：）=可：xi a? x2 ；aq xq x h al（3.12）求解如下的广义特征值

12、问题，即.* ,* , 一cik a= pm 】a,p = 62（3.13）其中*tik = x ik hx *tim 1 = ix i im ii x 1k和m为原结构离散化之刚度阵和质量阵，它们都是n阶方阵。求解式（3.13）,得到q个特征矢量，有i a1 ? = | a11a211 aq. a2= = |a122a22 t aq*”=| a2q aqq再按照ritz的变换，即式（3.12）,由特征矢量与,可计算出矢 量&,他，%即是q1i ；二，aji lx，i=1,2, ,q（3.14）j i现在用近h =地包来表示此变换阵，它就是我们要构造的假设模态矩阵。立4中心差分法（显示法）现在

13、开始讨论直接积分法，或称逐步积分法。前面讨论的模态迭加法，弁非总是有效的。当刚度矩阵ik,或质量矩阵im】，或阻尼矩阵匕】出现随时间变化时，或当外荷载激 起的振型太多，需要计算的特征对太大时， 就不宜于采用模态迭 加法，在这些情况下，采用逐步积分法是适宜的。中心差分法就 是其中的一种。这种方法的特点，是将动力方程在时间域上离散， 化成对时间的差分格式，然后根据初始条件，利用直接积分法逐 步求解出一系列时刻上的响应值。假定t=0时，位移、速度和加速度分别为已知的u0, u0和u；。再将求解的时间区间划分为n个等分，即削=工。我们要建立的积n分格式就是从已知的0,稣,2&,，t的解来计算下一个时间

14、 步的解。在中心差分法中，是按中心差分将速度和加速度矢量离散化 为.1. . ./、ut = yy(ut+at - ut_at )(3.15)1 ut =彳(uts - 2ut +u“)(3.16)于是上面二式，就将t时刻的速度和加速度用相邻时刻的位移来 表示了。考虑在t时刻的动力方程，有m ut + cut + kut = f t( 3.17)将式(3.15)和(3.16)代入式(3.17)中，得到fm 21t5 ut t=f 的佯皿ut-点叱点5 ut ta. .18)这样，上式就化为用相邻时刻的位移表示的代数方程组。由它可解出ut地。又由于它是利用t时刻的方程解得ut处的，所以，它 称为

15、显7k积分。弁且，还注意到，在求解uta时，需要用ut ,ut泊的值。于是，在计算开始时，即t=0时，要计算ua的值、 就需要u_t的值，他是未知的，因此，必须有一个启动的处理， 因而这种算法不是自起步的。由于u0,u0和u0是已知的，所以，由t=0时的式(3.15)和(3.16),可解得t2u_& = u0 - ahu。+ -2-uo(3.19)使用中心差分法的逐步求解过程如下：b. 初始计算(1)形成刚度矩阵k,质量矩阵m和阻尼矩阵c。(2)给定初始值u0,u0和u0。(3)选择时间步长at, 与&-弁计算积分常数：11c1a0 = , 2 , a1 = _ , , a2 = 2a0 ,

16、a3 =.,:t22：ta2(4)计算u0 =uo &uo+a3uo。(5)形成有效质量矩阵m = a0m十ajc。(6)三角分解m： m =ldlt。c. 对每个时间步计算(1)计算t时刻的有效载荷一 一 一 一 一 一 一一ft = ft -(k - a2m) ut - (a0m -ajc) uto(2)求解t +组时刻的位移tldltu/=ft。(3)如果需要计算t时刻的速度和加速度ut 二 a。ut t 2ut utut =a ut t u1应当指出，这种中央差分算法，左端的系数矩阵只与质量阵m和阻尼阵c有关，而与刚度阵k无关。如果质量阵和阻尼阵 是对角阵，那么在解方程时，就不需要对系

17、数阵进行三角分解， 即不需要解线性代数方程组，从第一步开始逐次直接求得各个时 刻ut也的值，这是中央差分格式就是一种显示的格式。此外，由于不求解代数方程组，也就不需要进行组集，它的右端项的形成 也只须在单元一级水平上，由每个单元对有效载荷矢量的贡献迭 加而成。因此，adina 程序规定，在用中心差分法时，必须使 用对角的质量阵和阻尼阵。 从计算稳定性角度来看，中心差分法 的缺点，在于它是条件稳定的，即当时间步长与太大 时，积分是不稳定的。所以，对步长的限制是t t tcr n 31这里，是临界步长值，tn是有限元系统的最小周期。这样， 当tn很小时，就限制了 &必须很小，所以求解所花的代价就很

18、大。立5线性加速度法和wilson-。法线性加速度法和wilson - 0法，都是属于逐步积分法。线性加 速度法是假定在t,t十4时间间隔内，即在步长m时间内，加速度 u（t+）呈线性变化，其表达式为ut+z = ut+ha（3.20）其中，a=（u匹-u）/at。但是，这个方法不是无条件稳定的， 所以在应用上受到限制。70年代初期，wilson推广了线性加速 度法，他假定在此步长 &更大的时间区间（t,t+e&）内，加速度仍 保持线性变化，经过证明，当 日之1.37时，这一方法是无条件稳定 的，这就是wilson - 0方法。这个方法的加速度表达式为ut =utai（3.21）式中aj= （

19、utrt ut）q t显然，对比式（3.20）和式（3.21）得知，线性加速度法是 wilson10-e法中，当8=1时的一个特例。所以，我们只讨论wilson -。法就够了。在0日与区间内，对式(3.21)进行积分，得到 2ut+t = ut 十u尸十 2(ut出-ut)(3.22)和12. 3u = ut+u户+彳ut 七 十7?(ut侬ut) (3.23) 261 t令7=阳，由上二式，有tut 梯=ut 十(ut+ait +ut)(3.24)和2 t2ut祕=ut 十日”ut + (ut他 + 2ut) (3.25)6从这二式，可将(t+e&)时刻的加速度和速度用位移来表示即66ut

20、琦=7tr(ut 期-ut)-ut-2ut(3.26)和，一3一一 -一 一- t,ut ,t =-t(ut -ut)-2ut-2-ut(3.27)于是，在t+e彳时刻的动力方程为mut5叱 t kut t =ft t(3.28)式中，ft t =ft 八ft t ft)将(3.26)和式(3.27)代入式(3.28),就得到关于u他的方程为63(k -i-t2m 何c1)ut-t =ft mf】t t-fj 66m(ut ut 2ut)(3.29)t t 1t 3tc(商ut 2ut /uj记k=k舟mc6.6.f-t =f(f _ft) m(hut:ut t t tj- -6 . . .

21、t.2ut) c(ut 2ut -ut) t2于是，式(3.29)可写为kut t =ft =t(3.30)求解方程(3.30),则得到ut他将求解得到的ut他,代入(3.26)中，就得到ut他。如在(3.21) 中，取弁将式(3.26)代入，有663ut t= 2(ut-ut) f ut (1-)ut(3.31)t t -t t将(3.21)代入式(3.22)和(3.23),并取工=t ,有tut - =ut (ut , ut)(3.32)2l tut t =uttut -6(ut t 2ut)(3.33)用wilson -9法逐步求解的过程如下：a .初始计算(1)形成刚度矩阵k,质量矩阵

22、m和阻尼矩阵c。(2)给出初始值uo, uo和uo。(3)选择时间步长at,取日=1.4,弁计算积分常数，ait：-.ta3 =t二1一3ta7 =2t2 6(4)形成有效刚度矩阵k* ： -k* =k +a0 m +a c q19(5)对k*作三角分解：k*卜ldhl】tb.对每个时间步计算(1)计算t +&时刻的有效载荷1f”f 与、-if% 1m 1 ao ：u)t a2：u)t 2：u tc 1 a1 ：u)t 2 ：u)t a3 ：u)t(2)计算t+8at时刻的位移l11d1l 屋u；_t二g(3)计算t+&时刻的位移，速度和加速度u t t =a4 u：t t - u：t a5

23、u：t a6 ：u：t；u；l：u)t a7 ,：u；ut,：u：；t ,t = ：u tt：u：t a8 ut . j 2，：u ；t与中心差分法相比较，wilson-日法是隐式积分，即每计算步，必须解一个线性代数方程组。当 日1.37时，它是无条件稳 定的。止匕外，这种算法是自起步的，t+&时刻的位移，速度和加 速度都可由t时刻的变量表示，不需要特别的起动处理。6 newmark 方法newmark在1959年提出的逐步积分格式， 故称为newmark 方法。它的基本假定是其中a和6是按积分的精度和稳定性要求可以调整的参数。当6 = 1, a=1时，它就是线性加速度法，所以，newmark

24、方法也 26可以理解为线性加速度法的一个小延伸。newmark法最初提出作为无条件稳定的一种积分格式是常平均加速度法，即假定从t到t+&时刻，加速度不变，取为常数工（总+%）。此时，取6,2 tt -t/20f =1 o常平均加速度法是应用得最广泛的逐步积分方法之一。4研究表明，当6 25,支之0.25（0.5+4时，newmark方法是无条件 稳定的。从式（3.34）和（3.35）可得到同小丹他用ut也及5、ut 和gt表示的表达式，即有），+=4（ u . u）ta u 巾 fu t （3.36） t t 0 i.m i - 1 i.c 1(5)对k作三角分解:；k=ilid *b.对每个

25、时间步计算(1)计算t+&时刻的有效载荷* f .t = f ,1ml 二 0 ： u 二 2 ： u：t i ： 3 u- i.c i ? 1 u u/t 上，4 t u：t，工5 t u t(2)求解t + &时刻的加速度和速度：uu)- ： u)-二 2： u -二 3 ： ut：a t 0t：. tt 2 t 3 t:u： 一 : u 6 u) ：u?t tt 6 t 7t , :t我们注意到 wilson-e法与newmark法的计算关系式，在形式上是相同的，只是其中的系数取不同的值而已。因此，它们可用同一计算机程序来实现。园7 houbolt方法这个差分格式是利用t+及、t、t-及

26、、t-2at四个时刻上位移 的三次插值多项式建立起来的。即假定1ut 也=tl(2ut 也-5ut +4ut-a-ut-24)(3.40)和1ut 也=61t(11ut 也-18ut +9uq-2ut-24)(3.41)这里认为ua，必和ut是已知的，而ut也是未知的。考 虑t+ &时刻的动力方程，有mut帖十cut也+kut站=ft也(3.42)将式(3.40)和(3.41)代入式(3.42)中，就得到求解ut+及时刻的方程为2 211 -m+c+k)ut地at6at/=什0+&网+京5加-$网+白5加、(3.43)1 1 一 . -2m c ut-2 t t 3 t由上式解得ut把后，代入式(3.40