地震碰撞分析中接觸單元模型的精細積分解法

張瑞杰, 李青寧, 尹俊紅

(西安建筑科技大學 土木工程學院，西安　710055)

?

地震碰撞分析中接觸單元模型的精細積分解法

張瑞杰, 李青寧, 尹俊紅

(西安建筑科技大學 土木工程學院，西安710055)

摘要：地震碰撞數值分析是結構抗震研究的重要內容，將精細積分法引入到地震碰撞數值分析中。結合5種常用接觸單元模型的特點，將碰撞力計算式表達成統一的形式，推導了碰撞階段的精細積分公式。將地震碰撞過程劃分為碰撞階段和非碰撞階段，對不同的階段采用不同的積分步長，并且提出一種基于線性加速度假定的碰撞時刻搜索法。用兩小球自由彈性碰撞模型對精細積分法程序進行驗證。并采用本文方法對Chau振動臺碰撞試驗進行了數值仿真。結論為：數值分析方法既保證了模擬地震碰撞的精度又具有較高的計算效率；在地震碰撞分析中應合理選擇接觸單元模型和模型參數；對于近似彈性碰撞情況Hertz模型和Hertz-damp模型數值仿真結果與試驗結果吻合較好。

關鍵詞：地震;結構碰撞;數值分析;精細積分法;接觸單元

震害調查表明，地震時相鄰結構或構件由于動力特性的差異導致的碰撞會引起局部破壞、落梁，甚至倒塌[1-3]。要揭示地震中橋梁結構的碰撞機理、進行碰撞研究及抗震設計，需要精確地模擬橋梁在地震中的碰撞行為。國內外眾多學者已相繼開展了地震碰撞試驗及數值分析研究[4-6]。目前，地震碰撞數值分析大多采用接觸單元法[7-8]，該方法是在普通有限元時程積分基礎上發展起來的，通過在可能碰撞的結構體之間設置接觸單元，當碰撞發生時激活接觸單元，實現結構碰撞行為的模擬。迄今，已發展了多種接觸單元模型[9-12]。

接觸單元法本質上是將接觸碰撞過程看作碰撞力的作用過程，在結構的動力平衡方程中加入碰撞力，得到考慮碰撞作用的動力平衡方程。動力平衡方程在時域內的數值求解有多種方法[13-14]。傳統的顯式積分法如中心差分法、隱式積分方法如Newmark-β法、Wilson-θ法都是將微分化成為差分算子，使得這些方法僅具有二階精度或三階精度，而且對時間步長非常敏感[14]。在碰撞動力平衡方程求解中，碰撞力是在當前階段的位移和速度響應基礎上計算得到，然后將其作為下一時間點的已知力進行求解。這造成了碰撞力施加時間的滯后。因而在地震碰撞模擬中，只能采用較小的時間步長來減輕這種不良影響，但這大幅增加了計算量。

鐘萬勰[15]提出了一種求解結構動力響應的精細積分法。鮑四元等[16]將其應用于彈性撞擊問題的動態響應求解，并指出精細積分法具有受時間步長限制小、精度高和無條件穩定的優點。但目前鮮有文獻報道將精細積分法應用到地震碰撞過程的模擬中。本文將精細積分法引入到地震碰撞過程模擬中，結合5種常用接觸單元模型的特點，將其碰撞力表達成統一形式，推導了精細積分公式。利用精細積分法對步長不敏感的特點，本文提出了一種變步長的措施，在沒有接觸碰撞時用較大步長，提高計算效率，在接觸碰撞階段用較小步長，以提高模擬地震碰撞的精度。編制了計算程序，通過對Chau振動臺試驗進行數值仿真驗證了本文方法的有效性和優越性。

1接觸單元模型

考察幾種常用的接觸單元模型，其模型假定和碰撞力特點見圖1。線彈簧模型、線彈簧阻尼(Kelvin-Voigt)模型、非線性彈簧(Hertz)模型、非線性彈簧阻尼(Hertz-damp)模型在碰撞階段的碰撞力表達式均與相對壓入量及速度差有關，可寫成如下統一形式：

(1)

當u1(t)-u2(t)-gp≥0

