結(jié)構(gòu)動力學(xué)中的常用數(shù)值方法

1、第五章 結(jié)構(gòu)動力學(xué)中的常用數(shù)值方法51結(jié)構(gòu)動力響應(yīng)的數(shù)值算法 當(dāng)c為比例阻尼、線性問題模態(tài)疊加最常用。但當(dāng)C無法解耦，有非線性存在，有沖擊作用（激起高階模態(tài)，此時模態(tài)疊加法中的高階模態(tài)不可以忽略）。此時就要借助數(shù)值積分方法，在結(jié)構(gòu)動力學(xué)問題中，有一類方法稱為直接積分方法最為常用。所識直接是為模態(tài)疊加法相對照來說，模態(tài)疊加法在求解之前，需要對原方程進(jìn)行解耦處理，而本節(jié)的方法不用作解耦的處理，直接求解。（由以力學(xué)，工程中的力學(xué)問題為主要研究對象的學(xué)者發(fā)展出來的）中心差分法的解題步驟1. 初始值計算（1） 形成剛度矩陣K，質(zhì)量矩陣M和阻尼矩陣C。（2） 定初始值，。（3） 選擇時間步長，使它滿足，并

2、計算 ，（4） 計算（5） 形成等效質(zhì)量陣（6） 對陣進(jìn)行三角分解2對每一時間步長（1） 計算時刻t的等效載荷 （2） 求解時刻的位移 （3） 如需要計算時刻t的速度和加速度值，則若系統(tǒng)的質(zhì)量矩陣和阻尼矩陣為對角陣時，則計算可進(jìn)一步簡化。紐馬克法的解題步驟1.初始值計算（1）形成系統(tǒng)剛度矩陣K，質(zhì)量矩陣M和阻尼矩陣C（2）定初始值，。（3）選擇時間步長，參數(shù)、。并計算積分常數(shù) ， ， ，（4）形成等效剛度矩陣 （5）矩陣進(jìn)行三角分解 2. 對第一時間步長（1）計算時刻的等效載荷 （2）求解時刻的位移 （3）計算時刻的加速度和速度 威爾遜-法的解題步驟1. 初始值計算 （1）形成系統(tǒng)剛度矩陣K，

3、質(zhì)量矩陣M和阻尼矩陣C（2）定初始值，。（3）選擇時間步長，并計算積分常數(shù) ， ， ，（4）形成等效剛度 （5）將等效剛度進(jìn)行三角分解 2.對每一個時間步長 （1）計算時刻的等效載荷 （2）求解時刻的位移 （3）計算在時刻的加速度、速度和位移 52 結(jié)構(gòu)動力響應(yīng)數(shù)值算法性能分析對公式（5.1）描述的線性系統(tǒng)結(jié)構(gòu)動力學(xué)問題，已經(jīng)有證明對整個多自由度的積分，等價于將模態(tài)分解后對單自由度的積分的結(jié)果進(jìn)行模態(tài)疊加，因此可以通過對單自由度問題的分析，來說明算法的特性，其中阻尼均假設(shè)為比例阻尼，這樣，模態(tài)分解后的單自由度結(jié)構(gòu)動力學(xué)方程為： （5-29）以下算法的性能分析，均將算法用于這個方程。521 算法

4、用于結(jié)構(gòu)動力學(xué)方程的有限差分表示將數(shù)值計算方法應(yīng)用于（5-29）, 即分別在相鄰的不同時刻應(yīng)用算法可得如下一般形式 （5-30）A為放大矩陣或稱逼近算子，為載荷逼近算子。 例如將Newmak方法應(yīng)用于方程（5-29）有： （5-33） ， 矩陣A的特征多項式為 （5-34）其中A1，A為該矩陣的兩個特征向量，分別為矩陣的跡的一半和矩陣的行列式 （5-35） （5-36） （5-37）對Newmak方法有：, （5-38）其中 h為時間步長，。Newmak方法放大矩陣的規(guī)模是二維的，因此特征值也只有兩個，可以根據(jù)它們進(jìn)行分析。有的算法放大矩陣是三維的，例如Wilson-方法,在無阻尼情況下放大矩

5、陣為： （5-39）放大矩陣A的特征多項式為： （540）其中A1，A2，A3為該矩陣的三個特征向量，分別為矩陣的跡的一半、各階主子式的和以及矩陣的行列式，對Wilson-方法有 （5-41）此外，在幾個不同時刻應(yīng)用數(shù)值算法，然后將方程中的速度和加速度項消去，可得數(shù)值算法關(guān)于位移的差分方程，例如Newmak方法，有 （5-42）很顯然，其特征方程與其放大矩陣A的特征方程是相同的，使用關(guān)于位移的線性多步方式和放大矩陣來說明算法性能是一樣的，只不過各有方便之處。522 算法的穩(wěn)定性分析設(shè)為放大矩陣A的特征值，則定義為A的譜半徑，若特征值互異，則的算法是穩(wěn)定的，但若有重特征根，則要求。如果算法的穩(wěn)定

