基于耐震時程法的特高拱壩極限抗震能力研究

馬天驍， 張燎軍， 張漢云， 史俊飛， 翟亞飛

(河海大學 水利水電學院， 南京 210098)

我國西南地區修建了大量的300 m級特高拱壩，這些高壩樞紐工程，規模和數量在世界上都是史無前例的，而處在地震動活躍地區的特高拱壩的抗震安全尤其值得關注[1-2]。在對特高拱壩進行抗震性能評價時，采用常規的時程分析法計算常常面臨著結構模型復雜、自由度數目較多、計算耗時長的問題。特別是當模型存在同時考慮壩體-基巖損傷的材料非線性和各類縫面的接觸非線性時，一次分析往往需要耗費大量的計算時間，效率較低。對于特高拱壩的極限抗震能力分析，則需要多次大量的時程分析計算才能最后確定，計算成本更高。為了解決這一問題，需要采用一種既能考慮時程分析中的動力效應，又能簡便高效的對特高拱壩抗震性能和失效模式進行評價的方法。Estekanchi 等[3-4]提出的耐震時程法，是一種基于時域的新型彈塑性動力分析方法。它的優點是可以通過一條加速度時程曲線來代表不同強度的地震系列，只需要計算一次耐震時程分析就能夠得到不同地震動強度下結構的動力響應。自該方法提出以來，在鋼框架及橋梁結構動力分析中得到了較多的應用[5-8]。在水工結構動力分析領域，徐強等[9-10]采用該方法對重力壩的動力損傷指標進行研究，Hariri-Ardebili等[11]采用耐震時程法分析了強震作用下拱壩損傷模式，并對高拱壩動力損傷進行安全評價。但對于300 m級特高拱壩，目前尚缺少采用耐震時程法全面考慮壩體及基巖的損傷、橫縫接觸非線性等方面的研究。

特高拱壩的極限抗震能力研究是地震時確保不發生庫水失控下泄的關鍵課題之一[12]。很多學者都對高拱壩的極限抗震能力進行研究。如Cheng等[13]提出位移突變、塑性區貫通及計算不收斂等判據來分析高拱壩的極限抗震能力；李德玉等[14]同時采用地震動超載和基巖結構面降強等方法分析了白鶴灘拱壩的極限抗震能力；張社榮等[15]從位移響應，塑性區貫通，壩體開裂等角度對重力壩的極限抗震能力進行研究；李靜等[16]采用地震動超載法，從壩體損傷及橫縫開度等角度研究了高拱壩的抗震安全性能。以上研究大多是采用地震動超載法，采用多次時程分析計算來確定高拱壩的極限抗震能力。但由于高拱壩計算模型較為復雜，大量的時程分析計算比較耗時，目前尚少有采用耐震時程法來研究高拱壩極限抗震能力的研究，對采用耐震時程法和常規多次時程分析法計算極限抗震能力的結果也缺少對比研究。

基于上述原因，本文采用耐震時程法研究高拱壩的動力損傷行為及極限抗震能力。以中國西南某300 m級特高拱壩為例，建立了可以同時考慮壩體混凝土及壩基巖體的損傷、橫縫的動態接觸非線性的有限元分析模型。基于耐震時程法，采用非線性最小二乘函數合成耐震時程曲線。分別采用耐震時程法和常規的時程分析法，對比分析高拱壩位移、損傷體積比、損傷耗能、橫縫最大開度等動力響應。結果表明，耐震時程法能夠有效的反映不同地震動強度下高拱壩的損傷情況，其結果與常規時程分析法的計算結果較為接近，計算效率較高，在節省計算時間和計算資源方面具有優勢。

1 基于耐震時程法的極限抗震能力分析方法

1.1 耐震時程法理論

耐震時程法(ETA法)是通過對特高拱壩計算模型輸入一條地震動強度隨時間逐漸增大的時程曲線，以分析壩體及基巖系統從線性階段到非線性階段的動力響應，直至最終破壞的全過程。該方法最主要的特點是隨著地震動持續時間的增加，地震動強度逐漸增大，且在某一個時間段內的目標加速度反應譜與該時段的持續時間t呈線性關系

