徐 強 姚文斌 陳健云

(1.大連理工大學 海岸與近海工程國家重點實驗室, 遼寧 大連 116024;2.大連理工大學 工程抗震研究所,遼寧 大連 116024)

我國水能資源豐富,尤其是西南地區(qū),具有較高的開發(fā)利用前景.因此,我國在該區(qū)域修建了數座百米級高拱壩,在灌溉、防洪、供水、發(fā)電、航運及水產養(yǎng)殖等領域發(fā)揮著重要作用.但是,大壩作為主要的擋水建筑物,一旦發(fā)生破壞,將直接影響下游人民的生命和財產安全.近年來,我國唐山、汶川等地發(fā)生了烈度較高的強地震動,超過了部分建筑物的設計地震烈度,使當地的各類建筑物遭受了不同程度的破壞.地震動作為大型結構損傷破壞的重要因素,研究拱壩在地震動作用下的損傷破壞機理是有必要的.

對于特高拱壩這種大型混凝土結構,由于其結構復雜,探究壩體在不同地震動強度下的非線性響應時,需要反復進行加速度時程的調幅計算.然而這一過程十分耗時,而且計算結果可能無法收斂.為了解決這一問題,Estekanchi[1]提出了一種耐震時程法(endurance time analysis,ETA),詳細闡述了耐震時程加速度曲線的合成方法,并驗證了該方法在保證計算精度的同時,能有效減少計算量.白久林、楊樂等[2]基于我國的抗震反應譜,合成了ETA 時程曲線,并將其與增量動力分析方法(IDA)進行了對比,結果驗證了ETA 方法的可行性.

地震的危險性分析、易損性分析和災害損失評估是預測地震風險的3個主要方面[3].其中,易損性分析作為評估結構安全性的核心內容,近年來在研究大壩的抗震性能方面被廣泛關注.大壩地震易損性分析是為了預測大壩受不同強度地震動作用時,能滿足各級性能水準要求的可能性.地震易損性分析按數據資料的來源可分為經驗法與解析法.經驗法就是利用現有地震記錄及震害資料分析易損性,解析法主要通過多組數值計算模擬地震中結構的反應.對高拱壩進行易損性分析可以發(fā)現結構在設計中存在的薄弱環(huán)節(jié),為后期的除險加固提供依據[4-6].

國內學者曾對拱壩的地震易損性進行了研究.金愛云[7]分別采用了云分析、多條帶分析(MSA)以及增量動力分析方法研究了高拱壩在地震動作用下的易損性曲線,結果表明多條帶分析方法的可靠度最高;鐘紅等[8]考慮了多種參數的不確定性,定義了拱壩的5級極限狀態(tài),分別為無損傷、輕微損傷、中度損傷、嚴重損傷和破壞.姚霄雯和蔣建群[3]等考慮了拱壩的材料非線性和橫縫接觸非線性,分析了塑性損傷、橫縫數量對拱壩動力響應的影響,以及15種地震強度指標與拱壩動力響應的相關性.最終,他們以壩頂位移和橫縫開度為指標,提出了基于性能的高拱壩抗震安全評價方法.王進廷等[9]以大崗山拱壩為例,考慮材料不確定性,定義了幾種極限狀態(tài),分別為輕微破壞、中等破壞和嚴重破壞,基于拱壩的損傷程度以及是否發(fā)生貫穿性損傷.章明旭[5]比較了直接擬合法和概率分布法用于擬合地震易損性曲線時的差異.范書立等[10-11]基于向量地震強度參數和多種工程需求參數,建立了以拱冠位移、橫縫開度和損傷體積比為指標的地震易損性曲線和曲面,采用了基于響應面的方法進行了拱壩的地震易損性分析,研究結果表明響應面分析比傳統(tǒng)的蒙特卡羅模擬法效率更高.這些研究成果為高拱壩抗震設計和評估提供了科學依據.

基于上述理論方法,本文根據水工結構設計規(guī)范[12]中的標準設計反應譜生成了20 組耐震時程加速度曲線;建立了白鶴灘拱壩-地基系統(tǒng)有限元模型,進行非線性時程分析,分析結構的損傷特點;選取相對位移、損傷體積比和最大橫縫開度為指標,提取各指標時程,并分別基于MSA 法、回歸曲線法、矩估計法及截斷最大似然估計法對結構的易損性進行了分析,得出易損性曲線.通過對高拱壩進行抗震安全評價及失效概率分析,可為結構的抗震設計提供可靠依據.

1 耐震時程法

