wilson法和newmark法地理論的過程

1、第三章離散化結(jié)構(gòu)動力方程的解法（2013.4.24）§3.1緒言對于一個實際結(jié)構(gòu)，由有限元法離散化處理后，應(yīng)用瞬時最小勢能原理可導出動力方程MH&+tcuku=（F（f）（3.1）這里，倔、個、“及f（t）分別表示加速度、速度、位移及所作用的外力矢量，他們都是與時間有關(guān)的。從數(shù)學的角度來看，式（3.1）是一個常系數(shù)的二階線性常微分方程組，對于它的求解原則上并無困難。但是，由于M、C和K的階數(shù)非常高，使得式（3.1）的求解必須花費很大的代價，便促使人們?nèi)で笠恍┬矢叩慕朴嬎惴椒āＤ壳?，用于求解式?.1）的方法，大致可分為兩大類。一是坐標變換法，它是對結(jié)構(gòu)動力方程式（3.1

2、），在求解之前，進行模態(tài)坐標變換，實際上就是一種Ritz變換，即把原物理空間的動力方程變換到模態(tài)空間中去求解。現(xiàn)在，普遍使用的方法是模態(tài)（振型）迭加法，即用結(jié)構(gòu)的前q階實際主模態(tài)集（主振型陣）構(gòu)成坐標變換陣進行變換。通過這一變換，實現(xiàn)降階，求較好的近似解，而且，還用解除耦合的辦法，簡化方程的計算。還有一種所謂假設(shè)模態(tài)法，即是用一組假設(shè)模態(tài)，構(gòu)成模態(tài)坐標變換陣進行變換，獲得一組降階的而不解耦的模態(tài)基坐標方程。顯然，這種方法的計算精度，取決于所假設(shè)的模態(tài)。用Ritz矢量法求解的近似模態(tài)作為假設(shè)模態(tài)，可得到滿足要求的精度。二是直接積分法，它是對式（3.1）在求解之前，不進行坐標變換，直接進行數(shù)值積分

3、計算。這種方法的特點是對時域進行離散，將式（3.1）分為各離散時刻的方程，然后，將該時刻的加速度和速度用相鄰時刻的各位移線性組合而成，于是，式（3.1）就化為一個由位移組成的該離散時刻上的響應(yīng)值，通常又稱為逐步積分法。線性代數(shù)方程組的解法與靜力時刻的位移來線性組合，就導致了各種不同的方法。主要有中央差分法，Houbolt方法，Wilson-0法和Newmark方法等。§3.2模態(tài)（振型）迭加法設(shè)有n個自由度的系統(tǒng)，在外力f（t）的作用下，常常被激起較低階的一部分模態(tài)（即振型），而絕大部分高階模態(tài)被激起的分量很小，一般可忽略不計。例如，在地震載荷作用下，通常，只有最低的二階，三階模態(tài)起

4、主要作用。所以，對于這樣的一些問題，采用模態(tài)迭加法是有效的。設(shè)有式（3.1）的n階動力方程，起主要作用的是其前q階模態(tài)，通常取q=n。按Ritz變換，則可將式（3.1）中的“用前q個模態(tài)的線性組合來表示，即u=Y他+Yb+.+Y他=£bY1122qqjjj=1=Y（3.2）其中，為結(jié)構(gòu)的已知的保留主模態(tài)矩陣，而丫是維的模nxqqx1態(tài)基坐標矢量，它形成了一個q維的模態(tài)空間。它表示在丫中，各階主模態(tài)所占有的成分的多少。假定已用第二章所述的某一方法解出，再將式（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階就小的多了，實現(xiàn)了降階，因而也就容易求解多了。若展開上述的M*的表達式，根據(jù)主模態(tài)（主振型）關(guān)于M的表達式，根據(jù)主模態(tài)的（主振型）關(guān)于M的正交性質(zhì)，可知m*=0（i主j）ij所以，M*是一個對角陣。同理可知K*也是一個對角陣。然而，在一般的情況下，C*是一個非對角陣，即在模態(tài)空間中，系統(tǒng)的的阻尼一般是耦合的。因此，式（3.3）是一個完全解耦的動力學方程。但是，它是一個已降階的q階的動力方程，可使用后面即將介紹的直接積分法求解。當系統(tǒng)的阻尼為比例阻尼時，即C可以表示為C*=aM*+PK*（3.4）則C*為對角陣。此

