一類非耗散的顯式時間積分方法

楊 超， 肖守訥， 陽光武， 朱 濤， 楊 冰

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

一類非耗散的顯式時間積分方法

楊 超， 肖守訥， 陽光武， 朱 濤， 楊 冰

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

針對大型復雜結構的動力學問題，為了得到高效的計算方法，以泰勒公式為基礎，通過調整參數改善算法的穩定性，得到了一類具有非耗散特性且以加速度為基本變量的顯式數值積分方法。對提出的方法進行了穩定性和精度分析，并通過多個算例對提出的方法、蛙跳式中心差分法、Newmark-β法和精確解進行比較。結果表明，無阻尼情況下，所提出的方法的穩定性條件與中心差分法相同；提出的方法的振幅衰減率為0，具有非耗散特性，其周期誤差約為隱式Newmark-β法的一半；所提出的方法給出了蛙跳式中心差分法和無耗散特性的翟方法的統一格式，并可以衍生出更多的顯式時間積分方法。

顯式方法； 非耗散； 差分格式； 精度； 穩定性

引 言

結構動力分析的數值解法可以分為模態疊加法和直接積分法，其中直接積分法按照求解方式還可以分為隱式法和顯式法。隱式法有Newmark-β法[1]和Wilson-θ法[2]等，通常應用于低頻響應占主導地位的長時間結構響應問題或準線性問題。顯式法中最具代表性的是中心差分法[3]，通常應用于高頻響應占主導地位的非線性問題或瞬態問題。對于大型非線性的動力學計算問題，顯式方法的應用更加廣泛。隨著結構動力響應分析向復雜非線性的方向發展，計算量和計算代價顯著提高，計算方法的速度和精度越來越受到重視。

文獻[4]提出了一種具有三階精度的顯式差分方法，適用于有阻尼系統的動力響應分析，高階精度需要增加中間變量而使整體計算量增加[5]。隱式算法也可以轉換為顯式算法，一般以位移為第一求解變量[6]。文獻[7]提出了Newmark更新精細積分法，具有良好的穩定性，但要求阻尼矩陣為對角矩陣，所以只能屬于半顯式的方法。很多新開發的算法一般都是基于泰勒展開或加權殘余法[8]，應用泰勒公式展開并截除高階小量會引入一定的誤差，而加權殘余法也同樣存在一定的誤差，所以在考慮計算速度的同時要盡量提高算法精度。

現有的研究大多以位移為基本變量，速度和加速度通過對位移的求導獲得，而以加速度為基本變量的方法很少受到重視。本文以加速度為基本變量，速度和位移通過對加速度的近似積分獲得，在泰勒展開的基礎上，通過調整積分參數改善算法的穩定性，提出了一類無耗散特性的顯式時間積分方法。文中提出的方法的穩定性條件與中心差分法相同，其周期誤差僅為Newmark-β法的一半。文中提出的顯式時間積分方法統一了蛙跳式中心差分法和翟方法(φ=φ=0.5)，并可以衍生出更多的無耗散特性的顯式積分格式。

1 顯式方法

顯式是指增量步結束時的狀態僅依賴于該增量步開始之前的位移、速度和加速度，計算時不需要進行迭代或矩陣求逆。結構的動力學運動方程為

(1)

1.1 中心差分法

文獻[9-10]中通常所描述的中心差分法(CDM)一般都是半顯式格式的，在阻尼矩陣為非對角的情況下會退化為隱式格式。中心差分法對運動方程直接進行積分，其主要公式如下

(2)

將式(2)帶入式(1)，得到以位移為第一求解變量的運動方程

(3)

中心差分法先利用式(3)求得tn+1時刻的位移，然后根據式(2)進行兩次求導得到tn時刻的速度和加速度。這種半顯式的中心差分法是以位移為基本變量，不存在算法阻尼，具有二階精度。

1.2 蛙跳式中心差分法

蛙跳式中心差分法是中心差分法的一種變異形式，是完全顯式格式的[9]。該方法的獨特之處在于：取單個時間步中間時刻的點作為普通中心差分法的第三點，每個循環內只涉及單個時間步長內的物理量，中間時刻的點只計算速度，整數時刻的點只計算位移和加速度，這種跳躍步進方式猶如蛙跳一般。在增量步開始時，首先計算當前時刻的加速度為：

(4)

(5)

(6)

1.3 翟方法