耐震時程法是一種動力推覆過程,其主要優(yōu)點在于保證預測結構響應精度的同時,可以減少計算量[13-14].該方法的核心思想是生成ETA 加速度時程曲線,使得不同時間段對應的反應譜呈現既定的倍數擬合關系.

ETA 加速度時程的生成過程是先生成平穩(wěn)加速度時程,然后通過反復迭代優(yōu)化,使得加速度反應譜與目標反應譜的擬合程度更佳.在地震動持續(xù)時間內,不同時間段對應的目標反應譜與目標反應譜的倍數關系為:

式中:tTarget為目標時間;SaC(T)為標準設計反應譜;SaT(T,t)為目標反應譜;t為任意時間.

根據加速度反應譜與位移反應譜的對應關系可以推導出位移譜的公式為:

式中:SuT(T,t)為時刻t的目標位移反應譜.

但是,想讓加速度時程上的任一時刻前對應的反應譜都能以目標譜完美擬合是不可能的,因此,可以采用式(3)對加速度時程進行無約束優(yōu)化,使其對應的反應譜與目標譜更加吻合:

式中:ag為初始生成的加速度時程曲線;Tmax是地震動時長;α是位移譜的權重系數,本文只考慮加速度反應譜的影響,因此α值取0;Sa(T,t)和Su(T,t)分別為周期T下0～t時刻的加速度反應譜和位移反應譜.

基于上述理論,根據標準設計反應譜,共生成20組隨機ETA 地震動時程,每組包含3條地震動時程,分別作為X、Y和Z向的動力輸入.其中,X向和Y向的最大峰值加速度調幅為0.6g,Z向最大峰值加速度調幅為0.4g.圖1展示了ETA 時程和IDA 時程的對比圖.從圖中可以觀察到,ETA 時程主要具有峰值加速度不斷增大的特點.分別對ETA 時程0～5 s、0～10 s、0～15 s和0～20 s各時段時程進行反應譜優(yōu)化擬合,擬合效果如圖2所示,可以看出,不同時段的反應譜與標準譜呈倍數關系.

圖1 ETA 與IDA 加速度時程對比

圖2 ETA 時程加速度反應譜

2 結構地震易損性分析方法

結構地震易損性是指結構在某一強度地震動作用下,結構達到極限狀態(tài)的概率[15],采用下式對其進行描述:

式中:F(im)為易損性概率;LS為指定的破壞指標極限狀態(tài),本文選取相對位移、最大橫縫開度和損傷體積比作為破壞指標;IM為表征地震動強度的指標,本文選取峰值加速度PGA 為地震動強度指標;P[LS|IM=im]為地震動參數IM=im時,結構達到極限狀態(tài)的概率.

2.1 基于回歸曲線的易損性分析方法

假設地震動需求(Dd)與IM呈對數相關性,采用直線擬合lnDd與lnIM之間的關系:

式中:a、b表示回歸參數;Dd表示地震動需求,也就是本文選取的損傷指標.

lnDd在每個地震動強度IM水平下,具有恒定的方差,其標準差可表示為:

2.2 基于矩估計的易損性分析方法

用含兩個參數的對數正態(tài)累計分布函數來定義易損性函數,其表達式為[17]:

式中:θ為在超越概率為50%的IM水平下易損性函數的中位數;β表示lnIM的標準差.通過對任意工況下結構開始發(fā)生破壞的IM值取對數來對易損性函數的參數進行估計,其均值和標準差分別為[18]:

式中:n為工況總數量;IMi為第i個工況下結構開始發(fā)生破壞的PGA 值.

2.3 基于截斷最大似然估計的易損性分析方法

最大似然估計法的原理是通過使似然函數達到最大值來對參數進行估計.通過設定所考慮的最大地震動強度IMmax,假設在進行的n次地震動分析中,有m次結構在低于IMmax的IM水平下發(fā)生破壞,通過式(12)對θ和β進行估計:

式 中:^θ和^β是 使 得 該 函 數 取 得 最 大 值 的θ和β值;IMi為第i個工況下結構開始發(fā)生破壞的PGA 值.該函數需要借助MATLAB 中的fmincon函數進行求解.將參數估計值代入式(9)中求解易損性曲線.

2.4基于多樣條分析(MSA)的易損性分析方法

MSA 分析方法通過分析每一個IM水平下,損傷指標超過設定閾值的工況所占總工況比例,基于最大似然估計對參數進行估計[19]:

式中:m表示不同地震動強度的個數;xj表示逐漸增大的第j個地震動強度;zj和nj分別表示工況總數和在地震動強度為xj條件下導致結構發(fā)生破壞的工況數.式(13)的求解過程同式(12).

