基于雙步長的顯式積分算法

基于雙步長的顯式積分算法

楊超，肖守訥，魯連濤

(西南交通大學牽引動力國家重點實驗室， 成都610031)

摘要：為了建立精度更高、穩定性更好的顯式積分算法，在雙步長線性加速度假設的條件下，通過積分推導，提出了雙步長顯式法和修正雙步長顯式法。對這兩種新方法進行了穩定性、精度和截斷誤差分析，并通過算例對兩種新方法、中心差分法和解析解進行比較。結果表明，無阻尼情況下，雙步長顯式法是不穩定的，而修正雙步長顯式法與中心差分法具有相同的穩定條件；大步長情況下修正雙步長顯式法和中心差分法的幅值誤差小于雙步長顯式法；小步長情況下修正雙步長顯式法的阻尼比誤差和周期誤差都比中心差分法小，具有良好的綜合性能。

關鍵詞：雙步長;顯式算法;穩定性;精度

中圖分類號:O327文獻標志碼：A

基金項目：國家自然科學

收稿日期：2013-09-10修改稿收到日期：2014-01-09

基金項目：國家自然科學基金資助項目(51078103,51278152)

收稿日期：2013-08-01修改稿收到日期：2014-01-09

Explicitintegrationalgorithmbasedondoubletimesteps

YANG Chao, XIAO Shou-ne, LU Lian-tao(TractionPowerStateKeyLaboratory,SouthwestJiaotongUniversity,Chengdu610031,China)

Abstract:In order to establish a dynamical explicit method with higher accuracy and better stability, an explicit algorithm of double time steps and a corrected explicit algorithm of double time steps were proposed through integration under the assumption of linear acceleration in two time steps. The stability, accuracy and truncation error analyses were performed for the two new algorithms. A numerical example was utilized for the comparison of the two new algorithms, the central difference method and the analytical solution. The results show that the explicit algorithm of double time steps is unstable under undamped condition, however, the corrected explicit algorithm of double time steps has the same stability condition as the central difference method. The amplitude error of the explicit algorithm of double time steps is bigger than the corrected explicit algorithm of double time steps and the central difference method when the time step is long. The corrected explicit algorithm of double time steps possesses excellent comprehensive properties with less damping ratio error and periodic error than the central difference method in the case of short time step.

Keywords:doubletimesteps;explicitalgorithm;stability;accuracy

逐步積分法是動力學問題數值解法中最有效的方法之一[1]。該方法可以分為隱式法和顯式法，隱式法有線性加速度法、Wilson-θ法[2]和Newmark-β法[3]，顯式法中最具代表性的是中心差分法[4-5]。對于大規模的動力學計算問題，由于無需求解耦聯方程組，顯式法在計算速度方面具有隱式法無可比擬的優越性[6]。

隨著結構動力分析向復雜、大型和非線性的方向發展，許多學者尋求建立精度更高和穩定性更好的顯式方法。王進廷等[6]介紹了四種與中心差分法相關的顯式差分方法，并分析了其穩定性。周正華等[7]提出了利用Lagrange插值公式推導出的顯式方法，該方法具有二階精度但穩定性略差。李常青等[8]還提出只需要位移項的三步二階精度顯式方法。

上述方法都是對顯式方法的有力推進，但是這些方法忽略了顯式方法一般是以加速度為基本變量，本文在雙步長線性加速度假設的基礎上通過積分推導出雙步長顯式法，并進一步提出修正雙步長顯式法。

1方法

根據D’Alembert原理，離散結構的動力學方程為：

(1)

對于非線性系統，實際的加速度時程曲線是非線性變化的，較長時間段內的加速度線性變化只是特例。如果將曲線的時間軸劃分為若干個足夠短的時間段，每個時間段就可以按照線性問題的方法求解，這樣就把非線性動力學問題轉化為若干個線性動力學問題來處理。

1.1雙步長顯式法

(2)

設加速度方程為

(3)

對上式連續積分兩次，得到速度、速度增量和位移增量公式為

(4)

(5)

(6)

把初始條件代入式(3)和式(4)，得到

(7)

(8)

將式(2)、式(7)、式(8)和tn+1代入式(5)和式(6)，得到

(9)

為了消除線性加速度假設造成的累積誤差，可用動力學方程求解相應時刻的加速度，將式(1)變換為顯式格式，得到

(10)

1.2雙步長顯式法的修正

雙步長顯式法是通過積分推導獲得的，而本文最關鍵的步驟就是將其拓展演化為精度更高、穩定性更好的一般格式。首先由式(1)得到

(11)

(12)

式中，η為時間步長控制參數，取值范圍為[0,1]。

由于對式(9)作了一系列修正，線性加速度假設便不再成立。對于自由度規模龐大的動力學問題，在不需要求逆矩陣的情況下，利用上述解耦線性方程，就可以由式(11)、(12)構成的循環快速地計算出不同時刻的位移、速度和加速度。下文中將該方法稱為修正雙步長顯式法。

