第十章-動力問題的有限元法34頁

1、第十章 動力問題的有限元法動力問題：結(jié)構(gòu)受載荷作用沒有達(dá)到平衡狀態(tài)，或由于結(jié)構(gòu)的彈性和慣性而圍繞平衡位置振動時，其位移、應(yīng)力、應(yīng)變都是時間的函數(shù)，除各點(diǎn)有位移外，還有速度和加速度。10-1 結(jié)構(gòu)的動力學(xué)方程在動態(tài)情況下，結(jié)構(gòu)承受的載荷可以隨時間而變化，是時間的函數(shù)。按有限元方法，將此載荷以作功等效的原則，移置到結(jié)點(diǎn)上，得到的結(jié)點(diǎn)載荷列陣也應(yīng)該是時間的函數(shù)，用表示。此外，結(jié)構(gòu)在運(yùn)動中，各點(diǎn)除位移外，尚有速度和加速度，分別用、和表示。顯然，、和不僅是坐標(biāo)的函數(shù)，還是時間的函數(shù)，即 采用有限元法對結(jié)構(gòu)進(jìn)行離散，對任一單元，若單元的結(jié)點(diǎn)位移為、單元的結(jié)點(diǎn)速度為、單元的結(jié)點(diǎn)加速度向量為，則單元內(nèi)的位移可

2、由結(jié)點(diǎn)位移插值得到，于是 （10-1） 按照達(dá)朗貝爾原理：一個運(yùn)動中的物體，若加上慣性力，則可看為靜力平衡狀態(tài)。設(shè)材料密度為，則結(jié)構(gòu)內(nèi)單位體積的慣性力為，這對結(jié)構(gòu)來說，相當(dāng)于受到另外一種體力分布載荷，大小與質(zhì)點(diǎn)的加速度成正比，而方向與加速度方向相反。另外，在介質(zhì)中運(yùn)動的質(zhì)點(diǎn)總是會受到阻力作用，如阻力與速度成正比，則單位體積的阻力為，表達(dá)式中為阻力系數(shù)，負(fù)號表明阻力方向與速度相反。其次，設(shè)結(jié)構(gòu)還受到體力和面積力作用，則結(jié)構(gòu)所受總的體力為： 根據(jù)達(dá)朗貝爾原理，結(jié)構(gòu)受到所有力（包括慣性力）作用處于靜力平衡，于是 （10-2）上式等式右邊的第一積分相當(dāng)于體力的等效結(jié)點(diǎn)力。將慣性力和阻力項代入上式，得

3、（10-3）上式右邊第一、第二項即為一般靜力分析時的等效結(jié)點(diǎn)力，令 （10-4） （10-5） （10-6）、分別為慣性力和阻力對應(yīng)的等效結(jié)點(diǎn)力，考慮到積分與結(jié)點(diǎn)位移無關(guān)，于是 （10-7） （10-8）式中 （10-9）稱為單元質(zhì)量矩陣，而 （10-10）稱為單位阻尼矩陣。于是，單元動力方程化為 （10-11）若結(jié)構(gòu)系統(tǒng)有個單元，則對單元動力方程擴(kuò)大后相加，得 （10-12）式中稱為整個結(jié)構(gòu)的質(zhì)量矩陣，而稱為整個結(jié)構(gòu)的阻尼矩陣。2質(zhì)量矩陣式（10-9）所表示的單元質(zhì)量矩陣稱為一致質(zhì)量矩陣，因為形成質(zhì)量矩陣所用的形函數(shù)與位移插值函數(shù)所用的形函數(shù)相同，單元的動能和位能是相互協(xié)調(diào)的。此外，還可采用

4、一種較為簡單的質(zhì)量矩陣，它假定單元的質(zhì)量團(tuán)聚在單元的結(jié)點(diǎn)處，因此稱為團(tuán)聚質(zhì)量矩陣。單元的一致質(zhì)量矩陣是一滿陣，而團(tuán)聚質(zhì)量矩陣為一對角陣。以常應(yīng)變?nèi)切螁卧獮槔?，一致質(zhì)量矩陣團(tuán)聚質(zhì)量矩陣為 3阻尼矩陣式（10-10）所示的單元的阻尼矩陣還可進(jìn)一步表示為由此可見，當(dāng)阻力正比于速度時，阻尼矩陣正比于單元的質(zhì)量矩陣，這種阻尼稱為粘滯阻尼。此外，還可假定阻尼力正比于應(yīng)變速度，此時可見此時阻尼矩陣正比于單元的剛度矩陣。由材料內(nèi)摩檫引起的阻尼通?？杉僭O(shè)為這種情況。但結(jié)構(gòu)的阻尼往往是多重因素造成的，為了能更一般地描述結(jié)構(gòu)的阻尼，有關(guān)文獻(xiàn)建議將單元阻尼表示為式中為常系數(shù)，實際計算中，通常取前兩項，于是得 （10

