基于擬力法的增量動力分析方法

郝潤霞 ， 楊作續， 李 鋼， 余丁浩， 賈 碩

(1. 內蒙古科技大學 土木工程學院，內蒙古 包頭 014010；2. 大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024)

增量動力分析(Incremental Dynamic Analysis,IDA)方法[1-2]是基于性態地震工程(Performance-based Earthquake Engineering,PBEE)中確定結構在不同水準地震動響應的參數化分析方法，可以準確的觀察結構從彈性、彈塑性直至動力失穩倒塌的全過程中結構性態的發展變化規律[3]，同時可以克服靜力推覆分析(Push-Over Analysis,POA)方法特別是在結構接近倒塌時因近似靜力分析而存在的問題[4]。然而其分析計算時運算量大、耗時長等計算效率問題，卻嚴重地影響了其在工程應用中的發展。

近年來，國內外學者圍繞IDA的計算效率問題開展了一系列的研究。在簡化分析方法方面，如基于靜力推覆分析的簡化近似IDA(SPO2IDA)方法[5]，該方法主要適用于基本振型控制的結構體系；基于模態推覆分析(Modal Pushover Analysis ,MPA)的簡化快速模態增量分析(Modal Incremental Dynamic Analysis, MIDA)方法[6-7]，該方法考慮了模態選取、結構延性和滯回模型對其計算精度的影響，當考慮的高階振型貢獻越多時[8]，其計算精度越高，且保證精度的前提下，能夠有效提高計算效率，但對于不規則結構體系適用性較差[9]；基于等效二自由度模型(Equivalent Dual Degree of Freedom ,EDDOF)的簡化IDA方法[10]，該方法以單向偏心結構的等效二自由度模型應用于IDA，對單向偏心結構進行抗震性能評估，在保證較高的精度前提下，計算效率得到了大幅提高。對比分析SPO2IDA方法與MIDA方法，SPO2IDA方法計算效率更為突出，但其計算精度較差[11]。在優化非線性求解方面，以快速非線性時程分析(Fast Nonlinear Analysis,FNA)方法[12]作為計算工具提出的快速增量動力分析(Fast Incremental Dynamic Analysis ,FIDA)方法[13]，與傳統IDA方法、MIDA方法相比，FIDA方法在計算效率與精度上優勢明顯。此外，如基于Nataf變換點估計法與單地震動記錄IDA相結合的單地震動記錄隨機IDA方法[14]，該方法可以使IDA分析的計算精度和計算效率達到平衡。以IDA曲線混合統計方法進行統計分析，并直接以結構離散IDA曲線進行地震易損性曲線及地震失效概率計算，可有效避免多余的統計與擬合[15]。IDA方法在結構非線性動力時程分析時常采用傳統變剛度方法，該方法在任一時刻非線性迭代步求解時，結構整體剛度矩陣需要實時合成與分解，計算量大，效率低。

擬力法(Force Analogy Method,FAM)[16-17]作為一種高效的結構非線性分析新方法，自提出以來，被廣泛應用到結構振動控制[18]、消能減震[19]、抗震性能評估[20]等領域。之后，李鋼等[21]將其形變分解思想擴展到材料層面，與有限元相結合，提出了基于FAM的精細化有限元非線性分析方法。在求解過程中結構整體剛度矩陣始終保持不變，非線性信息通過局部塑性矩陣予以體現，迭代計算時避免了傳統變剛度方法中結構整體剛度矩陣的實時合成與分解，且局部塑性矩陣相比傳統變剛度法的結構整體剛度矩陣規模小的多，可有效改善傳統變剛度法的計算效率問題[22]。

本文提出了基于擬力法的IDA計算方法，建立了基于FAM的IDA運動方程，考慮運用FAM非線性求解時整體剛度矩陣始終保持不變的特性以及迭代求解平衡方程時步步更新位移增量的特點，實現特定剛度矩陣存儲與調用，并給出本文方法與傳統IDA方法的時間復雜度函數，通過數值算例驗證本文方法的準確性和高效性。

1 基本理論

1.1 增量動力分析方法

IDA法是一種基于非線性動力時程分析的參數化分析法，其實質是將單一強度的動力時程分析擴展成為多重強度的動力時程分析，以獲取結構在不同強度地震動作用下的動力響應。