(1)

式中：tt為目標時間點；SaC(T)為事先確定的反應譜；T為結構的自振周期；t為任意時間點；SaT(T,t)為t時刻的目標加速度反應譜。從式(1)可知，在tt和SaT(T,t)確定的條件下，采用耐震時程法得到的加速度時程曲線中，從起始時間點到t時刻的加速度反應譜值與時間t成正比。因為位移反應與加速度反應存在相關性，故耐震時程法加速度時程的目標位移反應譜可以表示為

(2)

式中，SuT(T,t)為t時刻時的目標位移反應譜。在某一固定時刻t，同時滿足式(1)和式(2)的加速度時程在一定的精度許可范圍下是可以獲得的；為了使在加速度時程上的任何點均滿足式(1)和式(2)則需要求解滿足式(3)的無約束變量的優化問題[17]

α[Su(T,t)-SuT(T,t)]2}dtdT

(3)

1.2 基于耐震時程法的極限抗震能力分析

在對特高拱壩進行非線性耐震時程分析后，為了方便對比分析常規時程分析法和耐震時程法的計算結果，需要將天然地震動的強度換算成耐震時程法對應的時間，地震動的等效耐震時程換算關系可以采用式(4)

(4)

同理，耐震時程曲線上各個時間點與天然地震動強度換算關系可采用式(5)

(5)

式中：tET為不同強度地震動的等效耐震時間；tTarget為耐震時程曲線總時間；S1為地震動的調幅參數；PGAC為規范中時程分析所用地震加速度時程的最大值。

通過式(4)和(5)可以從耐震時程分析計算結果中的到在各級地震動強度下壩體及基巖的動力響應。

損傷體積比(DVR)能夠定量的反映大壩的損傷程度,損傷體積比采用下式計算[18]

(6)

式中：n為壩體單元數目；i為壩體單元編號；di(t)為i號單元在t時刻對應的損傷因子；Vi為i號單元的體積。損傷體積比的范圍為0～1，損傷體積比越大，表示壩體損傷情況越嚴重。當損傷體積比為0時，表示壩體無損傷，當損傷體積比為1時，表示壩體所有單元已完全損傷。

高拱壩在地震損傷的過程中常常伴隨著損傷耗能的增加，因此損傷耗能指標也可以在一定程度上反映高拱壩動力損傷情況，損傷耗能采用下式計算

(7)

式中：MD為壩體損傷耗能；σ為壩體單元應力；εck為開裂應變；V為單元體積。

本文綜合考慮在各級地震強度下拱冠梁的最大相對位移、橫縫最大開度、壩體損傷程度、損傷體積比、損傷耗能等因素，進而確定大壩的極限抗震能力。

2 特高拱壩有限元模型及耐震時程曲線的生成

2.1 工程概況及有限元模型

以中國西南某300 m級特高拱壩為例，建立能夠同時考慮高拱壩—基巖相互作用的非線性有限元分析模型(圖1)。壩址區地震基本烈度為Ⅷ度，特征周期Tg=0.35 s。該特高拱壩最大壩高285.5 m，拱冠梁壩底處厚60 m，壩頂處厚14 m。壩頂中心線弦長602.2 m，弦高比為2.17。正常蓄水位時水深為275.5 m，相應的下游水深為53.5 m。壩體共設置橫縫30條。模型共劃分單元31 453個，自由度數為115 149，橫縫30條。由于特高拱壩計算模型規模較大，全部采用精細網格劃分會大量提高計算時間，但網格過大會造成計算精度下降。為了平衡這個問題，本文在壩體和近壩部分基巖采用精細網格劃分，遠端基巖采用較為稀疏的網格劃分。

圖1 特高拱壩非線性有限元模型Fig.1 Nonlinear finite element model of ultra-high arch dam

2.2 材料參數及荷載

壩體混凝土及基巖具體材料參數如表1所示，按照規范規定，混凝土的動態彈性模量較靜態彈性模量提高50%，強度參數較靜態提高20%。基巖材料彈性模量不提高，其損傷本構關系，目前尚缺乏充分研究。本文采用文獻[19]提供的方法，將基巖的材料參數折算為混凝土的強度參數，采用下式計算

