趙鵬舉，于曉輝，陸新征

(1.哈爾濱工業(yè)大學結構工程災變與控制教育部重點實驗室，黑龍江，哈爾濱 150090;2.清華大學土木工程安全與耐久教育部重點實驗室，北京 100084)

汶川地震的震害調查發(fā)現(xiàn)，許多鋼筋混凝土(RC)框架結構，如居民住房和教學樓等，發(fā)生了較為嚴重的破壞甚至倒塌[1－2]。由于建筑結構發(fā)生倒塌是造成人員傷亡的主要原因，因此，合理評估建筑結構在地震作用下的抗倒塌能力一直是地震工程研究的熱點。近年來，針對我國量大面廣的RC框架結構的抗地震倒塌的設計[3]和能力評估[4]開展了一系列研究。

目前，結構抗地震倒塌能力評估的主要方法是增量動力分析(Incremental Dynamic Analysis, IDA)方法[5]。該方法選取一組地震動作為地震輸入，通過對每條地震動進行調幅，最終獲得結構發(fā)生倒塌所對應的地震動強度，并將其作為結構的抗地震倒塌能力。考慮地震動之間的不確定性，最終獲得的結構抗地震倒塌能力通常體現(xiàn)為地震倒塌易損性的形式。我國學者采用IDA方法對RC框架結構的抗地震倒塌能力開展了一系列研究。例如：范萍萍等[4]對不同抗震設防烈度和抗震等級的RC框架結構進行了地震倒塌易損性分析，比較了不同抗震等級下框架結構的抗地震倒塌能力的差異。呂大剛[6]采用IDA方法，比較了分別采用FEMA準則和結構動力失穩(wěn)點作為結構倒塌判據(jù)的RC框架結構地震倒塌易損性。于曉輝和呂大剛[7]將可靠度分析中的均值一次二階矩法與IDA方法相結合，提出了一種隨機IDA方法來考慮結構不確定性對地震倒塌易損性的影響。

要進行IDA分析，一個關鍵的問題是選擇合適的結構動力學數(shù)值計算方法。在已有的結構動力學數(shù)值計算方法中，逐步積分法由于其簡單高效的特點被廣泛采用。逐步積分法分為隱式算法與顯式算法。其中，顯式算法無須求解聯(lián)立方程組，對于多自由度結構以及非線性分析問題，與隱式算法相比，顯式算法有更高的計算效率[8－9]。目前，顯式算法已經(jīng)被一些商業(yè)化軟件所采用[10]，并在土木工程領域得到了一些應用。例如：陳國興等[11]針對地鐵地下結構比較了ABAQUS在并行計算時隱式算法與顯式算法的計算效率與精度，表明顯式算法具有遠高于隱式算法的計算效率；梁艷晨等[12]基于ANSYS顯式動力模塊模擬了碎石墊層地基對框架結構抗震能力的提高。近年來，作為一種免費開放式有限元模擬軟件平臺，OpenSees軟件得到了廣泛應用。張書豪[13]基于蛙跳法提出了一種新的顯式積分算法，并將該算法集成于OpenSees軟件中，同時對該算法的精確性進行了驗證。由于顯式方法不存在隱式算法常有的收斂性問題，其在非線性分析計算時能更好地模擬結構的極限狀態(tài)，也更加適合進行結構的抗地震倒塌能力的評估。

鑒于顯式算法在進行結構抗地震倒塌能力分析中的優(yōu)勢，本文以OpenSees軟件作為分析平臺，采用文獻[13]所提出的基于蛙跳法的顯式積分算法，通過IDA方法對按我國規(guī)范設計的RC框架結構的抗地震倒塌能力進行了評估。在分析過程中，考慮了4種不同倒塌判據(jù)的影響，包括：我國《建筑抗震設計規(guī)范》[14]推薦的結構層間位移角限值θmax=1/50(判據(jù)1)、美國FEMA 356[15]推薦的θmax=4%(判據(jù)2)、文獻[16]建議的IDA曲線切線剛度退化為初始剛度的20%(判據(jù)3)以及結構主要構件豎向倒塌位移超過1 m(判據(jù)4)。其中，當采用倒塌判據(jù)4時，結構已經(jīng)進入強非線性，此時用顯式算法可以很好地獲得結構在這一狀態(tài)附近的非線性狀態(tài)，而隱式算法則會由于收斂性的問題無法逼近這一極限狀態(tài)。

