考慮參數不確定性的無筋砌體結構地震易損性分析

蔣亦龐，蘇 亮，黃 鑫

(浙江大學空間結構研究中心，浙江，杭州 310058)

在歷次地震中，砌體結構的大量破壞與倒塌造成了非常嚴重的人員傷亡和財產損失[1-2]，而無筋砌體結構由于未設置合理的抗震構造措施，其地震破壞則尤為嚴重。由于無筋砌體結構仍是我國廣泛存在的建筑結構形式之一[3]，對其進行準確的地震易損性分析具有非常重要的工程實際意義。

結構地震易損性分析需要考慮不確定性的影響。以往的研究認為地震易損性分析的不確定性主要源于地震動作用，其影響大小則通常是通過對多條地震動輸入的增量動力分析(incremental dynamic analysis, IDA)進行計算考察[4-5]。然而，隨著近年來不確定性研究的不斷深入，越來越多的研究成果表明[6-10]，當結構參數不確定性較大或結構響應非線性程度較高時，結構參數的不確定性將會帶來無法忽視的影響。

目前對結構參數不確定性的研究主要集中于鋼筋混凝土結構、鋼結構以及橋梁結構等[6-10]，國內外對砌體結構參數不確定性影響的研究則尚處于起步階段。Parisi等[11]利用靜力推覆法研究了8個材料參數不確定性對一2層無筋砌體結構的抗震能力影響，分析結果顯示結構變形能力產生了較大的離散性，表明材料參數不確定性對砌體結構抗震能力的影響較大，但文章沒有考慮地震動不確定性的影響。Lu等[12]基于IDA分析，采用簡化的多自由度層剪切模型，利用蒙特卡羅法和一次二階矩(first-order second-moment, FOSM)方法對加筋砌體結構的區域震害預測進行了結構參數敏感性研究。研究表明，FOSM方法能夠在較少的計算量下，得到與蒙特卡羅法相似精度的結果；結構參數的不確定性對單個結構的地震易損性較大，而在區域震害分析時則可以被忽略。當然由于文章僅采用了簡化的多自由度層剪切模型，其精確程度值得進一步計算驗證。

鑒于無筋砌體結構參數的高離散性及其地震響應的強非線性[11]，在其地震易損性分析中考慮結構參數的不確定性顯得尤為必要。正是基于這樣的工程背景，本文以4幢不同層數的無筋砌體結構為研究對象，在OpenSees中建立其等效框架有限元模型，采用IDA法和FOSM法考察了地震動以及包括荷載、材料和阻尼在內的9個參數不確定性對結構地震易損性的影響，并給出了相應的工程建議。

1 無筋砌體結構模型

1.1 無筋砌體結構參數

采用層數為3層~6層的4幢典型無筋砌體結構為研究對象[13]，其中3層結構模型示意圖如圖1(a)所示。除層數不同之外，其余3幢結構的參數均與3層結構相同。不考慮結構平面外破壞及扭轉效應，在此僅取一榀平面無筋砌體墻作簡化分析。單層單跨墻的尺寸為層高3 m、跨長3.6 m，平面外結構跨度3.9 m。單層單跨墻的門窗洞口尺寸為寬1.2 m、高1.2 m，位于墻體中心。墻厚370 mm，墻材料取MU20的燒結普通磚和M10的砂漿。樓板假定為鋼筋混凝土材料，厚120 mm，計算時僅考慮其荷載作用。

1.2 結構建模方法

結構模型采用等效框架法建模，在OpenSees上進行有限元分析。等效框架法是目前被廣泛采用的砌體結構建模方法，用于砌體結構的整體抗震分析，兼具高計算效率和高計算精度的優點[14-15]。等效框架法的基本思想是將墻肢和窗群梁等效成框架柱和框架梁單元，墻肢與窗群梁相交的部分由于在地震中一般不會破壞則等效成剛域，如圖1(b)所示。

圖1 3層無筋砌體結構及其等效框架模型Fig.1 3-storey unreinforced masonry structure and its equivalent frame model

等效框架模型在OpenSees中實現如下[14]：

1) 單元：墻肢和窗群梁選用基于柔度法的梁柱單元建模，剛域處采用足夠大的彈性模量本構來模擬；

2) 等效高度：墻肢和窗群梁的等效高度選用Dolce[16]提出的準則，即外側墻肢的等效高度為相鄰洞口兩頂點延伸出來的30°傾斜線所包圍的高度，而內側墻肢和所有窗群梁的等效高度與相鄰洞口尺寸相同；

3) 截面：墻肢和窗群梁截面采用纖維模型，各個截面的單個纖維尺寸約為40 mm×40 mm，沿截面長、寬方向分別等分；