圖1為IDA方法的基本流程圖，即按照一定的調幅法則對待施加的一條(或多條)地震動記錄的地震動強度(Intensity Measure,IM)進行調幅，使之成為涵蓋結構從彈性、彈塑性直至動力失穩倒塌所需的多重強度；然后對給定的結構模型在多重強度地震動作用下分別進行非線性動力時程分析，得到每一強度下的地震最大響應結果即結構的工程需求參數(Engineering Demand Parameters,EDP)；將每一條地震動記錄下一系列的離散點(EDP,IM)點繪制于二維坐標系中，并選用合適的插值方法擬合EDP-IM關系曲線，生成一條(或多條)IDA曲線；通過在IDA曲線定義結構不同性態水準的極限狀態如：正常使用(Normal Occupancy,NO)極限狀態，立即使用(Immediate Occupancy,IO)極限狀態，生命安全(Life Safe,LS)極限狀態，防止倒塌(Collapse Prevention,CP)極限狀態，并對其趨勢及離散狀況進行統計分析，可以較為準確的觀察結構在彈塑性響應歷程中強度、剛度及變形的發展變化規律，能夠反映出結構在不同水準地震動作用下地震需求能力和整體抗倒塌能力。

圖1 IDA方法流程圖Fig.1 Flow chart of IDA method

1.2 擬力法

FAM作為結構非線性分析新方法，其核心思想是在結構非線性分析過程中引入局部塑形機制，將結構的變形進行彈塑性分解。通過在結構中設置塑性插值點來獲取塑性應變分布形式，使結構的非線性狀態可由塑性插值點處塑性變形予以體現；結構控制方程則由虛功原理及塑性插值點處的內力平衡條件給出。同時在給出結構控制方程時僅考慮塑性插值點處進入非線性狀態時塑性插值點(對于處于彈性狀態的塑性插值點處的塑性自由度為零，可進行消元處理)。對于含有n個結點位移自由度、d個塑性自由度(塑性插值點處產生的塑性自由度數)的結構體系，其非線性求解時每一迭代步中，FAM的結構控制方程可表示為

(1)

結合Newton-Raphson迭代方法，對于式(1)進一步整理，得到每一非線性迭代步的求解方程可表示為

(2)

其中，

(3)

式中：Kf為局部塑性剛度矩陣，反映結構中進入塑性狀態的局部非線性信息，其規模與結構中參與計算的塑性自由度數量有關。

對于傳統變剛度法每一非線性迭代步的求解方程可表示為

(4)

對比式(2)和式(4)可知，式(2)中系數矩陣項實為式(4)中結構切線剛度矩陣的逆矩陣的等效表達。由此可知，FAM通過引入Kf，不僅避免了傳統變剛度法對切線剛度矩陣實時合成與分解時計算效率低的問題，而且保留了Newton-Raphson迭代收斂速度快的特性。式(2)即為Woodbury求逆公式，求解過程中僅需對Kf進行實時的合成與分解即可。

2 基于擬力法的增量動力分析

2.1 運動方程

結構在地震動作用下的運動方程可表示為

(5)

對于式(1)進一步展開可寫成

(6)

(7)

對于式(6)以全量形式可進一步表示為

(8)

將式(8)代入式(5)，并考慮引入地震動強度調幅系數λi(i=1,2,...,n)，用以表示不同的調幅強度，即可得到多重強度下的基于FAM的IDA運動方程為

(9)

2.2 方程求解

這里以Newmark-β方法[23]為例，以增量形式對式(9)運動方程進行逐步積分求解，得到λi調幅強度下時間步長為Δtj(j=1,2,...,n表示不同的時間步長)的位移增量向量的等效平衡方程可表示為

(10)

其中，

(11)

(12)

(13)

式(13)為Woodbury求逆公式，對于式(13)可進一步表示為

ΔX=ΔX1+ΔX2

(14)

其中，

(15)

(16)

(17)

2.3 剛度矩陣存儲與調用

本文方法的特定剛度矩陣存儲與調用流程示意如圖2所示。

圖2 本文方法的特定剛度矩陣存儲與調用流程圖Fig.2 Flow chart of storage and call of a particular stiffness matrix of this method

3 計算效率評價