1 基于蛙跳法的顯式算法基本原理

顯式算法與隱式算法的主要差別在于顯式算法是由已知數(shù)據(jù)來表示下一積分步的狀態(tài)量(如式(1))，而隱式算法對下一積分步的狀態(tài)量不能直接給出公式，需要通過解方程來求解(如式(2))[17]。

式中：ωi、ti、ωi+1與ti+1分別為第i個積分步與第i+1個積分步的狀態(tài)量；y0為初值；h為步長；f表示導函數(shù)。

中心差分法[18]是基本的顯式積分算法，其以有限差分在數(shù)值上代替求導。對于運動方程：

將其位移項進行泰勒展開可得：

將式(4)和式(5)相加和相減并舍去高階小量O(Δt2)可得：

將式(6)和式(7)代入式(3)中，即可得：

由式(8)、式(9)和式(10)可見，對于顯式算法，只要給定初值x0以及x-Δt即可依次求解之后任一時刻的狀態(tài)量。其中x-Δt可由式(5)舍去高階小量得到：

式中：x0和為初始條件；可由式(3)求得(t=0)。

蛙跳法在1970年由Hockney[19]提出，在分子動力學領域有廣泛應用[20－21]。蛙跳法可以認為是中心差分法的一種變異格式[18]，其與中心差分法的差別在于蛙跳法是在半個時間步上求解速度，如式(12)和式(13)所示。

由式(12)與式(13)可見，速度與位移和加速度均未在同一時刻定義，故在將式(12)與式(13)代入運動方程式(3)時，可將代替，獲得：

根據(jù)式(14)，只需已知x0和即可依次求解之后任意時刻的狀態(tài)量。在式(14)中將代替，采用了速度近似，會導致算法精度低于傳統(tǒng)的中心差分法(只有一階精度)[13]。此外，使用該方法會造成系統(tǒng)的勢能與動能不能在同一時刻定義[22]。為解決上述問題，張書豪[13]提出了一種修正的蛙跳格式顯式算法。該方法首先計算xt+Δt：

計算得到xt+Δt后，繼而利用中心差分法可計算得t時刻的速度與加速度(見式(18)和式(19))，再利用蛙跳法基本步驟式(12)和式(13)，即可得修正后的速度和位移。

該顯式算法對單自由度彈性體系的穩(wěn)定性條件如式(20)所示[13,18]。其中，顯式積分極限步長Δt需滿足：

式中：ζ為單自由度體系的阻尼比；ω為單自由度體系的自振頻率。

對于有限單元體系，顯式積分極限步長Δt應取所有單元極限步長的最小值。單元對應的極限步長與單元的特征長度與單元材料波速有關，梁柱單元以及桁架單元的極限步長為[23]：

式中：L為梁柱單元長度；V為單元材料波速，取壓縮波速；E為材料彈性模量；ρ為單元密度。

與基于蛙跳法的顯式算法相配合使用的阻尼格式為質量阻尼與振型阻尼的組合形式，采用這種阻尼形式可以有效消耗大量的計算資源，能夠確保方程解耦并且有效抑制不合理的高階振型的影響[13, 24]。

2 RC框架結構的設計和建模