5、-13）即認(rèn)為阻尼矩陣是質(zhì)量矩陣和剛度矩陣的線性組合，這種阻尼通常稱為比例阻尼。10-2 無阻尼自由振動為了計算動力方程，首先討論無阻尼自由振動，即考慮 的情況于是，動力學(xué)方程（10-12）變?yōu)?（10-14）這便是無阻尼自由振動方程，該方程組具有下列形式的解： （10-15）將上式代入（10-14）式，得 （10-16）記 可將式（10-16）改寫為 （10-17）求解方程組（10-17）的問題稱為廣義特征值問題，滿足方程（10-17）的解及其相應(yīng)的向量分別稱為特征值和特征向量。由此求得的就是結(jié)構(gòu)的固有頻率，而為相應(yīng)的固有振型。我們已經(jīng)知道，當(dāng)結(jié)構(gòu)上的約束足以消除剛體位移時，其總體剛度矩陣是

6、對稱正定的。由于材料的質(zhì)量密度，不難證明，一致質(zhì)量矩陣是正定的，而團(tuán)聚質(zhì)量矩陣是正定或半正定的對角陣，例如，對梁、板單元，當(dāng)考慮轉(zhuǎn)動自由度而忽略轉(zhuǎn)動慣量時，團(tuán)聚質(zhì)量矩陣因與轉(zhuǎn)動自由度相對應(yīng)的對角元為零而成為半正定（不為零的元素全為正值）。自由振動時，結(jié)構(gòu)中各點(diǎn)的振幅不可能全為零，因而齊次方程組必有其系數(shù)行列式等于零，即它稱為廣義特征方程，如果矩陣的階數(shù)為，則廣義特征方程是的次代數(shù)方程，由此可解出個廣義特征值，從而可得到個固有頻率。對每個固有頻率，由式（10-17）可確定維列向量，它代表個質(zhì)點(diǎn)的振幅構(gòu)成的振型。個特征值和相應(yīng)的個特征響量常稱為個特征對。求標(biāo)量和非零向量，使?jié)M足方程組的問題，稱為標(biāo)

7、準(zhǔn)特征值問題，方程（10-19）的解和稱為矩陣的特征值和特征向量。廣義特征值問題可以化為標(biāo)準(zhǔn)特征值問題，如果質(zhì)量矩陣是正定的，可以將進(jìn)行三角分解其中為下三角陣，方程（10-17）可寫為 用前乘上式，并注意到，則上式變?yōu)?若記 ， （10-20）就得到了式（9-17）所對應(yīng)的標(biāo)準(zhǔn)特征值問題。如果是正定的，則也是對稱正定的，且當(dāng)為對角陣時，是一個與相同的帶狀矩陣，在求得標(biāo)準(zhǔn)特征值問題的特征值和特征向量后，即可求得固有頻率，并利用（10-20）求出相應(yīng)的振型。 （10-21） （10-22）在總體剛度矩陣是正定的情況下，可以將進(jìn)行三角分解采用同樣的方法，也可將（10-15）化為如下的標(biāo)準(zhǔn)特征值問題

8、（10-23）其中 若是對稱正定的，則也是對稱正定的，但它一般不是帶狀的。在求得特征值和特征向量后，由式（10-24）可求得固有頻率和和相應(yīng)的振型。 （10-25） （10-26）目前標(biāo)準(zhǔn)特征值問題已有充分的研究，有許多現(xiàn)成的算法和程序可供使用。10-3 特征值和特征向量的性質(zhì)在研究特征值問題求解之前，首先介紹特征值和特征向量的性質(zhì)。1 實對稱陣的特征值問題 可以證明：對于和均為實對稱陣的廣義特征值問題，其所有的特征值都為實數(shù)，所有的特征向量也為實向量。 事實上，階對稱陣和的元素都是實數(shù)，且有 ，（）對廣義特征值問題，應(yīng)有即 用的共軛數(shù)乘上式的兩端，并對求和，得 （10-28）由于為對稱陣，因