算法時間復雜度分析[24]作為一種行之有效的評價算法效率的分析方法，其將運行的算法從計算機中抽象出來，使其量化結果不受計算條件的影響，依據算法時間復雜度與算法計算效率之間的關系(算法時間復雜度越高，算法計算效率越低)，可直接用于定量的評價算法的計算效率。文獻[25]中采用算法時間復雜度理論對比分析了FAM與傳統變剛度方法的計算效率。基于文獻[25]已有理論，本文將算法時間復雜度理論應用于IDA分析的效率評價中，考慮特定矩陣的調用特性，忽略矩陣調用等次要運算的時間復雜度，給出本文方法與傳統IDA方法的算法時間復雜度函數，從數學理論的角度用以揭示本文方法較傳統IDA方法在IDA分析時所具有的計算高效性。本節從IDA的實質出發，考慮非線性動力時程分析中，結點位移的求解最為耗時，且本文方法與傳統IDA方法主要區別也在于此，故本節僅對結構結點位移求解部分進行算法時間復雜度分析討論。

本節以單一地震動記錄的IDA分析為例，分別對本文方法與傳統IDA方法的算法時間復雜度進行分析比較。對于本文1.2節給出的結點位移自由度數為n，塑性自由度數為d的結構體系，考慮因算法的迭代次數對復雜度分析的影響，假定兩種方法在計算結點位移時均采用Newton-Raphson迭代法，即保證相同的收斂速度，記IDA分析每一強度下的非線性動力時程分析所需迭代次數為Ni(i對應于IDA分析時調幅系數的下標，用以表征不同調幅的強度)。

3.1 算法時間復雜度

(1)本文方法

考慮本文方法在求解結點位移時，各系數矩陣項的特殊性質(如稀疏性、帶狀性、對稱性、可存儲性等)以及計算次序對算法時間復雜度的影響，忽略矩陣調用等次要運算的復雜度，以最合適的計算方式進行算法時間復雜度分析。則本文方法的一次非線性動力時程分析的時間復雜度，如表1所示。

表1 本文方法的時間復雜度

本文方法在單一地震動記錄下總時間復雜度可表示為

(18)

(2)傳統IDA方法

傳統IDA方法在單一地震動記錄下總時間復雜度可表示為

(19)

3.2 效率對比

(20)

從上述分析可知，di值的大小是影響兩種方法的時間復雜度差異的重要參數。取ηmax=dimax/n為塑性自由度數理論最大值dimax與結點位移自由度數比值，以結點位移自由度n=500，n=1 000，n=2 000為例，分別給出如圖3所示兩種方法的時間復雜度隨塑性自由度比值ηmax的變化的關系曲線。

從圖3中可以看出，本文方法的時間復雜度隨ηmax的增大而增大；時間復雜度增長速率隨著結點位移自由度數的增大而逐漸增強；達到dimax時的塑性自由度比例隨著結點位移自由度數的增大而減小。對于不同結點位移自由度的結構來說，當ηmax較小時，本文方法的效率優勢明顯；隨著ηmax的逐漸增大，優勢逐漸減弱；當d≥dimax，本文方法的效率將劣于傳統方法。

圖3 塑性自由度比值與時間復雜度關系曲線Fig.3 The relationship between the ratio of plastic degree of freedom and time complexity

事實上，隨著λi不斷增大，結構實現了從線彈性、彈塑性直至動力失穩倒塌的全過程，在這一過程中結構往往只有少部分截面或區域(薄弱部位)產生了非線性變形即此部位塑性插值點處的塑性自由度被激活，而這部位相對于整個結構所占比例較小，如框架結構的非線性變形一般集中于梁或柱的端部，構件的中間大部分區域并不會進入非線性狀態。因此，通常情況下，結構出現di≥dimax情況的概率較小，不會影響整個IDA分析過程的計算效率評價。

4 數值算例

4.1 工程概況

本文所選結構模型為8層鋼筋混凝土框架結構，抗震設防烈度Ⅷ度(0.30g),設計地震分組為第一組，場地類別為Ⅱ類，結構阻尼比為5%。結構平面布置如圖4所示，樓面恒荷載標準值為5.0 kN/m2， 活荷載標準值為2.0 kN/m2，屋面恒荷載為6.0 kN/m2，屋面活荷載為0.5 kN/m2。樓層層重按1.0×恒載+0.5×活載折算。首層層高3.9 m，其余各層層高3.3 m，結構總高度27 m。結構梁柱尺寸及配筋見表2。混凝土強度等級，柱、梁、樓板均為C30；梁柱縱向鋼筋為HRB400，箍筋為HPB300。由于結構平面對稱，選取一榀橫向框架(3軸線)進行建模與分析，結構立面模型如圖5所示。

