大規(guī)模動力系統(tǒng)改進(jìn)的快速精細(xì)積分方法

1、大規(guī)模動力系統(tǒng)改進(jìn)的快速精細(xì)積分方法高強(qiáng)，吳鋒，張洪武*自然科學(xué)基金(10902020, 10632030, 10721062, 2009AA044501)，遼寧省博士啟動基金(20081091)，遼寧省重點實驗室項目(2009S018)，鐵道部科技研究開發(fā)課題(2010T001-C)，大連理工大學(xué)青年教師培養(yǎng)基金。通訊作者：張洪武，大連理工大學(xué)工程力學(xué)系，電話：86 411 84706249; E-mail address: zhanghw，林家浩，鐘萬勰(大連理工大學(xué)，工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，工程力學(xué)系，遼寧大連 116024）摘要：本文提出一種針對大規(guī)模動力系統(tǒng)的改進(jìn)的快速精細(xì)積

2、分方法(FPIM)。以精細(xì)積分方法為基礎(chǔ)，利用大規(guī)模動力系統(tǒng)矩陣的稀疏性和動力問題的物理特性，分析了矩陣指數(shù)的特殊結(jié)構(gòu)，并基于此給出一種計算大規(guī)模動力系統(tǒng)矩陣指數(shù)及其動力響應(yīng)的高效率方法。關(guān)鍵詞：動力系統(tǒng)；稀疏矩陣；精細(xì)積分；矩陣指數(shù)；快速算法1 引 言對于線彈性結(jié)構(gòu)的動力學(xué)問題，比較成熟和常用的時域積分方法是Newmark方法1和Runge-Kutta方法2-4，由于計算穩(wěn)定性和精度的要求，這兩種方法的積分步長一般都取得比較小，計算量比較大。鐘萬勰提出了精細(xì)積分方法5-7，這種方法計算精度非常高，穩(wěn)定性好，允許采用很大的積分步長，特別是在處理剛性問題時具有明顯優(yōu)勢。精細(xì)積分方法提出后，得到了

3、很多發(fā)展8-10，但是這種方法的一個缺點是在處理規(guī)模很大的系統(tǒng)時，由于計算矩陣指數(shù)的計算量比較大，效率是一個主要問題。本文針對大規(guī)模動力系統(tǒng)提出一種計算其動力響應(yīng)的高效率算法。本文以精細(xì)積分方法為基礎(chǔ)5-7，利用動力問題的物理特性，利用一個關(guān)鍵思想，也就是結(jié)構(gòu)動力問題能量傳遞的有限性，來降低矩陣指數(shù)的計算量。利用上述思想，本文分析了大規(guī)模動力結(jié)構(gòu)對應(yīng)矩陣指數(shù)的稀疏性質(zhì)，并基于此給出一種計算矩陣指數(shù)的高效率方法。在高效和精確計算矩陣指數(shù)的基礎(chǔ)上，給出了大規(guī)模動力系統(tǒng)響應(yīng)的高效率和高精度算法。2動力系統(tǒng)的精細(xì)積分方法 假設(shè)系統(tǒng)的剛度矩陣、阻尼矩陣和質(zhì)量矩陣分別為，和，那么結(jié)構(gòu)動力學(xué)方程為其中為外力

4、。方程可以寫為狀態(tài)空間形式，即其中其中為動量。 數(shù)值積分時，將積分區(qū)間等分成步長為的積分區(qū)間，即若記則方程的解可以表示為其中 通過方程計算結(jié)構(gòu)的動力響應(yīng)，需要解決兩個主要問題，一是要準(zhǔn)確高效的計算方程定義的矩陣指數(shù)，二是要計算方程中的積分。對于常見的外載荷，方程中積分大多可解析積分，例如(1) 如果外力在一個積分時間步長內(nèi)為多項式，即，則其中上式中和分別為矩陣指數(shù)對應(yīng)的分塊矩陣，即(2) 如果外力在一個積分時間步長內(nèi)為簡諧變化，即，則 因此，上述計算過程的關(guān)鍵是計算矩陣指數(shù)。目前，計算矩陣指數(shù)最好的方法是精細(xì)積分方法5-7。但是，精細(xì)積分法計算矩陣指數(shù)的計算量比較大，即使對于為稀疏矩陣的情況，

5、也很難利用其稀疏性，得到的矩陣指數(shù)可能是滿陣。3 改進(jìn)的快速精細(xì)積分方法 精細(xì)積分方法基于兩個要點，一個是加法定理，另一個是增量計算和存儲。對于給定矩陣，它對應(yīng)的矩陣指數(shù)有如下性質(zhì)，即如果足夠大，則矩陣比較小，那么矩陣的矩陣指數(shù)可用如下的階Taylor級數(shù)近似，即精細(xì)積分方法5-7將的矩陣指數(shù)分為兩部分，即然后對增量部分應(yīng)用加法定理，即通過執(zhí)行次方程，然后將加上單位矩陣，則得到對應(yīng)的矩陣指數(shù)。上面簡要地介紹了精細(xì)積分方法，此方法編寫程序，計算精度非常高。但是，對于大規(guī)模系統(tǒng)，由于系統(tǒng)的自由度數(shù)目很大，計算矩陣指數(shù)將非常耗費時間和內(nèi)存。雖然，對于大規(guī)模動力系統(tǒng)，其剛度、阻尼和質(zhì)量矩陣是稀疏矩陣，

6、從而也是稀疏矩陣，但是直接通過以上的精細(xì)積分方法計算其矩陣指數(shù)，在計算過程，特別是方程的加法定理的合并過程中，矩陣將逐漸變?yōu)闈M矩陣或非常稠密的矩陣，很難利用矩陣的稀疏性質(zhì)，從而計算量很大。本文利用結(jié)構(gòu)動力響應(yīng)的物理特點，從物理上說明，大規(guī)模動力系統(tǒng)對應(yīng)的矩陣指數(shù)理論上應(yīng)該是稀疏矩陣，這樣在計算過程中可大大節(jié)約計算量，從而給出一種計算矩陣指數(shù)的高效算法，且可節(jié)省存儲要求。而原始精細(xì)積分方法之所以得到非常稠密的矩陣是因為計算誤差造成的。眾所周知的事實是能量在一維桿中的傳播速度是有限值，同理，在離散的結(jié)構(gòu)中，雖然其能量的傳播速度很難確切得到，但其動力學(xué)響應(yīng)的能量傳遞速度肯定是有限的，這將對矩陣指數(shù)的

