基于雙步長(zhǎng)的顯式積分算法

基于雙步長(zhǎng)的顯式積分算法

楊超，肖守訥，魯連濤

(西南交通大學(xué)牽引動(dòng)力國(guó)家重點(diǎn)實(shí)驗(yàn)室， 成都610031)

摘要：為了建立精度更高、穩(wěn)定性更好的顯式積分算法，在雙步長(zhǎng)線性加速度假設(shè)的條件下，通過(guò)積分推導(dǎo)，提出了雙步長(zhǎng)顯式法和修正雙步長(zhǎng)顯式法。對(duì)這兩種新方法進(jìn)行了穩(wěn)定性、精度和截?cái)嗾`差分析，并通過(guò)算例對(duì)兩種新方法、中心差分法和解析解進(jìn)行比較。結(jié)果表明，無(wú)阻尼情況下，雙步長(zhǎng)顯式法是不穩(wěn)定的，而修正雙步長(zhǎng)顯式法與中心差分法具有相同的穩(wěn)定條件；大步長(zhǎng)情況下修正雙步長(zhǎng)顯式法和中心差分法的幅值誤差小于雙步長(zhǎng)顯式法；小步長(zhǎng)情況下修正雙步長(zhǎng)顯式法的阻尼比誤差和周期誤差都比中心差分法小，具有良好的綜合性能。

關(guān)鍵詞：雙步長(zhǎng);顯式算法;穩(wěn)定性;精度

中圖分類號(hào):O327文獻(xiàn)標(biāo)志碼：A

基金項(xiàng)目：國(guó)家自然科學(xué)

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

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目(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

逐步積分法是動(dòng)力學(xué)問(wèn)題數(shù)值解法中最有效的方法之一[1]。該方法可以分為隱式法和顯式法，隱式法有線性加速度法、Wilson-θ法[2]和Newmark-β法[3]，顯式法中最具代表性的是中心差分法[4-5]。對(duì)于大規(guī)模的動(dòng)力學(xué)計(jì)算問(wèn)題，由于無(wú)需求解耦聯(lián)方程組，顯式法在計(jì)算速度方面具有隱式法無(wú)可比擬的優(yōu)越性[6]。

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

上述方法都是對(duì)顯式方法的有力推進(jìn)，但是這些方法忽略了顯式方法一般是以加速度為基本變量，本文在雙步長(zhǎng)線性加速度假設(shè)的基礎(chǔ)上通過(guò)積分推導(dǎo)出雙步長(zhǎng)顯式法，并進(jìn)一步提出修正雙步長(zhǎng)顯式法。

1方法

根據(jù)D’Alembert原理，離散結(jié)構(gòu)的動(dòng)力學(xué)方程為：

(1)

對(duì)于非線性系統(tǒng)，實(shí)際的加速度時(shí)程曲線是非線性變化的，較長(zhǎng)時(shí)間段內(nèi)的加速度線性變化只是特例。如果將曲線的時(shí)間軸劃分為若干個(gè)足夠短的時(shí)間段，每個(gè)時(shí)間段就可以按照線性問(wèn)題的方法求解，這樣就把非線性動(dòng)力學(xué)問(wèn)題轉(zhuǎn)化為若干個(gè)線性動(dòng)力學(xué)問(wèn)題來(lái)處理。

1.1雙步長(zhǎng)顯式法

(2)

設(shè)加速度方程為

(3)

對(duì)上式連續(xù)積分兩次，得到速度、速度增量和位移增量公式為

(4)

(5)

(6)

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

(7)

(8)

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

(9)

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

(10)

1.2雙步長(zhǎng)顯式法的修正

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

(11)

(12)

式中，η為時(shí)間步長(zhǎng)控制參數(shù)，取值范圍為[0,1]。

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

2穩(wěn)定性、精度和截?cái)嗾`差分析

目前，多自由度系統(tǒng)數(shù)值算法的穩(wěn)定性和精度，可以通過(guò)研究算法對(duì)單自由度系統(tǒng)的響應(yīng)性質(zhì)來(lái)實(shí)現(xiàn)。單自由度振動(dòng)系統(tǒng)的動(dòng)力學(xué)方程為：

(13)

式中:ζ、ω和fn分別為振動(dòng)系統(tǒng)的阻尼比、自振圓頻率和激勵(lì)載荷。

2.1穩(wěn)定性分析

為了導(dǎo)出穩(wěn)定性條件，聯(lián)立方程組，得到雙步長(zhǎng)顯式法和修正雙步長(zhǎng)顯式法的遞推公式[10]分別為

(14)

(15)

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

A1=

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

0<Ω<2

(18)

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