圖4 結構平面布置簡圖(單位：mm)Fig.4 Sketch of structural plane layout(unit:mm)

圖5 結構立面布置簡圖(單位：mm)Fig.5 Sketch of structural vertical layout(unit:mm)

樓層梁截面/mm配筋/mm2中柱截面/mm配筋/mm2邊柱截面/mm配筋/mm21300×6004 688(BL)6 541(ZL)700×7008 758700×70011 7842300×6004 688(BL)6 541(ZL)700×7005 024700×7005 0243300×6004 688(BL)6 541(ZL)700×7005 024700×7005 0244300×6004 688(BL)6 541(ZL)700×7005 024700×7005 0245300×6004 112600×6003 786600×6003 7866300×6004 112600×6003 786600×6003 7867300×6004 112600×6003 786600×6003 7868300×6002 518600×6003 786600×6003 786注:BL為邊梁;ZL為中梁

4.2 有限元模型

分別采用有限元分析軟件ABAQUS和FAM程序建立分析模型進行結構的動力非線性分析，采用迭代法，不考慮P-delta效應，兩種方法的梁柱均選用纖維梁單元予以模擬，單元類型為B32。其中每個構件劃分3個單元，選取的一榀橫向框架結構共計168個單元，864個結點位移自由度，1 008個塑性自由度，梁柱截面纖維劃分均為10×10，每根鋼筋相當于1根纖維；鋼筋選用彈塑性隨動硬化單軸本構模型；混凝土本構的骨架曲線選用Kent-Park[28]模型，采用Yassin[29]提出的滯回規則，其中ABAQUS中材料本構模型調用清華大學開發的子程序PQ-Fiber進行定義。結合兩種方法對此框架結構進行單一地震下的IDA分析。

4.3 IDA分析及對比

為便于量化地震動強度等級，本文選用PGA作為地震動強度指標，地震動如圖5所示方向輸入，考慮場地類型、震級、峰值加速度等因素，以FEMA P695[30]推薦的22組遠場地震動記錄(共44條水平分量)中隨機選取一條地震動記錄Northridge波(NORTHR-MUL279)并歸一化處理后作為IDA分析的地震動輸入，持時取20 s。根據Vamvatsikos等的建議，依據hunt & fill調幅法則進行PGA調幅，給定PGA初值為0.035g，調幅步長取0.2g，步長增量取0.05g。考慮到層間位移角與樓層層間的變形能力以及節點的轉動能力有關，能夠反映結構的層間位移延性和整體位移延性，選取最大層間位移角θmax作為工程需求參數以及倒塌點的判定指標。結合文獻[31]給出鋼筋混凝土框架結構的不同性態水準下的最大層間位移角限值，如表3所示，當給定最小Δtj下數值算法仍不收斂或結構的最大層間位移角θmax>4%時，可認為結構倒塌并結束計算。

表3 RC框架結構不同性態水準下層間位移角限值

通過IDA分析，得到該結構的工程需求參數最大層間位移角θmax與地震動強度指標PGA的關系曲線，即IDA曲線，如圖6所示。

由圖6可以看出，對于結構的變形而言，本文方法與傳統方法的計算結果基本一致，誤差0.8%～5.5%，兩者的IDA曲線基本吻合；在PGA處于0～0.2g內時，此時調幅強度較小，結構處于線彈性狀態，結構的動力響應與地震動強度基本成線性關系，計算結果誤差在1%左右，兩種方法的IDA曲線基本一致；隨著調幅強度的逐漸增大，結構達到屈服，非線性狀態逐步發展，在PGA處于0.2～0.6g內時，計算結果誤差在3.5%左右；當繼續調幅地震動強度(PGA≈0.625g后)，曲線開始表現出了較小的差異，計算結果誤差在5.5%左右，究其原因是兩種方法在計算過程中引入的計算累積誤差的差異性以及結構動力響應與地震動強度之間顯著的非線性關系等造成的；在PGA處于0.8～0.9g內時，結構的性能變化特征表現出強化現象，使得曲線的斜率出現上升；當PGA>0.9g后，結構在微小的地震動強度激勵下即可引起較大的地震動力響應，此時結構的最大層間位移角θmax趨于無窮大，結構體系已達到整體動力失穩狀態。結合IDA方法的最終目的是通過IDA曲線簇進行統計分析，從概率的意義上用以評估結構的抗震性能。綜上所述，本文方法在計算精度方面是可予以保證的，對結構的IDA分析是合理準確的。