圖1　接觸單元模型Fig.1 Model of contact elements

(2)

式中：m1、m2分別為碰撞體的質量；ζ為與恢復系數e有關的碰撞阻尼常數，當e=0時表示兩碰撞體之間發生完全塑性碰撞，e=1時表示兩碰撞體之間發生完全彈性碰撞。

線彈簧阻尼接觸單元假設在碰撞期間均存在能量耗散。由于阻尼系數ck為一常數，因而在接觸開始時刻會出現突加的阻尼力，而在恢復過程中，阻尼力可能超過彈性恢復力，引起碰撞力為負值的情況，與實際情況不符。

(3)

(4)

(5)

2精細積分法求解碰撞動力方程

對于圖2所示的相鄰的兩個單自由度系統碰撞模型，其動力平衡方程表達式為：

(6)

圖2　碰撞結構模型Fig.2 Model of colliding structures

方程(6)分兩種情況進行討論，其一在非接觸階段，Fc(t)=0，即為兩個獨立的動力平衡方程；其二，在接觸碰撞階段Fc(t)采用公式(1)計算。下面詳細推導接觸碰撞階段的精細積分法公式。

將公式(1)代入方程(6)，將碰撞力項中的彈性恢復力在左端進行分解并與結構剛度矩陣合并，而阻尼力仍作為非齊次項處理。可寫成如下形式：

(7)

對方程(7)，引入變量P，令

(8)

則有，

(9)

(10)

方程(7)可轉化為Hamilton體系下的一階常微分方程。

(11)

其中系統的轉換矩陣H為：

(12)

式(11)為結構碰撞階段動態響應的狀態方程。其一般解為：

(13)

假定接觸開始的時間點是t0時刻，Vt0是非接觸結束時刻的狀態向量。將接觸碰撞階段分成時間步長為Δt的若干時間間隔， 則tk=t0+kΔt(k=0,1,2,…)，而tk+1=tk+Δt。可以認為Δt時段內荷載按線性變化，滿足：F(τ)=F(tk)+[F(tk+1)-F(tk)](τ-tk)/Δt=F0+F1(τ-tk)，根據式(11)的解答可得V(tk)、V(tk+1)之間的轉換關系。

V(tk+1)=T[V(tk)+H-1(F0+H-1F1)]-

H-1[F0+H-1F1+F1Δt]

(14)

式中T為指數矩陣：

(15)

Ta矩陣采用泰勒級數展開求解：

(16)

式中：l為截斷階數，在碰撞問題求解中可取l=5,N=20[14-15]。求出Ta矩陣后，由于Ta矩陣中的元素非常小，若直接代入式(15)求矩陣T，則由于計算機的舍入誤差會使得指數矩陣T的精度喪失殆盡。利用指數函數的加法原理，指數矩陣T的求解最終采用下式進行。

T=I+Ta

(17)

對于非接觸碰撞階段的動力方程，同樣可以通過“增元降階”的方式轉化成如(11)式的狀態方程。此階段，系統的轉換矩陣H為：

非齊次項：

(19)

利用式(15)~式(17)可求得非接觸階段的指數矩陣T，代入式 (14)可遞推任意非碰撞時刻的狀態向量，進而求解速度和加速度。值得注意的是，除時刻為0時初始狀態為已知條件，其余非接觸階段的初始狀態是接觸階段結束時的狀態。

綜上所述，將接觸單元模型的彈性恢復力作為未知量處理，通過分解合并得到接觸碰撞階段的動力平衡方程。之所以將阻尼力作為已知量處理，采用前一時刻的值近似代替，是因為在Kelvin-Voigt模型和Jan-Hertz-damp模型中阻尼力是突變的，將其作為未知量在左端分解，會引起數值計算的不穩定。對于線彈簧模型，本文解法為精確解；對于Kelvin-Voigt模型、Hertz模型、Hertz-damp模型、Jan-Hertz-damp模型，碰撞彈簧剛度以及阻尼力是時間相關函數，其tk+1時刻值由tk時刻的值近似代替。由于Δt很小，這種近似引起的誤差也很小。精細積分法的求解分為兩個階段進行。對于不同的階段可以采用不同的積分步長。為進一步提高計算精度，需要精準的找到開始碰撞的時刻，為此，采用了線性加速度的假定。

