基于精細(xì)積分法的相鄰結(jié)構(gòu)彈塑性地震碰撞分析

張瑞杰, 李青寧, 王天利, 葉 毅, 孫建鵬

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

罕遇地震發(fā)生時，結(jié)構(gòu)可能因較大的變形而產(chǎn)生開裂、屈服和相對滑動，從而使結(jié)構(gòu)進(jìn)入塑性狀態(tài)。同時，因動力特性的差異或行波效應(yīng)的影響，相鄰結(jié)構(gòu)的位移差值可能超過結(jié)構(gòu)預(yù)留間隙，從而產(chǎn)生碰撞現(xiàn)象。這種相鄰結(jié)構(gòu)的碰撞，往往會產(chǎn)生巨大的破壞作用，如房屋建筑的垮塌、橋梁結(jié)構(gòu)的局部損壞甚至落梁[1-3]。結(jié)構(gòu)進(jìn)入塑性狀態(tài)是材料非線性問題，相鄰結(jié)構(gòu)碰撞是狀態(tài)非線性問題。對于罕遇地震下相鄰結(jié)構(gòu)的地震時程反應(yīng)分析，應(yīng)同時考慮結(jié)構(gòu)的彈塑性和碰撞。時程反應(yīng)分析中的相鄰結(jié)構(gòu)碰撞模擬方法多采用接觸單元法[4]。該方法是在兩可能接觸點之間添加具有一定剛度和阻尼特性的接觸單元，當(dāng)接觸發(fā)生時激活接觸單元。實質(zhì)上該方法是將碰撞接觸過程處理為碰撞力的作用過程，根據(jù)碰撞力計算假定的不同，發(fā)展了多種接觸單元模型。

對于結(jié)構(gòu)地震響應(yīng)過程，區(qū)分為分離狀態(tài)和碰撞狀態(tài)，并且，對不同的狀態(tài)可分別采用結(jié)構(gòu)切線剛度矩陣建立增量動力平衡方程或割線剛度矩陣建立全量動力平衡方程。動力平衡方程的常用數(shù)值積分方法，有中心差分法、Newmark-β法、Wilson-θ法等方法。近年來，求解結(jié)構(gòu)動力響應(yīng)的精細(xì)積分法[5]的提出，為動力平衡方程的求解增添了新的思路，因而，一經(jīng)提出就得到了眾多學(xué)者的重視和發(fā)展，并很快推廣到非線性動力平衡方程的求解[6]。申永康等[7]應(yīng)用精細(xì)積分法對建筑結(jié)構(gòu)的爆破地震反應(yīng)進(jìn)行彈塑性分析，郭澤英等[8]對精細(xì)積分進(jìn)行改進(jìn)，并將精細(xì)積分法應(yīng)用于短肢剪力墻結(jié)構(gòu)的彈塑性地震反應(yīng)分析。張瑞杰等[9～12]將精細(xì)積分法引入到地震碰撞動力平衡方程的數(shù)值求解，但只適用兩個單自由度的彈性碰撞問題。兩個單自由度體系的彈性地震碰撞僅有一個接觸單元，算法設(shè)計相對較簡單。對于多個相鄰結(jié)構(gòu)，地震發(fā)生時，有多個可能的接觸碰撞點，且結(jié)構(gòu)可能進(jìn)入塑性狀態(tài)。該類結(jié)構(gòu)可以簡化為多個單自由度結(jié)構(gòu)彈塑性多點碰撞模型，疊加了狀態(tài)非線性和材料非線性問題。在前人研究的基礎(chǔ)上，進(jìn)一步將精細(xì)積分法推廣到多個單自由度結(jié)構(gòu)的彈塑性多點碰撞動力平衡方程的求解。

1 彈塑性碰撞動力平衡方程及其精細(xì)積分解法

1.1 彈塑性動力平衡方程

(1)

同時，結(jié)構(gòu)反應(yīng)滿足全量動力平衡方程：

(2)

若t時刻為接觸碰撞狀態(tài)，全量動力平衡方程為：

(3)

