用于非線性動力分析的一種高效精細積分單步法

王海波， 陳 晉， 李少毅

(中南大學 土木工程學院，長沙 410075)

用于非線性動力分析的一種高效精細積分單步法

王海波， 陳 晉， 李少毅

(中南大學 土木工程學院，長沙 410075)

非線性動力方程；精細積分法；Romberg數值積分；龍格庫塔法；單步法

傳統數值積分方法，如中心差分法、Wilson法、Newmark法、Houbolt法等，在求解非線性動力方程中有著廣泛的應用，但上述方法存在能量耗散及相位誤差。基于上述原因，為更好滿足工程實際的需求，求解非線性動力問題的高精度算法得到了發展。1994年前后，鐘萬勰等[1-2]提出精細積分方法，為非線性系統的時程分析開辟了嶄新的途徑。趙秋玲等[3-4]提出了一種非線性系統線性化的精細積分法，但計算結果精度不高；張洵安等[5]提出了一種精度較高的線性化迭代精細積分方法；裘春航等[6]將非線性部分用j次多項式來近似，然后借助于分段直接積分法進行求解，通過確定j值，可獲得一系列具有不同精度的近似解，但需對狀態矩陣求逆。

為避免對狀態矩陣求逆，克服病態矩陣或逆矩陣不存在的缺點，張素英等[7]采用增維的方法，將非齊次動力方程化為齊次動力方程，然后假設各時間段狀態矩陣為常矩陣，再利用精細積分求解；葛根等[8]對文獻[7]進行了改進，提高了增維精細積分法的計算效率，但計算結果精度不高；李金橋等[9]先將非線性部分在tk時刻展開成泰勒級數形式，然后將非線性積分項泰勒展開法或高斯-勒讓德數值積分法與精細積分法結合進行求解，一定程度上提高了計算結果的精度；任傳波等[10]提出了一種不論是線性激勵還是非線性激勵都具有較高求解速度和精度的計算格式，但僅討論了線性動力系統；向宇等[11]提出了一種較高精度的齊次擴容精細積分法，但顯著增加了計算量；王海波等[12]將精細積分法和Adams-Bashforth-Moulton多步法相結合，提出了一種高精度、高效率的精細積分多步法，但需要使用單步法計算表頭；譚述君等[13]基于Duhamel項精細積分方法，構造了幾種求解非線性微分方程的數值算法；江小燕等[14]根據Hamilton變作用定律構造時空有限元矩陣，并結合傳遞矩陣原理和精細積分思想提出了精細時空有限元方法；李煒華等[15]利用改進的歐拉法對狀態未知量進行預估及校正，然后結合辛時間子域法對非線性方程進行求解。

本文將精細積分法和Romberg數值積分結合，對計算過程中待求的vk+j/m(j=1,2,…，m)，利用當前時刻vk，通過二階龍格庫塔法進行預估，提出了一種避免狀態矩陣求逆的高效精細積分單步法。數值算例驗證了算法的有效性。

1 高效精細積分單步法

討論非線性動力方程：

(1)

利用指數矩陣將其轉化為如下同解積分方程：

(2)

式中，vk+1、vk分別表示所求解向量在時刻tk+1和tk的值。

令

T=eHΔt

(3)

G(τ)=eH(tk+1-τ)f(v,τ)

(4)

式(2)變為：

vk+1=Tvk+v*(tk+1)

(5)

其中

(6)

這樣，將vk+1的表達式分解成了2項，第1項可直接采用精細積分精確得到，下面討論第2項的數值求解方法。

本文采用Romberg積分對第2項進行數值直接積分。

由文獻[10]可知，式(6)的Romberg積分格式為：

(7)

(j=1,2,…,m)(m=1,2,3,…)

其中

(8)

(9)

式(7)中，m的取值要視問題的求解精度而定，m取值越大精度越高。本文取m=2對式(6)進行計算，由式(4)、式(7)～式(9)可以得到如下計算步驟：

(10)

(11)

(12)

(13)

(14)

(15)

T1/2=T1/4·T1/4，T=T1/2·T1/2

(16)

(17)

(18)

(19)

(20)

2 算例分析

各種方法的前5 s分析結果及前40 s位移絕對誤差(誤差的絕對值，用R()表示，在下文中除另加說明外，誤差均采用絕對誤差表示)時程曲線分別見表1、圖1。

表1 v1的數值結果比較(Δt=0.1 s)

從本例Δt=0.1 s計算結果(表1和圖1)可知：本文單步法計算誤差明顯小于一次預-校法[12]和單步法1[6](按文獻[6]中公式(23)計算)，隨著時間的增大，這種優勢越來越明顯。

圖1 三種方法解的誤差R(v1)

算例2 考慮帶參數a和b的Duffing方程[6]