4) 單元的軸向與彎曲行為：截面纖維的單軸本構模型采用圖2(a)所示的Kent-Park模型[17]，以模擬單元的軸向及彎曲行為。圖中極限抗壓強度fu和極限壓應變εu取值分別為0.88倍抗壓強度fm及1.6倍抗壓強度處壓應變εm[18]；

5) 單元剪切強度：單元剪切強度Vu取代表對角破壞的剪切強度Vu1和代表剪切滑移破壞的剪切強度Vu2的較小值。Vu1參考Turn?ek等[19]提出的計算公式：

式中：ft為砌體材料抗拉強度；σ0為截面平均壓應力(σ0=N/lt)，l和t分別為單元截面的寬度和厚度，系數b由單元的長寬比確定，如下式：

式中，h為單元長度。Vu2的計算公式如下[20]：

式中，fv0為砌體抗剪強度。

圖2 模型所使用的兩種本構模型Fig.2 Two constitutive laws used in model

6) 單元剪切行為：選用如圖2(b)所示的Pinching4單軸本構模型[21]來描述單元的剪切行為，并以截面組合的形式加到纖維截面上，組合后單元響應為彎曲與剪切本構的共同作用結果。模型骨架曲線可由點1~4的坐標位置(γi，Vi)確定，骨架曲線正向與反向對稱，其中γi與Vi分別代表第i點對應的位移與荷載，其計算詳見式(4)；模型卸載與再加載曲線則由a、b、c、d、e與f的坐標位置確定，其中b點坐標(γm,Vm)代表歷史加載的最大位移γm及對應的荷載Vm；a點坐標為(0.5γm, 0.25Vm)；c點橫坐標由卸載剛度確定，其縱坐標為-0.05V3。d、e與f點分別與a、b、c點關于原點對稱。其余控制滯回反應中逐漸增加的卸載剛度退化、再加載剛度退化以及強度退化的16個參數均使用Lowes等[22]的建議值。

式中：G為砌體材料剪切模量；AS為單元截面積。

1.3 結構不確定性參數

選取9個無筋砌體結構不確定性參數，分別為：砌體材料密度ρm、鋼筋混凝土板材料重度γ、活荷載q、砌體抗壓強度fm、砌體抗拉強度ft、砌體抗剪強度fv0、砌體彈性模量E、砌體剪切模量G、阻尼比ξ。由于目前針對砌體結構不確定性參數之間的相關性研究仍不充分，現有文獻不足以獲得各參數間的相關系數，因此本文作簡化分析，假設各參數之間相互獨立。同時，本文也不考慮模型參數在結構空間分布上的不確定性。各參數詳細統計特征見表1。

表1 結構不確定性參數詳細統計特征Table 1 Detailed statistical characteristics of structural uncertainty parameters

2 地震動參數及結構極限狀態

2.1 地震動的選擇

既有研究表明，在IDA分析中選取20條地震動記錄足以考慮地震動的不確定性影響[4,9]。考慮到IDA分析比較耗時，本文從文獻[24]建議的100條遠場地震動記錄中選取20條作為地震動輸入考慮地震動的不確定性。所選地震動均來自美國太平洋地震研究中心PEER的強震數據庫[26]，如表2所示。所選地震波的震級分布在5.65級~7.6級，峰值加速度(peak ground acceleration, PGA)分布在0.05g~0.81g。同時，來自同一場地震事件的地震動記錄不超過兩條。

2.2 地震動強度指標的選取

地震易損性分析中常用的地震動強度參數(intensity measure, IM)包括地面峰值加速度PGA和結構基本周期對應的加速度譜值Sa(T1, 5%)，T1為結構基本周期。由于結構層數和結構不確定性參數的改變都會使得結構的基本周期成為一個不確定值，如采用Sa(T1, 5%)作為IM參數會使分析變得更為復雜，故選取PGA作為IM參數。在進行IDA分析時，分別將20條地震動記錄的PGA分別調幅為0.05g~0.8g(間隔0.05g)，對結構進行彈塑性時程分析。

表2 20條天然地震動記錄Table 2 20 natural ground motion records

2.3 結構極限狀態的定義

選取最大層間位移角為結構抗震性能指標，將無筋砌體結構的極限狀態(limit state, LS)分為三種，即輕微破壞(LS1)、中等破壞(LS2)、嚴重破壞(LS3)，各極限狀態的層間位移角限值如表3所示[27]。

表3 層間位移角限值Table 3 Allowance for story drift angle

3 考慮結構不確定性的地震易損性分析

