wilson法和newmark法地理论的过程

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

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

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

4、主要作用。所以，对于这样的一些问题，采用模态迭加法是有效的。设有式（3.1）的n阶动力方程，起主要作用的是其前q阶模态，通常取q=n。按Ritz变换，则可将式（3.1）中的“用前q个模态的线性组合来表示，即u=Y他+Yb+.+Y他=£bY1122qqjjj=1=Y（3.2）其中，为结构的已知的保留主模态矩阵，而丫是维的模nxqqx1态基坐标矢量，它形成了一个q维的模态空间。它表示在丫中，各阶主模态所占有的成分的多少。假定已用第二章所述的某一方法解出，再将式（3.2）代入（3.1）,并左乘以t，可得M叫舞+C*Y+K*Y=F*（3.3）式中M*=tMK*=tKC*=tCF卜=tF显然，

5、式（3.3）是一个q阶的微分方程组。由于q<n，所以，它比式（3.1）的n阶就小的多了，实现了降阶，因而也就容易求解多了。若展开上述的M*的表达式，根据主模态（主振型）关于M的表达式，根据主模态的（主振型）关于M的正交性质，可知m*=0（i主j）ij所以，M*是一个对角阵。同理可知K*也是一个对角阵。然而，在一般的情况下，C*是一个非对角阵，即在模态空间中，系统的的阻尼一般是耦合的。因此，式（3.3）是一个完全解耦的动力学方程。但是，它是一个已降阶的q阶的动力方程，可使用后面即将介绍的直接积分法求解。当系统的阻尼为比例阻尼时，即C可以表示为C*=aM*+PK*（3.4）则C*为对角阵。此

6、外，若系统的阻尼是一般的的线性阻尼，并非比例阻尼，但是只要结构的固有频率不相等，而且不十分接近，则可用舍去C*阵中的非对角元来实现C*的对角阵，也不会引起太大的误差。在上述两种情况下，可以获得对于模态坐标的完全解耦的动力学方程。即式（3.3）是q个独立的方程，每个方程只包含一个未知量，相互之间不耦合。因而式（3.3）可按单自由度的动力学方程写为m*輿+c*&+k*y=F*(t)(i=1,2,q)(3.5)iiiiiiiiii或y&+2o&+®2y=f*(t)(i=1,2,q)(3.6)iiiiiii其中20=c*/m*,f(t)=F*(t)/m*。式(3.6)

7、可用直接积iiiiiiiiii分法计算，或用Duhamel积分求得其解为1-y(t)=iftf(t)e-訥(t-t)sinOdx+0'1(3.7)e-%o/asinO.t+bcosO.t(i=1,2,q)iiiiii式中，O.=叫、,1_勺2，而a,b由初始条件y=(M*)_1TMu00(3.8&=(M*)-itMU000得出的yio与&o决定。由于有阻尼的存在，由初始条件所激发的振动，随时间的增长而衰减以致消失。因此，常可不计式(3.7)中的第二项，即是由初始条件激发的自由衰减振动。计算出y(t)后，便可利用式(3.2)，计算出物理坐标的响应iu(t)。数学计算步骤可

8、归纳如下:第一步：根据结构的离散化模型，建立系统的M,K以及F(t)，并进行结构的固有特性分析，即求解特征值问题(K-32M)©=0求出前q阶特征对(,帥),(i二1,2,.,q)ii第二步：形成模态阵二©©.©，并建立模态基坐nxq12q标下的动力方程少+疋3y+32y=f(t)(i=1,2,.,q)iiiiiii其中f(t)二丄0tF(t)，而m二©TM(H。根据实验结果或经验imiiiii数据确定各阶主振动中的比例阻尼g。i第三步：求解主模态基坐标的动力方程，有1y(t)=Jtf(T)ePi(t-T)sinB(t-t)dz，其中，iB0i

9、iiB=3Ji-g2。iii第四步：进行坐标变换后，求得动力响应u=y§3.3模态假设法上节所述的模态迭加法，是用系统的真实主模态组成的模态矩阵，再对系统的物理坐标进行模态坐标变换，从而在主模态空间中得到降阶并解耦的动力学方程，这样来实现简化计算。而这里提出的假设模态法，则是用一组假设模态矩阵，对系统的物理坐标进行模态坐标转换，从而在模态空间中得到一组只降阶的动力学方程。若令假设模态矩阵为，而m<<n，进行坐标变换，即nxm&Q)=q(t)(3.10)nx1nxmmx1把它代入式(3.1)，并左乘t，则可得到降阶的动力学方程为(3.11)M&+C&