取初值v1(0)=0.3，v2(0)=0.5，采用Matlab編程計算，圖2(a)、圖2(b)分別給出了當線性阻尼系數a較小，b=2、b=45時，前20 s內3種方法的位移誤差時程曲線(所取點為1 s的整數倍)。從圖2中可知：隨著激勵幅值b的增加，本文單步法的誤差有增加趨勢，但誤差值均較小；本文單步法與單步法2[6](按文獻[6]中公式(16)計算)相比，誤差稍大，但避免了狀態矩陣求逆且無需對非齊次項進行求導運算，計算效率較高，兩者計算時間比較見表2；本文單步法的誤差明顯小于一次預-校法[12]，計算精度較高。

(a) a=0.3,b=2,Δt=0.2 s

(b) a=0.3,b=45,Δt=0.05 s

表2 單步法2[6]與本文單步法計算時間比較

Tab.2 Comparisons of computation time for the single-step method 2 in Ref.[6]and proposed algorithm

時間步數執行10次的平均CPU時間/s單步法2[6]本文單步法1000.18410.132610000.60060.3744100003.81731.6770注:兩種算法的計算條件:a=0.3,b=2,Δt=0.2s

以前述初值為起點，計算到200 s為止，軌跡線從100 s開始繪制，如圖3所示。

由圖3(a)和圖3(b)可見，當暫態結束后，所有軌跡線都集中在一起，形成周期吸引子圖形；而圖3(c)則顯示了軌跡線仍分散的混沌吸引子。一般地說，當線性阻尼系數a較小時，隨著激勵幅值b的增加，周期解可能分支為對稱解，再經過一系列倍周期過程而進入混沌解。本例的收斂性和數值穩定性都較好，本文單步法在Δt=0.1 s和Δt=0.02 s算出的結果相近，對圖3(a)的問題，即使采用步長Δt=0.5 s，本文單步法和單步法2[6]也能算出誤差較小的結果，而一次預-校法[12]的結果已發散。

(a) a=0.3,b=2

(b) a=0.3,b=20

(c) a=0.3,b=45

算例3 受迫Van der Pol方程[15]為

取α=1，β=0.2，ω=1.2，f=1.0，時間步長Δt=0.1 s，初始條件為：v1(0)=1.5，v2(0)=0.5，采用預估校正-辛時間子域法[15]及本文單步法，通過Matlab編程計算，前5 s分析結果及兩種方法計算時間比較分別見表3、表4。從表3及表4可知：本文單步法與預估校正-辛時間子域法的計算結果均具有較高的精度，雖然前者精度略低于后者，但前者計算效率卻遠高于后者，且隨著計算步數的增加，這種優勢更加明顯。

表3 v1的數值結果比較(Δt=0.1 s)

算例4 考慮平方非線性的二自由度動力學方程[8]

取初值v1(0)=v2(0)=0.1，v3(0)=v4(0)=0，本文單步法在時間步長Δt=0.1 s、0.5 s、1.0 s前15 s分析結果見表5。

表4 預估校正-辛時間子域法[15]與本文單步法計算時間比較

Tab.4 Comparisons of computation time for the predictor-corrector symplectic time-subdomain algorithm in Ref.[15] and the proposed algorithm

時間步數執行10次的平均CPU時間/s預估校正-辛時間子域法[15]本文單步法10013.75340.105220025.34020.149350059.79210.2109

由表5可知：隨著時間步長的增大，本文單步法的精度逐漸降低；本文單步法在較大步長(Δt=1.0 s)情況下仍具有較理想的計算結果。

3 結 論

表5 v1、v2的數值結果比較

(2) 本文單步法是一種高精度、高效率和穩定性較好的方法。數值計算結果表明：取m=2時，本文單步法比單步法1[6]及一次預-校法[12]計算精度要高；與單步法2[6]、預估校正-辛時間子域法[15]相比，本文單步法在保持相同高精度計算結果的基礎上，計算效率更高，在求解多自由度、強非線性動力響應問題中具有較大優勢。

[1] 鐘萬勰. 結構動力方程的精細時程積分法[J]. 大連理工大學學報, 1994, 34(2): 131-136.

ZHONG Wanxie. On precise time-integration method for structural dynamics[J]. Journal of Dalian University of Technology, 1994, 34(2): 131-136.

[2] LIN Jiahao, SHEN Weiping, WILLIAMS F W. A high precision direct integration scheme for structures subjected to transient dynamic loading[J]. Computer & Structures, 1995, 6(1): 120-130.

[3] 趙秋玲. 非線性動力學方程的精細積分法[J]. 力學與實踐, 1998, 20(6): 24-26.

ZHAO Qiuling. An accurate integration method for solving nonlinear dynamic problems[J]. Mechanics in Engineering, 1998, 20(6):24-26.

[4] 呂和祥,蔡志勤,裘春航. 非線性動力問題的一個顯式精細積分算法[J]. 應用力學學報, 2001, 18(2): 34-40.

Lü Hexiang, CAI Zhiqin, QIU Chunhang. An explicit precise integration algorithm for nonlinear dynamics problems[J]. Chinese Journal of Applied Mechanics, 2001, 18(2): 34-40.

