張繼鋒鄧子辰，2?徐方暖張凱

（1.西北工業(yè)大學力學與土木建筑學院，西安 710072）

（2.大連理工大學工業(yè)裝備結構分析國家重點實驗室，大連 116023）

一種新的改進精細直接積分法*

張繼鋒1鄧子辰1，2?徐方暖1張凱1

（1.西北工業(yè)大學力學與土木建筑學院，西安 710072）

（2.大連理工大學工業(yè)裝備結構分析國家重點實驗室，大連 116023）

針對結構動力方程轉化為狀態(tài)空間方程后矩陣維數(shù)增加而導致計算量增大的問題，考慮狀態(tài)空間方程中所含外部荷載的特點，提出了一種新的改進精細直接積分法.給出了利用梯形公式、復化梯形公式、辛普生公式、復化辛普生公式、科特斯公式、高斯公式計算杜哈姆積分時的計算格式，分析了不同計算格式下的計算精度和計算效率.數(shù)值算例表明本文改進方法的正確性.

結構動力方程，直接積分，分塊計算，精細積分，改進方法

引言

結構動力方程廣泛應用于航空、航天、航海等多個領域，其求解一直備受關注.由鐘萬勰提出的精細積分法［1-2］，為結構動力方程的求解提供了一種高精度方法，并在各個領域得到了廣泛應用［3-5］.由于矩陣求逆有計算量大和穩(wěn)定性差的缺點，為了避免矩陣求逆，張森文等［5］利用辛普生積分公式計算杜哈姆積分，提出了狀態(tài)方程直接積分法.儲德文等［6］進一步討論了積分方法的選擇，指出科特斯積分和高斯積分是精度的較高，是較好的積分方法.后來一些學者嘗試用不同的數(shù)值積分方法來求解杜哈姆積分［7-9］，取得了比較好的結果.由于結構動力方程轉化為狀態(tài)空間方程后矩陣維數(shù)增加，并且數(shù)值積分都是利用插值方法，因而存在插值點過多導致計算量增大的問題，計算時間較長，效率需要提高.為了提高計算效率，有很多學者做了許多非常有益工作［10-12］.

本文在文獻［10］的基礎上，針對結構動力方程轉化為狀態(tài)空間方程后矩陣維數(shù)增加而導致計算量增大的問題，考慮狀態(tài)空間方程中外部荷載的特點，提出一種新的改進精細直接積分法.給出利用梯形公式、復化梯形公式、辛普生公式、復化辛普生公式、科特斯公式、高斯公式計算杜哈姆積分時的計算格式，同時分析不同計算格式下的計算精度和計算效率.

1　精細直接積分法的簡化計算

1.1精細直接積分法

采用集中質量法或有限元法進行動力學分析時，有如下所示的結構動力方程：

其中A和F（t）分別為：

利用常微分方程理論，對式（2）的非齊次方程求解可以采用如下迭代格式：

直接積分法一般是對式（3）的杜哈姆項數(shù)值積分，其間出現(xiàn)的exp（AΔt）可用精細計算［2］得到.

1.2精細直接積分法簡化計算［10］

式（3）的計算過程中將會出現(xiàn)如下矩陣運算：

對T和F（t）進行如下分塊：

下面通過式（6）進行簡化計算：

（1）當用梯形公式對杜哈姆項進行計算時，可得：

（3）當用辛普生公式對杜哈姆項進行計算時，可得：

（5）當用科特斯公式對杜哈姆項進行計算時，可得：

2　改進精細直接積分法

當n=3時上式中的參數(shù)為：w1=8/9，y1=0，

2.1改進精細直接積分法及計算格式

一般進行動力學分析時不一定每個質點上都作用有外荷載，所以外荷載可能出現(xiàn)有零項，因而可以通過對矩陣分塊，去掉對應的零項，將外荷載是零的部分行不參加矩陣運算.

在矩陣運算以前將外荷載是零的行對應的部分重新分塊，根據(jù)外荷載有零的特點，假設外荷載不為零的部分為m行，可以在式（6）的基礎上繼續(xù)分塊如下：

去掉零項對應的部分，則相應的式（6）可以改寫為：

下面分別給出梯形公式、復化梯形公式、辛普生公式、復化辛普生公式、科特斯公式、高斯公式計算時的改進計算格式，這些計算格式都是以式（13）和式（14）為基礎進行計算，可以不同程度的提高計算效率.

（1）當用梯形公式對杜哈姆項進行計算時，利用式（14）分塊計算可得：

（2）當用復化梯形公式對杜哈姆項進行計算時，利用式（14）分塊計算可得：

