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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

17、0.923252.15319.2結(jié) 論本文提出了一種針對(duì)大規(guī)模動(dòng)力系統(tǒng)的改進(jìn)的快速精細(xì)積分方法(FPIM)，利用大規(guī)模動(dòng)力系統(tǒng)矩陣的稀疏性和動(dòng)力問題能量傳遞有限的物理特性，表明對(duì)于合理的積分步長，矩陣指數(shù)具有稀疏結(jié)構(gòu)，并基于此給出一種計(jì)算大規(guī)模動(dòng)力系統(tǒng)動(dò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)動(dòng)力方程的精細(xì)時(shí)程積分法。大連理工大學(xué)學(xué)報(bào)，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é)報(bào)，41(6)： 961-966，2009。 10. 孫雁。奇異系統(tǒng)矩陣的精細(xì)積分。上海交通大學(xué)學(xué)報(bào)，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