圖1　雙步長(zhǎng)顯式法的穩(wěn)定性 Fig.1 Stability of explicit algorithm of double time steps

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

2.2精度和截?cái)嗾`差分析

要使振動(dòng)方程的位移解表示一個(gè)振動(dòng)響應(yīng)，積分逼近算子矩陣必須具有一對(duì)共軛復(fù)根，即

(19)

(20)

(21)

考慮無(wú)阻尼和阻尼比為0.2的兩種情況，將雙步長(zhǎng)顯式法、修正雙步長(zhǎng)顯式法(η=0.25)和中心差分法的阻尼比之差及周期延長(zhǎng)率進(jìn)行比較。由圖3可見(jiàn)，無(wú)阻尼情況下，中心差分法和修正雙步長(zhǎng)顯式法的算法阻尼都為0，說(shuō)明幅值誤差為0，雙步長(zhǎng)顯式法的算法阻尼為負(fù)值，考慮到位移解中的-ξ，幅值振蕩發(fā)散，這與穩(wěn)定性結(jié)果是一致的；當(dāng)振動(dòng)系統(tǒng)阻尼比為0.2時(shí)，在Ω=(0,0.4)范圍內(nèi)，修正雙步長(zhǎng)顯式法的阻尼比誤差比中心差分法小，隨著Ω進(jìn)一步增大時(shí)，修正雙步長(zhǎng)顯式法的阻尼比誤差將大于中心差分法的誤差。

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

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

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

雙步長(zhǎng)顯式法和修正雙步長(zhǎng)顯式法必然伴隨著一定的截?cái)嗾`差。根據(jù)具有拉格朗日余項(xiàng)的泰勒公式，分別對(duì)這兩種方法進(jìn)行截?cái)嗾`差分析。

(22)

(23)

式中,Rn1和Rn2分別為雙步長(zhǎng)顯式法和修正雙步長(zhǎng)顯式法的位移截?cái)嗾`差；ξ1和ξ2都在tn與tn+1之間。

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

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

3算例

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

圖5　振動(dòng)系統(tǒng)位移響應(yīng) Fig.5 Displacement response of vibration system

為了能夠明顯區(qū)別不同方法的優(yōu)缺點(diǎn)，取較大的時(shí)間步長(zhǎng)10ms，即Ω=1，分別用解析解、雙步長(zhǎng)顯式法、修正雙步長(zhǎng)顯式法和中心差分法計(jì)算位移響應(yīng)，如圖5所示。由波峰波谷值可以看出，雙步長(zhǎng)顯式法的幅值誤差最大，考慮到數(shù)值截?cái)?，修正雙步長(zhǎng)顯式法的算法阻尼比大于中心差分法，但修正雙步長(zhǎng)顯式法的周期誤差要小于中心差分法，這與2.2節(jié)的結(jié)論相同。

4結(jié)論

本文在雙步長(zhǎng)線性加速度假設(shè)的基礎(chǔ)上，通過(guò)積分推導(dǎo)得到雙步長(zhǎng)顯式法，并進(jìn)一步提出修正雙步長(zhǎng)顯式法，通過(guò)穩(wěn)定性分析、精度分析、截?cái)嗾`差分析和數(shù)值算例，得到以下結(jié)論：

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

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

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

參考文獻(xiàn)

[1]張相庭. 結(jié)構(gòu)振動(dòng)力學(xué)[M]. 上海:同濟(jì)大學(xué)出版社,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動(dòng)力分析方法與工程實(shí)例[M]. 北京:中國(guó)水利水電出版社,2008: 129-135.

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

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

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

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

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

[8]李常青,樓夢(mèng)麟. 中心-偏心差分法[J]. 同濟(jì)大學(xué)學(xué)報(bào), 2011(2): 179-186.

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

[9]朱昌允,秦國(guó)良,徐忠. 譜元方法求解波動(dòng)方程時(shí)顯式與隱式差分方法的比較[J]. 西安交通大學(xué)學(xué)報(bào), 2008(9): 1142-1145.

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

[10]李常青,樓夢(mèng)麟,蔣麗忠. 結(jié)構(gòu)動(dòng)力方程求解中隱式格式向顯式格式的轉(zhuǎn)換[J]. 振動(dòng)與沖擊, 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]翟婉明. 大型結(jié)構(gòu)動(dòng)力分析的Newmark顯式算法[J]. 重慶交通學(xué)院學(xué)報(bào), 1991(2): 33-41.

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

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

通信作者王志偉男，博士，教授，博士生導(dǎo)師，1963年生

郵箱：wangzw@jnu.edu.cn

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

通信作者支旭東男，博士，教授，博士生導(dǎo)師，1977年生