(8)

表1 壩體及基巖材料參數表Tab.1 Material parameter of dam body and bedrock

式中：ft為折算后基巖的抗拉強度；c和φ分別為地基巖體的黏聚力和內摩擦角。本工程壩基巖體c取2.5 MPa，φ取53.47°，計算得到折算后的抗拉強度ft為1.65 MPa。

本文采用Abaqus中的混凝土塑性損傷模型(CDP模型)來模擬壩體混凝土及基巖的損傷行為。混凝土受拉應力應變關系曲線是基于GB50010—2010《混凝土結構設計規范》推薦的混凝土應力應變關系得出的。損傷因子的計算是基于能量損失理論，將脆性固體材料的損傷因子D定義為[20]

(9)

式中：W0為無損材料的應變能密度；Wε為損傷材料的應變能密度；E0和E分別為無損材料和損傷材料的四階彈性系數張量；ε為相應的二階應變張量。

混凝土及基巖的受拉應力應變關系曲線如圖2所示，相對應的損傷因子與應變關系曲線如圖3所示。

圖2 混凝土及基巖受拉應力應變曲線Fig.2 Tensile stress-strain curves of concrete and bedrock

圖3 混凝土及基巖受拉損傷因子與應變關系曲線Fig.3 Relationship between tensile damage factor and strain of concrete and bedrock

施加的靜力荷載主要有壩體自重、上下游的靜水壓力，淤沙壓力以及溫度荷載。正常蓄水位工況下，上游水位高程為600 m，淤沙高程為490 m。

對于高拱壩而言，采用修正的Westergaard法模擬壩體與庫水的動力相互作用[21]

(10)

式中：p為水深h處的動水壓力；a為水深h處的加速度；H0為庫水深度；h為計算點的水深；ρ為水體密度。

壩體采用阻尼比為5%的瑞利阻尼計算。壩基四周設置黏彈性邊界來模擬地基輻射阻尼，通過在邊界施加彈簧和阻尼單元來實現。采用波動輸入法，按文獻[22]提出的計算公式，通過編制將底邊界的入射波場和側邊界的自由波場轉化為對應的解析應力場的程序，對Abaqus進行二次開發，通過在邊界施加等效節點力的形式實現地震動輸入。

2.3 耐震時程曲線的生成

基于2.1節的式(3)，本文采用NB 35047—2015《水電工程水工建筑物抗震設計規范》[23]中的標準反應譜為目標譜，取特征周期Tg=0.35 s，阻尼比為5%，反應譜最大代表值βmax=2.50。參照規范第5.5.8條的要求：“采用時程分析法計算地震作用效應時，應以5%的設計反應譜為目標譜，生成至少3套人工模擬地震加速度時程作為基巖的輸入地震動加速度時程”。為了避免偶然性，本文采用非線性最小二乘法生成3條(ETA-1～ETA-3)時長為40 s的耐震時程曲線。如圖4(a)為第一條耐震時程曲線ETA-1的加速度時程曲線，從圖中可以看出，隨著時間的增加，地震動強度逐漸增大。圖4(b)給出了40 s的耐震時程曲線在不同時間段的反應譜，分別為：0～10 s，0～20 s，0～30 s，0～40 s。可以看出4個時間段的反應譜能與目標譜能夠很好的吻合，類似的，圖5(a)、(b)以及圖6(a)、(b)分別表示第二條(ETA-2)和第三條(ETA-3)耐震時程曲線及對應的反應譜。經過計算，三套耐震時程曲線的平均誤差值在20%以內，擬合誤差較小。

(a) 耐震時程曲線(ETA-1)

(b) 加速度反應譜(ETA-1)圖4 耐震時程曲線及對應的加速度反應譜(ETA-1)Fig.4 Endurance time accelerogram and corresponding acceleration response spectrum(ETA-1)

(a) 耐震時程曲線(ETA-2)

(b) 加速度反應譜(ETA-2)圖5 耐震時程曲線及對應的加速度反應譜(ETA-2)Fig.5 Endurance time accelerogram and corresponding acceleration response spectrum(ETA-2)