9、此即式（10-28）的左端是實數(shù)，同理，式（10-28）的右端也是實數(shù)。因此，由式（10-28）可知必為實數(shù)，從而由（10-27）解得的特征向量也比為實向量。顯然，若為特征向量，則也必為特征向量。2 特征向量的正則化由上述可知，若為特征向量，則也必為特征向量，即相差一常數(shù)因子的一系列平行線均為特征向量，為確定起見，我們規(guī)定： （10-29）由此確定的特征向量稱為正則化特征向量。式（10-29）就為正則化條件。3 特征向量的正交性設(shè)和為廣義特征值問題的兩個特征對，則有分別用和前乘上兩式，得 （） （）（）式轉(zhuǎn)置減（），考慮到和均為實對稱陣，得若，且，則 （10-30）這就是特征向量關(guān)于的正交性，

10、將式（10-29）和（10-30）合并，可寫成 （10-31）其中為kronecker符號。將式（10-31）代入（）得 這說明特征向量對矩陣也具有正交性。假設(shè)是由個特征向量組成的特征向量矩陣，該矩陣稱為模態(tài)矩陣，即 （10-33）設(shè)是由特征值組成的對角陣，即 （10-34）則特征向量的正交性可表示為4 雷利（Rayleigh）商關(guān)于特征值的性質(zhì)還可以從雷利商出發(fā)去研究它，對于廣義特征值問題，雷利商定義為 （10-36）式中為任一向量，設(shè)它可以由特征向量的線性組合表示，即 （10-37）式中是由個系數(shù)組成的向量。將（10-37）代入（10-36），注意到特征向量的正交性，得展開后可以寫成 （1

11、0-38）將特征值由小到大排列：可見，雷利商在最小特征值和最大特征值之間，即并且，當(dāng)時這說明，當(dāng)取第階特征向量時，雷利商達(dá)到它的一個極值，該極值就是對應(yīng)的特征值，特別是這就是雷利商的極小值原理，次原理常被工程界作為預(yù)估振型的依據(jù)，用以計算特征值，從而求得基頻的近似值。如果，即向量與前階特征向量正交，則由式（10-38）看出：且當(dāng)時這就是說，是雷利商在向量與前階特征向量正交條件的極小值。5 移位 對于廣義特征值問題，可以使作一移位10-4 逆迭代法求解代數(shù)特征值問題的方法很多，這里介紹迭代法。對于階廣義特征值問題當(dāng)和均為正定陣時，我們已經(jīng)證明了它有個特征對。由個特征向量可以構(gòu)成模態(tài)矩陣同時，為了

12、說明迭代過程，令， （10-39）注意式（10-39）與式（10-34）中的區(qū)別，此時廣義特征方程變?yōu)橛媚B(tài)矩陣表示為 （10-40）式（10-40）兩邊右乘任意向量得 （10-41）若已知，將看為某一未知量，則由代入試式（10-41），求得 （10-42）即，經(jīng)過一次迭代后，對應(yīng)的特征向量的分量擴(kuò)大倍。如此經(jīng)過次迭代后得 （10-43）由于，當(dāng)?shù)螖?shù)較大時，有 （10-44）而前后兩次向量的模之比的平方為 （10-45）由此可見，迭代格式使向量逐漸與平行，前后兩次向量的模之比趨于，向量的收斂速率為，特征值的收斂速率為。實際上，這樣的迭代是無法進(jìn)行的，因為和正是需要求解的。為了克服這一矛盾，

13、并考慮到初向量的任意性，只要把作為初向量，把作為迭代后的向量，上述迭代過程就可進(jìn)行了。此時 （10-46）實際進(jìn)行迭代時還需要解決一個問題，從（10-43）可以看到，若，則的模趨于無窮大；若，則的模趨于零。不加處理在計算上是無法實施上述迭代的。為解決這一問題，只要在迭代過程中不斷使初向量的模取為1即可。逆迭代法的實施步驟概述如下：（1） 選取初向量，計算；（2） 由方程 解出；（3） 計算 （4） 計算近似特征值注意：此處，若滿足精度要求則轉(zhuǎn)向第五步，計算特征向量，否則計算下一次迭代的初向量：上式表明，不斷使迭代向量取正則化，然后轉(zhuǎn)向第（2）步進(jìn)行下一次迭代。（5） 計算特征向量取特征值為結(jié)束