6、性要求對步長的選取有限制，稱算法是有條件穩(wěn)定的，反之為無條件穩(wěn)定的。放大矩陣的譜半徑小于等于1成立的充分條件是 （5-43）對的放大矩陣 （5-44）上兩式是關(guān)于算法自由參數(shù)的不等式，由它可以判斷算法是否無條件穩(wěn)定，若不是，將給出穩(wěn)定條件。 例51分析Newmak方法、Wilson-方法的穩(wěn)定性解： 將（5-38）代入（5-43）有顯然，當(dāng) （5-45）算法無條件穩(wěn)定。當(dāng)且 （5-46）算法穩(wěn)定，但為條件穩(wěn)定，其中為臨界采樣頻率。由于（5-43）式僅僅是充分條件，所以可進(jìn)一步按照穩(wěn)定性的定義得到節(jié)敘述的無條件穩(wěn)定條件。對Wilson-方法，將（5-41）代入（5-44）得 （5-47）容易看出

7、，其中第一，二，五不等式恒成立，對第三，四不等式若希望對任意的均成立，則有：求解上述不等式得 （5-48）實(shí)際使用中通常選取=1.4 算法的相容性和收斂性直接積分算法的相容性、收斂性分析同樣要使用其位移型的差分方程，或?qū)?yīng)的單步多值形式。在算法（5-30）式中，用精確解代替近似解，即可得到局部截斷誤差表達(dá)式，用符號表示 （5-49）局部截斷誤差表達(dá)式用放大矩陣的特征量以最常用的線性三步法為例可表示為 （5-50）其中分別為對應(yīng)的的放大矩陣的三個特征向量，然后將 在點(diǎn)進(jìn)行泰勒展開，然后利用運(yùn)動平衡方程化簡即可。若局部截斷誤差表達(dá)式為步長的s（s>0）階小量，則稱算法是s階相容的。對的放大矩

8、陣，可仿照上述步驟，來驗證算法的相容性。在經(jīng)典的數(shù)值算法收斂性分析理論中，一個重要的結(jié)論就是相容加穩(wěn)定等于收斂，其相容的階數(shù)就是算法的精度階。收斂性的含義也是當(dāng)時間步長趨于零，算法的數(shù)值解趨于精確解。對直接積分算法該定理同樣可以證明是成立的。例52分析Newmak法的相容性和精度解：其局部誤差仿照(550)式得： (5-51)也可以由(5-42)是直接求得，即將在點(diǎn)泰勒展開，并注意到在時刻的運(yùn)動方程有： （5-52）顯然，當(dāng)物理阻尼為零時，選擇算法是二階的，即截斷誤差是步長的二階小量。物理阻尼的存在，使算法精度降了一階，但若同時選擇，算法精度仍然是二階的，一般稱為Newmak線加速度法。顯然N

9、ewmak方法中有兩個參數(shù)待定，每種特定的選取都是一個特定的算法，最常用的幾個算法見表5-1表5-1： 常用的Newmak族直接積分算法方法名稱穩(wěn)定條件無阻尼問題精度階有阻尼問題精度階類型1/21/4平均加速度方法（梯形法）無條件21隱式1/21/6線性加速度方法22隱式1/20中心差分方法21顯式如果在一個時間步內(nèi)需要求解一個隱式的方程組，則稱算法是隱式的，反之不需要求解方程，直接計算即可得到下一時刻的值，則稱算法是顯示的。從5.1節(jié)的Newmak方法的計算步驟可以看出，這類方法是隱式的，但對于中心差分方法，若質(zhì)量矩陣和阻尼矩陣都是對角矩陣就可以顯示地計算。顯然顯示方法計算量要小得多。讀者可

10、自行分析Wilson-方法的精度，不難分析，無論是無阻尼還是有阻尼其精度都是2階的，它也是隱式方法。 算法耗散和彌散特性算法的精度，在小步長的情況下可以通過局部截斷誤差分析來說明比較，但是，在實(shí)際計算過程中，步長的選取可能不是很小，此時如何來度量算法的計算精度，當(dāng)然可以針對有解析解的問題進(jìn)行大量的數(shù)值計算，將數(shù)值解與解析解進(jìn)行比較來分析算法的計算精度。理論上還可以通過數(shù)值耗散(disspation)和彌散（dispersion）來輔助度量與分析，為引出這兩個概念的含義，我們?nèi)匀灰詥巫杂啥扔凶枘嶙杂烧駝訂栴}為例，該問題的解析解為： (5-53)當(dāng)直接積分算法用于這樣的問題，前小節(jié)已經(jīng)講述過它可以