2穩定性、精度和截斷誤差分析

目前，多自由度系統數值算法的穩定性和精度，可以通過研究算法對單自由度系統的響應性質來實現。單自由度振動系統的動力學方程為：

(13)

式中:ζ、ω和fn分別為振動系統的阻尼比、自振圓頻率和激勵載荷。

2.1穩定性分析

為了導出穩定性條件，聯立方程組，得到雙步長顯式法和修正雙步長顯式法的遞推公式[10]分別為

(14)

(15)

式中，A1、A2分別為兩種方法的積分逼近算子矩陣，令Ω=Δtω，其公式如下

A1=

值得一提的是，當式(17)取ζ=0，通過求其特征方程的解析解發現，不論η取何值，其積分步長穩定區間都與中心差分法一樣，穩定條件如式(18)所示，這也是修正雙步長顯式法的特別之處。

0<Ω<2

(18)

由于式(16)和式(17)的特征方程都是高次方程，利用數值方法編程分析其穩定性是比較直觀的方法。圖1是對無阻尼和四種小阻尼情況，雙步長顯式法的積分逼近算子矩陣的譜半徑隨Ω的變化情況，在無阻尼情況下，矩陣的譜半徑都大于1，算法是不穩定的；只有在有阻尼情況下，該積分格式才穩定，所以該方法只能應用于有阻尼振動系統。修正雙步長顯式法在無阻尼和小阻尼情況下，其積分逼近算子矩陣的譜半徑隨Ω的變化情況如圖2所示，有阻尼情況下，穩定時間步長的區間隨阻尼比的增大而減??；無阻尼情況下，譜半徑在Ω=(0,2)范圍內恒等于1，表示該算法的振幅既不衰減也不發散，處于臨界穩定狀態。

圖1　雙步長顯式法的穩定性 Fig.1 Stability of explicit algorithm of double time steps

圖2　修正雙步長顯式法的穩定性(η=0.25) Fig.2 Stability of corrected explicit algorithm of double time steps (n=0.25)

2.2精度和截斷誤差分析

要使振動方程的位移解表示一個振動響應，積分逼近算子矩陣必須具有一對共軛復根，即

(19)

(20)

(21)

考慮無阻尼和阻尼比為0.2的兩種情況，將雙步長顯式法、修正雙步長顯式法(η=0.25)和中心差分法的阻尼比之差及周期延長率進行比較。由圖3可見，無阻尼情況下，中心差分法和修正雙步長顯式法的算法阻尼都為0，說明幅值誤差為0，雙步長顯式法的算法阻尼為負值，考慮到位移解中的-ξ，幅值振蕩發散，這與穩定性結果是一致的；當振動系統阻尼比為0.2時，在Ω=(0,0.4)范圍內，修正雙步長顯式法的阻尼比誤差比中心差分法小，隨著Ω進一步增大時，修正雙步長顯式法的阻尼比誤差將大于中心差分法的誤差。

圖3　阻尼比之差 Fig.3 The difference of damping ratio

圖4中，無阻尼情況下，中心差分法和修正雙步長顯式法的周期延長率相同；振動系統阻尼比為0.2時，在Ω=(0,1.2)范圍內，修正雙步長顯式法的周期誤差比中心差分法的周期誤差小，并且隨著Ω進一步增大到最大穩定極限時，中心差分法的周期誤差急劇增大。在穩定步長范圍內，修正雙步長顯式法的周期誤差比雙步長顯式法的小。

圖4　三種方法的周期延長率 Fig.4 The period elongation of three methods

雙步長顯式法和修正雙步長顯式法必然伴隨著一定的截斷誤差。根據具有拉格朗日余項的泰勒公式，分別對這兩種方法進行截斷誤差分析。

(22)

(23)

式中,Rn1和Rn2分別為雙步長顯式法和修正雙步長顯式法的位移截斷誤差；ξ1和ξ2都在tn與tn+1之間。

由式(22)可知，雙步長顯式法的位移解具有3階精度；當η=0.5時修正雙步長顯式法的位移解具有2階精度。中心差分法的位移解具有2階精度，故修正雙步長顯式法最高階精度與中心差分法的相同。

雖然三種方法中，雙步長顯式法的精度階數最高，但是在無阻尼情況下是不穩定的，并且在有阻尼情況下，其幅值誤差大于中心差分法和修正雙步長顯式法，其周期誤差大于修正雙步長顯式法。

3算例

單自由度的彈簧質量阻尼系統，質量m為10kg，彈簧剛度k為100N/mm，粘性阻尼c為0.4N·s/mm，外力為0，初始條件：初始位移為0，初始速度v0為10mm/s。

圖5　振動系統位移響應 Fig.5 Displacement response of vibration system

