金愛云， 王進廷， 潘堅文

(清華大學 水沙科學與水利水電工程國家重點實驗室， 北京 100084)

我國西南地區(qū)是許多大江大河的發(fā)源地，高山峻嶺密布，河流有著極大的天然落差，高混凝土壩多修建于此。然而西南地區(qū)也是我國的地震重災區(qū)，這些高壩的庫容都在幾十甚至上百億立方米，一旦在地震作用下發(fā)生潰決，將引發(fā)不可估量的生命財產損失和生態(tài)災難。因此，大壩在地震作用下的安全性問題成為大壩設計的重中之重。眾所周知，地震一般由前震，主震和余震構成。統(tǒng)計表明，約90%的強主震伴隨著強余震的發(fā)生[1]，主震強度越高則伴隨的余震強度越大[2]。近年來一些大震如2008年汶川地震、2013年蘆山地震和2015年尼泊爾Gorkha地震等，均伴隨發(fā)生了大量強余震[3-5]。大量震害資料表明強余震也具有較大的破壞能力。當主震導致結構損傷后，在強余震或多次余震的作用下，結構的損傷將會進一步累積，并可能最終導致結構倒塌。因此，有必要在目前高壩的抗震設計中進一步考慮余震的影響。

國內外已有一些學者對主余震序列作用下結構的動力響應規(guī)律進行了研究。楊佑發(fā)等[6]研究了余震對掉層結構易損性分析結果的影響，結果表明余震會提高掉層結構的超越概率。韓建平等[7]考慮了不同的主余震序列構造方法，研究了余震方向性和余震次數的影響。柳春光等[8]將兩次地震首尾相連構建主余震序列，研究了橋梁結構在兩次地震作用下的動力響應，表明連續(xù)兩次地震作用下結構的地震危險性增加。Jalayer等[9]以抗彎鋼架為例，考慮余震發(fā)生的時間依賴性，研究了余震引起的損傷累積效應和對極限狀態(tài)超越概率的影響。Wang等[10]以某32層框架芯管結構為例，研究了主余震序列作用下高層結構的地震脆弱性，結果表明隨著主震引起的破壞狀態(tài)變得越來越嚴重，主震下受損的建筑物對于余震激發(fā)變得越來越脆弱，并且對余震強度越來越敏感。Wen等[11]以鋼筋混凝土框架結構為研究對象，提出了在主余震序列作用下研究結構脆弱性的框架。在高壩領域，熊瑜等[12]將數值模擬與模型試驗結合，研究了復雜地基情況下主余震序列對重力壩穩(wěn)定性的影響，結果表明壩體變形相對單獨主震情況下大幅增加。翟亞飛等[13]結合NGA地震動衰減模型構建主余震序列，結果表明在主震導致壩體受損的情況下余震能引起明顯的殘余變形。王高輝等[14]研究了主余震序列作用下重力壩的損傷累積破壞效應，表明余震自身地震特征對壩體位移的重要影響。Zhang等[15]建立了局部損傷指標和整體損傷指標，對比了實測主余震序列和人工主余震序列，表明實測主余震序列能引起更大的損傷累積。Alliar等[16]研究了在考慮余震和排水系統(tǒng)效率降低的情況下重力壩的地震安全性。Pang 等[17-18]以某200 m高的混凝土面板堆石壩(concrete face rockfill dam，CFRD)為例，采用彈塑性分析的方法研究了主余震序列對壩體抗震性能的影響，結果表明CFRD在經歷余震時會遭受額外的地震損傷，且脆弱性顯著增加。

目前，高壩領域對主余震序列的研究尚缺乏對拱壩的研究成果。此外，已有的研究著重于分析余震是否會帶來影響，表明余震可能引起結構危險性增加，然而缺乏對余震作用下結構的破壞行為進行分析，也并未量化主震受損帶來的影響。從工程角度而言，人們更關心的是主震受損的結構還能承受多大強度的余震，也即極限抗震能力的變化程度，以便對后續(xù)修復工作進行指導。對該問題的研究需要在不同主震受損狀態(tài)下進行多次余震試算直至結構失效，同時應該考慮地震的不確定性。因此，為兼顧計算成本與計算結果的合理性，本文基于耐震時程法(endurance time analysis, ETA)[19]構建主震-ETA余震序列以進行高拱壩動力分析，并與單獨主震的情況進行比較，以量化主震受損對高拱壩極限抗震能力的影響，并給出不同主震強度下高拱壩極限抗震能力的變化情況。