（3）當用辛普生公式對杜哈姆項進行計算時，利用式（14）分塊計算可得：

（5）當用科特斯公式對杜哈姆項進行計算時，利用式（14）分塊計算可得：

2.2改進算法精度分析

對于計算精度來講，對比計算式（6）和（14）可

當n=3時上式中的參數(shù)為：w1=8/9，y1=0，以發(fā)現(xiàn)二者的差別在于對矩陣分塊后去掉了零項，因此式（15）～式（20）的改進過程中不會有計算精度的變化，所以改進方法仍能保持原有算法的精度.

2.3改進算法效率分析

通過式（14）的改進可以將2n×n次乘法變?yōu)?n×m次乘法.本改進方法節(jié)約的時間會隨著m的變化而變化，其極限情況就是外荷載只有一行不為零，只有一行參加計算，即最大能將矩陣與矩陣的相乘變?yōu)榫仃嚭蛿?shù)字相乘，計算效率可以提高.

只考慮矩陣乘法計算，定量分析一個迭代步內(nèi)的各計算格式下的計算效率，以下分別給出各計算格式減少時間的預估公式：

（1）梯形公式對杜哈姆項進行計算

式（7）一個迭代步內(nèi)乘法計算量為：

利用式（15）改進后一個迭代步內(nèi)乘法計算量為：

改進后和改進前二者比值為：

上式中m的范圍為1～n，故改進方法最多可節(jié)約大約33%的時間.

（2）辛普生公式對杜哈姆項進行計算

式（9）一個迭代步內(nèi)乘法計算量為：

利用式（17）改進后一個迭代步內(nèi)乘法計算量為：

改進后和改進前二者比值為：

上式中m的范圍為1～n，故改進方法最多可節(jié)約大約50%的時間.

（3）科特斯公式對杜哈姆項進行計算

式（11）一個迭代步內(nèi)乘法計算量為：

利用式（19）改進后一個迭代步內(nèi)乘法計算量為：

改進后和改進前二者比值為：

上式中m的范圍為1～n，故改進方法最多可節(jié)約大約67%的時間.

（4）高斯積分公式對杜哈姆項進行計算

式（12）一個迭代步內(nèi)乘法計算量為：

利用式（20）改進后一個迭代步內(nèi)乘法計算量為：

改進后和改進前二者比值為：

上式中m的范圍為1～n，故改進方法最多可節(jié)約大約60%的時間.

3　算例分析

文獻［10］已經(jīng)說明了改進方法保持了計算原有算法的計算精度，為了驗證本文改進方法的計算效率，選取如下結構動力方程：

其中：

荷載列向量：

取n=1000，在相同步長情況下，利用本文方法和文獻［10］方法對算例進行了計算，計算十次平均所需時間的比較結果列于表1.

表1　數(shù)值計算效率的比較（單位：秒）Table 1 Comparison of numerical calculation efficiency（unit：second）

從表1可得本文方法提高了計算效率，針對本算例，改進梯形公式可以節(jié)約大約32%的時間，改進辛普生公式可以節(jié)約大約49%的時間，改進科特斯公式可以節(jié)約大約66%的時間，改進高斯公式可以節(jié)約大約59%的時間.算例計算和前文的時間預估公式保持一致，也證明了時間預估公式的有效性.

4　結論

本文考慮狀態(tài)空間方程外荷載的特點，提出了一種新的改進精細直接積分法.推導了利用梯形公式、復化梯形公式、辛普生公式、復化辛普生公式、科特斯公式、高斯公式計算杜哈姆積分時的計算格式，給出了改進梯形公式、改進辛普生公式、改進科特斯公式、改進高斯公式的節(jié)約時間公式，證明了改進梯形公式最大可以節(jié)約33%的時間，改進辛普生公式最大可以節(jié)約50%的時間，改進科特斯公式最大可以節(jié)約67%的時間，改進高斯公式最大可以節(jié)約60%的時間.數(shù)值算例表明本文改進方法提高了計算效率.

1 Zhong W X，Williams.A precise time step integration method.JournalofMechanicalEngineeringScience，1994，208：427～430

2 鐘萬勰.應用力學的辛數(shù)學方法.北京：高等教育出版社，2006（Zhong W X.Symplectic solution methodology in applied mechanics.Beijing：Higher Education Press，2006（in Chinese））

3 鐘萬勰，林家浩，高強.分層介質中非平穩(wěn)隨隨機波的精細求解.動力學與控制學報，2003，1（1）：1～8（Zhong W X，Lin J H，Gao Q.Precise computation of non-stationary random waves in stratified materials.Journal ofDynamics and Control，2013，1（1）：1～8（in Chinese））