圖1 單自由度結(jié)構(gòu)多點碰撞模型Fig.1 Multi point collision model of SDOF structures

圖2 不考慮剛度退化雙直線模型Fig.2 Bi-linear model without considering stiffness degradation

圖3 考慮剛度退化雙直線模型Fig.3 Bi-linear model considering stiffness degradation

1.2 分離狀態(tài)增量動力平衡方程的精細(xì)積分法

首先討論分離狀態(tài)增量動力平衡方程的精細(xì)積分法。對方程(1)進(jìn)行降階，引入向量p，令

(4)

則有，

(5)

方程(1)可變換為：

(6)

(7)

式中:t時刻的系統(tǒng)矩陣H為：

(8)

在分離階段，將積分過程分成時間步長為Δt的若干時間間隔，則Δtk=tk-tk-1、Δtk+1=tk+1-tk是相鄰的兩個時間步，已知Δtk時間步的增量狀態(tài)時，Δtk+1時間步的增量狀態(tài)為：

(9)

Δt時段內(nèi)非齊次項按線性變化，滿足：f(τ)=f(tk)+[f(tk+1)-f(tk)](τ-tk)/Δt=f(tk)+Δf(τ-tk)，Δf表示非齊次項tk至tk+1時段的變化率。根據(jù)式(9)的解答可得vΔtk、vΔtk+1之間的轉(zhuǎn)換關(guān)系。

vΔtk+1=T(tk)[vΔtk+H(tk)-1(f(tk)+H(tk)-1Δf)]-
H(tk)-1[f(tk)+H(tk)-1Δf+Δf·Δt]

(10)

以上討論針對增量動力平衡方程，對于全量動力平衡方程精細(xì)積分法求解原理是相同的。應(yīng)注意，對于理想彈塑性模型切線剛度為0的屈服階段，增量法是不適用的。精細(xì)積分法是一種顯式的遞推方法，對于增量動力平衡方程，是從前一個時間步的結(jié)構(gòu)位移、速度、加速度增量，遞推其后一時間步的位移、速度、加速度增量。這就要求，采用增量法的時間步長應(yīng)該相等，且前一個時間步的增量狀態(tài)已知。當(dāng)前時刻的位移、速度、加速度是之前積分步長增量累加結(jié)果，應(yīng)同時滿足結(jié)構(gòu)全量動力平衡方程(2)。同樣的，通過全量動力平衡方程精細(xì)積分法遞推出結(jié)構(gòu)位移、速度、加速度后，也可計算出時間步內(nèi)結(jié)構(gòu)位移、速度、加速度的增量，滿足增量動力平衡方程(1)。這是全量法和增量法能夠混合使用的理論基礎(chǔ)。

1.3 碰撞動力平衡方程精細(xì)積分法

接觸單元模型的碰撞力可分解為彈性恢復(fù)力和阻尼力兩部分，彈性恢復(fù)力表達(dá)為碰撞彈簧剛度與侵徹位移的乘積，阻尼力表達(dá)為阻尼系數(shù)與速度差的乘積。根據(jù)碰撞彈簧剛度和阻尼系數(shù)的不同假設(shè)，誕生了不同的接觸單元模型，如線性彈簧模型[14]、線性彈簧-阻尼(Kelvin-Voigt)模型[15]、Hertz模型[16]、非線性彈簧阻尼(Hertz-damp)模型[17]、Jan-Hertz-damp模型[18-19]。這些接觸單元的碰撞力可統(tǒng)一表達(dá)為：

(11)

(12)

基于大量的碰撞試驗，Jankowski認(rèn)為接觸碰撞階段可劃分為接近階段和恢復(fù)階段，能量損失主要發(fā)生在接近階段。基于此，對Hertz-damp模型的阻尼力項進(jìn)行了改進(jìn)，即Jan-Hertz-damp模型。該模型仍可表達(dá)成式(11)形式。式中阻尼系數(shù)cp為：

(13)

cp=0

(14)

假設(shè)時刻圖3模型中所有接觸點均為接觸碰撞狀態(tài)，方程(3)中碰撞力列向量的表達(dá)式為：

(15)