7、結(jié)構(gòu)產(chǎn)生很大影響。根據(jù)方程，如果外力為零，則方程可寫為如下分塊形式，即由上述方程可以得到矩陣指數(shù)元素的物理意義，即：的第行第列元素表示第個自由度上給定單位位移，而其他自由度位移為零且所有自由度動量為零時，經(jīng)過一個積分步長后，第個自由度的位移響應(yīng)；的第行第列元素表示第個自由度上給定單位動量，而其他自由度動量為零且所有自由度位移為零時，經(jīng)過一個積分步長后，第個自由度的位移響應(yīng)；的第行第列元素表示第個自由度上給定單位位移，而其他自由度位移為零且所有自由度動量為零時，經(jīng)過一個積分步長后，第個自由度的動量響應(yīng)；的第行第列元素表示第個自由度上給定單位動量，而其他自由度動量為零且所有自由度位移為零時，經(jīng)過一

8、個積分步長后，第個自由度的動量響應(yīng)。由于能量在結(jié)構(gòu)中的傳遞速度是有限的，假設(shè)某個自由度上有擾動，在較小的時間內(nèi)，必然只能傳播到有限的自由度，而不可能傳播到所有自由度。根據(jù)上面給出的的物理含義，則它們一定是稀疏矩陣。這樣，既可以將矩陣指數(shù)作為稀疏矩陣計算和存儲，從而節(jié)約計算量和存儲空間。對于給定矩陣和積分步長，原始的精細(xì)積分方法按如下過程計算矩陣指數(shù)：首先確定，然后令，其次按照方程計算，然后執(zhí)行N次方程，最后將與單位矩陣相加得到矩陣指數(shù)。若采用上述的計算過程，雖然矩陣為稀疏矩陣，但是如果觀察根據(jù)方程計算得到的矩陣，可以發(fā)現(xiàn)它可能是一個滿矩陣或很不稀疏的矩陣，但仔細(xì)觀察會發(fā)現(xiàn)，此矩陣有很多很小的元

9、素。另一方面，根據(jù)上面的分析，此時的矩陣相當(dāng)于很小的時間步長對應(yīng)的矩陣指數(shù)減去一個單位矩陣，因此按照上面分析的能量傳播的性質(zhì)，它一定是非常稀疏的矩陣。因此，此時矩陣中非常小的非零元素是計算中的誤差，它們理論上應(yīng)該為零，應(yīng)該將它們判斷出來，并令它們?yōu)榱?，從而將轉(zhuǎn)換為稀疏矩陣存儲。具體過程如下：從上面分析給出的矩陣指數(shù)的物理意義可知，矩陣的四塊子矩陣的物理意義不同，因此我們將分為四塊，假設(shè)為矩陣中絕對值最大的元素，并給定一個誤差標(biāo)準(zhǔn)，如，則中的元素如果滿足其絕對值小于，則認(rèn)為此元素應(yīng)該為零，根據(jù)此原則可將稀疏存儲。同樣，可將稀疏化。以上過程定義為一個矩陣的稀疏化。同理，方程給出的加法定理同樣有上述