11、寫成形如(5-42)的關(guān)于位移的有限差分形式，就是可以得到一個關(guān)于位移的有限差分方程，對于一個收斂的且有一定精度的算法，這個差分方程通常有一對共扼復(fù)根， (5-54)其中，該兩根稱為主根，其它根稱為寄生根(spurious roots)。解的一般形式可寫為 (5-55)式中的稱為算法阻尼比，當(dāng)有物理阻尼存在時，它還包括了物理阻尼的影響，稱為算法頻率，對應(yīng)的稱為算法周期?？梢钥吹缴鲜角皟身椗c(5-53)式形式是相同的，這給了我們與精確解進(jìn)行比較的可能，如果 1） 寄生根的影響較小，即 。2） 解表達(dá)式中得常數(shù)差別不太大。這樣，我們就可以通過比較不同算法的算法阻尼比和相對周期誤差，來比較算法的計算

12、精度。此時，對無阻尼問題，可以很明顯看到算法阻尼比會使得數(shù)值解曲線的幅值與解析解相比要降低而產(chǎn)生振幅衰減，這就是所謂的算法的數(shù)值耗散。同時，不同的算法的算法周期與精確的周期會有一定的誤差，這個誤差一般用相對周期誤差來表示，它會使得數(shù)值解曲線上產(chǎn)生周期的延長或縮短，即所謂的數(shù)值彌散。實(shí)際分析時，可以首先通過求解放大矩陣的特征方程得到特征方程的主根和寄生根，若主根可以表示為： (5-56)并注意到式(5-54)有 (5-57) (5-58) (5-59) (5-60)對結(jié)構(gòu)動力學(xué)問題，一般總希望算法在低頻段有較小的耗散和彌散，而且在低頻段，當(dāng)時，可以很方便地獲得它們的近似解析表達(dá)式。在高頻段算法的

13、耗散特性，用譜半徑來說明更適合，一般用來度量算法對高頻分量的耗散特性，特別地，當(dāng)時，稱算法具有高頻漸進(jìn)消去特性，即當(dāng)算法計算一步以后高頻極限完全地被耗散掉，而其它高頻分量由高到低漸進(jìn)地被耗散。由前面的敘述可以看到，算法放大矩陣的特征向量決定了算法對應(yīng)特征方程的根，也就決定了算法的穩(wěn)定性，同時確定了譜半徑、以及算法耗散和彌散特性。這些特性有時也稱算法的譜特性。放大矩陣相同的不同算法，稱為互相相同的，放大矩陣不同但特征值相同，稱算法互相相似，或稱算法頻譜等價的。對于算法的耗散特性，應(yīng)該說明的是高頻耗散特性對實(shí)際的結(jié)構(gòu)動響應(yīng)求解是有益的，因為實(shí)際結(jié)構(gòu)進(jìn)行有限元離散以后計算出的高頻行為并不真正代表系統(tǒng)

14、的物理行為，它是結(jié)構(gòu)系統(tǒng)在空間進(jìn)行有限元離散的結(jié)果，是虛假的行為，而不具備高頻耗散特性的算法，是將所有頻率上的響應(yīng)全部進(jìn)行了積分，盡管步長取得相對較大時，高頻的積分不準(zhǔn)確，這樣的計算結(jié)果顯然與系統(tǒng)實(shí)際的反應(yīng)不一致。但是同時要注意，它在低頻也不同程度地引入了數(shù)值耗散，這樣這些算法就不適合進(jìn)行長時間的計算，因為長時間以后應(yīng)該精確計算的低頻響應(yīng)，由于耗散特性的存在，已經(jīng)被耗散得面目全非，因此，有耗散特性的直接積分方法只適合計算瞬態(tài)的、短時間內(nèi)的低頻動力響應(yīng)。例5-3 分析Newmak族算法頻譜特性解：對Newmak族算法來說，當(dāng)時，算法才可能有二階精度，我們僅討論這一種情況。此時算法放大矩陣的兩個特

