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

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

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

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

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

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

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

顯式方法； 非耗散； 差分格式； 精度； 穩(wěn)定性

引 言

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

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

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

1 顯式方法

顯式是指增量步結(jié)束時的狀態(tài)僅依賴于該增量步開始之前的位移、速度和加速度，計算時不需要進(jìn)行迭代或矩陣求逆。結(jié)構(gòu)的動力學(xué)運動方程為

(1)

1.1 中心差分法

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

(2)

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

(3)

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

1.2 蛙跳式中心差分法

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

(4)

(5)

(6)

1.3 翟方法

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

(7)

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

(8)

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

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

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

(9)

假設(shè)加速度的導(dǎo)數(shù)滿足向后差分公式，則有

(10)

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

(11)

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

(12)

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

1.計算起步

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

2.循環(huán)開始

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

3 穩(wěn)定性與精度分析

多自由度系統(tǒng)中各類積分算法的穩(wěn)定性和精度一般可以通過研究算法對單自由度系統(tǒng)的響應(yīng)性質(zhì)來實現(xiàn)。單自由度振動系統(tǒng)的動力學(xué)方程為

(13)

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

3.1 穩(wěn)定性分析

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

(14)

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

(15)

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

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

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

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

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

3.2 精度分析

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

(16)

(17)

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

(18)

(19)

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

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

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

4 算 例

采用文獻(xiàn)[14]的算例，兩自由度振動系統(tǒng)單端受恒定力激勵，振動方程為

(20)

(21)

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

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

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

表1 點1的位移結(jié)果比較

表2 點2的位移結(jié)果比較

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

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

5 工程應(yīng)用

5.1 鐵路貨車縱向振動

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

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

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

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

5.2 鐵路客車垂向振動

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

車輛模型參數(shù)為：車體質(zhì)量和點頭慣量Mc=39 500 kg，Jc=2 312 000 kg·m2，構(gòu)架質(zhì)量和點頭慣量為Mt=2 200 kg，Jt=2 200 kg·m2，輪對質(zhì)量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，車輛定距和轉(zhuǎn)向架軸距為lc=18 m，lt=2.4 m。車輛的4個輪對采用正弦同相激勵，假設(shè)激勵力為ft=1 000sin(50t) N。該系統(tǒng)最小和最大圓頻率分別為5.4和61.2 rad/s，由于阻尼的影響需要取時間步長為Δt=0.004 s，約為最小周期的1/25，外載激勵周期的1/32，最大周期的1/290，分別采用文中提出的加速度顯式法(AEM)和Newmark-β法計算該模型。以小步長的線性加速度法(LAM)的結(jié)果作為參考標(biāo)準(zhǔn)值，該小步長取Δt=0.000 4 s。

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

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

6 結(jié) 論

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

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

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

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

[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] 張石磊,王煥定,張曉志. 結(jié)構(gòu)地震響應(yīng)分析新方法[J]. 國際地震動態(tài), 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] 郭澤英,李青寧. 動力反應(yīng)分析的顯式積分方法及其穩(wěn)定性[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] 張雄, 王天舒. 計算動力學(xué)[M]. 北京:清華大學(xué)出版社, 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]翟婉明. 車輛-軌道耦合動力學(xué)[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]. 應(yīng)用力學(xué)學(xué)報, 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]張志超,趙巖,林家浩. 車橋耦合系統(tǒng)非平穩(wěn)隨機振動分析[J]. 振動工程學(xué)報, 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

國家自然科學(xué)基金資助項目(51275432，51405402)；鐵道部科技研究開發(fā)計劃資助項目(2011J022-A)；四川省科技廳應(yīng)用基礎(chǔ)研究項目(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