10、+Kq=g(t)其中M=0M®,k=otkbc=C0'Q(t)=0tf(t)。它们分别对应于假设模态坐标q的质量矩阵、阻尼矩阵、刚度矩阵与广义力列阵。因为矩阵0中的各列都是假设模态，它们般不具有正交性，所以，m、c和k*都不是对角阵。于是，方程(3.11)是不能解耦的方程组，但它却是比式(3.1)的阶数要低得多了。显然，对式(3.11)采用直接积分法求解，将比对式(3.1)求解要简便得多。这是假设模态法的优点。假设模态法的计算精度，很显然地是取决于假设模态阵中模态假设的好坏与质量。因此，应用假设模态法能否成功的关键在于确定出个适宜的假设模态矩阵。在第五章中，我们介绍了几种构造

11、假设模态的方法。实际上，在§2.9中介绍的Rayleigh-Ritz分析，可认为是种假设模态法。它的作用，在于降低方程的阶数，简化计算。它的基本思想是，事先假定出若干近似的特征矢量，然后按照这些特征矢量的最佳线性组合，而算得前若干阶特征值的近似值。显然，运用这种方法时，其计算精度与事先假定的特征矢量的近似程度和数量有关。按照Ritz变换的思想，找到了近似的特征矢量X后，i(=0,q)，即有)=aX+aX+aX=XA(3.12)1122qq(3.13)求解如下的广义特征值问题，即KA=pMa,p=rn2其中k=xPkxM=xtmxk和m为原结构离散化之刚度阵和质量阵，它们都是n阶方阵。

12、求解式（3.13），得到q个特征矢量，有A=a1ai112A=a2a2a1qa2q<Ai=faqaqq12再按照Ritz的变换，即式(3.12)矢量e，1aqTq由特征矢量A，可计算出（），£，即是2qe=faixijij=1212精彩文档(3.14)现在用nxq=q好来表示此变换阵，它就是我们要构造的假设模态矩阵。§3.4中心差分法（显示法）现在开始讨论直接积分法，或称逐步积分法。前面讨论的模态迭加法，并非总是有效的。当刚度矩阵k，或质量矩阵m，或阻尼矩阵c出现随时间变化时，或当外荷载激起的振型太多，需要计算的特征对太大时，就不宜于采用模态迭加法，在这些情况下，采用

13、逐步积分法是适宜的。中心差分法就是其中的一种。这种方法的特点，是将动力方程在时间域上离散，化成对时间的差分格式，然后根据初始条件，利用直接积分法逐步求解出一系列时刻上的响应值。假定t=0时，位移、速度和加速度分别为已知的u,&和鈕。000再将求解的时间区间划分为n个等分，即At=T。我们要建立的积n分格式就是从已知的0，At，2At，t的解来计算下一个时间步的解。在中心差分法中，是按中心差分将速度和加速度矢量离散化为&=-（u-u）（3.15）t2Att+Att-At&-（ua-2u+ua）（3.16）tAt2t+Attt-At于是上面二式，就将t时刻的速度和加速度用相

14、邻时刻的位移来表示了。考虑在t时刻的动力方程，有M&+Ct&+Ku=F（3.17）tttt将式（3.15）和（3.16）代入式（3.17）中，得到吕M+丄CuA=F-fK-ALiMu-fALM-CuaIAt22At丿t+AttIAt2丿tI、At22At丿t-At（3.18）这样，上式就化为用相邻时刻的位移表示的代数方程组。由它可解出u。又由于它是利用t时刻的方程解得u的，所以，它t+Att+At称为显示积分。并且，还注意到，在求解u时，需要用u，t+Attu的值。于是，在计算开始时，即t二0时，要计算u的值、t-AtAt就需要u的值，他是未知的，因此，必须有一个启动的处理，-

15、At因而这种算法不是自起步的。由于u,谒和粵是已知的，所以，000由t二0时的式（3.15）和（3.16），可解得A12UA=U-At&&(3.19)-At0020使用中心差分法的逐步求解过程如下：A初始计算(1)形成刚度矩阵K，质量矩阵M和阻尼矩阵Co(2) 给定初始值u,U和Uo000(3) 选择时间步长At,At<At，并计算积分常数：cr1，1，1oa=a=a2a,a=°0At212At203a2(4) 计算uu-Atu+auo-At0030(5) 形成有效质量矩阵雨aM+aCo01(6) 二角分解M:应LDLToB对每个时间步计算(1) 计算t时刻的有