15、征量為： (5-61)則譜半徑不考慮物理阻尼時，對任意的有 (5-62) (5-63)就是無阻尼時，Newmak族算法不存在數(shù)值耗散，但有一定的相對周期誤差。對有阻尼問題， (5-64) 算法的超調(diào)特性譜半徑這個指標(biāo)對算法性能的影響還需要進(jìn)一步說明的是，它只決定算法的長期特性，即可以保證隨著算法計算步數(shù)的增加計算過程是數(shù)值穩(wěn)定的。但對無條件穩(wěn)定的算法，由于步長大小選擇沒有限制，一般在滿足指定精度的條件下，盡可能取較大的時間步長，對于非零初始條件問題，在計算開始的幾步可能會出現(xiàn)初始數(shù)據(jù)及其誤差（如初始位移，速度的測量誤差，初始加速度的計算誤差）被放大的現(xiàn)象，這稱為超調(diào)(overshoot)。這種

16、現(xiàn)象是放大矩陣A病態(tài)，有較大的條件數(shù)而產(chǎn)生的。實(shí)際應(yīng)用時，由于當(dāng)算法是收斂的，不會出現(xiàn)超調(diào)。一般為簡單起見，只分析當(dāng)時，在計算的第一步是否會出現(xiàn)超調(diào)。例54 分析Newmak平均加速度法的超調(diào)特性為分析簡便起見，將Newmak平均加速度法用于無阻尼自由振動問題，此時其放大矩陣為： 在時，可得近似等式 (5-65)其中分別表示關(guān)于的零次和一次關(guān)系式。由此可知算法在位移上無超調(diào)，但由于初位移的影響，在速度上有關(guān)于線性超調(diào)現(xiàn)象。前面提過Wilson-方法有很強(qiáng)的超調(diào)現(xiàn)象，對無阻尼問題，從放大矩陣的各元素的表達(dá)式中，很容易得到在時，可得近似等式 (5-66) (5-67)其中分別表示關(guān)于的二次和一次關(guān)

17、系式。由此可知Wilson-方法在位移上關(guān)于初位移有二次超調(diào)，同時關(guān)于初始速度有一次超調(diào)。在速度上有關(guān)于初位移一次超調(diào)。另外，也可以用數(shù)值計算的方法對指定的初始條件，計算出近似解，然后與精確解比較，或計算系統(tǒng)能量范數(shù)： (5-68)然后將比較。顯然，由于直接積分方法適合于短時間的瞬態(tài)問題計算，因此超調(diào)現(xiàn)象也是必須加以注意的。綜上所述，對一個數(shù)值積分算法理論上要分析其相容性，穩(wěn)定性，數(shù)值的耗散與彌散特性，對無條件穩(wěn)定的算法還要分析其超調(diào)特性。這樣才可能對算法的本質(zhì)有深入的了解，進(jìn)而指導(dǎo)數(shù)值計算結(jié)果的解釋與分析。此外，由于直接積分方法對結(jié)構(gòu)運(yùn)動平衡方程進(jìn)行數(shù)值積分的目的在于估計結(jié)構(gòu)真實(shí)的動力響應(yīng)。

18、為了精確地預(yù)計結(jié)構(gòu)的動力響應(yīng)，要求模態(tài)分解以后所有的單自由度平衡方程都必須被精確地積分，但在直接積分法中，對所有的方程積分都相當(dāng)于采用相同的步長h，所以時間步長的選擇必須針對系統(tǒng)的最小周期。如果是系統(tǒng)的最小周期的話，則h選為，其中一般n10。對條件穩(wěn)定的算法，當(dāng)然還要同時考慮這個選取是否滿足算法穩(wěn)定性的要求。對有條件穩(wěn)定的算法，要求，若n取10，多數(shù)的條件穩(wěn)定算法的穩(wěn)定條件都滿足，不難驗證表51中的條件穩(wěn)定算法全都滿足。但需要注意的是，對于大型、復(fù)雜的實(shí)際結(jié)構(gòu)，經(jīng)過有限元離散以后通常都有上萬，甚至幾十萬個自由度，其最大固有頻率通常都很大，也就是系統(tǒng)的最小周期非常小，此時，按來選取步長就非常小，

19、這會大大增加計算量。而實(shí)際工程上只關(guān)心較低階的固有頻率，同時結(jié)構(gòu)的響應(yīng)也主要由若干較低階的響應(yīng)構(gòu)成，因此在計算時高頻可以不用精確積分，就積分出那些主要的，感興趣的低頻響應(yīng)就可以了。也就是步長可選擇為，比大倍。由于實(shí)際情況中可能會非常大，這樣條件穩(wěn)定算法的穩(wěn)定條件就可能無法滿足，而無條件穩(wěn)定算法對步長的選取就沒有穩(wěn)定性的限制，因此對于實(shí)際的結(jié)構(gòu)動力響應(yīng)計算，多數(shù)都使用無條件穩(wěn)定算法。5.3 矩陣特征值問題及解法531 .問題分類 這原本是廣義特征值問題，但可以化為標(biāo)準(zhǔn)特征值問題，前乘，得 令，則有 這是代數(shù)里的標(biāo)準(zhǔn)的矩陣特征值問題，對結(jié)構(gòu)動力學(xué)問題而言 M對角 對稱 M對稱非對角 不一定對稱 非