本文采用文獻[25]中6個RC框架結構作為結構對象，包括2個5層結構，2個8層結構，2個10層結構。其中，每一高度結構采用2種地震設防烈度：Ⅶ度(0.1g)和Ⅷ度(0.2g)。為方便以后論述，將上述結構編號為RC5-1、RC5-2、RC8-1、RC8-2、RC10-1、RC10-2。其中，RC后第一個數(shù)字代表樓層數(shù)，RC后第2個數(shù)字中，1代表設防烈度為Ⅶ度(0.1g)，2代表設防烈度為Ⅷ度(0.2g)。所有結構的底層層高為4.4 m，其余層高3.3 m。所有結構均按照設計地震動分組為第1組，場地類型為Ⅱ類(場地特征周期0.35 s)進行設計。結構標準層活荷載為2.0 kN/m2，標準層恒荷載為4.5 kN/m2，屋面恒荷載為7.0 kN/m2，基本雪壓取0.30 kN/m2。

結構平面和立面布置如圖1所示。各框架梁柱截面配筋如圖2所示。其中5層框架的B-1和B-2分別表示邊跨梁和中間跨梁，8層框架和10層框架的B-1和B-2分別表示結構下4層或下5層的邊跨梁和中間跨梁，B-3表示結構上4層或上5層的梁。

圖1 結構平面和立面圖 /mmFig.1 Plan view and elevation view of structures

采用OpenSees進行RC框架結構的有限元建模，動力分析時采用基于蛙跳法的顯式算法。由于dispBeamColumn單元為基于有限單元剛度法的分布塑性單元，無須單獨定義塑性鉸，可以很好地用于框架結構梁柱模擬[26－27]。而且張書豪[13]整合于OpenSees的顯式算法是基于該單元的，因此本研究中RC梁、柱采用dispBeamColumn單元進行模擬，每一桿件采用3個單元，每個單元采用3個高斯積分點。結構分析中考慮P-Δ效應影響。混凝土采用Concrete01本構模型，鋼筋采用Steel02本構模型。材料強度取平均值，材料強度的標準差按文獻[26]和文獻[28]取值?？紤]箍筋約束對混凝土軸向抗壓能力的提升，采用Saatcioglu ＆ Razvi模型[29]進行描述。本文以圖1中間一榀框架為研究對象進行有限元建模，不考慮樓板的影響。

圖2 各框架截面配筋圖 /mmFig.2 Section reinforcement diagrams of each frame

3 顯式算法與隱式算法的對比分析

選取本文第2節(jié)中RC5-1結構作為示例，對比研究顯式算法和隱式算法對結構地震反應分析的影響。其中，隱式算法使用瑞利阻尼，顯式算法使用質量阻尼與前20階振型阻尼進行組合。選取El-Centro(1940)地震動記錄作為地震輸入，分別將地震動的峰值加速度(PGA)調幅至0.1g和3.0g，對RC框架頂點的位移時程及運算耗時進行對比。圖3和圖4分別給出了PGA=0.1g和PGA=3.0g所對應的RC框架頂點(左上角頂點)位移時程。

由圖3和圖4中可以看出，在PGA=0.1g時，顯式算法與隱式算法的結構響應幾乎一致。當PGA=3.0g時，隱式算法由于無法收斂，在6.57 s計算終止，而顯式算法無收斂性問題。換言之，采用顯式算法可以更加接近結構倒塌狀態(tài)。

在配備Inter? Core? i5-7300U CPU @2.60 GHz的計算機上，隱式算法和顯式算法的耗時情況如表1所示。從表1可以看出，當PGA=0.1g時，隱式算法因易收斂，而耗時較小，平均單步耗時僅為0.0150 s。當PGA=3.0g時，隱式算法由于存在收斂性的問題，平均單步耗時為0.1766 s，較PGA=0.1g所對應的平均單步耗時增大了約 12倍。對于顯式算法，無論PGA大小，平均單步耗時保持穩(wěn)定，均為0.008 s左右。當PGA= 3.0g時，顯式算法平均每秒分析耗時約為隱式 算法的2倍，較PGA=0.1g時(約21倍)，顯式算法與隱式算法的總體計算速度的差異有了顯著縮小。