10、類似的現(xiàn)象，因此，若直接應(yīng)用方程，幾次合并后，矩陣將變?yōu)闈M陣或稠密矩陣，同樣可進(jìn)行上述的矩陣稀疏化。對于大規(guī)模動力系統(tǒng)，對于一個合理的時間步長，某一處的擾動經(jīng)過一個時間步長后，其影響只是局部化的，不會傳播到整個結(jié)構(gòu)，因此系統(tǒng)對應(yīng)的矩陣指數(shù)一定是稀疏矩陣，則可將精細(xì)積分方法作如下改進(jìn)，即對于給定矩陣和積分步長，我們可給出如下的快速精細(xì)積分算法(FPIM)，來計算矩陣指數(shù)。1. 由于是稀疏矩陣，將按照稀疏矩陣存儲；2. 確定8,9；3. 計算；4. 計算；5. 將矩陣稀疏化； 6. 執(zhí)行如下語句，即For iter=1:NR=R*( R + 2*I );將R稀疏化;End7. 得到矩陣后，將其與單

11、位矩陣相加，即得到和對應(yīng)的矩陣指數(shù)。上述計算過程與原始精細(xì)積分方法相比，只是增加了矩陣的稀疏化過程，但是這樣的處理將極大地提高計算效率，具體比較請見數(shù)值算例。得到矩陣指數(shù)后，還須考慮外力作用的影響。根據(jù)上面分析矩陣指數(shù)結(jié)構(gòu)的相同思想，可以對外力部分采用完全相同的處理方法，由于基本思想完全相同，本文不作詳細(xì)介紹。4 數(shù)值算例算例1：考慮如圖1所示的彈簧質(zhì)量系統(tǒng)，若取，則系統(tǒng)包含2000個質(zhì)點，本文隨機(jī)選擇每個質(zhì)點的質(zhì)量和彈簧的剛度，結(jié)果分別如圖2和3。此結(jié)構(gòu)很容易寫出其剛度矩陣和質(zhì)量矩陣，而阻尼矩陣取為。假設(shè)每個質(zhì)點上都作用相同的外力，則系統(tǒng)的運動方程為其中 分別采用本文方法(FPIM)、原始精

12、細(xì)積分方法(PIM)、4階R-K方法和Newmark方法積分此問題，積分區(qū)間為。本文方法和原始精細(xì)積分方法的積分步長為，4階R-K方法采用Matlab的ODE45函數(shù)計算，并將絕對誤差和相對誤差均設(shè)置為，Newmark方法的積分步長為。以R-K方法為參考解，分別計算本文方法、原始精細(xì)積分方法和Newmark方法的相對誤差。位移和動量的相對誤差分別定義為其中和分別表示R-K方法給出的位移和動量，而和可取本文方法、原始精細(xì)積分方法或Newmark方法積分得到的位移和動量。本文方法積分到時，各個質(zhì)點的位移和動量如圖4所示。Newmark方法、本文方法和原始精細(xì)積分方法的相對誤差如圖5中。圖5表明，本

13、文方法和原始精細(xì)積分方法的精度都非常高，而本文方法和原始積分方法給出的結(jié)果非常接近，并沒有損失計算精度，這說明了本文方法的正確性。圖5同時表明Newmark方法采用步長積分，精度也沒有達(dá)到本文方法的精度。四種方法的計算時間如表1所示，可以看到，本文方法的計算效率要大大優(yōu)于原始的精細(xì)積分方法、Matlab的ODE45和Newmark方法。對于Matlab的ODE45和Newmark方法要達(dá)到較高的計算精度，必須采用小的多的積分步長，因此效率降低，Newmark方法要達(dá)到更高的精度，還必須減小步長，從而效率還要降低。按照前文的分析，由于矩陣指數(shù)中存在大量的零元素，原始精細(xì)積分方法效率較低的原因在于

14、浪費了大量零元素的乘法運算。圖 1 彈簧質(zhì)量系統(tǒng)圖 2 各質(zhì)點的質(zhì)量圖 3 各彈簧的剛度圖 4 本文方法給出的時各質(zhì)點的位移和動量圖 5算例1相對誤差比較表 1 算例1計算效率比較FPIMODE45PIMNewmarkCPU time(s)18.31162.7412.5440.4算例2：考慮如圖6所示的平面應(yīng)力的動力學(xué)問題，結(jié)構(gòu)由兩種材料組成，兩種材料的密度和泊松比相同，而楊氏模量不同，它們分別為有限元網(wǎng)格和節(jié)點編號如圖7所示，采用三節(jié)點三角形單元，質(zhì)量矩陣采用集中質(zhì)量矩陣，共有3200個單元，1681個節(jié)點，施加邊界條件后，共有3280個自由度。初始位移為零，各自由度上具有初始動量1，無外力