3 數值模擬

3.1 有限元模型

以白鶴灘拱壩為研究對象,其壩高為289 m,壩頂寬度13 m,最大壩底寬度72 m,壩頂高程為834 m,壩頂中心線弧長為709 m.壩體-地基系統(tǒng)有限元模型如圖3所示;壩體共設有30條橫縫,橫縫分布如圖4所示.

圖3 白鶴灘拱壩-地基系統(tǒng)有限元模型

圖4 白鶴灘拱壩橫縫分布

分析中考慮的靜態(tài)荷載有壩體自重和靜水壓力.壩體混凝土密度為2 400 kg/m3,地基密度為2 800 kg/m3.正常蓄水情況下,壩體上游靜水壓力水位為825 m,下游靜水壓力水位為604 m.壩體動水壓力按照Westergaard附加質量添加,采用附加質量的形式模擬不可壓縮水體對結構的動水壓力.地震動采用ETA 生成的耐震加速度時程輸入.

采用接觸邊界模型模擬混凝土拱壩的橫縫,每個橫縫接觸面包含主面和從面,接觸關系[18]如圖5所示,法向接觸壓力的計算為:

圖5 主從面接觸模型

式中:初始接觸距離c取100 mm,初始壓力p0取0.3 MPa.

為模擬地基的彈性恢復能力,采用黏彈性人工邊界[20]來實現有限地基對無限地基的模擬作用,該人工邊界還能吸收邊界輻射能量,模擬效果良好.示意圖如圖6所示.人工邊界上,節(jié)點上的參數由式(15)計算給出:

圖6 黏彈性人工邊界示意圖

式中:Kn和Ks分別表示彈簧法向和切向的彈性剛度;Cn和Cs分別表示法向和切向的阻尼;Cp和Cs表示地震波的傳播速度;ρ為地基密度;r為震源距;λ和G為拉梅常數;a和b為修正系數,a取0.8,b取1.1.

壩體混凝土材料和地基巖石的其他力學性能參數分別見表1和表2.

表1 壩體混凝土材料參數

表2 地基基本參數

3.2 動力響應分析

對白鶴灘拱壩進行動力學分析.圖7展示了壩體結構在峰值加速度水平(PGA)分別為0.2g、0.4g和0.6g時,拱冠梁、下游面和上游面的拉應力損傷分布情況.從圖中可以觀察出,隨著PGA 的增大,損傷成明顯的加重趨勢.在PGA 為0.2g時,只有壩踵和下游面下部存在局部輕微拉應力損傷;PGA 達到0.4g時,壩踵處損傷加重,下游面多處出現大面積損傷,但損傷值均為超過0.5,上游面損傷加重不明顯;當PGA 達到0.6g時,壩踵處損傷進一步加重,下游面出現大面積嚴重拉應力損傷,上游面也在多處呈現不同程度損傷.拱冠梁的拉應力損傷云圖也呈現明顯加重的趨勢.

圖7 不同PGA 水平下拉應力損傷云圖

圖8 各指標隨PGA 的變化時程

通過對損傷云圖的分析,在遭遇地震動作用時,最先被破壞的是壩踵部位,若PGA 繼續(xù)增大,下游面會出現損傷,上游面需要在更大的PGA 水平地震動作用下才會出現損傷.此時,壩體在中上部出現貫穿損傷,結構存在失效風險.

3.3 損傷指標分析

采用最大橫縫開度、損傷體積比和相對位移作為損傷評價指標.計算結果顯示,隨著PGA 值增大,各指標呈現明顯的增大趨勢.最大橫縫開度是取各時刻30條橫縫的最大開度作為該時刻的開度值.最大橫縫開度隨PGA 的增大平穩(wěn)增長,PGA 達到0.6g時,各工況最大橫縫開度最大值為9 cm,最小值約4.7 cm,平均值約為7 cm;損傷體積比[21]可以宏觀地衡量損傷單元所占的比例,損傷體積在PGA 較小時增速較為平緩,但隨著PGA 的增大旗增速逐漸加快,峰值加速度0.6g時,各工況損傷體積比最大值為0.28,最小值約0.09,平均值約為0.17;順河向的相對位移是指拱冠梁底部和頂部之間的順河向相對位移,與損傷體積比相反,在PGA 較小時增速較快,但隨著PGA 的增大旗增速逐漸減緩,峰值加速度0.6g時,各工況順河向相對位移最大值為0.34 m,最小值約0.18 m,平均值約為0.26 m.