(a) 耐震時程曲線(ETA-3)

(b) 加速度反應譜(ETA-3)圖6 耐震時程曲線及對應的加速度反應譜(ETA-3)Fig.6 Endurance time accelerogram and corresponding acceleration response spectrum(ETA-3)

3 計算結果分析

本小節同時采用耐震時程法和常規時程分析法對高拱壩進行非線性動力計算。首先，基于標準反應譜生成3組人工地震波(GM-1～GM-3)，以峰值加速度(PGA)作為強度指標，將每條地震動在0.1g～1.0g之間進行調幅計算，間隔為0.1g，共進行30次非線性動力分析。其次，采用2.3節生成的三條耐震時程曲線，進行三次耐震時程分析。為了全面反映高拱壩在地震過程中的動力響應，本小節采用壩體最大相對位移、壩體動態損傷分布、壩體損傷耗能、橫縫最大開度等指標來對高拱壩的極限抗震能力進行評價。

3.1 壩體最大相對位移及橫縫開度分析

拱冠梁最大位移和橫縫最大開度是評價特高拱壩地震易損性的重要動力響應參數。為了方便分析對比耐震時程法和常規時程分析法的計算結果，對兩種方法得到的最大相對位移及橫縫最大開度計算結果進行匯總，如圖7和8所示。

圖7 拱冠梁最大相對位移Fig.7 Maximum relative displacement of crown cantilever

圖8 橫縫最大開度Fig.8 Maximum opening of contraction joint

從圖7結果可以看出，由于選取的是從起始時刻至該時刻間最大的動力響應。因此得到的數值曲線隨時間呈逐漸增大的趨勢。從整體上來看，耐震時程法計算結果與常規時程分析法結果較為接近，三條ETA曲線的位移值集中分布在常規時程分析法結果的均值附近，結果吻合較好。

從圖8的結果可以看出，對于壩體橫縫最大開度而言，一般而言超過50 mm就很可能會造成止水撕裂，威脅大壩的安全。采用ETA法計算得到的壩體橫縫最大開度整體上處于時程分析法的均值水平，耐震時程法和常規時程分析法的計算結果具有很好的一致性。當ETA時間在0～30 s時，對應最大PGA約為0.75g，三條ETA曲線值均小于50 mm，說明在罕遇地震作用下，橫縫不會發生較大的張開而導致止水被破壞。

3.2 壩體—基巖損傷分析

為了能夠對大壩的損傷情況進行定量的分析，本小節通過采用損傷體積比(DVR)和損傷耗能這兩個指標來反映大壩的損傷程度。繪制它們隨時間變化曲線。從圖9、10中可以看出，兩種方法得到的壩體損傷體積比很接近，結果離散性較小。從損傷耗能曲線上來看(圖10)，也可以得到相似的規律。

圖9 壩體受拉損傷體積比Fig.9 Tensile damage volume ratio of dam

圖10 壩體損傷耗散能Fig.10 Damage dissipation energy of dam

由于三條耐震時程計算得到的壩體損傷分布規律非常相近，本小節以其中一條耐震時程曲線的損傷計算結果為例，分析特高拱壩的動力損傷規律。

當ETA曲線處在0～10 s時，對應的最大峰值加速度約為0.25g。從圖11(a)～(c)可以發現，此時壩體主要處于線彈性狀態，在壩體上下游面及拱冠梁均未出現明顯的損傷區域。由于考慮了基巖的損傷，使得壩基附近壩體的應力得到釋放，因此在拱冠梁的壩踵壩趾處未出現受拉損傷。

(a) 拱冠梁(t=10 s)

(b) 下游壩面(t=10 s)

(c) 上游壩面(t=10 s)

(d) 拱冠梁(t=20 s)

(e) 下游壩面(t=20 s)

(f) 上游壩面(t=20 s)

(g) 拱冠梁(t=30 s)

(h) 下游壩面(t=30 s)

(i) 上游壩面(t=30 s)

(j) 拱冠梁(t=40 s)

(k) 下游壩面(t=40 s)

(l) 上游壩面(t=40 s)圖11 各時刻壩體受拉損傷分布圖Fig.11 Tensile damage distribution of dam at each time