采用積分步長為Δt，假定ti時刻u1-u2-gp<0(沒有接觸)，而ti+1時刻u1-u2-gp>0(已經接觸)，說明接觸開始時刻ti+τ介于ti時刻與ti+1時刻之間。通過精細積分法已解出ti時刻和ti+1時刻的位移、速度、加速度。假定ti時刻到ti+1時刻加速度按線性變化，則位移變化為三次曲線形式。ti+τ時刻的位移表達式為：

(20)

假定ti+τ開始接觸，即滿足u1-u2-gp=0條件：

(21)

3算例驗證

3.1理論驗證

通過對兩小球自由彈性碰撞的數值模擬，驗證精細積分法在求解碰撞問題時的精確性和高效率。兩小球質量均為2 kg，剛度均為210.25 N/m，左側小球初始位移0.1 m，釋放后自左向右運動與右側小球發生碰撞，不考慮系統的能量損失，碰撞模型[5]見圖3(a)。該模型滿足動量守恒的理論位移時程見圖3(b)。在用精細積分法求解時，接觸單元采用線性彈簧模型，碰撞彈簧剛度參數取1.547 4×106N/m。

圖4~圖6給出了Matlab平臺上的精細積分法、商業軟件Ansys 平臺上的Newmark法、Matlab平臺上Willson-θ法計算的前2 s位移響應時程和加速度響應時程，三種算法的分析過程及結果對比列于表1。

圖3　兩質點自由彈性碰撞模型及其理論解[5]Fig.3 An illustration of two point masses on free vibration and its theoretical solution

圖4　精細積分法求解位移時程及加速度時程Fig.4 Time history of displacement and acceleration with PIM

圖5 Newmark法求解位移時程及加速度時程Fig.5 Time history of displacement and acceleration with Ansys method

圖6　Willson-θ法求解位移時程及加速度時程Fig.6 Time history of displacement and acceleration with Willson-θ method

算法步長/s碰撞非碰撞計算時間/s加速度峰值/(m·s-2)質點2位移峰值/m質點2精細積分0.00010.0011.2638.40.100Newmark0.00010.001313637.40.100Willson-θ0.000010.000011201709.30.106

與理論解相比，精細積分法位移時程的周期，幅值均與理論解一致，但是碰撞的位置與理論解有微小的差異，并且隨著時間推移差異有增加的趨勢。這是因為理論解基于動量守恒定律，碰撞是在瞬間完成的，而采用線性彈簧接觸模型的碰撞力是有作用過程的，本例接觸過程約為0.002 5 s，這一過程隨著接觸彈簧剛度的增大而減小，其位移解也更接近圖5理論解。

總體上，精細積分法和Ansys 軟件Newmark法計算的位移時程響應與滿足動量守恒的理論值基本吻合，表明精細積分法和Ansys軟件Newmark法均具有較高的計算精度。Wilson-θ法計算的前2 s位移峰值與理論解相差6%，加速度峰值相差11%，且位移峰值和加速度峰值有隨時間增加的趨勢，說明該方法計算精度相對較差。

精細積分法在非碰撞階段的步長為0.001 s，碰撞階段的步長為0.000 1 s，2 s內總積分點2145個；Ansys 軟件Newmark法的積分步長與精細積分法相同，2 s內積分點數同樣為2145個，但精細積分法總計算時間僅1.2 s，而在Ansys軟件Newmark法需313 s，可知在保持相同計算精度的情況下，采用精細積分法效率得到很大提高。Wilson-θ法計算程序由于未采用變步長算法，當積分步長為0.000 01 s時，2 s內積分點數達到20 000個，計算時間為1 201 s，若積分步長取0.000 1 s，得到的位移時程及加速度時程是發散的，說明Wilson-θ法對積分步長非常敏感，不適用于接觸過程短暫的結構碰撞分析。

3.2試驗驗證

Chau進行了地震激勵下兩鋼塔的振動臺碰撞試驗[6]，文獻[13]應用基于FENAP平臺的隱式積分功能對該試驗進行了數值仿真并對比了不同接觸單元計算結果。本文基于精細積分法在Matlab平臺上編制了地震碰撞分析程序，并針對Chau試驗進行數值仿真。