20、對稱矩陣的結(jié)構(gòu)特征值問題計算量很大，此時作如下處理：1 對M進(jìn)行Cholesky分解（因為M對稱正定） L（下三角陣）（上三角陣）代入方程 前乘得 提出U得 令 ， 則有 可以征明與A有相同的特征值，因為它們是相似的矩陣，但是對稱的因為 所以對稱。 只不過特征向量有所變化，求出以后 至于Cholesky分解，可直接由M（或K）的元素計算得到 的各非零元素 i=j+1,j+2,.,n上式依次取j=1,2,.,n，即可求得L的下三角部分各列元素。又若記 則其下三角之i=1,2,.,n行元素，依次為 j=1,2,i-15.3.2 特征值、特征向量的一些特性1.對角陣、三角陣、塊對角陣、塊三角陣對角陣

21、和三角陣的特征值就是這些矩陣對角元素的數(shù)值。塊對角陣和塊三角陣的特征值就是這些矩陣的對角線上各個子塊的特征值。2.幾種特殊矩陣（1）實(shí)對稱矩陣的特征值必為實(shí)數(shù)，其特征向量也可選為實(shí)向量。（2）反對稱矩陣的特征值或者為純虛數(shù)，或者為零。（3）正定對稱矩陣的特征值全部大于零。（4）對稱矩陣對應(yīng)于不同特征值的特征向量彼此正交。3.與原矩陣A相關(guān)聯(lián)的幾種矩陣設(shè)矩陣A的特征值是，其對應(yīng)的特征向量是x，則（1）陣的特征值是，特征向量是x；（2）陣的特征值是（），特征向量x；（3）是A的轉(zhuǎn)置矩陣，則的特征值就是；（4）陣的特征值是，特征向量仍為x,m為整數(shù)；（5）若A非奇異，則A的逆矩陣存在，記為，是的特征

22、值為，特征向量為x；（6）若矩陣B與A相似，即有可逆陣P存在，使 則B的特征值也是，特征向量是4.特征值的積與積若矩陣A的特征值為、.，則有 它們可作為校核、估計甚至計算特征值的手段，讀者可以自行驗證。533 特征值問題的基本計算方法 由于高于四次的一元代權(quán)方程無法求精確解，所以必然要采用迭代方法求計算求近似解。目前的迭代方法，只不過迭代的方法和技巧不同。已有很多成熟方法。 方法選擇時主要取決于 1) 所要求的特征對數(shù)目 2）K、M的階，3）K、M的帶寬以及是否帶狀。 根據(jù)用到的基本關(guān)系大致可分四類 1.向量迭代法（冪迭代法Power iteration method） 基本關(guān)系 根據(jù)迭代格式

23、不同又分： 正向迭代 （冪迭代）逆迭代 （反冪迭代）主要用于求特征向量。 2.變換法 基本關(guān)系： 質(zhì)量歸化為的模態(tài)矩陣 雅可比迭代（Jacobi）（小模型實(shí)對稱陣標(biāo)準(zhǔn)特征值問題的全部特征對方法） 廣義雅可比迭代（求全部特征對，小模型的廣義特征值問題，或者有大量非對角線零元素和少量對角線零元素問題） 豪斯霍爾德（Householder）3.多項式迭代法（不單獨(dú)使用）基本關(guān)系： 顯式多項式迭代 隱式多項式迭代4.斯圖姆（Sturm）序列法 利用的Sturm序列性質(zhì)來求解。 實(shí)際工程結(jié)構(gòu)的動力學(xué)問題絕大多數(shù)使用有限元法方法離散求解。隨著計算機(jī)存儲和速度能力的成倍提高，所建模型越來越精確，幾十萬、上百萬個自由度的計算成為可能。為此發(fā)展了一些針對大型特征值問題的方法，它們綜合上述典型方法和技巧，常用的有：多項式迭代 行列式探索法（帶寬較窄的低階特征值問題。 子空間迭代法（Subspace iteration）(與行列式法相似，但計算量小) 蘭索斯（Lanczos）法向量迭代 HQRI(HouseholderQR 迭代方法)（針對大型滿陣的大多數(shù)或全部的特征值和向量）逆迭