1 分析模型

1.1 拱壩-庫水-地基有限元模型

選取大崗山混凝土雙曲拱壩為分析案例。大崗山拱壩最大壩高210 m，壩頂高程1 135 m，正常蓄水位高程1 130 m。有限元模型如圖1所示，由37 120個八結點六面體實體單元和53 817個結點組成。壩體由28條橫縫分為29個壩段，橫縫采用接觸邊界模型[21]進行模擬，使用Westergaard附加質量[22]模擬大壩與庫水之間的相互作用，采用黏彈性人工邊界[23]來模擬半無限地基的輻射阻尼效應。計算中考慮了壩體混凝土的非線性行為，并采用Lee等[24]提出的混凝土塑性損傷本構模型。壩體阻尼采用Rayleigh阻尼，取第一階和第五階的阻尼比為5%。壩體混凝土密度2 400 kg/m3，初始彈性模量31.2 GPa，泊松比0.17，抗拉強度ft為3.25 MPa，斷裂能Gf為280 N/m，極限拉應變εf為4×10-4。地基彈性模量26.0 GPa，密度2 625 kg/m3，泊松比為0.25。

(a) 壩體-地基有限元模型

1.2 耐震時程法

耐震時程法由Estekanchi等提出，是一種動態(tài)的Push-Over方法，通過擬合出滿足特定條件的人工地震動開展結構動力分析，其重要特征之一是地震強度隨時間持續(xù)增長。具體來說，耐震時程法中人工地震動在某時段(0～t)內計算得到的加速度反應譜(本文定義為t時刻對應的加速度反應譜)與該時段的持時t呈線性關系

(1)

式中：T為結構自振周期；t為地震動任意時間點；tR為預先指定的目標時刻；Sa，Tar(T,t)為地震動在t時刻的目標加速度譜；Sa，R(T)為預先指定的tR時刻對應的加速度譜。通過指定tR和Sa，R(T)，便能得到任意時刻的目標加速度譜，并以此為目標進行地震動擬合。

耐震時程法中，人工地震動在0～t時段內計算得到的位移反應譜為

(2)

式中，Sd，Tar(T,t)為t時刻的目標速度譜。對某一時刻而言，可以使擬合得到的地震動在一定精度條件下同時滿足式(1)和式(2)。然而，耐震時程法在理論上要求式(1)和式(2)對任意時刻都成立，這顯然是沒有必要，也無法實現的。通常，簡化為在某些時刻及一定精度條件下滿足式(1)和式(2)，此時耐震時程法中的地震擬合實際上是一個無邊界的優(yōu)化問題

Sa,Tar(T,t)]2+αD[Sd(T,t)-

Sd，Tar(T,t)]2}dtdT

(3)

式中：ag為擬合得到的地震動；Sa(T,t)和Sd(T,t)分別為在t時刻基于ag計算得到的加速度譜和位移譜；αD為位移譜的權重系數，根據Estekanchi等的建議，αD越小則加速度反應譜在短周期范圍內擬合得越好，因此本文中取αD=0，Tmax和tM分別為擬合中所考慮的最大周期和最長地震持續(xù)時間，分別為3 s與30 s。

采用以上步驟所擬合得到的地震動的示意圖如圖2所示。通過人為控制，可以使耐震時程法中的最大地震強度值(本文用基頻對應的譜加速度Sa(T1)表示地震動強度值，T1為基頻對應的自振周期。)盡可能大，從而通過一次ETA計算便可以得到結構從完好到破壞的過程中不同地震強度所對應的結構動力響應，相當于一次完整的增量動力分析法(incremental dynamic analysis ，IDA)計算，能顯著減少計算工作量。

圖2 耐震時程法擬合的地震動示意圖

1.3 地震記錄

為研究余震作用下高拱壩的破壞行為，量化主震受損對高拱壩極限抗震能力的影響，本文提出將高拱壩極限抗震能力損失程度定義為

(4)