對于更一般的接觸碰撞情況，假設(shè)n-1個接觸單元中的第i個接觸點未產(chǎn)生接觸碰撞，只需將式(15)中kpi、cpi置為0即可。以往情況下式(15)作為方程2的非齊次項處理。由于求解時刻的位移、速度未知，采用前一時間步的位移、速度來近似計算碰撞彈簧剛度、阻尼系數(shù)以及碰撞力。這里介紹另一種處理方式，將式(15)代入方程(3)，其中碰撞力的彈性恢復(fù)力在左端進(jìn)行分解并與結(jié)構(gòu)恢復(fù)力項合并，而阻尼力仍作為非齊次項處理，由前一時間步位移、速度近似計算。對于線性彈簧模型和Kelvin-Voigt模型完全消除了彈性恢復(fù)力項的迭代誤差。之所以將阻尼力作為非齊次項處理，是因為對于Kelvin-Voigt模型，阻尼力是突變的，而Jan-Hertz-damp阻尼模型阻尼力是位移差和速度差的函數(shù)，在碰撞開始時刻迅速增大，隨后減小。實踐證明，這兩種模型將阻尼力作為未知量在左端分解，會引起數(shù)值計算的不穩(wěn)定。經(jīng)過處理后，碰撞動力平衡方程(3)進(jìn)一步寫為如下形式：

(16)

對碰撞動力平衡方程(16)的精細(xì)積分法推導(dǎo)過程與方程(1)類似，不再贅述。碰撞狀態(tài)的動力平衡方程也可以采用結(jié)構(gòu)切線剛度矩陣建立成增量形式，但要求在第一個時間步采用全量法精細(xì)積分，獲得增量方程精細(xì)積分的啟動條件。由于碰撞持續(xù)過程非常短暫，結(jié)構(gòu)割線剛度變化較小，所以體現(xiàn)不出增量法的優(yōu)勢。

1.4 界點的處理

對于n個單自由度彈塑性結(jié)構(gòu)，均可定義獨立的彈塑性模型。以圖1所示不考慮剛度退化的雙直線模型為例說明，在載荷過程中，存在A、B、C、D、E界點。A界點是彈性階段正向加載時屈服點；B界點是正向屈服至反向卸載階段界點，判斷條件是速度由正變負(fù)；C界點是反向卸載屈服界點；D界點是反向卸載至正向加載界點，判斷條件是速度由負(fù)變正；E界點是正向加載屈服界點。在界點處結(jié)構(gòu)反應(yīng)滿足連續(xù)性條件。當(dāng)荷載處于OA彈性狀態(tài)時，求解可以采用全量法進(jìn)行。當(dāng)荷載通過A界點進(jìn)入塑性狀態(tài)時，為避免迭代可采用增量法求解，但應(yīng)注意理想彈塑性模型切線剛度為0的屈服階段不適用。在當(dāng)前時刻，根據(jù)精細(xì)積分求解的結(jié)構(gòu)位移、速度，可對各結(jié)構(gòu)彈塑性階段進(jìn)行判斷。當(dāng)結(jié)構(gòu)彈塑性階段改變，即意味著界點產(chǎn)生。為了精確吻合滯回曲線，需要對界點產(chǎn)生的時刻以及界點處的結(jié)構(gòu)反應(yīng)進(jìn)行精確分析。

(17)

而B、D界點，判定條件是ti+τ時刻速度為0，可建立如下方程：

(18)

方程式(17)為一元三次方程，可采用二分法求解[0 Δt]內(nèi)的解。方程式(18)為一元二次方程，可直接求解[0 Δt]內(nèi)的解。ti+τ時刻即為精準(zhǔn)界點產(chǎn)生時刻，進(jìn)而可按線性加速度假定求解該時刻的位移、速度、加速度。舍棄ti+1時刻求解結(jié)果，存儲ti+τ時刻結(jié)果。采用增量法時，由于通過界點后的第一個時間步的增量狀態(tài)未知，應(yīng)通過全量法迭代求解一次，為增量法的啟動創(chuàng)造初始條件。

1.5 分離碰撞狀態(tài)轉(zhuǎn)換

