潘子悅 陳燈紅 趙藝園 劉云輝

(1.防災減災湖北省重點實驗室(三峽大學),湖北 宜昌 443002;2.三峽大學 土木與建筑學院,湖北 宜昌 443002)

目前,我國一批特大型混凝土重力壩擬建、正在建設或已建成,這些高混凝土重力壩工程,除了具有規(guī)模宏大、效益顯著的特點外,另一重要的共同點是它們均處于地震活動頻繁的西南地區(qū),大壩抗震設防水平高,抗震設計難度大,國內外少有先例.

國內外對于高壩抗震的研究主要有四個方面[1-3],包括壩址地震動的確定、筑壩材料動態(tài)力學特性和破壞機理研究、大壩地震響應分析以及大壩動力安全評價.其中,大壩-地基-庫水系統(tǒng)的地震響應分析一直是國內外學者研究的熱點.在數(shù)值建模中,該系統(tǒng)分為近場部分和遠場部分,其中如何精確描述遠場截斷邊界的輻射條件是關鍵.

在大壩-地基動力相互作用方面,現(xiàn)有求解方法可分為兩類,即局部的近似解法和全局的精確解法[4].近似方法包括黏性邊界、黏彈性邊界、多次透射邊界、無限元及其他的高階邊界等,精確方法有邊界元法、比例邊界有限元法等.目前,用于混凝土壩-庫水-地基系統(tǒng)地震響應分析中的主要方法有黏性邊界、黏彈性邊界和多次透射邊界.Chopra[5]采用黏性邊界計算了二維、三維混凝土壩-庫水-地基系統(tǒng)的地震響應.杜修力等[6]將黏彈性人工邊界結合顯式有限元的時域波動求解方法用于小灣拱壩-地基開放系統(tǒng)的地震響應分析.劉云賀等[7]建立了拱壩地震自由場輸入模型,比較、驗證了黏性邊界和黏彈性邊界的吸能效果.張楚漢等[8]考慮了無限地基輻射阻尼、混凝土損傷開裂非線性和橫縫接觸非線性等關鍵影響因素,研究了大崗山拱壩的地震響應.何建濤等[9]綜合考慮壩體混凝土和壩基巖體的材料非線性、無限地基輻射阻尼等因素,對Koyna重力壩進行了地震響應分析.王海波等[10]采用全壩段三維有限元模型、計入地基輻射阻尼以及壩段間動態(tài)接觸這些因素,研究了其對混凝土重力壩地震響應的影響.

在大壩-庫水動力相互作用分析方面,自Westergaard開創(chuàng)性研究以來,國內外研究者對于庫水可壓縮性、庫底邊界的吸收作用等方面的研究一直不斷.林皋等[11]推導了綜合考慮庫水可壓縮性和庫底吸收邊界的壩面動水壓力方程,提出了一種求解壩面動水壓力的半解析方法.王翔等[12]基于標量波動方程建立了求解混凝土壩動水壓力波的高階雙漸近透射邊界,發(fā)展了混凝土壩-庫水動力相互作用時域分析的耦合模型.高毅超等[13]將高階雙漸近透射邊界直接嵌入到近場有限元方程中,建立了大壩-庫水動力相互作用的直接耦合分析模型.王俊等[14]通過對比流固耦合模型與附加質量模型的壩體損傷來分析結果,發(fā)現(xiàn)流固耦合模型分析大壩-水庫系統(tǒng)相互作用更貼近實際情況.陳燈紅等[15]研究了庫水壓縮性及分縫布置對高拱壩非線性地震響應的影響.劉明志等[16]探討了庫水可壓縮性對重力壩動力特性和地震響應的影響.

在這些研究中,學者們建立了多種模型來模擬大壩-地基、大壩-庫水及大壩-庫水-地基耦合系統(tǒng)的動力相互作用,但在有限元數(shù)值分析中,對于基礎與庫水的模擬方法以及基礎與庫水截斷邊界的范圍尚未達成一致,大多數(shù)情況據(jù)經驗選取.本文將以黃登重力壩12號壩段為研究對象,考慮不同地基尺寸、地基輻射阻尼和庫水可壓縮性這些因素,分析大壩地震動力響應作用,并分別與無質量地基及附加質量模型的分析結果進行對比,以此來分析這些因素對混凝土重力壩地震響應的定量影響,為大壩抗震設計提供參考.