本文采用FOSM方法考慮結構的不確定性，首先對FOSM方法作簡單介紹。假設隨機變量y是隨機向量的函數y=f(X)，FOSM方法將該函數在隨機向量X的均值μX處近似展開為一階泰勒級數形式，從而估計隨機變量y的均值μy和方差分別如式(5)及式(6)所示：

式中：ρij為xi和xj的相關系數；σxi為xi的標準差。由下式計算：

上述公式詳細的推導過程可見文獻[12]及其他相關文獻。此外，值得一提的是，上述方法實際上是平均值一次二階矩法(mean value first-order second-moment, MVFOSM)，屬于FOSM方法的一種。然而在很多文獻中通常將該方法直接稱為FOSM方法[12,28]，為簡化起見，本文也直接將其稱為FOSM方法。

在本文中，隨機向量X由9個結構不確定性參數確定，如下式所示：

結構的地震易損性[24]定義為給定地震強度水平時結構達到或超過某種極限狀態的條件概率：

式中：D≥C表示結構達到或超過某種極限狀態，其中D為地震需求，C為抗震能力；FR(x)為稱為地震易損性函數，一般符合對數正態分布。采用基于地震動強度(IM-based)的地震易損性函數來展開研究，針對某一特定極限狀態LSi，其表達式為：

式中：Φ[·]是標準正態概率分布函數；Ci為極限狀態LSi時結構抗震能力限值；mR和βR分別為結構地震易損性函數的中位值和對數標準差。

針對某一特定極限狀態LSi的結構地震易損性曲線，可由下述步驟獲得[29]：

1) 選擇NRTR條地震動記錄對結構進行IDA分析；

2) 統計當地震動強度IM=im時，NRTR條地震動記錄中導致結構超過該極限狀態的地震動數量Ni，采用基于頻率的統計方法獲得結構超過該極限狀態的失效概率為：

3) 將獲得的失效概率用式(10)擬合，可獲得結構對應該極限狀態的易損性曲線，如圖3所示。

將隨機變量y定義為結構達到某一極限狀態的平均地震動強度mR的對數，即：

圖3 地震易損性曲線的擬合Fig.3 Fitting of a fragility curve

由文獻[12, 30]可知，當只考慮地震動的不確定性時，式(10)的對數標準差為βR=βRTR，表征地震動不確定性的影響；當同時考慮地震動和結構參數的不確定性時，由結構參數不確定性貢獻的對數標準差βMDL可由FOSM方法通過式(6)獲得，此時定義新易損性函數的總對數標準差βTOT可由式(13)獲得[30]：

由上述原理可知，當進行考慮結構參數不確定性的地震易損性分析時，若結構數量為n個，結構不確定性參數為m個，則需要進行的IDA分析次數為：NIDA=n×(2×m+1)。

4 結果分析

以表2中的20條地震動記錄作為輸入，按照第3節所介紹的計算方法，對層數為3層~6層的4幢無筋砌體結構進行4×(2×9+1)=76次IDA分析和地震易損性分析。其中IDA分析在主頻為2.8 GHz的Intel i9-7960X處理器配置電腦中進行，共耗時約180 h。

圖4為4幢無筋砌體結構在不同極限狀態下的地震易損性曲線，其中實線表示只考慮地震動不確定性時的計算結果，虛線為同時考慮地震動和結構參數不確定性時的計算結果。為考察地震動和結構參數不確定性對無筋砌體結構地震易損性的影響程度，表4進一步列出了βRTR、βMDL和βTOT三個分別考慮不同不確定性參數組合的易損性曲線對數標準差值。

由地震易損性曲線特性可知，曲線的傾斜程度表征了不確定性對地震易損性的影響程度。以圖4(c)中的6號易損性曲線為例，5層無筋砌體結構在PGA處于0.35g~0.7g時，其達到LS3極限狀態的概率為20%~80%，這表明不確定性導致了結構地震響應存在很大離散性，從而使得結構達到LS3極限狀態所需要的地震強度也變得較為離散。因此，圖4的結果表明不確定性對無筋砌體結構的地震易損性影響是不可忽視的。而這種影響在結構極限狀態LS3時表現的更為明顯，體現在圖4中即LS3對應的易損性曲線傾斜度明顯大于LS1和LS2。究其原因，是因為結構極限狀態LS3所對應的地震動輸入增強，使得結構的非線性程度增大，放大了不確定性對無筋砌體結構地震易損性的影響。與圖4的分析結果相呼應，表4中每幢結構極限狀態LS3所對應的βTOT值均明顯大于極限狀態LS1和LS2所對應的βTOT值。

圖4 考慮和不考慮結構參數不確定性的易損性曲線Fig.4 Fragility curves with and without structural parameter uncertainty