圖6 結構的IDA曲線Fig.6 IDA curve of the structure

4.4 計算效率分析

通過4.3節IDA分析，可獲取該結構在每一強度下每一非線性迭代步產生的塑性自由度數的曲線。由于篇幅有限，以調幅強度為0.635g為例，現給出該強度下每一時間增量步的塑性自由度數曲線，如圖7所示。

圖7 塑性自由度數曲線Fig.7 Number of plastic degrees of freedom curve

由圖8可以看出，在整個時程分析的過程中，傳統方法的時間復雜度均高于本文方法；本文方法產生的時間復雜度的最大值約為1.94×107，但僅僅發生在最大迭代次數所在的時間增量步內，其每一時間增量步的時間復雜度的平均值約為2.16×106，將每一時間增量步下的時間復雜度匯總可得到本文方法的總時間復雜度約為4.32×109；對于傳統IDA方法產生的時間復雜度的最大值約為6.65×107，其每一時間增量步的時間復雜度的平均值約為9.29×106。兩種方法時間復雜度對比可知，本文方法的時間復雜度的平均值(約為2.16×106)僅為傳統IDA方法時間復雜度的平均值(約為9.29×106)的23.25%，即該強度下本文方法的總時間復雜度僅為傳統IDA方法的總時間復雜度的23.25%。

圖8 時間復雜度曲線Fig.8 Time complexity curve

現給出地震動記錄Northridge波下的IDA的不同調幅強度的時間復雜度以及復雜度比，如表4所示。

表4 本文方法與傳統IDA方法的時間復雜度對比

由表4中可知，兩種方法的時間復雜度隨著調幅強度的不斷增大，呈現增大的趨勢。當調幅強度較小時，如0.035g，0.135g，0.235g，0.36g，此時結構處于線彈性狀態或僅有極少部分進入塑性狀態，結構中被激活的塑性自由度數遠小于結構的結點位移自由度數，此時本文方法的時間復雜度相對于傳統方法的優勢較為明顯，僅為傳統方法時間復雜度的10%以下。隨著調幅強度的不斷增大，結構進入塑性狀態的程度加大，此時塑性自由度數不斷增加，使得本文方法的時間復雜度相對于傳統方法的優勢逐漸縮減。當調幅強度為0.035g時，本文方法的時間復雜度量級僅為108，復雜度比約為6%；而當調幅系數為0.916 25g時，此時結構處于倒塌極限狀態，時間復雜度量級為1011，復雜度比約為56%。兩強度之間的時間復雜度量級相差近103之多，復雜度比相差近50%。將各個調幅強度下的時間復雜度累積求和，本文方法的總的時間復雜度(1.597 8×1011)僅約為傳統IDA方法(3.854 8×1011)的41.45%。由上述分析可知，本文方法較傳統IDA方法具有較明顯的計算效率優勢且隨著調幅強度的逐漸增大而減少。

5 結 論

本文針對傳統IDA方法在結構非線性動力時程分析時計算量大、耗時長的問題，從優化非線性求解方面出發，以FAM作為結構非線性分析工具，提出了一種基于FAM的IDA方法，并通過8層鋼筋混凝土框架結構IDA分析，得到如下結論：

(1) 通過引入附加虛擬荷載項，使得結構非線性求解過程中整體剛度矩陣始終保持等效彈性剛度矩陣不變(時間步長一定時)，避免了傳統IDA方法等效切線剛度矩陣實時分解與更新，從而提高了計算效率。

(2) 通過建立特定剛度矩陣存儲與調用機制，實現了IDA分析全過程中特定矩陣的重復調用回代求解，以減少計算工作量，進而提高了計算效率。

(3) 計算精度上，本文方法與傳統IDA方法的計算結果較為吻合，驗證了本文方法的準確性。

(4) 計算效率上，對非線性動力時程分析最為耗時的結點位移求解部分，采用時間復雜度評價方法進行了定量地分析，針對本文中算例，基于擬力法的IDA方法的時間復雜度約為傳統IDA方法的41.45%，驗證了本文方法較傳統IDA方法具有明顯的計算效率優勢。