14、迭代。這是基本迭代格式，當(dāng)時，收斂非?？臁＿@種迭代過程稱為逆迭代法。如果迭代方向反過來，即 也可求得最大特征值，這種迭代稱為冪法。逆迭代法除求最小特征值和特征向量外，還可以擴(kuò)大運(yùn)用范圍。和克萊姆-史密特（cram-schmidt）正交化過程相結(jié)合，可以用來求取最低的幾階特征對；也可求得重特征值對應(yīng)的特征向量；和移頻法相結(jié)合可以對奇異的剛度矩陣進(jìn)行迭代，可以改善某一近似特征值的精度和求取對應(yīng)的特征向量。現(xiàn)逐一敘述如下：由式（10-43）可見，假如，而，那么迭代結(jié)果應(yīng)該得到第階特征對。當(dāng)求得前階特征向量以后，要做到這一點(diǎn)很容易，只要選擇初向量與前階特征向量正交。雖然任意初向量不一定與正交，但可以構(gòu)

15、造一個新的初向量： （10-47）與前階特征向量正交，只要取 （10-48）這就是克萊姆-史密特（cram-schmidt）正交化過程。從迭代公式看，似乎這種正交化過程可以一勞永逸，一旦初向量與前階特征向量正交，則所有的迭代向量都與前階特征向量正交。但實際上并非如此，由于計算機(jī)誤差的原因，不可避免產(chǎn)生低階特征向量，迭代到最后還是得到最低階的特征向量。為了克服這一點(diǎn)，必須每迭代一次都進(jìn)行正交化處理，不斷把前階特征向量從迭代向量中清除掉。對于重特征值，用帶正交化過程的迭代，可得出重特征值對應(yīng)的一組正交特征向量。用逆迭代法求得第一特征對后，利用正交化過程的迭代法，可以依次求取各階特征對。為保證計算精

16、度，前面的特征向量必須計算精確些。逆迭代一般用來計算少數(shù)前幾階特征對。移頻逆迭代可以用來改善某幾階近似特征值的精度和求取相應(yīng)的特征向量。設(shè)已經(jīng)知道第階近似特征值，則對廣義特征值問題作以下移頻 （10-49）得到新的特征值問題 （10-50）式中新特征值問題與原特征值問題的特征向量相同，各階特征值移動一個距離。新特征值問題絕對值最小的特征值是，與它鄰近的兩個特征值是和。由此可知，對（10-50）進(jìn)行逆迭代，將得到和，從而得到。其收斂率為和中較大的一個，即和中較大的一個，所以只要比較接近收斂便很快。最后還要提一下初向量怎么選取的問題。初向量選取得好壞，對迭代效率影響較大，顯然。要求初向量選得應(yīng)盡量

17、接近所求的振型。一個簡單的選取辦法是，把結(jié)構(gòu)簡單地考慮為各個自由度是互不耦合的，認(rèn)為特征方程便為非耦合形式所以，將從小到大排序，得各階近似特征向量。設(shè)排在第位，記，。一般取第一階初向量為 取第階初向量為10-5 振型疊加法在10-5中我們已經(jīng)導(dǎo)出了系統(tǒng)的動力學(xué)方程（10-12）在一定的初始條件下求解上述方程，得出任一時刻的位移、速度和加速度響量，這就是結(jié)構(gòu)動力響應(yīng)的計算。數(shù)學(xué)上方程（10-51）是一個二階線性微分方程，原則上可用解常系數(shù)微分方程組的標(biāo)準(zhǔn)方法求界。但是如果矩陣的階數(shù)很高，求解將要付出高的代價，而有限元分析中，有的方法是行之有效的。這里我們將介紹兩種動力響應(yīng)的計算方法：振型疊加法和

18、逐步積分法。本節(jié)先介紹振型疊加法。振型疊加法的基本思想是，利用系統(tǒng)無阻尼自由振動的振型將動力學(xué)方程轉(zhuǎn)換為一組互不耦合的微分方程，然后分別求解這些方程，并將結(jié)果疊加而得到方程（10-51）的解。首先解出方程（10-51）的所對應(yīng)的無阻尼自由振動的固有頻率和振型。設(shè)結(jié)構(gòu)被激發(fā)起個振型，其相應(yīng)的特征對為。其次，將位移矢量看成是各振型的線性組合，即引入變換 （10-52）式中為振型矩陣，將此變換方程代入（10-51），兩端前乘，由的正交性，得 （10-53）式中 當(dāng)阻尼取比例阻尼時，即則為方便起見，設(shè) （10-54）此時，方程（10-53）就變?yōu)榛ゲ获詈系亩A常微分方程，即 （10-56）此時，初始條