當ETA曲線處在0～20 s時，對應的最大峰值加速度約為0.5g。此時壩體下游面首先出現損傷開裂，從拱冠梁的損傷分布圖中可以看出，損傷區域主要集中出現在壩體的中上部位，并從下游面逐漸向上游面擴展。壩體上游面整體上未出現明顯的損傷區域，兩側岸坡壩段的壩踵處開始出現輕微的損傷，河床壩段的壩踵處完好，未出現明顯的損傷區域，但壩基損傷區域進一步擴大，并向橫向發展。

當ETA曲線處在0～30 s時，對應的最大峰值加速度約為0.75g。此時的地震動強度已經達到了罕遇地震以上，從損傷分布圖來看，壩體下游面的損傷區域進一步擴大，特別是中上部區域損傷較為嚴重。從拱冠梁及上游面的損傷分布來看，此時損傷已經從下游貫穿到上游，并在上游面的中上部位形成了明顯的損傷區域。兩岸岸坡壩段壩踵處的損傷區域進一步擴大，河床壩段壩踵處未出現明顯的損傷。但壩基出現了較大范圍的損傷，可能會危及帷幕的安全，會給壩體的穩定性造成安全隱患。

當ETA曲線處在0～40 s時，對應的最大峰值加速度約為1.0g。從圖中可以看出，此時大部分壩體均處于嚴重損傷狀態，損傷范圍已經貫穿壩體，但損傷未必一定會造成壩體開裂，也有學者提出損傷因子超過0.7時可以認為混凝土發生宏觀裂縫，該判據也可作為判斷混凝土開裂范圍的參考。考慮到此時大壩的損傷程度，可以判定大壩已失去承載能力。

需要說明的是，單元尺寸的大小對于損傷結果具有較大的影響。本研究同時采用耐震時程法和常規方法來對比分析特高拱壩的極限抗震能力，由于兩種方法采用的是同一套網格，因此網格尺寸不影響兩者的對比分析。

3.3 大壩極限抗震能力分析

為了全面的評價大壩極限抗震能力，需要綜合考慮壩體的最大位移、橫縫的最大開度、壩體及基巖的損傷程度等因素。從3.1～3.2節分析可知，采用耐震時程法得到的壩體—基巖動力響應與采用常規的時程分析法得到的結果非常相近。耐震時程法能夠全面的反映出大壩和基巖從完好到逐漸破壞直至最終失效的過程。在耐震時程分析的前10 s時間段，僅基巖出現輕微損傷，壩體完好，此時大壩處于安全狀態。當處在第20 s時刻時，盡管壩體的局部區域發生了損傷，但仍未貫穿壩體，橫縫最大開度仍達不到破壞止水的程度，此時大壩仍能夠承載。當處在30 s的時刻時，壩體的損傷區域進一步擴大，并且貫穿了拱冠梁的中上部，基巖損傷區域進一步橫向發展，可能會破壞帷幕進而威脅大壩的穩定性。此時大壩承受的最大峰值加速度PGA為0.75g，綜合上面的分析，可以判定此時大壩處于極限狀態。

4 結 論

本文基于耐震時程法來研究特高拱壩—壩基的動力損傷行為，并進行極限抗震能力分析。計算中綜合考慮了壩體混凝土及壩基巖體的損傷等材料非線性問題以及橫縫開合等接觸非線性問題，得到以下結論：

(1) 耐震時程法具有較高的計算效率，計算得到的壩體位移、損傷、橫縫最大開度整體上與常規時程分析法接近，且計算結果離散性小。

(2) 耐震時程法能夠全面的反映出大壩和基巖從完好到逐漸破壞直至最終失效的過程。通過一次耐震時程計算，即可大致判斷出大壩的極限抗震能力，可以大大簡化極限抗震分析的計算工作量。

(3) 采用耐震時程法對某300 m級特高拱壩進行極限抗震能力分析，綜合位移、橫縫開度、壩體及基巖的損傷情況等方面，確定其極限抗震能力為0.75g，為相關工程的極限抗震能力分析提供一種新思路。