兩鋼塔塔高均為2 m，1號鋼塔質量m1=204.0 kg，基頻f1=2.3 Hz，阻尼比ξ1=1.4%；2號鋼塔質量為m1=146.5 kg，基頻f1=2.9 Hz，阻尼比ξ1=1.6%；兩鋼塔間初始間隙為0.011 m。兩鋼塔固定在振動臺上，沿縱向輸入1940年El-Centro地震動記錄NS分量。該地震動記錄持時53.73 s，峰值加速度為0.348 g。分析中將鋼塔簡化為兩個單自由度體系，剛度分別為k1=4 297.5 N/m，k2=4 864.0 N/m，兩鋼塔的碰撞剛度參數及阻尼參數見表2。

表2　兩鋼塔接觸模型的參數[6,13]

數值分析中，非接觸階段精細積分時間步長為0.01 s，接觸碰撞階段的積分時間步長為0.000 01 s。

圖8為Chau試驗所測得的1號鋼塔速度響應時程。圖9(a)~(e)為基于本文精細積分法采用5種接觸單元模型計算的1號鋼塔速度響應時程。對比圖8、圖9發現，5種接觸單元模型計算的1號鋼塔速度響應時程與試驗結果趨勢吻合，但數值分析結果的速度響應峰值偏小約25%。文獻[6,13]數值分析結果也比試驗結果明顯偏小。這是由于振動臺試驗雖然輸入地震動記錄是峰值0.348 g的El-Centro波，但臺面實際加速度會產生變異。若在數值分析中采用實測振動臺臺面加速度時程，數值分析結果與試驗結果吻合程度會更高。

圖7　兩鋼塔碰撞振動臺試驗[6,13]Fig.7 Shaking table test on pounding of two steel towers

圖8　試驗測得1號鋼塔速度響應時程[6]Fig.8 Tested time history of velocity of tower 1

圖9　精細積分法計算1號鋼塔速度響應時程Fig.9 Calculated time history of velocity of tower 1 by PIM

圖10(a)~(e)給出5種接觸單元模型計算的相對壓入量-碰撞力關系曲線。圖10(a)碰撞力與相對壓入量呈線性關系，符合線彈簧接觸單元模型假定，斜率為接觸彈簧剛度；圖10(b)碰撞力包含彈性恢復力和阻尼力兩項，在接近階段碰撞力大于彈性恢復力，且在接觸瞬間出現突加碰撞力，而恢復階段碰撞力小于彈性恢復力，且在碰撞結束時出現負值，這符合Kelvin-Voigt模型特點；圖10(c)碰撞力與相對壓入量關系為1.5次拋物線，接近階段與恢復階段碰撞力路徑相同，符合Hertz模型假定；圖10(d)、(e)分別符合Hertz-damp和Jan-Hertz-damp模型假定。兩種模型均是在Hertz模型的基礎上引入阻尼項，因而接近階段與恢復階段碰撞力路徑不再重合。但不同的阻尼項假定所呈現的阻尼力特點有差別，Hertz-damp模型阻尼力隨著相對壓入量逐漸增加，在接近最大壓入量時迅速減小為0，在恢復過程阻尼力為負值；而Jan-Hertz-damp模型阻尼力呈現隨著相對壓入量快速增大，然后基本保持穩定，在接近最大壓入量時逐漸減小為0，在恢復過程不出現阻尼力的特點。圖10(a)~(e)計算結果與理論假定吻合良好，說明采用本文精細積分法能精確的模擬符合不同接觸單元模型假定的碰撞過程。

計算發現，采用線彈簧、Kelvin-Voigt模型計算的碰撞時間約為0.77 ms，采用Hertz模型、Hertz-damp模型、Jan-Hertz-damp模型碰撞時間為1.89~1.98 ms。對于線彈簧、Kelvin-Voigt模型本文選取的碰撞階段的積分步長約為碰撞歷程的1/77；對于Hertz、Hertz-damp、Jan-Hertz-damp模型積分步長約為碰撞歷程的1/200。在碰撞階段采用很小的積分步長保證了碰撞過程模擬的精度，而在非碰撞階段采用的積分步長是碰撞階段積分步長的1 000倍，極大減少了計算工作量，保證了計算效率。