6、外，若系統(tǒng)的阻尼是一般的的線性阻尼，并非比例阻尼，但是只要結(jié)構(gòu)的固有頻率不相等，而且不十分接近，則可用舍去C*陣中的非對角元來實現(xiàn)C*的對角陣，也不會引起太大的誤差。在上述兩種情況下，可以獲得對于模態(tài)坐標的完全解耦的動力學方程。即式（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決定。由于有阻尼的存在，由初始條件所激發(fā)的振動，隨時間的增長而衰減以致消失。因此，常可不計式(3.7)中的第二項，即是由初始條件激發(fā)的自由衰減振動。計算出y(t)后，便可利用式(3.2)，計算出物理坐標的響應(yīng)iu(t)。數(shù)學計算步驟可

8、歸納如下:第一步：根據(jù)結(jié)構(gòu)的離散化模型，建立系統(tǒng)的M,K以及F(t)，并進行結(jié)構(gòu)的固有特性分析，即求解特征值問題(K-32M)©=0求出前q階特征對(,帥),(i二1,2,.,q)ii第二步：形成模態(tài)陣二©©.©，并建立模態(tài)基坐nxq12q標下的動力方程少+疋3y+32y=f(t)(i=1,2,.,q)iiiiiii其中f(t)二丄0tF(t)，而m二©TM(H。根據(jù)實驗結(jié)果或經(jīng)驗imiiiii數(shù)據(jù)確定各階主振動中的比例阻尼g。i第三步：求解主模態(tài)基坐標的動力方程，有1y(t)=Jtf(T)ePi(t-T)sinB(t-t)dz，其中，iB0i

9、iiB=3Ji-g2。iii第四步：進行坐標變換后，求得動力響應(yīng)u=y§3.3模態(tài)假設(shè)法上節(jié)所述的模態(tài)迭加法，是用系統(tǒng)的真實主模態(tài)組成的模態(tài)矩陣，再對系統(tǒng)的物理坐標進行模態(tài)坐標變換，從而在主模態(tài)空間中得到降階并解耦的動力學方程，這樣來實現(xiàn)簡化計算。而這里提出的假設(shè)模態(tài)法，則是用一組假設(shè)模態(tài)矩陣，對系統(tǒng)的物理坐標進行模態(tài)坐標轉(zhuǎn)換，從而在模態(tài)空間中得到一組只降階的動力學方程。若令假設(shè)模態(tài)矩陣為，而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)。它們分別對應(yīng)于假設(shè)模態(tài)坐標q的質(zhì)量矩陣、阻尼矩陣、剛度矩陣與廣義力列陣。因為矩陣0中的各列都是假設(shè)模態(tài)，它們般不具有正交性，所以，m、c和k*都不是對角陣。于是，方程(3.11)是不能解耦的方程組，但它卻是比式(3.1)的階數(shù)要低得多了。顯然，對式(3.11)采用直接積分法求解，將比對式(3.1)求解要簡便得多。這是假設(shè)模態(tài)法的優(yōu)點。假設(shè)模態(tài)法的計算精度，很顯然地是取決于假設(shè)模態(tài)陣中模態(tài)假設(shè)的好壞與質(zhì)量。因此，應(yīng)用假設(shè)模態(tài)法能否成功的關(guān)鍵在于確定出個適宜的假設(shè)模態(tài)矩陣。在第五章中，我們介紹了幾種構(gòu)造

11、假設(shè)模態(tài)的方法。實際上，在§2.9中介紹的Rayleigh-Ritz分析，可認為是種假設(shè)模態(tài)法。它的作用，在于降低方程的階數(shù)，簡化計算。它的基本思想是，事先假定出若干近似的特征矢量，然后按照這些特征矢量的最佳線性組合，而算得前若干階特征值的近似值。顯然，運用這種方法時，其計算精度與事先假定的特征矢量的近似程度和數(shù)量有關(guān)。按照Ritz變換的思想，找到了近似的特征矢量X后，i(=0,q)，即有)=aX+aX+aX=XA(3.12)1122qq(3.13)求解如下的廣義特征值問題，即KA=pMa,p=rn2其中k=xPkxM=xtmxk和m為原結(jié)構(gòu)離散化之剛度陣和質(zhì)量陣，它們都是n階方陣。

12、求解式（3.13），得到q個特征矢量，有A=a1ai112A=a2a2a1qa2q<Ai=faqaqq12再按照Ritz的變換，即式(3.12)矢量e，1aqTq由特征矢量A，可計算出（），£，即是2qe=faixijij=1212精彩文檔(3.14)現(xiàn)在用nxq=q好來表示此變換陣，它就是我們要構(gòu)造的假設(shè)模態(tài)矩陣。§3.4中心差分法（顯示法）現(xiàn)在開始討論直接積分法，或稱逐步積分法。前面討論的模態(tài)迭加法，并非總是有效的。當剛度矩陣k，或質(zhì)量矩陣m，或阻尼矩陣c出現(xiàn)隨時間變化時，或當外荷載激起的振型太多，需要計算的特征對太大時，就不宜于采用模態(tài)迭加法，在這些情況下，采用