1 大壩-地基-庫水系統(tǒng)模型建立

1.1 大壩-地基系統(tǒng)相互作用模型

在大壩-地基相互作用的研究中,工程上最常用的是無質量地基模型,即只截取結構周圍一定范圍內的地基,并假定地基是線彈性、無質量的.但實際上地基是有質量的半無限介質,若使用無質量地基模型則無法考慮遠域地基的輻射阻尼效應.為了模擬遠域地基的輻射阻尼效應,更加真實地反映壩體地震響應,需要在其邊界處施加人為處理的邊界條件.黏彈性邊界是由Deeks、Randolph[17]在黏性邊界基礎上提出的,其基本思想是在截斷邊界上設置并聯(lián)的彈簧-阻尼器物理元件,劉晶波等[18]基于球面波動方程將黏彈性邊界推廣應用到三維情況.黏彈性邊界具有能同時模擬散射波輻射和半無限地基的彈性恢復能力的優(yōu)點,且能克服黏性邊界引起的低頻漂移問題,穩(wěn)定性好,并且物理意義明確,容易與有限元程序結合,因此在實際工程中得到了一些應用.

在黏彈性邊界的施加過程中,對于人工邊界上的彈簧、阻尼器的彈簧剛度和阻尼系數(shù),不同學者給出了相應建議.本文二維黏彈性人工邊界等效物理系統(tǒng)的法向與切向彈簧系數(shù)、阻尼系數(shù)[19]分別取為：

法向邊界：

切向邊界：

式中：KBT、KBN分別為彈簧切向與法向剛度系數(shù)；R為近場結構幾何中心至人工邊界點的距離；cs、cp分別為SV 波和P波波速；G 為介質剪切模量；ρ為介質密度；αN、αT分別為黏彈性人工邊界的法向、切向參數(shù),本文αN取為1.0,αT取為0.5.

在施加二維人工邊界的彈簧-阻尼器元件時,只需將式(1)和(2)中的系數(shù)乘以邊界上結點影響面積施加在截斷邊界結點上即可,其二維示意圖如圖1所示.

圖1 二維黏彈性邊界物理意義示意圖

對于黏彈性人工邊界上的波動輸入問題,采用將波動散射問題轉化為波源問題的方法來實現(xiàn),即通過在人工邊界節(jié)點上施加等效荷載來實現(xiàn),施加在黏彈性邊界節(jié)點的切向與法向等效荷載[19]可表示為：

式中：n為人工邊界外法線方向余弦向量；K B和C B分別為以黏彈性人工邊界彈簧系數(shù)、阻尼系數(shù)為元素的對角矩陣分別為自由場位移、速度和應力矢量.

1.2 大壩-庫水相互作用模型

水體對壩體的作用可以分為兩種,一是靜水壓力,二是地震作用下的慣性力,大壩在地震荷載的持續(xù)作用下,水體會隨著大壩不斷地做往復運動,反之水體產生的慣性力又反作用于壩體上,即為動水壓力.在水體對壩體表面產生的動水壓力的研究中,常有兩種計算方法：

1)附加質量法

該方法由Westergaard于1933 年提出,假定壩面是直立剛性的,庫水上游方向無窮遠,庫底為剛性水平面.該方法可忽略庫水的可壓縮性以及壩體的變形因素,其表達式為：

式中：p為壩面某點受到的動水壓力；ρ為庫水質量密度；H為壩前庫水深度；Z為該點在壩基面以上的高 度為壩面結點加速度.

2)流固耦合法

該方法首先將水體假設為均勻、可壓縮的、無旋度無黏性的理想流體,基于這種假定,可以得到以壓力p為未知量的波動方程[20]：

其中：▽為拉普拉斯算子；p為動水壓力；c為壓聲波波速；K為流體體積模量；ρ為流體密度.

在設定邊界條件后,考慮結構位移和流體速度勢組成的混合未知量的有限元方程[20]為：

式中：M、C、K分別為質量矩陣、阻尼矩陣和剛度矩陣；Q為流固交界面上的耦合矩陣；u和p分別為位移和動水壓力；f為外力矢量；r為近場水作用于水庫的相互作用力矢量.下標s為壩體上的自由度；下標f表示近場水域的自由度；下標r表示截斷邊界上的自由度.

2 算例驗證