圖10　碰撞力-相對壓入量關系曲線Fig.6 Curve of pounding force varying with penetration displacement

表3給出基于精細積分法的5種接觸單元模型數值計算結果與試驗結果的對比。從表3分析發現，線彈簧模型和Kelvin-Voigt模型碰撞力峰值比試驗結果明顯偏高，計算結果不合理；Hertz、Hertz-damp、Jan-Hertz-damp模型碰撞力峰值與試驗結果偏差分別為18.7%、2.3%、-10.6%，Hertz-damp模型計算結果最為接近；Hertz模型和Hertz-damp模型計算的碰撞次數與實測結果最為接近；從超過60 kN碰撞次數看，線彈簧、Kelvin-Voigt模型計算明顯不合理，而Hertz模型與試驗結果最為接近。這是因為Chau所做的振動臺碰撞試驗接近彈性碰撞，因而與Hertz模型計算結果吻合較好。

以上分析表明，本文提出的精細積分法應用于地震碰撞模擬是可行的，且具有較高的精度。為準確模擬地震碰撞，合理選擇接觸單元模型并且確定模型的參數是關鍵。在近似彈性碰撞時Hertz、Hertz-damp模型具有較好的效果。在相同的剛度參數和恢復系數下，Jan-Hertz-damp模型計算的碰撞力最小，說明其能量耗散大于Hertz-damp模型。

表3　數值計算結果與試驗結果對比

4結論

本文將5種常用接觸單元模型碰撞力表達成統一形式，并與精細積分法相結合，推導了碰撞階段的精細積分公式。在Matlab環境下編制了計算程序，進行了理論解驗證和試驗驗證，主要結論如下：

(1)精細積分法是一種顯式積分方法，對積分時間步長不敏感，在求解時將結構運動過程分為碰撞階段和非碰撞階段，在非碰撞階段求解時采用較大的步長，而在碰撞階段采用較小的步長，不但保證了碰撞模擬的精度且獲得了較高的計算效率，適于結果碰撞過程的模擬。

(2)基于線性加速度假設的碰撞時刻搜索法，可以精確求解碰撞時刻，從而進行碰撞階段和非碰撞階段的求解轉換。

(3)在地震碰撞模擬中，應合理選擇接觸單元模型和模型參數。本文采用Hertz模型及Hertz-damp模型計算結果與試驗結果較為吻合，說明在近似彈性碰撞時，這兩種模型具有較高的模擬精度。

參 考 文 獻

[1] Rosenblueth E, Meli R. The 1985 earthquake: causes and effects in Mexico City[C]//Concrete International, 1986, 8: 23-34.

[2] Priestley M J N, Seible F, Calvi G M. Seismic design and retrofit of bridge [M]. USA, New York: John Wiley and Sons INC, 1996: 8-20.

[3] Kasai K, Maison B F. Building pounding damage during the 1989 Loma Prieta earthquake. Engineering Structures 1997;19:195-207.

[4] Jankowski R. Pounding force response spectrum under earthquake excitation.[J]. Engineering Structures, 2006, 28(8): 1149-1161.

[5] Zhu Ping, Abe M, Fujino Y. Modelling three-dimensional non-linear seismic performance of elevated bridges with emphasis on pounding of girders[J].Earthquake Engineering and Structural Dynamics, 2002, 31:1891-1913.

[6] Chau K T, Wei X X, Guo X. Experimental and theoretical simulations of seismic poundings between two adjacent structures [J]. Earthquake Engineering and Structural Dynamics, 2003, 32:537-554.

[7] Chau K T, Wei X X. Pounding of structures modeled as non-linear impacts of two oscillators [J]. Earthquake Engineering and Structural Dynamics, 2001, 30:633-651.

[8] Jankowski R. Non-linear viscoelastic modeling of earthquake-induced structural pounding [J]. Earthquake Engineering and Structural Dynamics, 2005, 34(6):595-611.