16、效载荷FF-(K-aM)u-(aM-aC)u。tt2t01t-At。(2) 求解t+At时刻的位移LDLTuF。t+Att(3) 如果需要计算t时刻的速度和加速度&=a(u-2u+u)t0t+Attt-AtU&a(U-U)t1t+Att-At应当指出，这种中央差分算法，左端的系数矩阵只与质量阵M和阻尼阵C有关，而与刚度阵K无关。如果质量阵和阻尼阵是对角阵，那么在解方程时，就不需要对系数阵进行二角分解，即不需要解线性代数方程组，从第一步开始逐次直接求得各个时刻u的值，这是中央差分格式就是一种显示的格式。此外，由t+At于不求解代数方程组，也就不需要进行组集，它的右端项的形成也只须

17、在单元一级水平上，由每个单元对有效载荷矢量的贡献迭加而成。因此，ADINA程序规定，在用中心差分法时，必须使用对角的质量阵和阻尼阵。从计算稳定性角度来看，中心差分法的缺点，在于它是条件稳定的，即当时间步长At太大时，积分是不稳定的。所以，对步长的限制是TAt<At=ncr冗这里，At是临界步长值，T是有限元系统的最小周期。这样，crn当T很小时，就限制了At必须很小，所以求解所花的代价就很大。n§3.5线性加速度法和Wilson-°法线性加速度法和Wilson-0法，都是属于逐步积分法。线性加速度法是假定在t,t+At时间间隔内，即在步长At时间内，加速度U&

18、t+t)呈线性变化，其表达式为&=&+TA(3.20)t+Tt其中，A二(U&U)/At。但是，这个方法不是无条件稳定的，t+T所以在应用上受到限制。70年代初期，Wilson推广了线性加速度法，他假定在此步长At更大的时间区间(t,t+0At)内，加速度仍保持线性变化，经过证明，当0、1.37时，这一方法是无条件稳定的，这就是Wilson-0方法。这个方法的加速度表达式为&=&+t-A(3.21)t+Tt1式中A=(&0-&)/0At1t+0Att显然，对比式(3.20)和式(3.21)得知，线性加速度法是Wilson-0法中，当9=1

19、时的一个特例。所以，我们只讨论Wilson-0法就够了。在0<i<0At区间内，对式(3.21)进行积分，得到T2&=&+&1+话(&-&)(3.22)t+ttt29Att+9Att和,6A(&t+9at-殍)3.23)3.24)3.25)1u=u+T+&*T2+t+Ttt2令T=0At，由上二式，有u&t+9At=u&t+芈(&a+&)2t+9Att92At2他舷=伦+9At&t+6(ut+6At+2&t)从这二式，可将(t+9At)时刻的加速度和速度用位移来表示即&=

20、qt(u-u)-：U-2&t+9At92At2t+9Att9Attt(3.26)和斶a=a1(u-u)-2&-呼&t+9At9Att+9Attt2t(3.27)于是，在t+eVt时刻的动力方程为Mu&+Cu&+Ku=F%t+aAtt+aAtt+aAtt+aAt(3.28)式中，F%t+脑=+9(Ft+At-Ft)将(3.26)和式(3.27)代入式(3.28)，就得到关于町的方程为t+9Vt(K+30AtC)叫低=F+0(F-F)tt+Att+M(Lu+&+2U)02At2tUAttt+C(蟲u+2&+0At&)0Attt2t(3

21、.29)记%=K+蟲M+島CF=F+0(F-F)+M(二u+16-&t+0Attt+0Att02At2t0Att+2娉)+CGu+2U+0At&)t0Attt2t于是，式(3.29)可写为%u=F(3.30)t+0Att+0At求解方程(3.30)，则得到ut+0At将求解得到的u，代入(3.26)中，就得到&。如在(3.21)t+0Att+0At中，取T=At，并将式(3.26)代入，有&=(u-u)+U+(1-3)&(331)t+At03At2t+0Att02Att0t将(3.21)代入式(3.22)和(3.23)，并取心t，有d=u+At(&

22、;+&)(3.32)t+Att2t+Attu=u+Att&+竺(娉+2娉)(3.33)t+Attt6t+Att用Wilson-0法逐步求解的过程如下：A初始计算(1) 形成刚度矩阵K，质量矩阵M和阻尼矩阵C。(2) 给出初始值u，&u和&u&。000(3) 选择时间步长At，取。二1.4，并计算积分常数，6,3，a=7a=a=2a，0(OAt)210At21OAt，a，a，a=5a=0-ja=2j324050,3AtAt2a=1,a=-ja=607286(4)形成有效刚度矩阵K*:K*=k+am+ac对k*作三角分解:k*=ldltB对每个时间步计算1)