19、件化為 其分量形式為由式（10-56）及其初始條件，利用杜哈美爾(Duhamel)積分得式中，而可由初始條件確定。于是方程（10-51）的解為 （10-59）關(guān)于比例阻尼系數(shù)的確定：待定系數(shù)的一般由實驗確定，設(shè)和分別為第和第個固有頻率，和分別為第和第個振型對應(yīng)的阻尼比，即實際阻尼與臨界阻尼的比值，可由實驗測定。根據(jù)式（10-54）得 由此解得 （10-60a） （10-60b）上式中一般取第一和第二階特征值和相應(yīng)的阻尼比進(jìn)行計算。在求得后，其它阻尼比可由下式求得 由上式可以看出，比例阻尼中高階阻尼很高，這實際上就等于在計算中將高階振型的響應(yīng)消除了。事實上，高階振型的響應(yīng)也是很小的，忽略它是合理

20、的。所以振型疊加法不必考慮全部振型，只需取少數(shù)最低幾階振型即可。 振型疊加法適用于象地震等只激發(fā)起較少振型，所需計算時間較長這類問題。例：考慮一個兩個自由度系統(tǒng)，其運(yùn)動方程為 （a）初始條件是 當(dāng)，試用振型疊加法求系統(tǒng)的位移響應(yīng)。解：該系統(tǒng)對應(yīng)的廣義特征值問題是 （b）方程存在非零解的條件為系數(shù)行列式等于零，即 （c）上式是關(guān)于特征值的特征方程，令，則特征方程化為 (d)由上式求得 (e) (f)將(e)代入（b）得顯然方程的系數(shù)行列式為零，設(shè)代入，則得，于是為系數(shù)因子，無論為何值，方程均滿足，為了唯一確定特征響量，根據(jù)正則化條件：于是，得 (g)將(f)代入（b），并由特征向量的正則化條件，

21、同理可得 (h)利用(g)和(h)將原問題轉(zhuǎn)換為互不耦合的運(yùn)動方程 (i) (j)初始條件化為 應(yīng)用杜哈美爾積分式，可求得方程(i)和(j)的解為最后將兩個振型疊加，10-6 逐步積分法逐步積分法根據(jù)動力學(xué)方程（10-51），引進(jìn)某些假設(shè)：任一時刻動力學(xué)方程都成立；于是建立由時刻結(jié)構(gòu)狀態(tài)響量到時刻狀態(tài)響量的遞推關(guān)系，從而從時刻的初始狀態(tài)向量出發(fā)，一步一步第求出各時刻的狀態(tài)向量。在時刻之間的位移、速度和加速度采用某種假設(shè)，根據(jù)假設(shè)的不同，可以有各種個樣的方法。我們這里介紹紐馬克(Newmark)法和Wilson-法。1 Newmark法 建立逐步積分法的關(guān)鍵是建立由時刻的狀態(tài)向量的遞推關(guān)系。時刻

22、有三組未知量，它們滿足動力學(xué)方程（10-51），即 （10-62）方程有三個未知量，顯然還應(yīng)補(bǔ)充二組方程才能求解。這兩個補(bǔ)充方程由速度和位移的泰勒級數(shù)展開采用某種近似得到。Newmark法取速度的一次展開式式中是在區(qū)間中某點(diǎn)的值。Newmark法取近似假設(shè) 于是有 （10-63）位移的二次展開式并對取類似的假設(shè) 于是又有 （10-64）式（10-62）、（10-63）、（10-64）便是Newmark法的基本公式，由此可以建立從時刻到時刻狀態(tài)向量的遞推關(guān)系。由（10-64）可得 （10-65）將（10-63）和（10-65）代入（10-62）得到關(guān)于的線性方程組，從而可以解出，然后將其代入（1

23、0-64）求出加速度向量；將代入（10-63）求出。參數(shù)的選擇對算法的影響很大，算法穩(wěn)定性分析指出，當(dāng)時，Newmark法是無條件穩(wěn)定的，這時可只根據(jù)精度的要求選擇步長。取時，Newmark法即為平均加速度法。根據(jù)以上分析，綜合計算步驟如下：1 初步計算（1） 形成剛度矩陣，質(zhì)量矩陣和阻尼矩陣；（2） 獲得初始狀態(tài)向量；（3） 選擇步長，以及參數(shù)，（），并計算下列常數(shù)， ，（4） 形成“剛度”矩陣 （5） 分解矩陣 2 對每個時間步計算（1） 計算時刻的“載荷”向量 （2）求時刻的位移（2） 計算時刻的加速度、速度 二、Wilson-法Wilson-法是線性加速度法的推廣，它把線性加速度假定的范圍從擴(kuò)充到。Wilson-法作如下的假設(shè) （10-66）在此假設(shè)下，積分便有速度和位移的表達(dá)式 （10-67） （10