有限元第六章-動(dòng)力問(wèn)題的有限元法

PAGEPAGE26第六章動(dòng)力問(wèn)題的有限元法6.1概述前面幾章所研究的問(wèn)題都屬于靜力問(wèn)題，其特點(diǎn)是施加到結(jié)構(gòu)上的外載荷不會(huì)使結(jié)構(gòu)產(chǎn)生加速度，且外載荷的大小和方向不隨時(shí)間變化，因而結(jié)構(gòu)所產(chǎn)生的位移和應(yīng)力也不隨時(shí)間變化。本章將要研究結(jié)構(gòu)分析中另一類(lèi)重要問(wèn)題的有限元解法，即動(dòng)力問(wèn)題的有限元解法。動(dòng)力學(xué)問(wèn)題的特點(diǎn)是，載荷是隨時(shí)間變化的，因而結(jié)構(gòu)所產(chǎn)生的位移和應(yīng)力是時(shí)間的函數(shù)，結(jié)構(gòu)會(huì)產(chǎn)生速度和加速度。由于結(jié)構(gòu)本身的彈性和慣性，結(jié)構(gòu)在動(dòng)力載荷的作用下，往往呈現(xiàn)出振動(dòng)的運(yùn)動(dòng)形態(tài)。結(jié)構(gòu)振動(dòng)是工程中一個(gè)很普遍很重要的問(wèn)題。有些振動(dòng)對(duì)我們有利，例如，振動(dòng)打樁，振動(dòng)選料，有些振動(dòng)對(duì)我們有害，例如，機(jī)床的振動(dòng)，儀器與儀表的振動(dòng)，橋梁、水壩及高層建筑在地震作用下的振動(dòng)等。因此，我們必須對(duì)振動(dòng)體本身的振動(dòng)特性以及它對(duì)外部激振力的響應(yīng)有一個(gè)明確的認(rèn)識(shí)，才能更好地利用它有利的一面，而避免它有害的一面，設(shè)計(jì)出更好的機(jī)械和結(jié)構(gòu)。振動(dòng)問(wèn)題主要解決兩方面的問(wèn)題。尋求結(jié)構(gòu)的固有頻率和主振型，從而了解結(jié)構(gòu)的固有振動(dòng)特性，以便更好地利用或減少振動(dòng)。分析結(jié)構(gòu)的動(dòng)力響應(yīng)特性，以計(jì)算結(jié)構(gòu)振動(dòng)時(shí)動(dòng)應(yīng)力和動(dòng)位移的大小及其變化規(guī)律。6.2結(jié)構(gòu)的振動(dòng)方程 結(jié)構(gòu)的振動(dòng)方程可用多種方法建立，這里我們使用達(dá)朗伯原理（動(dòng)靜法），仿照前幾章建立靜力有限元方程的方法，來(lái)建立動(dòng)力問(wèn)題的有限元方程。在靜力問(wèn)題中用有限元法建立的平衡方程是在振動(dòng)問(wèn)題中，對(duì)結(jié)構(gòu)的各節(jié)點(diǎn)應(yīng)用達(dá)郎伯原理所建立的振動(dòng)方程仍然具有與上式相同的形式，只不過(guò)節(jié)點(diǎn)位移是動(dòng)位移，節(jié)點(diǎn)載荷是動(dòng)載荷，它們都是時(shí)間的函數(shù)。上面的方程成為(6.1)上式中為節(jié)點(diǎn)的動(dòng)位移，它是時(shí)間的函數(shù)，是時(shí)刻的節(jié)點(diǎn)位移產(chǎn)生的彈性恢復(fù)力，它與該時(shí)刻的節(jié)點(diǎn)外力構(gòu)成動(dòng)態(tài)平衡。在動(dòng)態(tài)情況下，結(jié)構(gòu)承受的載荷（集中載荷，分布載荷）可隨時(shí)間而變化，是時(shí)間的函數(shù)。按有限元方法將此種載荷移置到節(jié)點(diǎn)上，得到的節(jié)點(diǎn)載荷向量也是時(shí)間的函數(shù)。此外，結(jié)構(gòu)在運(yùn)動(dòng)中，各點(diǎn)除位移以外，還有速度及加速度。按照達(dá)郎佰原理，有加速度的質(zhì)量應(yīng)附加有慣性力載荷。如材料的密度為，則結(jié)構(gòu)單位體積的慣性力為。這對(duì)結(jié)構(gòu)來(lái)說(shuō)，相當(dāng)于又受有另一種體積力，大小與點(diǎn)的加速度成比例，而方向與加速度方向相反。另外，在結(jié)構(gòu)運(yùn)動(dòng)過(guò)程中，還會(huì)受到周?chē)橘|(zhì)和來(lái)自內(nèi)部的阻力。精確地描述這種阻力的變化規(guī)律是很困難的，一般采用阻力與速度成比例的近似線性假定，如阻力系數(shù)為μ，則單位體積的阻力為。這對(duì)結(jié)構(gòu)來(lái)說(shuō)相當(dāng)于另一種體積力，大小與點(diǎn)的速度成比例，方向與速度方向相反。按有限元方法，用單元節(jié)點(diǎn)位移進(jìn)行插值表示單元內(nèi)部位移。(6.2)此處形函數(shù)仍只是位置的插值函數(shù)，與時(shí)間無(wú)關(guān)，則單元內(nèi)的速度和加速度分別為(6.3)以及(6.4)其中、為單元節(jié)點(diǎn)的速度及加速度向量。將單元慣性力與阻力作為體積力，按式(3.28)移置到單元各節(jié)點(diǎn)，就得到相應(yīng)的單元等效節(jié)點(diǎn)載荷向量，記為和，則有將式（6.3）、（6.4）代入上式，有式中(6.5)稱為單元質(zhì)量矩陣。由于推導(dǎo)式（6.5）時(shí)采用了與推導(dǎo)單元?jiǎng)偠染仃嚂r(shí)相一致的形函數(shù)，故式（6.5）所表示的質(zhì)量矩陣也稱為一致質(zhì)量矩陣。而(6.6)稱為單元阻尼矩陣，由式（6.5）和（6.6）可見(jiàn)單元質(zhì)量矩陣和單元阻尼矩陣是對(duì)稱的。將移置到節(jié)點(diǎn)上的動(dòng)載荷、慣性力、阻力作為載荷，按單元疊加，得到有限元節(jié)點(diǎn)位移方程：或(6.7)其中(6.8)稱為結(jié)構(gòu)質(zhì)量矩陣或總質(zhì)量矩陣，而(6.9)稱結(jié)構(gòu)阻尼矩陣。可見(jiàn)結(jié)構(gòu)質(zhì)量矩陣和結(jié)構(gòu)阻尼矩陣分別為單元質(zhì)量矩陣和單元阻尼矩陣的疊加，其疊加方法與結(jié)構(gòu)剛度矩陣的形成完全一樣，借助于單元定位向量，用單元集成法完成疊加過(guò)程。由于和是對(duì)稱的，因而疊加合成的和也是對(duì)稱的。 式（6.7）是節(jié)點(diǎn)位移的二階微分方程，稱為結(jié)構(gòu)的動(dòng)力方程式。對(duì)于不同的結(jié)構(gòu)，可以選用不同的單元，有不同的形函數(shù)矩陣，但動(dòng)力方程（6.7）的建立過(guò)程都是一樣的。當(dāng)結(jié)構(gòu)不受外載荷時(shí)，，如果再忽略阻尼，則動(dòng)力方程（6.7）式成為(6.10)這是系統(tǒng)的自由振動(dòng)方程。 彈性結(jié)構(gòu)的振動(dòng)實(shí)際上是連續(xù)體的振動(dòng)，位移是連續(xù)的，具有無(wú)限多個(gè)自由度。經(jīng)有限元離散化處理后，單元內(nèi)的位移按假定的位移形式來(lái)變化，可用節(jié)點(diǎn)位移插值表示。這樣，連續(xù)系統(tǒng)的運(yùn)動(dòng)就離散化為有限個(gè)自由度系統(tǒng)的運(yùn)動(dòng)了。如果全部節(jié)點(diǎn)有N個(gè)自由度，則式（6.10）就是N階的自由振動(dòng)微分方程了。6.3一致質(zhì)量矩陣與集中質(zhì)量矩陣 按照式（6.5），當(dāng)單元的位移插值形函數(shù)矩陣確定后，就可用此式算出它的一致質(zhì)量矩陣。例如，對(duì)平面三節(jié)點(diǎn)三角形單元，按照2.3節(jié)的敘述，形函數(shù)矩陣可用面積坐標(biāo)表示為其中I為二階單位陣。將上式代入式(6.5)可得平面三節(jié)點(diǎn)三角形單元的一致質(zhì)量矩陣?yán)梅e分公式(2.56)，可由上式求得(6.11)用式(6.5)可以計(jì)算出其它類(lèi)型單元的質(zhì)量矩陣。例如平面桁架單元的質(zhì)量矩陣為(6.12)而平面剛架單元的質(zhì)量矩陣為(6.13)式(6.11)、(6.12)及(6.13)中的Ｍ為單元質(zhì)量。由單元一致質(zhì)量矩陣疊加形成的結(jié)構(gòu)質(zhì)量矩陣，一般都是稀疏、帶狀的，但都有相當(dāng)?shù)陌霂?，如果將單元質(zhì)量矩陣近似作為對(duì)角型矩陣，單元疊加后的結(jié)構(gòu)質(zhì)量矩陣也是對(duì)角陣。這種對(duì)角型的質(zhì)量矩陣稱為集中質(zhì)量矩陣，而對(duì)角型的集中質(zhì)量矩陣對(duì)結(jié)構(gòu)的動(dòng)力分析是非常有利的?？梢詫卧馁|(zhì)量以某種方式分配在單元的節(jié)點(diǎn)上而得到單元集中質(zhì)量矩陣。例如當(dāng)質(zhì)量均勻分布時(shí)，將質(zhì)量平均分配給各節(jié)點(diǎn)。按照這種方法平面三節(jié)點(diǎn)三角形單元的集中質(zhì)量矩陣為(6.14)其中為６階單位陣。Ｍ為單元質(zhì)量。而四節(jié)點(diǎn)四面體單元的集中質(zhì)量矩陣為(6.15)其中[I]為12階的單位陣。對(duì)于矩形彎曲板單元，在不計(jì)轉(zhuǎn)動(dòng)質(zhì)量影響的條件下，它的集中質(zhì)量矩陣為(6.16)對(duì)于平面剛架單元，在不計(jì)轉(zhuǎn)動(dòng)質(zhì)量影響的條件下，它的集中質(zhì)量矩陣為(6.17)6.4阻尼矩陣各種工程結(jié)構(gòu)的阻尼力及其產(chǎn)生的機(jī)理是非常復(fù)雜的。從宏觀上看，阻尼有兩種主要形態(tài)。一種是結(jié)構(gòu)周?chē)承越橘|(zhì)產(chǎn)生的阻尼，稱為粘性阻尼。粘性阻尼的阻尼力一般近似與運(yùn)動(dòng)速度成正比。另一種是結(jié)構(gòu)材料內(nèi)部磨擦產(chǎn)生的阻尼，稱為結(jié)構(gòu)阻尼或材料阻尼。結(jié)構(gòu)阻尼的阻尼應(yīng)力一般近似地認(rèn)為與彈性體的應(yīng)變速率成正比。如果假定阻尼力與運(yùn)動(dòng)速度成正比，那么在運(yùn)動(dòng)的彈性體中任意點(diǎn)處單位體積上作用的阻尼力為式中——比例常數(shù)；——材料密度；——單元形函數(shù)矩陣；——單元節(jié)點(diǎn)速度矢量。可以將阻尼力看成是一種體積力，其等效的單元節(jié)點(diǎn)阻尼力向量為或?qū)懗梢虼耍瑔卧枘峋仃?6.18)它正比于單元質(zhì)量矩陣。如果假定阻尼應(yīng)力與彈性體的應(yīng)變速率成正比，則阻尼應(yīng)力可表示為式中，為比例系數(shù)，為彈性矩陣，為應(yīng)變矩陣，為單元節(jié)點(diǎn)速度向量。下面推導(dǎo)阻尼應(yīng)力的單元等效節(jié)點(diǎn)阻尼力向量。設(shè)單元等效節(jié)點(diǎn)阻尼力向量仍用表示，設(shè)節(jié)點(diǎn)發(fā)生虛位移，單元內(nèi)各點(diǎn)產(chǎn)生的虛應(yīng)變?yōu)椋瑒t。在虛位移上做的虛功為單元的虛應(yīng)變能為由得到或?qū)懗梢虼?，單元阻尼矩?6.19)它正比于單元?jiǎng)偠染仃囋趯?shí)踐中，要精確地計(jì)算各種單元的阻尼矩陣是很困難的。通常在程序設(shè)計(jì)中，假定結(jié)構(gòu)總體阻尼矩陣是結(jié)構(gòu)總體剛度矩陣與總體質(zhì)量矩陣的線性組合。稱為瑞利阻尼，其表達(dá)式為(6.20)其中，比例系數(shù)及可通過(guò)實(shí)驗(yàn)確定。采用瑞利阻尼近似，可以使運(yùn)動(dòng)方程求解大大簡(jiǎn)化，并且在程序中不必單獨(dú)存貯總體阻尼矩陣。在實(shí)際工程問(wèn)題中，阻尼的作用對(duì)結(jié)構(gòu)的動(dòng)力響應(yīng)的影響并不大，這種近似處理具有實(shí)用價(jià)值。6.5結(jié)構(gòu)的自振特性和特征值問(wèn)題結(jié)構(gòu)的自振特性是指結(jié)構(gòu)的振動(dòng)頻率和振型，求結(jié)構(gòu)的自振頻率和振型也稱對(duì)結(jié)構(gòu)進(jìn)行模態(tài)分析，是結(jié)構(gòu)動(dòng)力計(jì)算的主要內(nèi)容之一。計(jì)算經(jīng)驗(yàn)指出，結(jié)構(gòu)的阻尼對(duì)結(jié)構(gòu)的頻率和振型的影響很小，所以求頻率振型時(shí)可以不考慮阻尼的影響。此時(shí)系統(tǒng)的自由振動(dòng)方程如式(6.10)，即(6.21)當(dāng)系統(tǒng)作自由振動(dòng)時(shí)，各質(zhì)點(diǎn)作簡(jiǎn)諧振動(dòng)，各節(jié)點(diǎn)的位移可表示為(6.22)將(6.22)代入(6.21)，并消去公因子得到或(6.23)因此，求解(6.21)式就是尋找滿足式(6.23)的值和非零向量，這種問(wèn)題稱為廣義特征值問(wèn)題。記，和分別稱為廣義特性值和廣義特征向量。式(6.23)可以寫(xiě)成(6.24)這是一個(gè)齊次的線性方程組，若要有的非零解，系數(shù)行列式必須等于零，即(6.25)展開(kāi)此式可得如果彈性結(jié)構(gòu)的總體剛度矩陣和總體質(zhì)量矩陣的階數(shù)是n，則上述行列式展開(kāi)后成為的n次代數(shù)方程式，由此可以求出n個(gè)根，即n個(gè)廣義特征值，從而求出結(jié)構(gòu)的n個(gè)自振頻率。求得廣義特征值后，就可以利用式(6.24)算得對(duì)應(yīng)的廣義特征向量，它代表n個(gè)質(zhì)點(diǎn)的振幅構(gòu)成的振型。由式(6.24)可知，只能被確定到相差一常數(shù)因子的程度。N個(gè)特征值和相應(yīng)的n個(gè)特征向量常稱為n個(gè)特征對(duì)。在彈性結(jié)構(gòu)的振動(dòng)分析中，結(jié)構(gòu)自由度總數(shù)n往往很大，因此無(wú)法直接從上述代數(shù)方程求解廣義特征值。現(xiàn)在已研究出多種有效的廣義特征值問(wèn)題求解方法和相應(yīng)的程序模塊，供求解上述問(wèn)題使用。6.6特征值和特征向量的性質(zhì)下面討論廣義特征方程中特征值和特征向量的性質(zhì)，這些性質(zhì)在結(jié)構(gòu)動(dòng)力分析中是很重要的。１．當(dāng)和矩陣是實(shí)系數(shù)對(duì)稱矩陣時(shí)，其特征值一定是實(shí)數(shù)。如果為正定矩陣時(shí)，則特征值一定是正實(shí)數(shù)，如果為半正定矩陣時(shí)，則特征值一定是非負(fù)實(shí)數(shù)。并且，特征向量也是實(shí)向量。２．當(dāng)和是對(duì)稱矩陣時(shí)，廣義特征方程的不同特征根所對(duì)應(yīng)的特征向量，具有正交性。當(dāng)對(duì)特征向量進(jìn)行正則化處理后，即令時(shí)，則特征向量與質(zhì)量矩陣和剛度矩陣的正交性，可表達(dá)為(6.26)式中為克羅內(nèi)克符號(hào)。式(6.26)第一式稱作正則化條件，而把滿足該式的特征向量稱作正則化特征向量。 設(shè)是由特征值組成的對(duì)角陣，即(6.27)又設(shè)是由n個(gè)特征向量組成的特征向量矩陣，即(6.28)則特征向量的正交性，即式(6.26)也可表示為(6.29) ３．瑞雷（Rayleigh）商 關(guān)于特征值的性質(zhì)，也可從瑞雷商出發(fā)去研究它。對(duì)于廣義特征值問(wèn)題，瑞雷商定義為(6.30) 將特征值由小到大排列成如下順序：設(shè)一向量是n個(gè)特征向量的線性組合，可表示成(6.31)其中為式(6.28)定義的特征向量矩陣，而是由n個(gè)系數(shù)組成的向量。 將式(6.31)代入(6.30)，并注意到特征向量的正交關(guān)系式(6.29),得展開(kāi)后可寫(xiě)成(6.32)可見(jiàn)，瑞雷商在最小特征值和最大特征值之間并且，當(dāng)時(shí)這說(shuō)明當(dāng)取第i階特征向量時(shí)瑞雷商達(dá)到它的一個(gè)極值，該極值就是與對(duì)應(yīng)的特征值，特別是這就是瑞雷商的極小值原理。此原理常被工程界作為預(yù)估振型的依據(jù)去計(jì)算特征值，從而求得基頻的近似值。 如果，即向量和前m-1階特征向量正交，則由式(6.32)看出且當(dāng)時(shí)這就是說(shuō)，是瑞雷商在向量與前m-1階特征向量正交條件下的極小值。 ４．移位 在特征值問(wèn)題的求解過(guò)程中，廣泛使用所謂移位去加速計(jì)算的收斂性。 對(duì)廣義特征值問(wèn)題，可以使作一移位，即(6.33)然后去求解下列特征值問(wèn)題為了確定原特征值問(wèn)題與上述特征值問(wèn)題的關(guān)系，我們將式(6.33)代入上式，得(6.34)與式(6.24)對(duì)比，顯然有(6.35)因此，的特征向量與的特征向量是相同的，但特征值減小一個(gè)值。6.7逆迭代法 求解特征值問(wèn)題的方法很多，但基體上可分為兩大類(lèi)：變換法和迭代法。熟知的雅可比法就是一種變換法，可用于求問(wèn)題的全部特征對(duì)，具有簡(jiǎn)單和穩(wěn)定的優(yōu)點(diǎn)。迭代法有所謂正迭代和逆迭代兩種。由于有限元分析中，系統(tǒng)的自由度往往很高，而對(duì)工程結(jié)構(gòu)進(jìn)行動(dòng)力分析時(shí)，通常只需了解少數(shù)幾個(gè)較低階的特征值和相應(yīng)的特征向量，逆迭代法正是適應(yīng)上述要求而被廣泛應(yīng)用的一種方法；在求解大型特征值問(wèn)題的子空間迭代法中也要用到逆迭代法。現(xiàn)對(duì)此法作一討論。 設(shè)是某特征向量的一個(gè)初始近似向量，并假定。于是可求出式(6.23)右端的慣性力由于是任意假定的，一般不能滿足平衡條件，即而按照靜力平衡的關(guān)系可解得另一近似向量。通常是一個(gè)比更好的近似特征向量。通過(guò)反復(fù)迭代就能得到滿足精度要求的特征向量。 在逆迭代法中，假定是正定的。迭代公式為(6.36)由于特征向量各分量的模是任意的，只是它們之間保持一定的比例。為確定起見(jiàn)，在每次迭代中，對(duì)求出的須進(jìn)行正則化處理。利用正則化條件(6.26),可得正則化了的，即(6.37) 如果所選初始迭代向量不與第一階特向量關(guān)于矩陣正交，即則當(dāng)時(shí) 關(guān)系式(6.36)和(6.37)是逆迭代的基本公式，然而在計(jì)算機(jī)中執(zhí)行時(shí)，按如下方式進(jìn)行更為有效。假設(shè)，按，進(jìn)行下列迭代計(jì)算：(6.38)式中，若，則當(dāng)時(shí) 在迭代過(guò)程中，由(6.38)的第三式可求得由瑞雷商給出的特征值的近似值。該近似值可用來(lái)確定迭代法所達(dá)到的精度。當(dāng)條件(6.39)滿足時(shí)，迭代就終止，令是最后一次迭代，則有(6.40)為了保證特征向量準(zhǔn)確到r位精度，則應(yīng)要求特征值準(zhǔn)確到2s位精度，正如式(6.39)所表示的。這個(gè)要求可用瑞雷商來(lái)證明。 現(xiàn)在來(lái)證明逆迭代的收斂性。設(shè)初始迭代向量是n個(gè)特征向量的線性組合，可表示為式中系數(shù)不能為零。應(yīng)用式(6.36)和(6.37)進(jìn)行迭代，并注意到關(guān)系式(6.26)和(6.24)，可將第k次迭代后的向量和寫(xiě)成(6.41)如果假設(shè)，將上式分子和分母同乘以，則得(6.42)顯然，當(dāng)。由上式也可看出，當(dāng)為特征值時(shí)，趨于它們對(duì)應(yīng)的特征向量的線性組合。由上式可以看出，迭代收斂的速度與比值有關(guān)，比值愈小收斂愈快。當(dāng)時(shí)，收斂可能很慢；當(dāng)收斂就可能很快。如果我們采用移位，將移至原點(diǎn)附近但仍保持正值，則收斂速度會(huì)很快。 移位不但可用來(lái)加快收斂速度，而且可用來(lái)求中間特征對(duì)。如果作一移位，依據(jù)式(6.34)，使為負(fù)值，而為接近于零的正值，則逆迭代的結(jié)果就可求出特征對(duì)。同樣可以用移位來(lái)求出其它特征對(duì)。 當(dāng)用移位法依次自小至大求特征對(duì)時(shí)，有時(shí)很難收斂到指定的特征對(duì)，由于計(jì)算中的舍入誤差很可能收斂到已求出的特征對(duì)，尤其當(dāng)存在重特征值時(shí)更是如此。 為了不收斂到已求出的特征對(duì)，可以利用特征向量的正交性，在開(kāi)始選擇初始迭代向量時(shí)就進(jìn)行正交化處理。具體說(shuō)，在求出m-1階特征向量后，可以構(gòu)造一新的迭代向量(6.43)其中是任意向量，是待定常數(shù)，它們由與關(guān)于質(zhì)量矩陣的正交條件確定。由此，應(yīng)取(6.44) 由于與前m-1階特征向量關(guān)于正交，則用它作為初始向量進(jìn)行逆迭代可指望求得第m階特征對(duì)。為了防止計(jì)算中的舍入誤差而使偏離正交方向，還應(yīng)在每一步開(kāi)始都對(duì)進(jìn)行正交化處理。以上迭代過(guò)程稱為克萊姆－史密特(Gram-Schmidt)正交化。例6.1 圖6.1示一三層剛架結(jié)構(gòu)，各層的樓面質(zhì)量分別為；各層的側(cè)移剛度分別為。求剛架的固有頻率和型。圖6.1 解該剛架的質(zhì)量矩陣和剛度矩陣為 我們先用行列式法直接求這一問(wèn)題的特征值。由式(6.25)(a)式中 令式(a)矩陣行列式展開(kāi)、化簡(jiǎn)并令其等于零、得到三次方程，解之得從而得到而由式(6.23)可得三個(gè)主振型向量，用振型矩陣表示為 現(xiàn)在用式(6.38)的逆迭代法來(lái)解同一問(wèn)題。為以后迭代計(jì)算方便先求出剛度矩陣的逆矩陣采用初始迭代向量則當(dāng)k=1時(shí)以下各次循環(huán)做法相同，現(xiàn)將結(jié)果載入下表6.1內(nèi)。 從上表看出，比收斂得快些，經(jīng)過(guò)5次迭代就已經(jīng)達(dá)到了要求，。因此由得到可見(jiàn)第一振型的結(jié)果與前面用行列式法求得的結(jié)果相符合。表6.1122.516.5922.524.7518197.9798——0.6857280.7542980.548580212.41168.297283.9772112.411612.44597.95442191.63930.03308530.7301550.7321730.467946312.62858.247543.8605512.628512.37137.72110191.36640.00142630.7398900.7248210.452317412.66768.22833.8341712.667612.34257.66834191.35370.00006660.7420430.7229990.449196512.67598.223613.8284412.675912.33547.65696191.35300.00000320.742520.7225750.448524612.67768.222533.8272412.677612.33387.65448191.35300.00000020.7426220.7224830.4483816.8用逐步積分法求結(jié)構(gòu)動(dòng)力響應(yīng) 6.8.1概述求結(jié)構(gòu)的動(dòng)力響應(yīng)，在數(shù)學(xué)上就是要求出運(yùn)動(dòng)方程(6.7)的解答。式(6.7)是一個(gè)二階常系數(shù)微分方程組，可以用數(shù)值積分的方法對(duì)方程直接求解，即按時(shí)間增量逐步求解運(yùn)動(dòng)微分方程，直至反應(yīng)終了，這一方法稱作逐步積分法，由于在數(shù)值積分前不對(duì)運(yùn)動(dòng)方程進(jìn)行變換，所以又稱作直接積分法。逐步積分法既可用于求解線性結(jié)構(gòu)體系問(wèn)題——在整個(gè)動(dòng)力反應(yīng)過(guò)程中，，矩陣保持不變的問(wèn)題；也可用于求解非線性結(jié)構(gòu)體系的問(wèn)題——，，矩陣隨動(dòng)力反應(yīng)的過(guò)程而變化的問(wèn)題。此處，我們只討論線性結(jié)構(gòu)體系的問(wèn)題。逐步積分法求解運(yùn)動(dòng)微分方程的基本思路是：１．把連續(xù)的時(shí)間過(guò)程離散為有限個(gè)點(diǎn)，對(duì)于運(yùn)動(dòng)微分方程只要求它們?cè)谏鲜雒總€(gè)時(shí)間離散點(diǎn)上得到滿足，也就是說(shuō)，最終求解得到的只是位移、速度和加速度在有限個(gè)時(shí)間離散點(diǎn)上的值，而不是連續(xù)函數(shù)。２．在每個(gè)時(shí)間間隔內(nèi)，假定位移、速度和加速度符合某一簡(jiǎn)單的關(guān)系。而的選擇，要求保證計(jì)算的穩(wěn)定性與精確性。從這樣一個(gè)基本思路出發(fā)，形成了多種逐步積分的方法，如中心差分法；線性加速度法；Wilson-法；Newmark法等等。逐步積分法的優(yōu)點(diǎn)在于，計(jì)算的精確度可隨的減少而提高，當(dāng)大型結(jié)構(gòu)體系受到歷時(shí)較短的脈沖荷載作用時(shí)，由于要考慮高階振型的反應(yīng)使計(jì)算振型的數(shù)目較多，而逐步積分法求解動(dòng)力反應(yīng)時(shí)，可省去特征方程的求解工作，在經(jīng)濟(jì)上、效率上可能是有利的。但是，逐步積分法的缺點(diǎn)，在于它不能同時(shí)給出結(jié)構(gòu)的振型和自振頻率，而是直接給出結(jié)構(gòu)的反應(yīng)過(guò)程，并且當(dāng)反應(yīng)的時(shí)間歷程較長(zhǎng)，取的較短時(shí)，計(jì)算工作量也是相當(dāng)可觀的。 6.8.2Wilson-法Wilson-法實(shí)質(zhì)上是線性加速度法的推廣，線性加速度法假定從時(shí)刻到的時(shí)間內(nèi)加速度是線性變化的。而Wilson-法則假定在時(shí)間區(qū)間到內(nèi)加速度是線性變化的，如圖6.1，其中。當(dāng)時(shí)，Wilson-法法就退化為線性加速度法。Wilson證明，要達(dá)到無(wú)條件穩(wěn)定，必須選用，通常采用。圖6.2Ｗilson-法按Wilson-法的假定，在時(shí)間區(qū)間到內(nèi)加速度的數(shù)學(xué)表達(dá)式為(6.45)式中，為時(shí)間的增量，將式(6.45)積分，得速度表達(dá)式(6.46)和位移表達(dá)式(6.47)現(xiàn)在的問(wèn)題是如何根據(jù)時(shí)刻的位移、速度和加速度和上述位移、速度和加速度關(guān)系式以及運(yùn)動(dòng)方程來(lái)求時(shí)刻的位移、速度和加速度。將代入式(6.46)及(6.47)，得到時(shí)刻的速度、位移計(jì)算式(6.48)(6.49)選擇為基本未知量。由式(6.49)和(6.48)，可以得到用表示的及的算式，即(6.50)(6.51)時(shí)刻的運(yùn)動(dòng)方程是(6.52)式中，是時(shí)刻的總體載荷向量，但我們只知道時(shí)刻和時(shí)刻的總體載荷向量，和?？捎镁€性外推法求出為(6.53)將(6.50)，(6.51)代入(6.52)式，得到只有一個(gè)未知量的方程組(6.54)式中(6.55)(6.56)求解(6.54)式，可以得到。再將它代入式(6.51)和(6.50)得到和。將代入式(6.45)，(6.46)和(6.47)。得到、和為(6.57)(6.58)(6.59) 為編程方便，可以進(jìn)一步簡(jiǎn)化運(yùn)算過(guò)程。將式(6.50)代入(6.57)得(6.60)由式(6.57)得將上式代入(6.58)和(6.59)式得(6.61)(6.62) Wilson-法的計(jì)算步驟歸納如下：１．初始計(jì)算（１）形成總剛度矩陣和總質(zhì)量矩陣，而總阻尼矩陣不存貯。（２）形成初始值，并計(jì)算。（３）選取時(shí)間步長(zhǎng)和（一般取），計(jì)算積分常數(shù)。得（４）形成有效剛度矩陣（５）對(duì)進(jìn)行三角分解２．對(duì)每一時(shí)間步長(zhǎng)，循環(huán)計(jì)算（１）計(jì)算時(shí)刻的有效載荷向量（２）求解時(shí)刻的位移向量（３）計(jì)算時(shí)刻的加速度、速度和位移向量 例6.2有一簡(jiǎn)單的結(jié)構(gòu)系統(tǒng)，不考慮阻尼影響，它的運(yùn)動(dòng)方程是試用Wilson-法計(jì)算此結(jié)構(gòu)系統(tǒng)從0到3.5的位移響應(yīng)。已知結(jié)構(gòu)初始狀態(tài)是靜止的，受突加不變載荷作用。解：選擇時(shí)間步長(zhǎng)。為保證計(jì)算精度，通常取，其中為結(jié)構(gòu)系統(tǒng)自由振動(dòng)的最小周期。對(duì)于復(fù)雜結(jié)構(gòu)系統(tǒng)，一般高頻部分對(duì)結(jié)構(gòu)動(dòng)力響應(yīng)的貢獻(xiàn)不大，可取低頻部分的最小周期。本例的最小周期，取,計(jì)算12個(gè)時(shí)間步長(zhǎng)。將初始條件代入運(yùn)動(dòng)方程，求得加速度初始條件。取，計(jì)算Wilson－法的常數(shù)，得形成有效剛度矩陣對(duì)于每一個(gè)時(shí)間步長(zhǎng)，需先求解式中，然后計(jì)算加速度、速度和位移12個(gè)時(shí)間步長(zhǎng)按以上諸式計(jì)算的結(jié)果如下表所示：t234567891011120.006050.05250.1960.4900.9521.542.162.672.922.822.331.540.3661.342.643.924.885.315.184.613.823.062.522.29 6.8.3Newmark法Newmark法也是一種線性加速度法的推廣。它采用如下兩個(gè)基本假定：(6.63)(6.64)其中和是參數(shù)，根據(jù)積分的精度和穩(wěn)定性要求確定。當(dāng)和時(shí)，關(guān)系式(6.63)、(6.64)相當(dāng)于線性加速度法。當(dāng)和時(shí)，相當(dāng)于常平均加速度法。已經(jīng)證明，當(dāng)和時(shí)，Newmark法是無(wú)條件穩(wěn)定的。根據(jù)上述兩個(gè)假定，以及時(shí)刻的運(yùn)動(dòng)方程，就可以求得時(shí)刻的位移、速度和加速度的解。將作為基本未知量，由式(6.64)求得(6.65)將式(6.65)代入式(6.63)得(6.66)將式(6.65)和(6.66)代入（）時(shí)刻的運(yùn)動(dòng)方程(6.67)得到只有一個(gè)基本未知量的方程(6.68)式中(6.69)(6.70)解方程(6.68)得到后，再代入式(6.65)和(6.63)求出和。為了便于編寫(xiě)計(jì)算機(jī)程序，將Newmark法的計(jì)算步驟綜述如下： １．初始計(jì)算（１）形成剛度矩陣，質(zhì)量矩陣，阻尼矩陣不存貯。（２）形成初始向量、，計(jì)算。（３）選擇時(shí)間步長(zhǎng)和參數(shù)。計(jì)算積分常數(shù)，得（４）形成有效剛度矩陣（５）對(duì)進(jìn)行三角分解。２．對(duì)每一個(gè)時(shí)間步長(zhǎng)循環(huán)（１）計(jì)算在時(shí)刻的有效載荷向量（２）求解時(shí)刻的位移向量（３）計(jì)算時(shí)刻的加速度和速度 可以看出，Wilson-法和Neumark法的計(jì)算步驟是相同的，計(jì)算公式的形式也是一樣的，只是系數(shù)有所不同，故可以在一個(gè)有限元程序中同時(shí)存在這兩種直接積分法供用戶選擇。例6.3用Newmark計(jì)算例題6.2的結(jié)構(gòu)系統(tǒng)的響應(yīng)。取，。 解：選擇。將代入運(yùn)動(dòng)方程，得。計(jì)算紐馬克法的常數(shù) 形成有效剛度矩陣為 對(duì)于每一個(gè)時(shí)間步長(zhǎng)，先求解下列方程 式中然后計(jì)算 用以上諸式，按12個(gè)時(shí)間步長(zhǎng)依次計(jì)算的結(jié)果如下時(shí)間456100.006730.05040.1890.4850.9611.582.232.763.002.852.281.400.3641.352.684.004.955.345.134.483.642.902.442.316.9用振型疊加法求動(dòng)力響應(yīng)在直接積分法中，將整個(gè)動(dòng)力響應(yīng)歷程劃分成若干個(gè)時(shí)間步長(zhǎng)。為保證計(jì)算精度，必須把每個(gè)時(shí)間步長(zhǎng)的時(shí)間間隔取得比較小，每個(gè)時(shí)間步長(zhǎng)都要求解一次線性方程組 或盡管對(duì)于線性彈性結(jié)構(gòu)，有效剛度矩陣可以在初始計(jì)算時(shí)就三角分解好，但每個(gè)時(shí)間步長(zhǎng)進(jìn)行前消回代仍然要花費(fèi)不少機(jī)時(shí)。一般來(lái)講，計(jì)算短時(shí)間的動(dòng)力響應(yīng)，使用直接積分法是很有效的，但是，要進(jìn)行長(zhǎng)時(shí)間的動(dòng)力響應(yīng)計(jì)算，直接積分法就變得不經(jīng)濟(jì)了。如果在逐步求解開(kāi)始之前，先將運(yùn)動(dòng)方程進(jìn)行變換，把有效剛度矩陣化成對(duì)角矩陣，那么每一時(shí)間步長(zhǎng)的計(jì)算就非常簡(jiǎn)單，使逐步求解變得比較經(jīng)濟(jì)。振型疊加法在逐步求解之前；先將運(yùn)動(dòng)方程進(jìn)行變換化簡(jiǎn)?？紤]到彈性結(jié)構(gòu)無(wú)阻尼自由振動(dòng)時(shí)的振型向量具有對(duì)質(zhì)量矩陣和剛度矩陣的正交性，即 (6.71)式中，因此，采用振型矩陣作為變換矩陣，總體位移向量可表示為(6.72)或