23、2)3)计算t+At时刻的有效载荷圖=f+6(f-f)+也(au+a&+2/&)t+OAttt+Att0t2tt+cau+2&+a&)1tt3t计算t+OAt时刻的位移LDLTu=%t+OAtt+OAt计算t+At时刻的位移，速度和加速度核=a(u-u)+a品+a&t+At4t+OAtt5t6t&=4+a融+個)t+Att7t+Attu=u+At+a(4+2i&)t+Attt8t+Att与中心差分法相比较，Wilson-O法是隐式积分，即每计算一步，必须解一个线性代数方程组。当1.37时，它是无条件稳定的。此外，这种算法是自起步的，t+A

24、t时刻的位移，速度和加速度都可由t时刻的变量表示，不需要特别的起动处理。§3.6Newmark方法Newmark在1959年提出的逐步积分格式，故称为Newmark方法。它的基本假定是&=&+|"(1-5)倔+5倔Att+Atttt+Atu=u+rfAt+&+a&t+Attttt+AtA/2(3.34)(3.35)其中a和§是按积分的精度和稳定性要求可以调整的参数。5=1,a=1时，它就是线性加速度法，所以，Newmark方法也26可以理解为线性加速度法的一个小延伸。Newmark法最初提出作为无条件稳定的一种积分格式是常平均加速度

25、法，即假定从t到t+At时刻，加速度不变，取为常数1（個+倔）。此时，取5=1,2tt+At2a=1。常平均加速度法是应用得最广泛的逐步积分方法之一。4研究表明，当§>0.5，a>0.25（0.5+5时，Newmark方法是无条件稳定的。从式（3.34）和（3.35）可得到，&用u及&、屈t+Att+Att+Attt和u表示的表达式，即有t佃=丄（一u）一丄t+ataAt2t+attaAtt丄一12a丿厨（3.36）tu=邑(“-u)+h-5&+(1-)At&(3.37)t+AtaAtt+attIa丿t2xt考虑t+At时刻的动力方程，有M

26、倔+c4+ku=Ft+Att+Att+Att+At将式(3.36)和(3.37)代入(3.38)，就得到关于/的方程t+Att+At3.38)39)r%|u=ft+Att+At其中r岡-k+丄M1+AClaAt2aAt%-F+Mut+Att+AtAt2t+OLu+I_L-1&)Att12a丿+C-8-&+f-1Attl+('-11At&12a丿t求解方程(3.39),就可得到&，然后，根据式(3.36)t+At(3.37)可解出&和&。t+Att+AtNewmark方法逐步求解的过程如下：和式A.初步计算1)形成刚度矩阵k,质量矩阵M和阻

27、尼矩阵c。2)3)给定初始值u,u和U&000选择时间步长At，参数a和8，并计算积分常数。8>0.50,a>0.25(0.5+S)218a；a；0aAt2iaAt1.8a-1；a-1；324aAt(1-8)；a8At;67(4)形成有效刚度矩阵%:1a=；2aAtAt(8JI2I；52la丿=k对K%作三角分解：K=ldltB对每个时间步计算+am+aC01(1)计算t+At时刻的有效载荷%=f+M(au+au&+1a&)t+Att+At02t3+cvxu+a&+aU)1t4t5t(2)求解t+At时刻的加速度和速度=at+At&=a(u-

28、u)-a品-a&&t+At=i+Au+a&3t+Att6t7t+At我们注意到Wilson-9法与Newmark法的计算关系式，在形式上是相同的，只是其中的系数取不同的值而已。因此，它们可用同一计算机程序来实现。§3.7Houbolt方法这个差分格式是利用t+Jt、八t-At、t-2Jt四个时刻上位移3.40)3.41)是未知的。考的三次插值多项式建立起来的。即假定&=(2u-5u+4u-u)t+At/2t+Attt-Att-2At1&(11u-18u+9u-2u)t+At6Att+Attt-Att-2At这里认为u,u和u是已知的，而ut-2Att-Attt+At3.42)虑t+Jt时刻的动力方程，有M&+C&+Ku=Ft+Att+Att+Att+At将式(3.40)和(3.41)代入式(3.42)中，就得到求解ut+Jt时刻的方程为(211、M+-C+K)u(t26