圖3 PGA=0.1 g時RC框架頂點位移時程曲線Fig.3 The time history curve of the top displacement of the RC frame at PGA=0.1 g

圖4 PGA=3.0 g時框架頂點位移時程曲線Fig.4 The time history curve of the top displacement of the RC frame at PGA=3.0 g

表1 隱式算法和顯式算法耗時對比Table 1 A comparison of the time consumption of implicit algorithm and explicit algorithm

4 基于顯式算法的結構抗地震倒塌能力分析

4.1 結構倒塌判據(jù)

1) 判據(jù)1：根據(jù)我國《建筑抗震設計規(guī)范》[14]，RC框架結構在地震作用下，其薄弱層彈塑性層間位移角不應大于1/50，以防止結構倒塌。為此，本文取結構最大層間位移角θmax=1/50作為結構倒塌判據(jù)1。

2) 判據(jù)2：美國FEMA 356[15]規(guī)定混凝土框架結構防止倒塌(S-5，Collapse Prevention，CP)這一性能水準的層間位移角限值為4%。為此，本文取θmax=4%為結構倒塌判據(jù)2。

3) 判據(jù)3：根據(jù)文獻[16]，認為當結構IDA曲線的切線斜率低于初始線彈性階段斜率的20%時，結構將會發(fā)生倒塌，本文將其作為結構倒塌判據(jù)3。

4) 判據(jù)4：本文以結構變形達到不足以維持“安全使用空間”作為倒塌的物理定義[30]。因此，在進行結構地震反應分析中，可取結構構件豎向倒塌位移超過某一數(shù)值作為結構倒塌判據(jù)[4,30－31]。本文以結構主要構件豎向倒塌位移超過1 m作為判據(jù)4。本文所述的4種判據(jù)如圖5所示。

圖5 4種倒塌判據(jù)Fig.5 Four collapse criteria

4.2 IDA分析

選取文獻[32]推薦的22條遠場地震動進行IDA分析，所選地震動記錄的相關信息如表2所示。獲得不同結構的IDA曲線，并在IDA曲線上按上述四種結構倒塌判據(jù)獲得結構倒塌點。僅以RC5-1為例，圖6給出了采用判據(jù)1～判據(jù)3所獲得的結構倒塌點。圖7給出了采用判據(jù)4所獲得的結構倒塌點。圖8進一步給出了結構RC5-1在表2中第4條地震動作用下，當?shù)卣饎覲GA達到判據(jù)4所定義的倒塌強度時，結構θmax和豎向倒塌位移的時程曲線，以及在結構在豎向倒塌位移達1 m時的θmax大小。

表2 遠場地震動記錄Table 2 Far-field ground motion records

圖6 RC5-1框架IDA曲線及前三種倒塌判據(jù)的對應點Fig.6 IDA curves of RC5-1 frame and corresponding points of the first three collapse criteria

由圖8可見，當結構達到判據(jù)4所定義的倒塌狀態(tài)時，結構的θmax響應在短時間內急劇增加，可認為結構已經(jīng)崩潰。急劇增加的θmax響應反映了結構在該PGA的地震動作用下已經(jīng)更為真實地達到了“倒塌”這一極限狀態(tài)。當?shù)卣饎覲GA達到判據(jù)4所定義的倒塌強度時，整個時程的θmax(1015)已不具有實際物理意義，故需要說明的是，為方便展示，圖7中所繪倒塌點的狀態(tài)為達到判據(jù)4所定義倒塌強度前一個IDA調幅步所對應θmax和PGA。

圖7 RC5-1框架IDA曲線及第四種倒塌判據(jù)的對應點Fig.7 IDA curves of RC5-1 frame and corresponding points of the fourth collapse criterion

4.3 地震倒塌易損性分析