翟方法[12]是由翟婉明院士提出的一類顯式二步數值積分方法，該方法是基于著名的Newmark隱式積分方法構造的一類顯式積分方法，所以兩者的格式較為接近，其積分格式為

(7)

式中φ,φ為積分參數，一般取為0.5；在得到tn+1時刻的位移和速度后，就可以由式(8)得到該時刻的加速度。

(8)

翟方法具有二階精度，在積分參數都為0.5時不存在算法阻尼，與Newmark-β法具有類似的算法阻尼特性。該顯式方法也是以加速度為基本變量，只要其質量陣是對角陣或經對角化而不論阻尼陣如何復雜，積分過程中都無須聯立求解耦合方程組，避免了組裝等效剛度矩陣和矩陣求逆，這就大大節省了計算內存和時間，因而用于分析工程中大型非線性結構動力問題時，可以顯著地提高計算速度和效率。

2 基于加速度的顯式積分方法

蛙跳式中心差分法和翟方法具有共同的特征：以加速度為基本變量，無數值阻尼，即無耗散特性。根據這些特點，本文以泰勒展開公式為基礎，通過研究可變參數對算法穩定性的影響，得到了一類非耗散的顯式數值積分方法。首先用泰勒公式表示位移和速度

(9)

假設加速度的導數滿足向后差分公式，則有

(10)

把式(10)帶入式(9)，并略去高階項，得到

(11)

將式(11)中值不為1的系數以變量代替

(12)

式(12)中，推薦取α=1,β=γ,γ∈[0,1]，γ的取值盡量接近0.5，此時該算法僅有唯一的可變參量γ，具有非耗散特性，即不存在算法阻尼，這個結果是通過下文的穩定性分析得到的。當3個系數變量取上述值時，結合式(8)，即為本文要提出的顯式時間積分方法，其計算過程歸納如下：

1.計算起步

(1)形成質量矩陣M、剛度矩陣K和阻尼矩陣C；

2.循環開始

(2)更新矩陣M，K和C；

3 穩定性與精度分析

多自由度系統中各類積分算法的穩定性和精度一般可以通過研究算法對單自由度系統的響應性質來實現。單自由度振動系統的動力學方程為

(13)

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

3.1 穩定性分析

為了導出穩定性條件，由式(12)和式(13)，得到加速度顯式法的遞推公式為

(14)

式中A為積分逼近算子矩陣，令采樣頻率Ω=Δtω，得到

(15)

式中P=-2(1+γ)ζΩ-αΩ2,Q=2γζΩ+γΩ2。

由穩定性原理可知，λi為矩陣A的第i個特征值，當矩陣A的譜半徑ρ(A)=max|λi|<1時，特征方程有共軛復根，積分格式是收斂的；當譜半徑ρ(A)=1時，積分格式是臨界穩定的[13]。通過求解特征方程的解析解發現，當阻尼比ζ=0時，加速度顯式法的時間步長穩定區間與中心差分法是一樣的，在(0,2/ω)之間。

如圖1所示，在無阻尼情況下，α<1時，譜半徑大于1，積分格式是發散的，不能符合要求；α>1時，譜半徑存在小于1的區間，積分格式是收斂的，存在算法阻尼，即振幅會隨著時間增長而不斷衰減；當α=1且β≠γ時，積分格式存在起始點不為0的穩定區間，這意味著步長取某值時算法是穩定的，但進一步減小時間步長時算法可能不穩定；只有在α=1并且β=γ時，通過分析高次方程的解，無論γ取何值，譜半徑在Ω∈(0,2)都是恒等于1的，振幅既不衰減也不發散而保持恒定，此時算法具有非耗散特性。圖2顯示了阻尼對加速度顯式法的穩定性影響，隨著阻尼比的增大，算法穩定區間逐漸減小。譜半徑越小，表示算法收斂越快。

圖1 無阻尼時加速度顯式法的穩定性Fig.1 The stability of the acceleration explicit methods under undamped condition

圖2 阻尼對穩定性的影響(γ=0.2)Fig.2 The influence of damping to the stability (γ=0.2)

3.2 精度分析

積分逼近算子矩陣具有一對共軛復根，可以用這對共軛復根表示出振動方程的解

(16)

(17)

振幅衰減率AD反映了幅值誤差，周期延長率TD反映了周期誤差，見式(18)和(19)，因此算法的精確度可以由振幅衰減率和周期延長率來度量。

(18)

(19)