以國際大壩委員會發(fā)布的第十五屆大壩數(shù)值分析基準算例研討會主題A 中的Pine Flat重力壩地震響應分析為例[21],驗證構建的重力壩-地基-庫水系統(tǒng)地震響應分析模型的正確性.混凝土與地基的彈性模量E=22.41 GPa,泊松比ν=0.2,密度ρ=2 483 kg/m3；水體體積模量K=2.07 GPa,密度ρw=1 000 kg/m3.采用基于有限元法的黏彈性邊界模擬無限壩基,采用等效荷載法實現(xiàn)地震動輸入,采用聲學單元及無反射邊界條件模擬壩-庫動力相互作用[20].有限元網格采用四結點四邊形平面應變單元離散,共劃分了6 417個單元6 720個結點,計算網格如圖2所示,輸入的Taft地震波如圖3所示.

圖2 Pine Flat重力壩-地基-庫水系統(tǒng)有限元網格圖

為了說明構建耦合模型的優(yōu)越性,采用文獻[21]中基于比例邊界有限元法的計算模型結果進行了對比,其中采用聲學單元及基于標量波波動方程的高階雙漸近透射邊界模擬壩-庫動力相互作用,采用位移單位脈沖響應函數(shù)算法及相適配的地震動輸入方法模擬壩-基動力相互作用.壩頂水平向位移、加速度結果及相對水平位移結果的比較如圖4～6所示,結果表明本文構建的耦合模型結果與基于比例邊界有限元法模型的結果[21]相比具有良好的計算精度和效率.

圖5 壩頂水平加速度結果的比較

圖6 相對水平位移結果的比較

3 工程應用

3.1 工程概況

黃登水電站位于云南省蘭坪縣境內,上游與托巴水電站,下游與大華橋水電站相銜接.壩址控制流域面積9.19×104km2,總庫容約為1.5×109m3,多年平均流量901 m3/s.黃登水利工程為混凝土重力壩,最大壩高203 m,屬Ⅰ等大(1)型工程.本次模擬選用12號擋水壩段進行數(shù)值模擬計算,壩頂高程1625m,壩基面高程1422 m,壩高203 m,壩頂寬16 m,正常蓄水位1 619 m,如圖7所示.

圖7 12號擋水壩段斷面(單位：m)

計算內容共分為A、B、C 3種工況,分別考慮地基尺寸大小、地基輻射阻尼以及庫水可壓縮性等各項影響因素探究其對重力壩動力響應的影響.并分別以A、B、C表示這3種工況,各工況匯總見表1.其中H表示大壩高度,L1表示地基范圍從壩踵往上游方向延伸尺寸,L2表示從壩趾往下游方向延伸尺寸,H0表示地基深度,L表示庫水長度.

表1 計算工況匯總表

3.2 計算模型

計算采用平面有限元分析方法,基于Abaqus有限元仿真軟件實現(xiàn).計算中考慮了壩體自重及以正常蓄水位下的靜水壓力等荷載,地震荷載選用Koyna地震動.壩體采用平面應力單元,地基采用平面應變單元,網格尺寸均為2 m×2 m.當動水壓力以Westergaard公式計算的附加質量模擬,施加在壩體的上游迎水面時,建立有限元模型如圖8所示；當庫水動水壓力采用聲學單元模擬,施加在壩體的上游迎水面時,其有限元模型如圖9所示.工況A 和工況C 采用無質量地基模型,邊界條件為：對側邊界施加法向約束,底部邊界全部約束；工況B采用有質量地基模型,截斷人工邊界采用黏彈性邊界來模擬無限地基的輻射阻尼效應,其中,在Abaqus進行結構整體靜動力分析中,通過一種單元“生死”技術[22],實現(xiàn)靜態(tài)約束邊界到動態(tài)黏彈性邊界的穩(wěn)定切換,從而實現(xiàn)靜動力統(tǒng)一計算.

圖8 工況C1計算模型

圖9 工況C2計算模型

3.3 模型參數(shù)

壩體不模擬分區(qū),只采用一種混凝土材料,并假定為線彈性,動彈模取為靜彈模的1.5倍[23]；地基巖體假定為均質線彈性介質,且忽略材料阻尼.具體材料參數(shù)見表2.

表2 材料參數(shù)

此外,在實際的動態(tài)分析中還需考慮混凝土的阻尼,本文采用瑞利阻尼的標準形式來計算,即兩參數(shù)計算模型如下：

其中：ω1、ω2分別為系統(tǒng)的前兩階自振頻率；ξ為阻尼系數(shù),本文取ξ=5%[24].