基于獲得的結構倒塌點，可進一步采用對數(shù)正態(tài)分布函數(shù)進行擬合，獲得結構的倒塌易損性曲線，如圖9所示。本文獲得6個RC算例框架結構對應不同倒塌判據(jù)的倒塌易損性曲線，如圖10(a)、圖10(b)、圖10(c)所示。從結果可以看出，判據(jù)1明顯低估RC框架結構的抗倒塌能力，進一步說明了如判據(jù)1這類結構設計規(guī)范中提供的供框架結構設計時參考的層間位移角限值只能作為“結構不倒塌”的判據(jù)，不能用來作為“結構倒塌”的判據(jù)[33]。從易損性曲線可以看出，判據(jù)4對RC框架結構的抗倒塌能力的評估要高于判據(jù)3，判據(jù)2對RC框架結構的抗倒塌能力的評估則低于判據(jù)3。

圖8 RC5-1結構最大層間位移角時程及最大豎向倒塌位移時程Fig.8 The time histories of the maximum drift ratio and maximum vertical collapse displacement of RC5-1

圖9 RC5-1 四種倒塌判據(jù)易損性曲線對比Fig.9 Comparison of vulnerability curves of four collapse criteria for RC5-1

圖10 算例結構4種倒塌判據(jù)易損性曲線對比Fig.10 Comparison of vulnerability curves of four collapse criteria for case frames

4.4 地震倒塌模式分析

本節(jié)以RC5-1為例，對不同判據(jù)所對應的結構地震倒塌模式進行了對比。選擇表2中第4條地震動為例，當其分別達到不同判據(jù)所對應的倒塌強度時(見表3)，所對應的結構倒塌模式，如圖11所示。圖11中，采用真實結構比例給出了結構在前3種判據(jù)所對應的倒塌強度下殘余變形模式，以及判據(jù)4對應倒塌強度下結構臨近倒塌狀態(tài)的變形模式。圖11中還分別給出了對應各判據(jù)的頂層側向變形數(shù)值。

由圖11可以看出，判據(jù)1～判據(jù)3所定義的結構倒塌狀態(tài)的殘余變形較小。因此，判據(jù)1～判據(jù)3顯然不能真實地反映結構的倒塌狀態(tài)。從判據(jù)4所定義的結構倒塌狀態(tài)的結構變形來看，判據(jù)4對倒塌的定義顯得更加真實，相較于其他3種判據(jù)，判據(jù)4更好地表現(xiàn)出了結構的倒塌極限狀態(tài)。

針對判據(jù)3和判據(jù)4，進一步繪制了結構頂層側向變形時程，以及不同時刻對應的結構變形模式，如圖12所示。

圖11 4種判據(jù)致倒塌PGA下的位移模式Fig.11 The displacement states under the PGA that causes the collapse corresponding to the four collapse criteria

圖12 致倒塌PGA下判據(jù)3與判據(jù)4對應頂層側向變形時程及不同時刻結構變形模式Fig.12 The time histories of the lateral displacement of the top floor and the displacement states at different times under the PGA that causes the collapse corresponding to collapse criteria 3 and 4

由圖12可見，判據(jù)3致倒塌PGA作用下結構位移時程響應在整個時程中間一時刻達到最大，但是隨后又在波動中有減少趨勢。而判據(jù)4致倒塌PGA作用下結構位移時程響應在整個時程中呈發(fā)散趨勢，直至結構崩潰倒塌。

4.5 地震倒塌裕度比分析

地震倒塌裕度比(Collapse Margin Ratio, CMR)也稱作結構抗倒塌儲備系數(shù)，是文獻[32]推薦使用的結構抗地震倒塌能力評估指標。本文以PGA為地震動強度指標，罕遇地震對應的CMR為：

式中：PGA50%倒塌為50%倒塌概率對應的地震動PGA；PGA罕遇地震為參考文獻[14]取值的地震動強度。