在無阻尼情況下，將加速度顯式法、Newmark-β法和中心差分法(CDM)的振幅衰減率和周期延長率進行比較。由圖3可見，無耗散特性的加速度顯式法、Newmark-β法和中心差分法的振幅衰減率都為0，說明幅值誤差為0，而α=1.2時振幅會迅速衰減，這與穩定性結果是一致的。如圖4所示，無阻尼情況下，隱式的Newmark-β法的周期是隨著時間步長增加而逐漸變長的；無耗散特性的加速度顯式法和中心差分法的周期延長率相同，周期都是縮小的，且周期誤差約為Newmark-β法的一半。

圖3 幅值誤差Fig.3 The amplitude error

圖4 周期誤差Fig.4 The periodic error

4 算 例

采用文獻[14]的算例，兩自由度振動系統單端受恒定力激勵，振動方程為

(20)

(21)

分別用文中提到的蛙跳式中心差分法(LCDM)、翟方法、Newmark-β法和本文方法(γ=0;γ=0.2;γ=0.5;γ=0.8)分別計算該振動系統，除了本文方法取了兩種時間步長0.28和0.14 s，即最小周期的1/10和1/20，其他方法的時間步長都取0.14 s，質量點的加速度結果見圖5和6，位移結果列于表1和2中。

圖5 點1的加速度結果比較(Δt=0.14 s)Fig.5 The comparison of accelerations of point 1 (Δt=0.14 s)

圖6 點2的加速度結果比較(Δt=0.14 s)Fig.6 The comparison of accelerations of point 2 (Δt=0.14 s)

表1 點1的位移結果比較

表2 點2的位移結果比較

本文方法(γ=0)和蛙跳式中心差分法的位移和加速度結果完全一致。γ<0.5時本文方法的結果相對于精確解有一定的相位延遲，誤差較大。當γ=0.5時誤差最小，Newmark-β法次之。所有方法中，本文方法(γ=0.5)的位移計算結果與精確解擬合度最高，其均方根誤差也最小，其精度超過了隱式Newmark-β法。計算精度與計算速度通常是相互制約的兩個指標，雖然蛙跳式中心差分法的誤差偏大，但正因為其速度方面的優勢才會被顯式動力學軟件采用，而通過循環公式可以看出本文方法與該方法在計算量上幾乎相同，且精度可以通過積分參數調整。

從位移結果可以看出，本文方法的位移均方根誤差隨著時間步長增大而增大，本文方法對積分參數γ的選擇是比較敏感的，隨著γ從0增加到0.5，均方根誤差逐漸減小，當γ從0.5逐漸增大時，均方根誤差又開始增大，所以γ的推薦選取條件為：盡量接近0.5。

5 工程應用

5.1 鐵路貨車縱向振動

重載貨物列車在貨場通常要進行車輛之間的連掛，已連掛的車輛組受到被連掛車輛的沖擊而產生振動。將車輛和鉤緩裝置組成的系統簡化為由n個剛體串聯組成的彈簧-質量系統[15]，車輛與鉤緩裝置的數目相等。為了分析方便，假設沖擊端的車輛由簡諧力激勵，不受沖擊最末端的鉤緩裝置一端被固定，鉤緩裝置考慮為漸進軟化的非線性彈簧。

簡化模型的參數為：車輛數為100，簡諧力fn=106sin(18t) N，質量mi=50 000 kg，非線性彈簧剛度為ki=109[1-(ui-ui-1)2] N/m，i=1,2,…,n。該系統的最小和最大圓頻率分別為2.21和282.81 rad/s，則本文方法的最大時間步長為0.007 s。為了便于比較，采用Newmark-β法步長Δt=0.000 5 s的結果作為“精確解”。根據上節中的算例確定的積分參數選取原則和時間步長Δt=0.007 s，約為最小周期的1/π，外載激勵周期的1/50，最大周期的1/406，采用本文方法(γ=0.45)和Newmark-β法分別計算該模型。

采用硬件平臺為戴爾Inspiron580微機，軟件平臺為64位MATLAB 2010a。受外力作用車輛的位移結果如圖7所示。本文方法和Newmark-β法的計算結果與精確解在初始的半個最大周期基本重合，經過半個周期振動以后，振動趨勢一致，幅值和相位與精確解相比存在誤差，但本文方法的結果更接近精確解。當時間步長為0.007 s時，上述方法分別耗時0.484和1.967 s；當時間步長為0.000 5 s時，本文方法與Newmark-β法分別耗時6.582和27.04 s，相同條件下本文方法的速度約為Newmark-β法的4倍。此外，計算速度與程序語言的選取也是有關系的，若Newmark-β法取較大步長可以抵消計算速度上的弱勢，但計算精度會降低。