將結(jié)構(gòu)反應(yīng)過程分為分離狀態(tài)和接觸碰撞狀態(tài)，兩種狀態(tài)下的積分步長分別為Δt1、Δt2。精細(xì)積分法對步長不敏感，對分離狀態(tài)積分步長設(shè)定的可以大些。而對于接觸碰撞階段，由于碰撞過程本身很短暫，要采用較小的積分步長，該步長取值可通過試算后確定。

為進(jìn)一步提高計算精度，需要精準(zhǔn)的找到開始接觸和分離的時刻。假定ti時刻第k個接觸單元uk-uk+1-gpk<0(沒有接觸)，而ti+1時刻uk-uk+1-gpk>0(已經(jīng)接觸)，說明接觸開始時刻ti+τ介于時刻ti至ti+1時刻之間。已解出ti時刻和ti+1時刻的位移、速度、加速度。根據(jù)線性加速度假定，ti+τ時刻的位移滿足：

(19)

方程(19)為一元三次方程，可以利用二分法求解τ在區(qū)間[0 Δt]內(nèi)的解，ti+τ時刻即為碰撞起始時刻。程序設(shè)計時，需要在每一步遞推前對n-1個接觸單元循環(huán)判定是否接觸。一旦判定有分離-接觸碰撞狀態(tài)轉(zhuǎn)換，首先求解精確的碰撞時刻ti+τ和該時刻的結(jié)構(gòu)反應(yīng)，舍棄ti+1的結(jié)果。從ti+τ時刻積分步長變?yōu)榻佑|碰撞階段的積分步長，進(jìn)入碰撞階段求解。接觸碰撞狀態(tài)至分離狀態(tài)轉(zhuǎn)換的分離時刻可以采用同樣的原理計算。

需要注意的是，由于是多點碰撞，在當(dāng)前積分點可能有多個接觸點同時轉(zhuǎn)為接觸碰撞狀態(tài)，也可能在當(dāng)前接觸點未分離的情況下有新的接觸點進(jìn)入接觸碰撞狀態(tài)。對多接觸點判定產(chǎn)生接觸情況，通過求解各自的精確碰撞時刻，以最小值作為狀態(tài)轉(zhuǎn)換時刻。接觸點分離時，每一接觸點分離均求解精確分離時刻。所有接觸點均分離后，步長變化為分離狀態(tài)的步長。同時，對于每一積分步接觸單元狀態(tài)的判定作為組集剛度矩陣的依據(jù)，未產(chǎn)生碰撞的接觸單元相應(yīng)的彈簧剛度系數(shù)、阻尼系數(shù)置為0。對于結(jié)構(gòu)割線剛度，應(yīng)在每一積分步迭代求解。

2 算法設(shè)計

單自由度多點彈塑性地震碰撞分析程序設(shè)計較為復(fù)雜，整體上劃分為分析準(zhǔn)備階段和積分計算階段。算法具體流程見圖4。

圖4 單自由度彈塑性結(jié)構(gòu)多點地震碰撞算法流程圖Fig.4 Flow chart of multi point seismic collision algorithm for elastic-plastic SDOF structures

3 算例驗證

某三跨簡支梁橋墩梁間采用盆式橡膠支座，一端縱向固結(jié)一端活動，力學(xué)圖式見圖5。該模型簡化為三個單自由度體系，有兩個可能的鄰梁碰撞點，采用Kelvin-Voigt接觸單元模型。橋墩的彈塑性模型假定為不考慮剛度退化的雙直線模型。該例的結(jié)構(gòu)參數(shù)、雙直線模型參數(shù)、接觸單元模型參數(shù)見表1。沿縱向輸入1940年El-Centro地震動記錄NS分量。該地震動記錄持時53.73 s，峰值加速度調(diào)整為0.5 g，地震波步長0.02 s。

圖5 某三跨簡支梁地震碰撞模型Fig.5 The earthquake collision model of three spans simply supported beam