式中:dC為高拱壩震后極限抗震能力的損失程度；Cb為高拱壩完好時的極限抗震能力；Cp為高拱壩震后的極限抗震能力。

需要說明的是，本文中極限抗震能力代表高拱壩破壞時對應的地震強度(破壞強度)。由于地震存在不確定性，則應通過破壞易損性分析[25]求出單獨主震情況下高拱壩不同極限抗震能力的概率水平，同時考慮多條主震-ETA余震序列用于計算。因此，本文中根據DL 5073—2000《水工建筑物抗震設計規(guī)范》[26]中指定的加速度反應譜，采用文獻[27]中的擬合方法，擬合了12組三向人工地震動作為主震，同時依據式(3)擬合了12條ETA地震動作為余震，將它們隨機組合成12條主震-ETA余震序列。主震時長20 s，主震震級考慮Sa(T1)=0.4～1.0g共7個等級，主余震之間間隔5 s，ETA余震時長30 s，最大地震強度為2.0g，因此一共構建了84條不同的主震-ETA余震序列，如圖3所示。另外，在計算Cb時，本文采用以ETA地震記錄作為單獨主震的形式進行分析，以避免IDA中所需的多次計算。需要說明的是，由于本文主要探究主震受損情況下拱壩抗震能力的變化情況，為了使單獨主震和主余震序列的計算結果具備可比性，因此并未考慮主震和余震的頻譜特征差異。

圖3 主震—ETA余震示意圖

2 高拱壩非線性動力分析

2.1 單獨ETA主震分析

首先將ETA地震記錄單獨作為主震進行動力分析，并求出拱壩的破壞強度Cb。分析中所選擇的工程需求參數(engineering demand parameters, EDP)為損傷體積比(Rdt>0)、開裂體積比(Rdt>0.8)和頂拱朝上下游最大動力相對位移(Δu，Δd),其中dt為拉伸損傷因子。圖4給出了ETA分析結果中一些EDP隨著譜加速度(時間)的變化情況(ETA曲線)。圖4中每條ETA曲線都可以等價為一條IDA曲線，與IDA方法類似，ETA分析中損傷與位移結果也隨著地震強度的增加不斷增大。計算結束時(Sa(T1)=2.0g)損傷體積比Rdt>0和開裂體積比Rdt>0.8都接近100%，上游位移Δu甚至達到了數米，說明此時拱壩損傷非常嚴重，實際上早在計算結束前就已經破壞。因此，可以根據ETA的計算結果求出每條ETA地震記錄對應的Cb。

(a) 損傷

為了驗證ETA計算結果的合理性，需要將其與IDA的計算結果進行對比。在基于主震-ETA余震序列的動力分析中，將各條主震的計算結果進行匯總便得到IDA的計算結果。圖5將單獨ETA的計算結果與IDA的計算結果進行了對比。為方便比較，圖5中ETA曲線為12個ETA計算工況的平均值。從圖5中可以看到，在IDA的強度范圍內IDA的損傷結果大致分布在ETA均值的兩側，Rdt>0的吻合情況較好，尤其是在地震強度不大的時候。對Rdt>0.8來說ETA均值處于IDA數據范圍中偏下的區(qū)域，因此ETA對開裂體積比的結果略低于IDA。位移的差異比損傷差異稍大，在地震強度較小時(約0.6g之前)，IDA中Δd與Δu幾乎相等，但ETA中Δd均值明顯比Δu更大。此外，地震強度較小時ETA與IDA在Δd上并無太大差異，隨地震強度增大ETA出現高估Δd的現象，對于Δu而言在IDA的計算范圍內ETA偏小。

(a) 損傷

綜合以上結果，表明在IDA的計算范圍內，ETA與IDA的損傷分析結果接近，但位移結果上存在一定區(qū)別。與IDA相比，ETA的Δd更大，但Δu更小。整體而言，在所計算的IDA強度范圍內ETA的計算結果略微小于IDA，但仍具有可比性。

2.2 極限抗震能力