圖7 受力車輛的位移結果Fig.7 Displacement of the forced vehicle

5.2 鐵路客車垂向振動

車輛振動試驗臺通過激勵輪對使整節車廂振動，可以用于測試鐵路客車的舒適度，而舒適度與車體加速度直接相關。為了模擬車輛振動試驗臺測試鐵路客車垂向振動性能，根據達朗貝爾原理建立車輛的垂向動力學模型，模型簡圖和公式參見文獻[12]和[16]。該模型為7剛體10自由度模型，包括1個車體、2個構架和4個輪對，車體和構架考慮沉浮自由度和點頭自由度，輪對只考慮沉浮自由度。

車輛模型參數為：車體質量和點頭慣量Mc=39 500 kg，Jc=2 312 000 kg·m2，構架質量和點頭慣量為Mt=2 200 kg，Jt=2 200 kg·m2，輪對質量Mw=1 900 kg,一系懸掛剛度和阻尼為k1=2 130 000 N/m，c1=120 000 N·s/m，二系懸掛剛度和阻尼為k2=800 000 N/m，c2=217 400 N·s/m，車輛定距和轉向架軸距為lc=18 m，lt=2.4 m。車輛的4個輪對采用正弦同相激勵，假設激勵力為ft=1 000sin(50t) N。該系統最小和最大圓頻率分別為5.4和61.2 rad/s，由于阻尼的影響需要取時間步長為Δt=0.004 s，約為最小周期的1/25，外載激勵周期的1/32，最大周期的1/290，分別采用文中提出的加速度顯式法(AEM)和Newmark-β法計算該模型。以小步長的線性加速度法(LAM)的結果作為參考標準值，該小步長取Δt=0.000 4 s。

由于一系和二系阻尼的作用，車體響應逐漸趨于穩定，穩定后振動頻率與激勵頻率一致。由于步長與激勵周期的比值較小，所以計算結果之間差異很小，為了便于區別，圖8中還給出車體穩定以后半個周期的車體垂向加速度時間歷程。車體加速度在0.18 s以后趨于穩定，加速度顯式法和Newmark-β法的結果的變化趨勢與參考結果基本一致，但幅值上有差異，每個時間點對應的標準值位于加速度顯式法和Newmark-β法的結果之間，加速度顯式法的結果偏大，而Newmark-β法的結果偏小。

圖8 車體加速度Fig.8 Acceleration of the carbody

6 結 論

本文提出了一類以加速度為基本變量的非耗散的顯式時間積分方法，對加速度顯式法進行了穩定性分析和精度分析，并通過一個線性算例和兩個工程應用實例將加速度顯式法和常用的經典方法進行了詳細的比較，得到以下結論：

1) 加速度顯式法是條件穩定的，與中心差分法具有相同的穩定區間；有阻尼情況下，加速度顯式法的穩定區間隨著阻尼的增大而減小。

2) 加速度顯式法的振幅衰減率為0，具有非耗散特性；其周期延長率和中心差分法相同，周期誤差約為隱式Newmark-β法的一半。算例表明加速度顯式法在相同步長情況下具有較高的計算速度和精度。

3) 加速度顯式法給出了蛙跳式中心差分法和無耗散特性的翟方法的統一格式，并可以衍生出更多的顯式積分格式。

[1] Newmark N M. A method of computation for structural dynamics[J]. Journal of the Engineering Mechanics Division (ASCE), 1959, 85(3):67—94.

[2] Wilson E L, Farhoomand I, Bathe K J. Nonlinear dynamic analysis of complex structures[J]. Earthquake Engineering and Structural Dynamics, 1973，1(3): 241—252.

[3] 尚曉江, 蘇建宇, 王化鋒. ANSYS/LS-DYNA動力分析方法與工程實例[M]. 二版. 北京:中國水利水電出版社, 2008: 129—135.

[4] 張石磊,王煥定,張曉志. 結構地震響應分析新方法[J]. 國際地震動態, 2009，(6): 8—14.

Zhang Shilei, Wang Huanding, Zhang Xiaozhi. A new method in analyzing seismic response of the structure[J]. Recent Developments in World Seismology, 2009，(6): 8—14.