結(jié)構(gòu)參數(shù)結(jié)構(gòu)質(zhì)量m/kg初始剛度k/(N·m-1)阻尼比/%雙直線模型參數(shù)屈服位移屈服后剛度接觸單元參數(shù)碰撞剛度kP/(N·m-1)恢復(fù)系數(shù)e間隙gP/m1221 0004.8×10750.044.8×1062221 0001.9×10750.071.9×1063221 0003.2×10750.063.2×1064.7×1090.650.084.7×1090.650.08

采用精細(xì)積分法分析時，設(shè)定三結(jié)構(gòu)分離狀態(tài)積分步長為0.002 s，接觸碰撞狀態(tài)為0.000 2 s。同時為驗證算法的正確性，采用Ansys軟件進(jìn)行彈塑性地震反應(yīng)分析，積分步長為0.001 s。其中單自由度結(jié)構(gòu)和接觸單元均采用Combin40非線性彈簧單元模擬。

圖6～圖9分別給出了兩種方法分析的三跨簡支梁的相對位移時程、加速度時程、結(jié)構(gòu)滯回曲線、碰撞力時程。對比兩種方法的計算結(jié)果，相對位移時程曲線、加速度時程曲線、碰撞力時程曲線均是重合的。兩方法分析的三跨結(jié)構(gòu)的位移極值相差小于1%，結(jié)構(gòu)1～結(jié)構(gòu)3最大加速度值相差4.2%、1.5%、1.5%，最小值相差1.7%、1.4%、3.2%，接觸點1和接觸點2最大碰撞力值相差分別為0.7%、1.9%。滯回曲線對比表明，兩種方法滯回曲線趨勢一致，精細(xì)積分法具有更清晰的界點，結(jié)構(gòu)滯回曲線完全符合不考慮剛度退化的彈塑性模型假定。通過Ansys驗證表明，彈塑性結(jié)構(gòu)多點碰撞程序計算結(jié)果準(zhǔn)確，且計算速度更快。

圖6 三跨簡支梁位移時程Fig.6 Displacement time history of the three spans simply supported beam

圖7 三跨簡支梁加速度時程Fig.7 Acceleration time history of the three spans simply supported beam

圖8 三結(jié)構(gòu)的滯回曲線Fig.8 Hysteretic curve of three structures

圖9 兩接觸點碰撞力時程Fig.9 Pounding force time history of the two contact points

地震結(jié)束后結(jié)構(gòu)并沒有恢復(fù)到最初位置。這是因為，當(dāng)結(jié)構(gòu)位移超過屈服位移時，即進(jìn)入塑性狀態(tài)，恢復(fù)力與位移呈現(xiàn)非線性關(guān)系，并且卸載后結(jié)構(gòu)有殘余變形。當(dāng)兩相鄰結(jié)構(gòu)的相對位移差值超出結(jié)構(gòu)間隙時，兩結(jié)構(gòu)即發(fā)生碰撞，接觸點1產(chǎn)生2次碰撞、接觸點2產(chǎn)生4次碰撞，6次碰撞的接觸持時0.015 5～0.016 4 s。本例接觸碰撞階段的積分步長約為接觸持時的1/80，能精確的模擬結(jié)構(gòu)的接觸碰撞過程。從加速度時程可以分析出，碰撞時結(jié)構(gòu)的加速度急劇變化，碰撞力遠(yuǎn)大于結(jié)構(gòu)恢復(fù)力，這對于易受沖擊損壞的伸縮縫、固定支座、梁端、橋臺背墻等是非常不利的。

4 接觸單元參數(shù)影響分析

Kelvin-Voigt接觸單元模型輸入?yún)?shù)包括碰撞彈簧剛度、間隙和恢復(fù)系數(shù)，其中恢復(fù)系數(shù)決定阻尼特性。為研究碰撞彈簧剛度、間隙和恢復(fù)系數(shù)對結(jié)構(gòu)碰撞的影響，對于圖5結(jié)構(gòu)假定了2組碰撞剛度、3組結(jié)構(gòu)間隙、3組恢復(fù)系數(shù)，組合了18組接觸單元模型參數(shù)，地震輸入和結(jié)構(gòu)參數(shù)不變。表2給出了18組不同接觸單元參數(shù)對應(yīng)的兩個接觸點的碰撞次數(shù)、碰撞力峰值。