為了能夠明顯區別不同方法的優缺點，取較大的時間步長10ms，即Ω=1，分別用解析解、雙步長顯式法、修正雙步長顯式法和中心差分法計算位移響應，如圖5所示。由波峰波谷值可以看出，雙步長顯式法的幅值誤差最大，考慮到數值截斷，修正雙步長顯式法的算法阻尼比大于中心差分法，但修正雙步長顯式法的周期誤差要小于中心差分法，這與2.2節的結論相同。

4結論

本文在雙步長線性加速度假設的基礎上，通過積分推導得到雙步長顯式法，并進一步提出修正雙步長顯式法，通過穩定性分析、精度分析、截斷誤差分析和數值算例，得到以下結論：

(1)在無阻尼情況下，雙步長顯式法不穩定，而修正雙步長顯式法與中心差分法具有相同的穩定條件。

(2)雙步長顯式法具有3階精度，修正雙步長顯式法最高具有2階精度，與中心差分法的相同。時間步長較大時，中心差分法和修正雙步長顯式法的幅值誤差小于雙步長顯式法。

(3)η=0.25時，在(0,0.4/ω)時間步長范圍內，修正雙步長顯式法的阻尼比誤差和周期誤差都比中心差分法小，該方法具有良好的綜合性能。

參考文獻

[1]張相庭. 結構振動力學[M]. 上海:同濟大學出版社,2005: 184-197.

[2]WilsonEL.Acomputerprogramforthedynamicstressanalysisofundergroundstructures[R].SESM,ReportNo.68-1,Divisionofstructuralengineeringandstructuralmechanics,UniversityofCalifornia,Berkely, 1968.

[3]NewmarkNM.Amethodforcomputationofstructuraldynamics[J].ASCE, 1959, 85(3):67-94.

[4]尚曉江.ANSYS/LS-DYNA動力分析方法與工程實例[M]. 北京:中國水利水電出版社,2008: 129-135.

[5]莊茁. 基于ABAQUS的有限元分析和應用[M]. 北京:清華大學出版社,2009: 189-195.

[6]王進廷,杜修力. 有阻尼體系動力分析的一種顯式差分法[J]. 工程力學, 2002(3): 109-112.

WANGJin-ting,DUXiu-li.Anexplicitdifferencemethodfordynamicanalysisofastructuresystemwithdamping[J].EngineeringMechanics, 2002(3): 109-112.

[7]周正華,周扣華. 有阻尼振動方程常用顯式積分格式穩定性分析[J]. 地震工程與工程振動, 2001(3): 22-28.

ZHOUZheng-hua,ZHOUKou-hua.Stabilityanalysisofexplicitintegralmethodsfordampedvibrationequation[J].EarthquakeEngineeringandEngineeringVibration, 2001(3): 22-28.

[8]李常青,樓夢麟. 中心-偏心差分法[J]. 同濟大學學報, 2011(2): 179-186.

LIChang-qing,LOUMeng-lin.Central-eccentricdifferencemethod[J].JournalofTongjiUniversity, 2011(2): 179-186.

[9]朱昌允,秦國良,徐忠. 譜元方法求解波動方程時顯式與隱式差分方法的比較[J]. 西安交通大學學報, 2008(9): 1142-1145.

ZHUChang-yun,QINGuo-liang,XUZhong.ComparisonofexplicitcentraldifferencemethodandimplicitnewmarkmethodusingChebyshevspectralelementmethodforwaveequations[J].JournalofXi’anJiaotongUniversity, 2008(9): 1142-1145.

[10]李常青,樓夢麟,蔣麗忠. 結構動力方程求解中隱式格式向顯式格式的轉換[J]. 振動與沖擊, 2012(13): 91-94.

LIChang-qing,LOUMeng-lin,JIANGLi-zhong.Transformationofimplicitmethodtoexplicitmethodforsolvingstructuraldynamicequation[J].JournalofVibrationandShock, 2012(13): 91-94.

[11]HilbertHM.Collocation,dissipationand‘overshoot’fortimeintegrationschemesinstructuraldynamics[J].EESD, 1978, 6(1): 116-124.

[12]ZhaiWM.Twosimplefastintegrationmethodsforlarge-scaledynamicproblemsinengineering[J].InternationalJournalforNumericalinEngineering, 1996(39): 4199-4214.

[13]翟婉明. 大型結構動力分析的Newmark顯式算法[J]. 重慶交通學院學報, 1991(2): 33-41.

ZHAIWan-ming.Theexplicitschemeofnewmark’sintegrationmethodforlargestructuraldynamicanalysis[J].JournalofChongqingJiaotongInstitute, 1991(2): 33-41.

第一作者張巖女，碩士生，1989年5月生

通信作者王志偉男，博士，教授，博士生導師，1963年生

郵箱：wangzw@jnu.edu.cn

第一作者尹弘峰男，碩士，助理工程師，1986年6月生

通信作者支旭東男，博士，教授，博士生導師，1977年生