4 朱寶，鐘萬勰.半解析高階有限譜元法及其在波導介質層PBG結構濾波器優(yōu)化設計中的應用.動力學與控制學報，2014，12（4）：289～294（Zhu B，Zhong W X.High order semi-analytical spectral element method and its application in optimal design of PBG structure in waveguide filter.Journal of Dynamics and Control，2014，12（4）：289～294（in Chinese））

5 張森文，曹開彬.計算結構動力響應的狀態(tài)方程直接積分法.計算力學學報，2000，17（1）：94～97（Zhang S W，Cao K B.Direct integration of state equation method for dynamic response of structure.Chinese Journal of Computational Mechanics，2000，17（1）：94～97（in Chinese））

6 儲德文，王元豐.精細直接積分法的積分方法選擇.工程力學，2002，19（6）：115～119（Chu D W，Wang Y F.Integration formula selection for precise direct integration method. Engineering Mechanics，2002，19（6）：115～119（in Chinese））

7 Wang M F，Zhou X Y.Modified precise time step integration method of structural dynamic analysis.Earthquake Engineering and Engineering Vibration，2005，4（2）：287～293

8 高小科，鄧子辰，黃永安.基于三次樣條插值的精細積分法.振動與沖擊，2007，26（9）：75～77（Gao X K，Deng Z C，Huang Y A.A high precise direct integration base on cubic spline interpolation.Journal of Vibration and Shock，2007，26（9）：75～77（in Chinese））

9 富明慧，廖子菊，劉祚秋.結構動力方程的樣條精細積分法.計算力學學報，2009，26（3）：379～384（Fu M H，Liao Z J，Liu Z Q.Spline precise time-integration of structural dyanmic analysis.Journal of Computational Mechanics，2009，26（3）：379～384（in Chinese））

10 張繼鋒，鄧子辰，胡偉鵬.結構動力方程精細直接積分的簡化計算.動力學與控制學報，2008，6（2）：107～111（Zhang J F，Deng Z C，Hu W P.Simplified computation of precise immediate integration method for structure dynamic equation.Journal of Dynamics and Control，2008，6（2）：107～111（in Chinese））

11 徐建新，郭巧榮，卿光輝.可分型指數(shù)矩陣的快速精細積分法.動力學與控制學報，2010，8（1）：24～28（Xu J X，Guo Q R，Qing G H.A fast precise integration method for the separable exponential matrix.Journal of Dynamics and Control，2010，8（2）：24～28（in Chinese））

12 高強，吳鋒，張洪武，林家浩，鐘萬勰.大規(guī)模動力系統(tǒng)改進的快速精細積分方法.計算力學學報，2011，28（4）：493～498（Gao Q，Wu F，Zhang H W，Lin J H，Zhong W X.A fast precise integration method for largescale dynamic structures.Chinese Journal of Computational Mechanics，2011，28（4）：493～498（in Chinese））

ANEW IMPROVED PRECISE DIRECT INTEGRATION METHOD*

Zhang Jifeng1Deng Zichen1，2?Xu Fangnuan1Zhang Kai1
（1.School of Mechanics，Civil Engineering&Architecture，Northwestern Polytechnical University，Xi′an 710072，China）
（2.State Key Laboratory of Structural Analysis of Industrial Equipment，Dalian University of Technology，Dalian 116023，China）

Considering the characteristics of non-homogeneous external loadings in the state space formula conversed from structural dynamic equation，a new improved direct precise integration method is presented in this paper.Meanwhile，different improved formats for Duhamel integration are proposed based on the Trapezium formula、compound Trapezium formula、Simpson formula、compound Simpson formula、Cotes formula and Gauss formula. The precision and the efficiency of the calculation are examined under the different formats.The results of the case study show the validity of the improved method.

structural dynamic equation，direct integration，partitioning calculation，precise integration，improved method

13 January 2015，revised 21 April 2015.

E-mail：dweifan@nwpu.edu.cn

10.6052/1672-6553-2015-028

2015-01-13收到第1稿，2015-04-21收到修改稿.

*國家自然科學基金重點項目（11432010）、國家基礎研究973項目（2011CB610300）、111引智計劃項目（B07050）、高校博士點基金（20126102110023）及西北工業(yè)大學基礎研究基金（310201401JCQ01001）資助

E-mail：dweifan@nwpu.edu.cn

*The project supported by the National Natural Science Foundation of China（11432010），the National Basic Research Program of China（2011CB610300），the Programme of Introducing Talents of Discipline to Universities（B07050），the Doctoral Program Foundation of Education Ministry of China（20126102110023）and the Fundamental Research Foundation of Northwestern Polytechnical University（310201401JCQ01001）