13、逐步積分法是適宜的。中心差分法就是其中的一種。這種方法的特點，是將動力方程在時間域上離散，化成對時間的差分格式，然后根據(jù)初始條件，利用直接積分法逐步求解出一系列時刻上的響應(yīng)值。假定t=0時，位移、速度和加速度分別為已知的u,&和鈕。000再將求解的時間區(qū)間劃分為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）這樣，上式就化為用相鄰時刻的位移表示的代數(shù)方程組。由它可解出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，質(zhì)量矩陣M和阻尼矩陣Co(2) 給定初始值u,U和Uo000(3) 選擇時間步長At,At<At，并計算積分常數(shù)：cr1，1，1oa=a=a2a,a=°0At212At203a2(4) 計算uu-Atu+auo-At0030(5) 形成有效質(zhì)量矩陣雨aM+aCo01(6) 二角分解M:應(yīng)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應(yīng)當指出，這種中央差分算法，左端的系數(shù)矩陣只與質(zhì)量陣M和阻尼陣C有關(guān)，而與剛度陣K無關(guān)。如果質(zhì)量陣和阻尼陣是對角陣，那么在解方程時，就不需要對系數(shù)陣進行二角分解，即不需要解線性代數(shù)方程組，從第一步開始逐次直接求得各個時刻u的值，這是中央差分格式就是一種顯示的格式。此外，由t+At于不求解代數(shù)方程組，也就不需要進行組集，它的右端項的形成也只須

17、在單元一級水平上，由每個單元對有效載荷矢量的貢獻迭加而成。因此，ADINA程序規(guī)定，在用中心差分法時，必須使用對角的質(zhì)量陣和阻尼陣。從計算穩(wěn)定性角度來看，中心差分法的缺點，在于它是條件穩(wěn)定的，即當時間步長At太大時，積分是不穩(wěn)定的。所以，對步長的限制是TAt<At=ncr冗這里，At是臨界步長值，T是有限元系統(tǒng)的最小周期。這樣，crn當T很小時，就限制了At必須很小，所以求解所花的代價就很大。n§3.5線性加速度法和Wilson-°法線性加速度法和Wilson-0法，都是屬于逐步積分法。線性加速度法是假定在t,t+At時間間隔內(nèi)，即在步長At時間內(nèi)，加速度U&

18、t+t)呈線性變化，其表達式為&=&+TA(3.20)t+Tt其中，A二(U&U)/At。但是，這個方法不是無條件穩(wěn)定的，t+T所以在應(yīng)用上受到限制。70年代初期，Wilson推廣了線性加速度法，他假定在此步長At更大的時間區(qū)間(t,t+0At)內(nèi)，加速度仍保持線性變化，經(jīng)過證明，當0、1.37時，這一方法是無條件穩(wěn)定的，這就是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區(qū)間內(nèi)，對式(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(xiàn)%t+腦=+9(Ft+At-Ft)將(3.26)和式(3.27)代入式(3.28)，就得到關(guān)于町的方程為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，質(zhì)量矩陣M和阻尼矩陣C。(2) 給出初始值u，&u和&u&。000(3) 選擇時間步長At，取。二1.4，并計算積分常數(shù)，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法是隱式積分，即每計算一步，必須解一個線性代數(shù)方程組。當1.37時，它是無條件穩(wěn)定的。此外，這種算法是自起步的，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和§是按積分的精度和穩(wěn)定性要求可以調(diào)整的參數(shù)。5=1,a=1時，它就是線性加速度法，所以，Newmark方法也26可以理解為線性加速度法的一個小延伸。Newmark法最初提出作為無條件穩(wěn)定的一種積分格式是常平均加速度

25、法，即假定從t到t+At時刻，加速度不變，取為常數(shù)1（個+倔）。此時，取5=1,2tt+At2a=1。常平均加速度法是應(yīng)用得最廣泛的逐步積分方法之一。4研究表明，當§>0.5，a>0.25（0.5+5時，Newmark方法是無條件穩(wěn)定的。從式（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)，就得到關(guān)于/的方程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),就可得到&，然后，根據(jù)式(3.36)t+At(3.37)可解出&和&。t+Att+AtNewmark方法逐步求解的過程如下：和式A.初步計算1)形成剛度矩陣k,質(zhì)量矩陣M和阻

27、尼矩陣c。2)3)給定初始值u,u和U&000選擇時間步長At，參數(shù)a和8，并計算積分常數(shù)。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法的計算關(guān)系式，在形式上是相同的，只是其中的系數(shù)取不同的值而已。因此，它們可用同一計算機程序來實現(xiàn)。§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