對于拱壩而言，由于資料有限，在進行數值模擬時對破壞的定義并沒有統(tǒng)一標準。目前，包括拱壩動力模型試驗在內的一些研究中[28-29]將上下游裂縫貫穿視為拱壩破壞的標志，該方法簡單可行，且一定程度上符合人們的直觀判斷。盡管目前對裂縫貫穿后梁頂自由塊體的動力穩(wěn)定性尚缺乏認識，本文從實踐性的角度出發(fā)，采用這種方法定義拱壩的破壞，認為當連續(xù)兩個或以上壩段出現貫穿裂縫時拱壩破壞，并求出對應時刻的地震強度。其中裂縫貫穿的標志是損傷因子dt>0.8的損傷區(qū)從下游貫穿至上游，如圖6所示。

基于此判定標準，對每條ETA地震記錄都能獲取對應的破壞強度Cb，并通過矩法[30]擬合出拱壩的破壞易損性曲線，如圖7所示，表示某個地震強度下拱壩出現破壞的概率，或拱壩的極限抗震能力等于某個特定值的概率。

圖7 基于ETA分析的拱壩破壞易損性曲線

2.3 主震-ETA余震序列分析

2.1節(jié)、2.2節(jié)給出的都是單一地震作用下的結果，下面對基于主震-ETA余震序列的動力響應結果進行分析。圖8以主震強度Sa(T1)=1.0g的情況為例，給出動力響應結果隨時間變化的過程，圖8中0～20 s為主震，25～55 s為ETA余震。

(a) 損傷

圖8中動力分析結果隨時間的增加過程與單一地震時有明顯區(qū)別。因為25～30 s為自由振動，余震開始前主震的影響已經衰減到足夠小，因此壩體的最大動力響應在該時間范圍出現一個平臺。這也說明5 s的時間間隔是足夠的，不會對余震響應帶來影響。進入ETA余震作用時間后，當余震強度較小的時候，壩體最大動力響應保持不變，直到ETA余震的強度超過某一水平后壩體最大動力響應重新增加。然而，需要注意的是ETA余震重新引起壩體最大動力響應增加時對應的時刻(或ETA余震強度)。表1給出了當損傷體積比Rdt>0在ETA余震作用下顯著增加(增加10%)時對應的余震強度。

表1 ETA余震引起Rdt>0顯著增加時對應的強度

由表1可知，當拱壩動力響應在余震作用下顯著增大時，對應的余震強度可能大于或小于主震強度，但平均而言兩者相當。這與我們的認知有一定出入，因為大多數情況下研究人員認為壩體一旦產生裂縫就可能在并不強的余震作用下導致裂縫進一步擴展，而表中數據卻顯示即使對于主震為1.0g的情況(Rdt>0約60%，壩體已明顯損傷)，余震強度也必須足夠大才能使損傷進一步擴展。對于該現象，本文認為一方面它表明拱壩具有良好的抗震性能；另一方面也可能是有限元分析模型的局限性造成的。本文采用的混凝土塑性損傷本構模型本質上是從能量的角度模擬裂縫的生成和擴展，無法從幾何上對裂縫進行模擬。從圖6的梁截面損傷圖也能看出，若以損傷因子dt>0.8定義裂縫，則梁截面上裂縫范圍很廣且上下連成一片，而現實中混凝土是脆性斷裂，其拉裂縫一旦形成則周圍混凝土的拉應力得到釋放，裂縫寬度應該比較窄。對于混凝土裂縫擴展的模擬一直是熱點問題，Pan等采用擴展有限元模擬Koyna重力壩的壩頸拉裂縫，得到較好的模擬效果，但受限于分析理論，該方法暫時無法用于拱壩的三維非線性動力計算。目前，模擬高拱壩壩體混凝土開裂的方法仍然以塑性損傷理論為主，更有效的模擬方式仍需進一步探究。

3 高拱壩極限抗震能力損失