表4 不確定性對結構地震易損性的影響Table 4 Effect of uncertainty on the structural seismic fragility

由前文可知，βMDL與βRTR值的大小分別表征了結構參數及地震動不確定性對地震易損性的影響大小，其值越大表示影響越大。因此，計算了表4中βMDL與βRTR的比值并列于其中以對比結構參數以及地震動不確定性對無筋砌體結構地震易損性的影響大小，其最小值為6層LS3時的0.466，最大值為3層LS1時的1.191。同時，將表4中的βMDL與βRTR值圖形化作成柱狀圖進行直觀比較，如圖5所示，左斜線柱狀圖表示βRTR值，右斜線柱狀圖表示βMDL值。可見，結構參數的不確定性相比于地震動的不確定性是不可忽略的，有時甚至會超過地震動不確定性的影響，如3層結構對應于極限狀態LS1、LS2以及4層結構對應于極限狀態LS1的易損性分析結果中，βMDL值均比與βRTR值大。從圖4中也可以看出，考慮與不考慮結構參數不確定性的兩條易損性曲線并不重合，考慮結構參數不確定性后易損性曲線傾斜度有明顯增加。

圖5 地震易損性曲線對數標準差對比Fig.5 Comparison of logarithmic standard deviation of seismic fragility curves

同時，從表4及圖5可看出，βRTR值隨著結構層數的增加總體上呈現了明顯增長的趨勢，例如極限狀態LS2對應的βRTR值從3層的0.158增長到了6層的0.277。說明地震動不確定性的影響在大部分情況下會隨著結構層數的增加而增大。而同時也能發現，βMDL值并沒有隨著結構層數增加而增加，反而有下降趨勢。除了3層LS3的0.237外，其余的βMDL值均分布在0.138~0.191內，總體來看變化不大。這說明結構參數不確定性對不同層數無筋砌體結構地震易損性的影響基本相同。而這就導致，隨著結構層數的增加，結構參數不確定性相比于地震動不確定性的影響呈現減小趨勢，體現在圖4中即隨著層數的增加，考慮與不考慮結構參數不確定性的兩條地震易損性曲線越來越趨近于重合。由此可見，結構參數不確定性在較低層無筋砌體結構地震易損性中的相對影響更大，在計算時尤其需要得到重視。

為了進一步研究無筋砌體結構地震易損性分析中9個結構不確定性參數的敏感性大小，以3層無筋砌體結構為例，在前文計算的基礎上，再對每個參數分別取增加和減少2倍標準差進行地震易損性分析。此時，每個參數共有5種計算工況：以易損性曲線的地震強度中位值mR為評價指標，對每個變量的5種計算工況結果進行歸一化處理：

式中：mR(·)為參數xi取相應值而其他參數取平均值時的地震強度中位值；(·)為對應的歸一化結果。(·)偏離1越多，表示該參數敏感性越高。

結果如圖6所示。從圖6可見，除阻尼比外的8個參數變化時計算的易損性曲線地震強度中位值歸一化結果在0.85~1.1之間變化，變化范圍為0.25；而阻尼比變化時的歸一化結果最小值與最大值分別達到0.4和1.4，變化范圍為1，是前者的4倍。表明阻尼比的敏感性遠高于其他8個參數，是結構地震易損性分析結果的控制性參數。在無筋砌體結構地震易損性分析當中，阻尼比的準確取值尤其需要被重視。同時，中位值mR隨著阻尼比的增大而增大。這是由于隨著阻尼比的增大，結構耗能增加，結構響應隨之減小，從而導致結構易損性曲線的中位值不斷增大。

圖6 地震強度中位值隨各結構參數的變化規律Fig.6 Variation of the median earthquake intensity with different structural parameters

5 結論

本文以4幢典型無筋砌體結構為例，采用IDA和FOSM方法同時考慮地震動及結構參數不確定性對其地震易損性的影響，得到以下結論：

(1) 在無筋砌體結構的地震易損性分析中需考慮地震動和結構參數的不確定性影響。對于同一結構，地震破壞程度越高這種影響程度越大；

(2) 結構參數的不確定性相比于地震動的不確定性是不可忽略的，有時甚至會超過地震動不確定性的影響，且結構的層數越小，結構參數不確定性相比于地震動不確定性的影響越大，因此有必要在無筋砌體結構地震易損性分析中同時考慮這兩種不確定性；

(3) 以3層結構為例，對無筋砌體結構地震易損性分析中9個結構不確定性參數進行了敏感性分析。結果表明，阻尼比變化對易損性曲線地震強度中位值的影響可達到其他單個參數的4倍，其敏感性最高，同時隨著阻尼比的增大，易損性曲線的地震強度中位值也隨之增大。