基于算例框架的地震倒塌易損性，本文進一步計算了6個算例框架對應4種倒塌判據(jù)的倒塌裕度比，如表4所示，將不同結構對應不同倒塌判據(jù)的倒塌裕度比進行對比，如圖13所示。僅以RC5-1為例，相較于判據(jù)3而言，判據(jù)1對結構的抗倒塌能力低估了約56%，判據(jù)2對結構的抗倒塌能力低估了約13%，而判據(jù)4對結構的抗倒塌能力高估了約47%。綜合6個算例結構的倒塌裕度比結果來看，判據(jù)3對應的CMR結果為判據(jù)1對應的CMR結果的2.5倍左右，判據(jù)2對應的CMR結果為判據(jù)1對應CMR結果的2倍左右。這一結果進一步說明了判據(jù)1僅適合用作結構設計時的倒塌控制，不適合作為判斷結構倒塌的判據(jù)[33]。換言之，當結構達到判據(jù)1所定義的倒塌狀態(tài)時，結構仍具有很大的抗倒塌能力冗余度。

判據(jù)2相對于判據(jù)3而言，其所對應的結構抗倒塌能力低估了13%～34%，且低估程度隨著結構層數(shù)增加而加深?？梢妼τ赗C框架結構，判據(jù)2對應的倒塌狀態(tài)仍有一定的抗倒塌能力冗余度，換言之，F(xiàn)EMA 356[15]規(guī)定的混凝土框架結構CP性能水準對應的4%層間位移角限值對于國內RC框架結構而言，作為倒塌判據(jù)略顯保守。

判據(jù)4相較于判據(jù)3而言，其所對應的結構的抗倒塌能力結果提高了約20%～70%。這一結果說明判據(jù)4所定義的結構倒塌狀態(tài)比其他判據(jù)定義的倒塌狀態(tài)更為危險。由前文可知，判據(jù)4是直接基于倒塌物理定義的，結構在達到判據(jù)4所定義的倒塌狀態(tài)時，迅速完全崩潰并真正達到了完全倒塌。因此，判據(jù)4定義的倒塌狀態(tài)更接近結構的真實倒塌狀態(tài)，對于RC框架結構，判據(jù)4定義的結構抗地震倒塌能力可視為結構抗倒塌能力的上限，也可視為結構的極限抗倒塌承載能力判據(jù)。

表4 RC框架結構4種倒塌判據(jù)的倒塌裕度比Table 4 CMR of RC frames corresponding to four collapse criteria

圖13 4種倒塌判據(jù)下RC框架結構CMR比較Fig.13 Comparison of CMR of RC frame structures under four collapse criteria

5 結論

本文針對不同層數(shù)、不同設防烈度的6個RC框架結構，基于改進蛙跳法的顯式算法，進行了RC框架結構的抗地震倒塌能力分析，獲得了結構的倒塌易損性曲線和倒塌裕度比，得到如下結論：

(1) 顯式算法與隱式算法在結構非線性發(fā)展程度不強時，分析結果基本完全吻合。與隱式算法相比，顯式算法非線性分析的收斂性更好，也更易于模擬結構在倒塌狀態(tài)附近的非線性行為。

(2) 判據(jù)1和判據(jù)2對結構抗地震倒塌能力的評估結果較為保守，用其來評估結構抗地震倒塌能力，結構仍有較大抗倒塌能力冗余度。與判據(jù)3對應的結果相比，判據(jù)4對應的結構倒塌裕度比要高20%～70%。

(3) 判據(jù)4基于結構倒塌的物理意義，在判據(jù)4所對應的結構倒塌狀態(tài)，結構在地震動作用下變形急劇增大乃至崩潰達到完全倒塌。相較于前3種倒塌判據(jù)，判據(jù)4能更好地體現(xiàn)結構倒塌狀態(tài)的真實非線性行為。由于顯式算法較強的非線性分析能力，因此其可以不受制于收斂性問題，而可以方便利用判據(jù)4來進行倒塌能力的評估，有利于結構真實抗倒塌能力的評估。