基于精細積分法的相鄰結構彈塑性地震碰撞分析

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

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

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

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

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

1.1 彈塑性動力平衡方程

(1)

同時，結構反應滿足全量動力平衡方程：

(2)

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

(3)

圖1 單自由度結構多點碰撞模型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 分離狀態增量動力平衡方程的精細積分法

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

(4)

則有，

(5)

方程(1)可變換為：

(6)

(7)

式中:t時刻的系統矩陣H為：

(8)

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

(9)

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

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)

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

1.3 碰撞動力平衡方程精細積分法

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

(11)

(12)

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

(13)

cp=0

(14)

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

(15)

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

(16)

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

1.4 界點的處理

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

(17)

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

(18)

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

1.5 分離碰撞狀態轉換

將結構反應過程分為分離狀態和接觸碰撞狀態，兩種狀態下的積分步長分別為Δt1、Δt2。精細積分法對步長不敏感，對分離狀態積分步長設定的可以大些。而對于接觸碰撞階段，由于碰撞過程本身很短暫，要采用較小的積分步長，該步長取值可通過試算后確定。

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

(19)

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

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

2 算法設計

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

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

3 算例驗證

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

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

結構參數結構質量m/kg初始剛度k/(N·m-1)阻尼比/%雙直線模型參數屈服位移屈服后剛度接觸單元參數碰撞剛度kP/(N·m-1)恢復系數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

采用精細積分法分析時，設定三結構分離狀態積分步長為0.002 s，接觸碰撞狀態為0.000 2 s。同時為驗證算法的正確性，采用Ansys軟件進行彈塑性地震反應分析，積分步長為0.001 s。其中單自由度結構和接觸單元均采用Combin40非線性彈簧單元模擬。

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

圖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 三結構的滯回曲線Fig.8 Hysteretic curve of three structures

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

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

4 接觸單元參數影響分析

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

表2 不同接觸單元參數地震碰撞分析結果對比

分析表2數據發現：剛度參數、恢復系數不變時，隨著結構間隙的逐漸增加，碰撞發生的次數逐漸減少，結構的碰撞力峰值沒有明確一致的變化規律。如間隙從0.06 m增加到0.08 m，接觸點1碰撞力峰值減小而接觸點2碰撞力峰值增加；增加到0.10 m，兩接觸點碰撞力峰值均減小。剛度參數、結構間隙不變時，隨著恢復系數的減小，意味著碰撞阻尼增加，碰撞發生的次數基本不變，碰撞力峰值均減小?；謴拖禂祻?.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%。結構間隙、恢復系數不變時，接觸單元碰撞剛度減小1個數量級，碰撞次數基本不變，而碰撞力峰值大幅度減小。接觸點1碰撞力峰值減小了67.8%～70.0%，接觸點2碰撞力峰值減小了68.2%～73.3%。

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

6 結 論

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

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

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

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

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