在基于主震-ETA余震序列的動力分析中，當拱壩在主震結束時破壞則認為其沒有后續(xù)承載能力，抗震能力完全喪失，即dC=100%。若主震結束后拱壩沒有破壞，則在ETA余震計算過程求出拱壩破壞對應的破壞強度，即震后的極限抗震能力Cp，然后基于式(4)計算dC。表2給出了主震未引起拱壩破壞時dC的統(tǒng)計值，其中缺省值表示主震導致破壞(dC=100%)，表2中平均值的計算并不考慮該種情況。從表2中可以看出，只有主震強度達到一定水平時才會引起拱壩極限抗震能力下降，隨著主震強度的增加，主震引起的損傷增大，拱壩極限抗震能力逐漸減小，dC不斷增大。圖9給出了表2中dC的平均值隨著主震強度的變化情況，并對其進行擬合。關于擬合曲線的分布形式有以下兩個方面的考慮：首先主震強度為正值，且當主震強度很小時dC接近于0；其次當主震強度很大導致拱壩接近破壞時，dC趨近于1。因此采用對數正態(tài)分布比較符合這兩種情況。需要指出的是，圖9中曲線的擬合是基于最小誤差平方和原理進行的，由于數據點只分布在部分強度范圍，因此得到的曲線只是對當前數據的估計，當更多數據得到補充時，曲線的形狀可能會有所變化。

表2 主震未引起拱壩破壞時dC的統(tǒng)計

圖9 主震未引起拱壩破壞時dC(均值)隨主震強度變化情況

考慮到主震有可能引起拱壩破壞，因此結合圖7的拱壩破壞易損性曲線，對每個主震等級下的dC進行加權處理

dC|IM=im=P[D|IM=im]×1 +(1-P[D|IM=im])×dC|ND,IM=im

(5)

式中：IM為主震強度；im為具體強度值，即Sa(T1)的值；P[D|IM=im]為給定主震強度下拱壩破壞的概率，可通過破壞易損性曲線求出；dC|ND,IM=im為主震未導致破壞時dC的期望值，可以通過圖9的擬合曲線得到；dC|IM=im為考慮主震引起破壞時給定主震強度下dC的期望值。

dC|IM=im隨主震強度變化情況(極限抗震能力損失曲線)如圖10所示，同時也給出了圖9中不考慮主震破壞的情況作為對比。很顯然，由于考慮了主震引起破壞的可能性，主震強度較大時dC的期望值顯著增大了。此外，當主震強度較小時，拱壩很可能保持完好無損，因此極限抗震能力幾乎不發(fā)生損失，dC的期望值接近0。隨著主震強度的增大，dC的期望值也迅速增大，當拱壩接近破壞時逐漸趨近100%。根據現行GB 51247—2018《水工建筑物抗震設計規(guī)范》[31]中拱壩的標準設計反應譜，大崗山拱壩在設計地震下的Sa(T1)為6.94 m/s2，由圖10可知設計地震引起極限抗震能力損失的期望值小于10%，因此壩體預期能經受設計地震并仍有較強的抗震能力。

圖10 拱壩極限抗震能力損失曲線

4 結 論

本文對拱壩在遭受主震損傷后的極限抗震能力進行了研究，提出了極限抗震能力損失的概念?；谀驼饡r程分析法(ETA)擬合了ETA地震記錄，并構建了主震-ETA余震序列對大崗山拱壩進行了動力分析，最終給出大崗山拱壩的極限抗震能力損失曲線。本文得出的主要結論如下。

(1) ETA與IDA的計算結果具備一定可比性，但兩者并不嚴格等效。當目標加速度譜一樣時，ETA略微低估了壩體損傷和頂拱上游位移，高估了頂拱朝下游的位移，總的來說低估了壩體的動力響應程度。

(2) 當采用塑性損傷模型模擬壩體混凝土開裂時，拱壩的損傷在主震-ETA余震序列計算過程中逐漸累積，并在主震結束后保持一段時間的恒定。ETA余震開始作用后拱壩的損傷并不會立即恢復增加，且只有當ETA余震強度和主震強度相當時才能使損傷結果出現顯著增加。

(3) 主震達到一定強度時會引起壩體損傷，這導致拱壩的極限抗震能力出現一定程度的損失，且隨著主震強度增大，拱壩極限抗震能力的損失程度也隨之增大。就本文的計算結果而言，設計地震作用下大崗山拱壩極限抗震能力損失的期望值小于10%，壩體仍有較強抗震能力。

極限抗震能力損失曲線可以對未知地震的危害做出定量評價，作為工程決策的參考手段。后續(xù)的研究可以考慮針對此概念進行完善，豐富極限抗震能力損失的評價手段，并采用更合理的方法構建主余震序列，以提高分析結果的可靠程度。