表2 不同接觸單元參數(shù)地震碰撞分析結(jié)果對比

分析表2數(shù)據(jù)發(fā)現(xiàn)：剛度參數(shù)、恢復(fù)系數(shù)不變時，隨著結(jié)構(gòu)間隙的逐漸增加，碰撞發(fā)生的次數(shù)逐漸減少，結(jié)構(gòu)的碰撞力峰值沒有明確一致的變化規(guī)律。如間隙從0.06 m增加到0.08 m，接觸點1碰撞力峰值減小而接觸點2碰撞力峰值增加；增加到0.10 m，兩接觸點碰撞力峰值均減小。剛度參數(shù)、結(jié)構(gòu)間隙不變時，隨著恢復(fù)系數(shù)的減小，意味著碰撞阻尼增加，碰撞發(fā)生的次數(shù)基本不變，碰撞力峰值均減小?；謴?fù)系數(shù)從0.85減小為0.75，接觸點1的碰撞力峰值減小了3.8%～5.2%，接觸點2的碰撞力峰值減小4.6%～14.7%；減小為0.65，接觸點1的碰撞力峰值減小了7.7%～10.1%，接觸點2的碰撞力峰值減小9.7%～28.6%。結(jié)構(gòu)間隙、恢復(fù)系數(shù)不變時，接觸單元碰撞剛度減小1個數(shù)量級，碰撞次數(shù)基本不變，而碰撞力峰值大幅度減小。接觸點1碰撞力峰值減小了67.8%～70.0%，接觸點2碰撞力峰值減小了68.2%～73.3%。

由上述分析可推知：增加結(jié)構(gòu)的間隙可以減小碰撞發(fā)生的可能性，甚至消除碰撞現(xiàn)象；而當(dāng)碰撞不可避免發(fā)生時，通過增加柔性緩沖裝置，可大幅度降低碰撞力，從而減小碰撞對結(jié)構(gòu)的影響；增加碰撞過程的阻尼對于消減碰撞次數(shù)碰撞力意義不大。

6 結(jié) 論

(1)精細(xì)積分法是顯式的遞推方法，適用于求解彈塑性結(jié)構(gòu)多點地震碰撞問題。結(jié)構(gòu)進(jìn)入塑性狀態(tài)后，增量動力方程精細(xì)積分可以避免全量方程的迭代求解和割線剛度為0或無窮大情況下的數(shù)值求解困難。在界點處，通過精確求解界點轉(zhuǎn)換時刻和對應(yīng)的狀態(tài)，使得結(jié)構(gòu)的滯回曲線完全符合不考慮剛度退化的雙直線彈塑性模型假定。

(2)結(jié)構(gòu)碰撞過程非常短暫，通過全量方程精細(xì)積分法求解。為分離狀態(tài)和接觸碰撞狀態(tài)設(shè)定了兩種積分步長，兼顧了整體計算效率和碰撞計算精度。并且，通過求解精確的分離碰撞狀態(tài)轉(zhuǎn)換時刻和對應(yīng)的狀態(tài)，從而能夠精確模擬結(jié)構(gòu)的碰撞行為，得到符合接觸單元模型假定的碰撞動力響應(yīng)。

(3)通過對三跨簡支梁算例進(jìn)行Ansys分析，結(jié)果驗證了多點彈塑性碰撞精細(xì)積分算法和程序的正確性。

(4)三跨簡支梁地震碰撞分析表明，碰撞使得結(jié)構(gòu)加速度在短時間內(nèi)急劇變化，橋梁伸縮縫、固定支座、梁端、橋臺背墻等易受沖擊損壞。

(5)對接觸單元模型參數(shù)分析表明：增加結(jié)構(gòu)的間隙可以減少碰撞發(fā)生的次數(shù)，甚至消除碰撞現(xiàn)象；而當(dāng)碰撞不可避免發(fā)生時，通過增加柔性緩沖裝置，可大幅度降低碰撞力，從而減小碰撞對結(jié)構(gòu)的影響。