[5] 張洵安,姜節勝. 結構非線性動力方程的精細積分算法[J]. 應用力學學報, 2000, 17(4): 164-168.

ZHANG Xunan, JIANG Jiesheng. The precise integration algorithm for nonlinear dynamics equations of structures[J]. Chinese Journal of Applied Mechanics, 2000, 17(4): 164-168.

[6] 裘春航,呂和祥,鐘萬勰. 求解非線性動力學方程的分段直接積分法[J]. 力學學報, 2002, 34(3): 369-378.

QIU Chunhang, Lü Hexiang, ZHONG Wanxie. On segmented-direct-integration method for nonlinear dynamics equations[J]. Acta Mechanica Sinica, 2002, 34(3): 369-378.

[7] 張素英,鄧子辰. 非線性動力方程的增維精細積分法[J]. 計算力學學報, 2003, 20(4): 423-426.

ZHANG Suying, DENG Zichen. Incremental-dimensional precise integration method for nonlinear dynamic equation[J]. Chinese Journal of Computational Mechanics, 2003, 20(4): 423-426.

[8] 葛根,王洪禮,譚建國. 多自由度非線性動力方程的改進增維精細積分法[J]. 天津大學學報, 2009, 42(2): 113-117.

GE Gen, WANG Hongli, TAN Jianguo. Improved increment-dimensional precise integration method for the nonlinear dynamic equation with multi-degree-of-freedom[J]. Journal of Tianjin University, 2009, 42(2): 113-117.

[9] 李金橋,于建華. 非線性動力方程避免狀態矩陣求逆的級數解[J]. 四川大學學報(工程科學版), 2004, 36(4): 26-30.

LI Jinqiao, YU Jianhua. Series solution of nonlinear dynamics equations avoiding calculating the inversion of the state matrix[J]. Journal of Sichuan University(Engineering Science Edition), 2004, 36(4): 26-30.

[10] 任傳波,賀光宗,李忠芳. 結構動力學精細積分的一種高精度通用計算格式[J]. 機械科學與技術, 2005, 24(12): 1507-1509.

REN Chuanbo, HE Guangzong, LI Zhongfang. A high precision and general computational scheme in precise integration of structural dynamics[J]. Mechanical Science and Technology, 2005, 24(12): 1507-1509.

[11] 向宇,黃玉盈,袁麗蕓,等. 非線性系統控制方程的齊次擴容精細積分法[J]. 振動與沖擊, 2007, 26(12): 40-44.

XIANG Yu, HUANG Yuying, YUAN Liyun,et al. Extended homogeneous capacity high precision integration method for control equation of nonlinear systems[J]. Journal of Vibration and Shock, 2007, 26(12): 40-44.

[12] 王海波,余志武,陳伯望. 非線性動力分析避免狀態矩陣求逆的精細積分多步法[J]. 振動與沖擊, 2008, 27(4): 105-108.

WANG Haibo, YU Zhiwu, CHEN Bowang. Precise integration multi-step method for nonlinear dynamic equations to avoid calculating inverse of state matrix[J]. Journal of Vibration and Shock, 2008, 27(4): 105-108.

[13] 譚述君,高強,鐘萬勰. Duhamel項的精細積分方法在非線性微分方程數值求解中的應用[J]. 計算力學學報, 2010, 27 (5): 752-758.

TAN Shujun, GAO Qiang, ZHONG Wanxie. Applications of Duhamel term’s precise integration method in solving nonlinear differential equtions[J]. Chinese Journal of Computational Mechanics, 2010,27(5):752-758.

[14] 江小燕,王建國. 非線性動力方程的精細時空有限元方法[J]. 工程力學, 2014, 31(1): 23-28.

JIANG Xiaoyan, WANG Jianguo. The precise space-time finite element method for nonlinear dynamic equation[J]. Engineering Mechanics, 2014, 31(1): 23-28.

[15] 李煒華,王堉,羅恩. 求解非線性結構動力方程的預估校正-辛時間子域法[J]. 計算力學學報, 2014, 31(4): 453-458.

LI Weihua, WANG Yu, LUO En. Predictor-corrector symplectic time-subdomain algorithm for nonlinear dynamic equations[J]. Chinese Journal of Computational Mechanics, 2014, 31(4): 453-458.

An efficient precise integration single-step method for nonlinear dynamic analysis

WANG Haibo, CHEN Jin, LI Shaoyi

(College of Civil Engineering, Central South University, Changsha 410075, China)

nonlinear dynamic equations; precise integration method; Romberg numerical integration; Runge-Kutta method; single-step method

國家自然科學基金(50908230)

2016-03-14 修改稿收到日期：2016-06-15

王海波 男，博士，副教授，1974年4月生

陳晉 男，碩士生，1993年7月生

O322; TU311.3

A

10.13465/j.cnki.jvs.2017.15.024