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

王海波， 陳 晉， 李少毅

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

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

王海波， 陳 晉， 李少毅

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

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

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

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

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

1 高效精細積分單步法

討論非線性動力方程：

(1)

利用指數(shù)矩陣將其轉(zhuǎn)化為如下同解積分方程：

(2)

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

令

T=eHΔt

(3)

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

(4)

式(2)變?yōu)椋?/p>

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

(5)

其中

(6)

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

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

由文獻[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分析結(jié)果及前40 s位移絕對誤差(誤差的絕對值，用R()表示，在下文中除另加說明外，誤差均采用絕對誤差表示)時程曲線分別見表1、圖1。

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

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

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

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

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

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

(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，采用預(yù)估校正-辛?xí)r間子域法[15]及本文單步法，通過Matlab編程計算，前5 s分析結(jié)果及兩種方法計算時間比較分別見表3、表4。從表3及表4可知：本文單步法與預(yù)估校正-辛?xí)r間子域法的計算結(jié)果均具有較高的精度，雖然前者精度略低于后者，但前者計算效率卻遠高于后者，且隨著計算步數(shù)的增加，這種優(yōu)勢更加明顯。

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

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

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

表4 預(yù)估校正-辛?xí)r間子域法[15]與本文單步法計算時間比較

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

時間步數(shù)執(zhí)行10次的平均CPU時間/s預(yù)估校正-辛?xí)r間子域法[15]本文單步法10013.75340.105220025.34020.149350059.79210.2109

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

3 結(jié) 論

表5 v1、v2的數(shù)值結(jié)果比較

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

[1] 鐘萬勰. 結(jié)構(gòu)動力方程的精細時程積分法[J]. 大連理工大學(xué)學(xué)報, 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] 趙秋玲. 非線性動力學(xué)方程的精細積分法[J]. 力學(xué)與實踐, 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]. 應(yīng)用力學(xué)學(xué)報, 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] 張洵安,姜節(jié)勝. 結(jié)構(gòu)非線性動力方程的精細積分算法[J]. 應(yīng)用力學(xué)學(xué)報, 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] 裘春航,呂和祥,鐘萬勰. 求解非線性動力學(xué)方程的分段直接積分法[J]. 力學(xué)學(xué)報, 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]. 計算力學(xué)學(xué)報, 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]. 天津大學(xué)學(xué)報, 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] 李金橋,于建華. 非線性動力方程避免狀態(tài)矩陣求逆的級數(shù)解[J]. 四川大學(xué)學(xué)報(工程科學(xué)版), 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] 任傳波,賀光宗,李忠芳. 結(jié)構(gòu)動力學(xué)精細積分的一種高精度通用計算格式[J]. 機械科學(xué)與技術(shù), 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] 向宇,黃玉盈,袁麗蕓,等. 非線性系統(tǒng)控制方程的齊次擴容精細積分法[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] 王海波,余志武,陳伯望. 非線性動力分析避免狀態(tài)矩陣求逆的精細積分多步法[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項的精細積分方法在非線性微分方程數(shù)值求解中的應(yīng)用[J]. 計算力學(xué)學(xué)報, 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]. 工程力學(xué), 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] 李煒華,王堉,羅恩. 求解非線性結(jié)構(gòu)動力方程的預(yù)估校正-辛?xí)r間子域法[J]. 計算力學(xué)學(xué)報, 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

國家自然科學(xué)基金(50908230)

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

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

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

O322; TU311.3

A

10.13465/j.cnki.jvs.2017.15.024