[5] Idesman A V, Schmidt M, Sierakowski R L. A new explicit predictor-multicorrector high-order accurate method for linear elastodynamics[J]. Journal of Sound and Vibration, 2008, 310(1/2): 217—229.

[6] Bonelli A, Bursi O S. Explicit predictor-multicorrector time discontinuous Galerkin methods for linear dynamics[J]. Journal of Sound and Vibration, 2001, 246(4):625—652.

[7] 郭澤英,李青寧. 動力反應分析的顯式積分方法及其穩定性[J]. 地震工程與工程振動, 2008，28(2): 15—19.

Guo Zeying, Li Qingning. An explicit integral method and the stability of dynamic response analysis[J]. Journal of Earthquake Engineering and Engineering Vibration, 2008，28(2): 15—19.

[8] Stolle D. A direct integration algorithm and the consequence of numerical stability[J]. Journal of Sound and Vibration, 1995, 180(3): 513—518.

[9] 張雄, 王天舒. 計算動力學[M]. 北京:清華大學出版社, 2007: 140—207.

Zhang Xiong, Wang Tianshu. Computational Dynamics[M]. Beijing: Tsinghua University Press, 2007: 140—207.

[10]Wiberg N E, Li X D. Adaptive finite element procedures for linear and non-linear dynamics[J]. International Journal for Numerical Methods in Engineering, 1999, 46(10): 1 781—1 802.

[11]Hallquist J O. LS-DYNA Theory Manual[M]. Livermore Software Technology Corporation, 2006.

[12]翟婉明. 車輛-軌道耦合動力學[M]. 第2版. 北京:中國鐵道出版社, 2002: 107—127.

Zhai Wanming. Vehicle-Track Coupling Dynamics[M]. 2nd ed. Beijing: China Railway Publishing House, 2002: 107—127.

[13]Hilbert H M. Collocation, dissipation and ‘Overshoot’ for time integration schemes in structural dynamics[J]. EESD, 1978, 6(1): 116—124.

[14]李常青,樓夢麟,余志武,等. 近似平衡多項式加速度動力顯式算法[J]. 應用力學學報, 2011，28(5): 475—479.

Li Changqing, Lou Menglin, Yu Zhiwu, et al. Psedo-balance polynomial acceleration explicit method[J]. Chinese Journal of Applied Mechanics, 2011，28(5): 475—479.

[15]Chang S Y. Improved explicit method for structural dynamics[J]. Journal of Engineering Mechanics (ASCE), 2007,133(7):748—760.

[16]張志超,趙巖,林家浩. 車橋耦合系統非平穩隨機振動分析[J]. 振動工程學報, 2007, 20(5): 339—446.

Zhang Zhichao, Zhao Yan, Li Jiahao. Non-stationary random vibration analysis for vehicle-bridge coupled systems[J]. Journal of Vibration Engineering, 2007, 20(5): 339—446.

Non-dissipative explicit time integration methods of the same class

YANGChao,XIAOShou-ne,YANGGuang-wu,ZHUTao,YANGBing

(Traction Power State Key Laboratory, Southwest Jiaotong University, Chengdu 610031, China)

In order to obtain efficient calculation methods for dynamical problems of large and complex structures. A class of non-dissipitive explicit numerical integration methods is proposed through adjusting parameters to improve the stability of algorithms. Accelerations are used as the basic variables in the proposed methods which are on the basis of the Taylor's formula. The analyses of stability and accuracy are performed for the proposed methods. Finally, some numerical examples are given for the comparison between the proposed methods,and some current methods including the leapfrog central difference method (LCDM), the Newmark-βmethod and the exact solution. The results show that the proposed methods have the same stability condition as the LCDM under undamped condition. The amplitude attenuation of the proposed methods is 0, which means they are non-dissipative. The periodic error is about half of that of the implicit Newmark-βmethod. The proposed methods provide a unified format of the LCDM and the Zhai method without dissipation, and more explicit integration algorithms can be derived directly.

explicit method; non-dissipative; difference scheme; accuracy; stability

2014-01-05；

2014-07-16

國家自然科學基金資助項目(51275432，51405402)；鐵道部科技研究開發計劃資助項目(2011J022-A)；四川省科技廳應用基礎研究項目(2014JY0242)

O241.8； O313

A

1004-4523(2015)03-0441-08

10.16385/j.cnki.issn.1004-4523.2015.03.014

楊超(1988—)，男，博士研究生。電話：15882317785 ； E-mail： yangchaosky@foxmail.com