3.4 計算荷載

本工況中選取黃登混凝土重力壩12 號擋水壩段,考慮上游正常蓄水位情況進行數(shù)值模擬計算,作用荷載包括自重、靜水壓力、動水壓力和地震荷載.

地震荷載的計算采用了Koyna地震時程,地震總歷時為10 s,選擇固定時間步長Δt=0.01 s進行計算.水平和垂直加速度時程如圖10所示.

圖10 輸入的地震加速度時程

4 地基尺寸的敏感度分析

工況A 為探究地基尺寸對重力壩動力響應的影響,地基范圍為從壩踵往上游方向、從壩趾往下游方向、壩基往深度方向均分別延伸1.5倍、2倍、3倍壩高,并分別以A1、A2、A3表示這3種工況.動水壓力以Westergaard 公式計算的動水附加質量模擬,施加在壩體的上游迎水面.

4.1 位移分析

上述模型的壩頂、壩踵的水平加速度峰值、水平位移峰值匯總見表3,壩踵點與壩頂點的水平相對位移時程如圖11所示.

表3 工況A壩頂、壩踵水平加速度與位移峰值

圖11 工況A 相對水平位移時程曲線

由表3可以看出,地基尺寸對大壩加速度響應有顯著影響,壩頂加速度在2 倍壩高時達到最大,為26.53 m/s2,而地基尺寸對位移影響較小,當?shù)鼗秶謩e延伸1.5倍、2倍、3倍壩高時,壩踵點與壩頂點的水平相對位移峰值呈現(xiàn)遞增規(guī)律,3倍壩高時水平相對位移峰值較1.5倍壩高增加了8.7%.

4.2 應力分析

選取壩體壩踵和壩趾部位作為特征點,各工況下的主拉應力和主壓應力結果如圖12～13所示.各工況的最大主拉應力均發(fā)生在壩踵處,拉應力峰值分別為14.79、13.48、10.31 MPa；最大主壓應力發(fā)生在壩趾處,壓應力峰值分別為-31.23、-29.15、-28.71 MPa.由此可見,在一定范圍內,隨著地基尺寸的增大,壩體的主拉應力、主壓應力峰值均呈現(xiàn)遞減趨勢,地基尺寸為3倍壩高的最大拉應力峰值較1.5倍壩高減小了30.3%.

圖12 工況A 壩踵主拉應力時程曲線

圖13 工況A 壩趾主壓應力時程曲線

5 地基輻射阻尼對重力壩地震響應的影響

在考慮地基質量的情況下,本節(jié)探究地基輻射阻尼對重力壩動力響應的影響,截斷地基人工邊界選用黏彈性人工邊界,動水壓力采用附加質量模型,采用和工況A 相同的地基尺寸,即從壩踵往上游方向、從壩趾往下游方向、壩基往深度方向均分別延伸1.5倍、2倍、3倍壩高,并分別以工況B1、B2、B3表示.

5.1 位移分析

壩頂、壩踵的水平加速度、水平位移峰值匯總見表4,壩踵點與壩頂點的相對位移時程如圖14所示.

圖14 工況B相對水平位移時程曲線

表4 工況B壩頂、壩踵水平加速度與位移峰值

對比表4與表3可以看出,黏彈性邊界模型考慮了地基的輻射阻尼效應后較無質量地基模型,壩頂或壩踵的加速度峰值降低了6%～63%,壩踵及壩頂水平位移峰值降低了60%左右,相對位移峰值降低了50%左右；并且,隨著地基尺寸的增加,壩踵及壩頂點相對位移峰值呈現(xiàn)上升趨勢.由此可見,無質量地基模型夸大了壩體的地震響應,考慮地基輻射阻尼對大壩進行動力響應分析有著重要的影響.

5.2 應力分析

各工況下的主拉應力和主壓應力結果如圖15～16所示.當?shù)鼗秶謩e延伸1.5倍、2倍、3倍壩高時,壩踵的主拉應力峰值分別為6.33、5.63、4.97 MPa；壩趾的主壓應力峰值分別為-18.9、-19.29、-19.37 MPa,在一定范圍內,隨著地基尺寸的增大,拉應力峰值則呈遞減趨勢.同時,與無質量地基模型的結果比較,考慮黏彈性邊界的壩體主拉、壓應力峰值也出現(xiàn)不同程度的減小.最大主拉應力發(fā)生在壩踵點,最大主壓應力發(fā)生在壩趾,與無質量地基模型的結果相比,主拉應力峰值降低了50%～60%,主壓應力峰值降低了30%～40%.