為分析各個工況在不同PGA 水平下的離散性,分別計算了各指標在不同峰值加速度水平下的標準差,計算結果如圖9所示.從圖中可以觀察到,標準差的整體走勢是隨著PGA 的增大逐漸增大.說明PGA越大,各工況離散性越高.但各指標的標準差大小和增減趨勢均有所差異,PGA 在小于0.4g情況下,損傷體積比和最大橫縫開度各工況的標準差較小,但在0.4g后,各工況損傷體積比標準差迅速增大.順河向相對位移各工況標準差始終處于較高水平,說明各工況計算結果有較大差異性.

圖9 各指標標準差隨PGA 變化曲線

從離散性角度出發(fā),最大橫縫開度始終保持較小的標準差,更適合作為評判拱壩結構損傷破壞程度的指標.

3.4 易損性分析

本節(jié)以上文介紹的易損性分析方法,基于最大橫縫開度、損傷體積比和順河向相對位移等指標,根據以往研究成果[5,10],結合實際分析結果設定損傷閾值,對拱壩進行易損性分析,并評價了各方法與ETA方法相結合的分析效果.

3.4.1 最大橫縫開度

以最大橫縫開度為損傷指標,設定其損傷閾值為3 cm,分別采用截斷最大似然估計法、MSA、矩估計法和回歸曲線擬合法分析指標在各PGA 水平下達到該閾值的概率,計算結果如圖10所示.

圖10 以最大橫縫開度為指標的易損性曲線

分析結果顯示,矩估計與MSA 方法具有較高的一致性;通過截斷最大似然估計計算的易損性曲線在PGA 小于0.3g時偏大,0.35g之后與其他方法計算結果較為接近;通過一次函數擬合方法計算的易損性曲線在0.2g～0.4g范圍內偏小.

3.4.2 損傷體積比

以損傷體積比為損傷指標,設定其損傷閾值為0.02,通過各分析方法計算的易損性曲線如圖11所示.其中,矩估計法和MSA 法依然保持了較好的一致性,一次函數擬合和截斷最大似然估計法計算結果偏大.

圖11 以損傷體積比為指標的易損性曲線

3.4.3 順河向相對位移

以順河向相對位移為損傷指標,設定其損傷閾值為14 cm.易損性計算結果如圖12所示.同樣,矩估計和MSA 方法的計算結果高度一致,截斷最大似然估計的計算結果最大,一次函數擬合方法的結果稍大,且曲線局部出現震蕩現象.

圖12 以順河向相對位移為指標的易損性曲線

通過對上述3個損傷指標進行易損性分析,發(fā)現矩估計和MSA 方法可以較好地與ETA 方法相結合,對拱壩的損傷和易損性進行評估.原因是上述兩種方法均依賴于給定PGA 水平下拱壩結構的響應,且給定的PGA 水平越多,曲線精度越高.本文采用的耐震時程法可以給定任意峰值加速度水平下結構的響應,因此,耐震時程法與矩估計和MSA 法相結合可以以理想的精度計算易損性曲線;雖然一次函數擬合法同樣依賴于給定PGA 水平的數量,但給定更多PGA 水平下的結構響應并不能有效提高易損性曲線的精度,因此,分析結果往往有一定偏差,且曲線局部易出現震蕩現象;截斷最大似然估計需要預先給定地震動強度閾值,這一閾值的選取沒有確定方法,且該閾值往往對易損性分析結果產生影響.

4 結 論

本文建立了白鶴灘拱壩-地基有限元模型,基于耐震時程法理論,生成了耐震時程加速度曲線,作為動力輸入,作用于白鶴灘拱壩,進行動力時程分析.對多個PGA 水平下拱壩結構的拉應力損傷分布進行分析,以最大橫縫開度、損傷體積比和順河向相對位移為損傷指標,分別基于截斷最大似然估計法、MSA、矩估計法和一次函數擬合法分析了拱壩的易損性,主要得出以下結論:

1)拱壩在地震動作用下的損傷破壞主要集中在壩踵和下游面中上部,隨著地震動強度的不斷增強,上游面出現損傷,并出現貫穿損傷.

2)標準差能反映各工況的離散性,在最大橫縫開度、損傷體積比和順河向相對位移3個損傷指標中,隨著PGA 的增大,標準差均呈現增大趨勢,但最大橫縫開度具有相對較小的標準差,更適合作為評價指標.

3)由于計算原理不同,MSA 法和矩估計法更適合與ETA 法相結合,計算結構在不同PGA 水平下的易損性概率.