[9] Zanardo G, Hao H, Modena C. Seismic response of multi-span simply supported bridges to a spatially varying earthquake ground motion [J]. Earthquake Engineering and Structural Dynamics, 2002, 31: 1325-1345.

[10] Anagnostopoulos S A. Equivalent viscous damping for modeling inelastic impacts in earthquake pounding problems [J]. Earthquake Engineering and Structural Dynamics, 2004, 33(8): 897-902.

[11] Van Mier J G M, Pruijssers A F, Reinhardt H W, et al. Load-time response of colliding concrete bodies [J].Journal of Structural Engineering (ASCE), 1991, 117(2):354-374.

[12] Muthukumar S, DesRoches R. A Hertz contact model with non-linear damping for pounding simulation [J].Earthquake Engineering and Structural Dynamics, 2006,35(7): 811-828.

[13] 禚一,李忠獻,王菲. 橋梁地震碰撞分析中不同接觸單元模型的對比分析[J].工程力學, 2014, 31(3): 11-17.

ZHUO Yi, LI Zhong-xian,WANG Fei. Comparative analysis of different contact element models in seismic pounding analysis of bridges. [J]. Engineering Mechanics, 2014, 31(2): 11-17.

[14] 郭澤英,李青寧,張守軍. 結構地震發應分析的一種新精細積分法[J]. 工程力學, 2007, 24(4): 35-40.

GUO Ze-ying,LI Qing-ning,ZHANG Shou-jun. A new precise integration method for structure seismic response analysis.[J]. Engineering Mechanics, 2007,24(4):35-40.

[15] 鐘萬勰.結構動力方程的精細時程積分法[J].大連理工大學學報,1994, 34 (2): 131-136.

ZHONG Wan-xie. One precise time integration method for structural dynamics[J]. Journal of Dalian University of Technology, 1994, 34 (2): 131-136.

[16] 鮑四元,鄧子辰. 結構撞擊響應的一種彈性模型及其精細求解[J].工程力學,2008,25 (60):14-17.

BAO Si-yuan, DENG Zi-chen. An elastic model and its precise solution for structure impact response[J]. Engineering Mechanics, 2008,25 (60):14-17.

[17] 鄒立華,郭潤,黃凱,等.帶預應力橡膠支座相鄰隔震結構碰撞分析[J].振動與沖擊,2014,33(9):131-136.

ZOU Li-hua, GUO Run, HUANG Kai, et al. Pounding of adjacent isolated-structures with prestressed rubber bearings[J]. Journal of Vibration and Shock,2014,33(9):131-136.

Precise integration method for the solution of contact element model in earthquake pounding analysis

ZHANGRui-jie,LIQing-ning,YINJun-hong

(School of Civil Engineering, Xi’an University of Architecture and Technology, Xi’an 710055, China)

Abstract:The numerical analysis of earthquake pounding is one of the important points in bridges’ seismic study. A precise integration method was introduced into the numerical analysis of the structures’ earthquake pounding. In consideration of the characteristics of five kinds of contact element models that are commonly used, a unified expression of pounding force was put forward, and then the formulas for the precise integration method were deduced. The structures’ response process in earthquake was divided into two parts, the one is collision process and the other is none collision process. And different integration time steps were applied in different processes as well as a method to search precise contact moment was developed. The precise integration method program has been tested by using a two lumped masses on free vibration model. And then the numerical simulation of two steel structures’ earthquake pounding experiments implemented by Chau was carried out based on the presented precise integration method. The conclusions are: the precise integration method not only can ensure the precision of the pounding simulation but also has high computing efficiency; the contact element model as well as the model’s parameters should be choosen carefully; the numerical analysis results by the Hertz model and Hertz-damp model are in good agreement with the experiment results in the case of approximate elastic pounding.

Key words:earthquake; structure pounding; numerical analysis; precise integration method; contact element model

中圖分類號:TU973;U442.5

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.03.019

收稿日期：2014-09-10修改稿收到日期：2014-12-19

基金項目：國家自然科學基金(51078306)；高等學校博士學科點專項科研基金(20106120110004)；西安建筑科技大學重大科技項目創新基金(ZX0901)；陜西省重點學科建設專項資金資助項目(E01004)

第一作者 張瑞杰 男，博士生，講師，1982年生