15、作用。分別采用本文方法(FPIM)、原始精細(xì)積分方法(PIM)、4階R-K方法和Newmark方法積分此問題，積分區(qū)間為。本文方法和原始精細(xì)積分方法的積分步長為，4階R-K方法采用Matlab的ODE45函數(shù)計算，并將絕對誤差和相對誤差均設(shè)置為，Newmark方法的積分步長為。以R-K方法為參考解，分別計算本文方法、原始精細(xì)積分方法和Newmark方法的相對誤差。位移和動量的相對誤差的定義與算例1相同。本文方法積分到時，各個質(zhì)點方向的位移和動量如圖8所示，而方向的位移和動量如圖9所示。Newmark方法、本文方法和原始精細(xì)積分方法的相對誤差如圖10。圖10表明，本文方法和原始精細(xì)積分方法的精度

16、都很高，并且本文方法和原始積分方法給出的結(jié)果非常接近，這再次說明了本文方法的正確性。圖10同時表明Newmark方法采用步長積分，相對誤差的精度僅達(dá)到約量級。四種方法的計算時間如表2所示，可以看到，本文方法的計算效率要優(yōu)于Matlab的ODE45和Newmark方法，比原始精細(xì)積分方法更是快了約220倍，與原始精細(xì)積分方法相比效率得到了巨大的提高。圖6 兩種材料組成的平面應(yīng)力問題圖 7 有限元網(wǎng)格和節(jié)點編號圖8 時刻各節(jié)點方向的位移和動量圖9 時刻各節(jié)點方向的位移和動量圖10算例2相對誤差比較表 2 算例2計算效率比較FPIMODE45PIMNewmarkCPU time(s)105.5107

17、0.923252.15319.2結(jié) 論本文提出了一種針對大規(guī)模動力系統(tǒng)的改進(jìn)的快速精細(xì)積分方法(FPIM)，利用大規(guī)模動力系統(tǒng)矩陣的稀疏性和動力問題能量傳遞有限的物理特性，表明對于合理的積分步長，矩陣指數(shù)具有稀疏結(jié)構(gòu)，并基于此給出一種計算大規(guī)模動力系統(tǒng)動力響應(yīng)的高效率方法。數(shù)值算例表明，本文方法僅需在原始精細(xì)積分方法基礎(chǔ)上進(jìn)行簡單修改，即可將原始精細(xì)積分方法的效率大幅度提高，并且效率也比常用的4階R-K方法和Newmark方法高。參考文獻(xiàn)1. N. M. Newmark. A method of computation for structural dynamics. ASCE Journal

18、 of engineering mechanics, 85(3): 67-94, 1959.2. Hairer E, Lubich C, Wanner G. Geometric Numerical Integration: Structure-Preserving Algorithm for Ordinary Differential Equations (Second Edition). New York: Springer, 2006.3. Hairer E, Nrsett S P, Wanner G. Solving Ordinary Differential Equations I-N

19、onstiff Problems (Second Edition). Berlin: Springer, 1993.4. Hairer E, Wanner G. Solving Ordinary Differential Equations II-Stiff and Differential- Algebraic Problems (Second Edition). Berlin: Springer, 1996.5. 鐘萬勰。結(jié)構(gòu)動力方程的精細(xì)時程積分法。大連理工大學(xué)學(xué)報，34(2)：131-136，1994。6. Zhong WX, Williams FW. A precise time s

20、tep integration method. Proceedings of the Institution of Mechanical Engineers. Part C. Journal of Mechanical Engineering Science 1994; 208(6): 427-430.7. Zhong WX. On precise integration method. Journal of Computational and Applied Mathematics 2004; 163(1): 59-78.8. Zhang HW, Chen BS, Gu YX. An ada

21、ptive algorithm of precise integration for transient analysis. ACTA Mechannica Solida Sinica 2001; 14(3): 215-224.9. 譚述君，吳志剛，鐘萬勰。矩陣指數(shù)精細(xì)積分方法中參數(shù)的自適應(yīng)選擇。力學(xué)學(xué)報，41(6)： 961-966，2009。 10. 孫雁。奇異系統(tǒng)矩陣的精細(xì)積分。上海交通大學(xué)學(xué)報，42(8)：1217-1225，2008。A Fast Precise Integration Method for Large-scale Dynamic StructuresGao Qiang, Wu Feng, Zhang Hongwu, Lin Jiahao, Zhong Wanxie (State Key Laboratory of Structural Analysis of Indu