圖15 工況B壩踵主拉應力時程曲線

圖16 工況B壩趾主壓應力時程曲線

6 庫水壓縮性對重力壩地震響應的影響

本節(jié)分析庫水的壓縮性對重力壩動力響應的影響.不可壓縮性庫水采用Westergaard 附加質量模型；可壓縮庫水可用有限元模擬,在庫水面上采用自由表面邊界條件,在庫水遠端采用無反射邊界條件,在庫水-大壩、庫水-地基交界處采用流固耦合邊界條件[20].

工況C采用無質量地基模型,地基尺寸從壩趾往下游、深度均延伸1.0倍壩高,從壩踵往上游分別延伸3、4、5倍壩高.在工況C1中,動水壓力以Westergaard 公式計算的附加質量單元模擬；工況C2、C3、C4中,動水壓力采用聲學單元模擬.

6.1 位移分析

表5為4種工況下的壩體動力響應極值匯總,圖17為壩踵與壩頂點相對位移時程曲線.對于相同地基尺寸而言,可壓縮庫水模型較Westergaard附加質量模型的水平相對位移降低了13.8%；當庫水長度分別為3、4及5倍壩高時,水平相對位移峰值呈現(xiàn)上升趨勢.

表5 工況C壩頂、壩踵的水平加速度與位移峰值

圖17 工況C相對水平位移時程曲線

6.2 應力分析

選取壩體壩踵和壩趾部位作為特征點,不同工況下得到的壩體-庫水界面的主拉應力和主壓應力結果如圖18～19 所示.C 工況中主拉應力峰值分別為15.26、9.32、9.68、10.28 MPa；主壓應力峰值分別為-32.80、-23.25、-24.41、-25.50 MPa.由此可見,與附加質量模型結果相比,流固耦合模型的主拉應力以及主壓應力的峰值均有所降低,降幅在30%左右；當庫水長度分別為3、4及5倍壩高時,主拉應力峰值呈現(xiàn)遞增趨勢.

圖18 工況C壩踵主拉應力時程曲線

圖19 工況C壩趾主壓應力時程曲線

7 結論

本文基于構建的重力壩-地基-庫水系統(tǒng)模型以黃登混凝土重力壩12 號壩段為研究對象,分別以A～C 10種工況探究了地基尺寸、無限地基輻射阻尼效應及庫水可壓縮性這些重要因素影響的重力壩-地基-庫水系統(tǒng)地震響應.得出主要結論如下：

1)地基尺寸對大壩加速度響應有顯著影響,對位移影響較小,壩頂加速度在2倍壩高時達到最大,為26.53 m/s2,當?shù)鼗秶謩e延伸1.5、2、3 倍壩高時,壩踵與壩頂?shù)乃较鄬ξ灰品逯党尸F(xiàn)遞增規(guī)律,地基尺寸為3倍壩高時水平相對位移峰值較1.5倍壩高增加了8.7%,壩體的最大拉應力峰值呈現(xiàn)遞減趨勢,地基尺寸為3倍壩高的最大拉應力峰值較1.5倍壩高降低了30.3%.

2)黏彈性邊界模型考慮了地基的輻射阻尼效應后較之無質量地基模型,動力響應均有不同程度的降低,壩頂和壩踵的加速度峰值降低了6%～63%,位移峰值降低了60%左右,相對位移峰值降低了50%左右.同時,壩體主拉、壓應力峰值也出現(xiàn)不同程度的減小,拉應力峰值降低了50%～60%,壓應力峰值降低了30%～40%.由此可見,壩體與地基的動力相互作用不容忽視,應考慮輻射阻尼效應,黏彈性人工邊界模型比無質量地基模型更能有效地模擬遠場地基的輻射阻尼效應.

3)對比流固耦合模型與Westergaard附加質量模型發(fā)現(xiàn),附加質量模型會夸大動水壓力作用.當?shù)鼗叽缦嗤瑫r,流固耦合模型較附加質量模型的水平相對位移降低13.8%,壩體應力峰值降低30%左右.因此,對比流固耦合模型比較與傳統(tǒng)Westergaard附加質量模型,前者模擬壩體-庫水之間的相互作用更為合理.