其物理意義是結(jié)構(gòu)的總體位移向量可以用結(jié)構(gòu)振型的線性組合來(lái)表示。由于振型矩陣在求解結(jié)構(gòu)自由振動(dòng)，即求解特征值問(wèn)題時(shí)已得到，因此計(jì)算位移向量的問(wèn)題變成了計(jì)算未知向量的問(wèn)題。需要特別指出的是，振型矩陣與時(shí)間無(wú)關(guān)。將(6.72)式代入運(yùn)動(dòng)方程，并左乘得 （6.73）考慮到振型的正交性，則 式中(6.74)—固有頻率，等于i=1,2,…,n;—階單位矩陣。如果采用瑞利阻尼，即那么振型對(duì)阻尼矩陣也是正交的，其表達(dá)式為 命(6.75) 式中，稱為振型阻尼系數(shù)i=1,2,…,n。 因此式(7.73)可寫(xiě)成 (6.76)由于都是對(duì)角陣，(6.76)式是一組相互獨(dú)立的二階微分方程，可寫(xiě)成： (6.77)式中(6.78) (6.77)式是二階常微分方程，可以用直接積分法來(lái)求解，也可以用杜哈梅(Duhamel)積分來(lái)求得，如 (6.79)式中，。和由初始條件確定。杜哈梅積分一般采用數(shù)值積分來(lái)計(jì)算。 向量的初始值，即時(shí)的值，可由下式求出