考慮滲流體力的某大壩地震響應仿真分析

韓 彰,蘇懷智,李 慧

(1.河海大學水文水資源與水利工程科學國家重點實驗室,江蘇 南京 210098;2.河海大學水利水電學院,江蘇 南京 210098)

近年來,國家規劃、建造了一批土石壩[1-2],這些土石壩多處于高地震烈度區,需進行地震動力響應分析,以確保在建和已建土石壩的安全[3-5]。目前,土石壩安全及動力響應分析的研究手段主要有極限平衡法、物理模型試驗、數值仿真、剪切楔法及集中質量法等。沈振中等[6]推導了滿足力和力矩平衡的壩坡穩定極限平衡水平條分法的計算公式;劉小生等[7]開展了面板堆石壩振動模型試驗及動力分析研究;竇興旺等[8]開展了深覆蓋層上高堆石壩振動臺試驗與動力數值分析驗證研究。由于數值仿真能直接反映力學演化過程,且計算成本低,因此,數值仿真已被廣泛應用于土石壩的動力響應分析,如竇興旺等[9]研究了人工邊界法在土石壩動力分析中的應用;楊秀竹等[10]開展了地震力作用下軟基土石壩的動力響應分析;丁陸軍[11]基于有限元法和土體的等價線性分析,運用三維有限元動力分析程序對某面板堆石壩進行了動力分析,研究了該壩的整體受力情況; 曹學興等[12]改進了傳統的Hardin-Drnevich模型,提出了可以考慮土石料動力特性參數圍壓依賴性的改進動力本構模型;龐林祥等[2]提出了適用于高土石壩的地震動力響應分析計算方法;孔憲京等[13-14]、Piao等[15]、趙劍明等[16]、楊昕光等[17]都開展了相關研究。但上述研究主要集中在土石壩地形條件、動力本構模型、抗震性能、計算方法等方面,而考慮滲流體力的土石壩動力仿真研究國內外還鮮有報道。針對某水庫大壩混凝土連接壩段、均質壤土擋水壩段靜動力計算中須考慮滲透體力的問題,本文采用等效線性黏彈性模型,提出了一種在靜動力計算中滲流體力的施加方法,并采用MSC.Marc有限元軟件,通過Fortran語言編寫用戶子程序，對該壩地震響應過程開展了仿真分析。

1 工程概況

某水庫工程(從左至右)主要由混凝土擋水壩段、混凝土取水壩段、泄洪排沙底孔壩段、表孔溢流壩段、混凝土連接壩段、均質壤土擋水壩段等組成。混凝土擋水壩段為現澆C20混凝土重力壩,壩體斷面上游為垂直面,下游壩坡1∶0.8,分為2個小壩段,長度分別為26 m和33 m,總長59.0 m；壩頂高程1 068.50 m,壩基坐落在弱風化巖體上,壩基最低建基面高程1037.0 m,向左岸側逐漸抬高至1061.0 m,最大壩高31.5 m。混凝土連接壩段的連接采用插入式,連接壩段分為4個壩段,每段壩段長25 m,其中半插入段長75 m,插入段長25 m,總長100 m；壩基坐落在弱風化基巖上,建基面高程1029.00 m,最大壩高39.50 m。均質壤土擋水壩段上游壩坡自上而下分別為1∶3.0和1∶3.5,下游壩坡自上而下分別為1∶2.75和1∶3.0,上游壩坡采用現澆C20混凝土板護面,混凝土板厚200～300 mm，混凝土板下鋪設粒徑20～40 mm的砂礫石墊層,厚800 mm;下游壩坡高程1045.0 m以上坡面采用預制C20混凝土網格草皮護坡。該區地震震動峰值加速度0.10g,對應地震基本烈度為Ⅶ度。根據該水庫壩址工程場地地震安全性評價報告,確定該壩址場地50年內超越概率為10%的地震烈度值為7.23,超越概率為2%的地震烈度值為7.88,地震持續時間取30 s。

2 數值模擬方法

采用動力計算本構模型對該水庫均質壤土擋水壩段及其混凝土連接壩段非線性有限元動力特性進行數值模擬。

2.1 動力計算本構模型

a. 均質壤土擋水壩段壩體土料和棱體排水料的動力計算模型。地震荷載是一種非等幅等周期的不規則荷載,在一次地震中,土料將經歷數十次甚至上百次卸載和再加載的過程,并且它們之間無規律可循,為解決此問題,比較常用的方法是應用Masing規則,制定一個應力應變關系的骨架曲線,在此基礎上,可采用雙線性、黏彈性和彈塑性等多種本構模型[18-20]。由于等效線性黏彈性模型簡單,結果精確,動力計算分析采用等效線性黏彈性模型,即假定壩體土料為黏彈性體,采用等效剪切模量G和等效阻尼比λ這兩個參數來反映土體動應力應變關系的非線性和滯后性兩個基本特征,并表示為剪切模量和阻尼比與動剪應變之間的關系。

b. 混凝土連接壩段計算模型。混凝土連接壩段動力計算分析時采用線性彈性模型。

c. 接觸面的動力計算模型。接觸面單元的動力參數采用河海大學的試驗成果[18]。接觸面單元的動力模型剪切勁度Kc與動剪應變γd、阻尼比λc的關系分別為

(1)

(2)

其中Kcmax=Γσn0.7τf=σntanδ

式中:σn為接觸面單元的法向應力;δ為接觸面的摩擦角;λcmax為最大阻尼比;Β、Γ為試驗參數。

2.2 壩體結構網格、邊界條件與荷載

圖1 水庫壩體三維有限元模型

圖2 樁號位置示意圖

根據初步設計圖、地質勘測報告等資料,考慮壩體分區、施工程序及加載過程,對壩體進行剖分,建立三維有限元模型。網格采用六面體和四面體結構化網格,總共剖分6 560個單元,7 417個結點,其中混凝土連接壩段716個單元,均質壤土擋水壩段5 228個單元,排水棱體616個單元。壩體的三維有限元模型見圖1,樁號位置見圖2。動力計算中忽略壩體與基巖的動力相互作用,將壩體與基巖的接觸面設為三向約束。將重力荷載、地震荷載及正常蓄水位下穩定滲流場所形成的壩體滲透壓力施加于壩體。借助Abaqus商業軟件進行滲流計算,為了完成Abaqus與Marc的對接,編寫3個Fortran子程序處理滲流計算得到的流速,將流速轉化為坡降數據,作為滲流力施加的輸入數據,在Marc中以加體力邊界的方式完成滲流壓力的施加。編寫3個子程序的目的為:①找出Abaqus與Marc結點坐標對應關系(Marc退化單元篩分),進行結點坐標轉換以及Marc中8、6結點單元篩選;②找出Abaqus與Marc中單元對應關系,進行8、6、4結點單元轉換;③提取滲透流速,計算出滲透坡降,導出dat輸入文件。

壩體0+158斷面滲透壓力分布以及壩軸線斷面滲透壓力分布見圖3和圖4。蓄水期0+158斷面滲透壓力最大值為0.28 MPa,壩軸線斷面滲透壓力最大值為0.08 MPa。

圖3 壩體0+158斷面滲透壓力分布(單位:MPa)

圖4 壩軸線斷面滲透壓力分布(單位:MPa)

2.3 動力計算參數

通過試驗測得動剪切模量比Gd/Gdmax和動阻尼比λd與動剪應變γd的關系曲線。動力計算時輸入相應關系曲線的控制數據,根據應力應變值進行內插和外延取值。本工程壩料試驗結果如表1所示。壩料模量系數K′=696.0,模量指數n′=0.5,固結比kc=1.0。

混凝土連接壩段動力計算分析時采用的線性彈性模型,靜力彈性模量Es提高30%即為動力彈性模量Ed。接觸面單元的動力模型采用河海大學試驗成果得出的接觸動本構模型[18]。接觸面的摩擦角取32°，最大阻尼比取0.2，Β、Γ分別取2.0和22.0。

表1 不同動剪應變下壩料的動剪切模量比和動阻尼比

2.4 壩體動力非線性有限元計算步驟

在進行動力計算分析之前,必須先進行靜力計算分析,以獲得動力分析壩體的初始應力狀態。靜力分析方法比較簡單,這里不再贅述。計算中采用MSC.Marc有限元軟件進行建模與仿真,Marc具有較強的處理非線性問題的能力,并具有良好的二次開發接口,通過Fortran語言自編用戶子程序，對該壩地震響應過程開展了仿真分析，計算的主要步驟如下:

步驟2求出土體單元的初始動剪切模量Gdmax,0,土體單元的初始阻尼比根據經驗取為5%。

步驟3整個地震歷程劃分為若干個時段。

步驟4對每個時段的動剪切模量迭代求解。

步驟5用Willson-θ法建議的放大時間間隔ΔT=θΔt代替實際時間間隔Δt,對每個時段進行時程分析。

步驟6計算各單元的質量矩陣和剛度矩陣,形成總體質量矩陣M和剛度矩陣K,采用空間迭代法求出壩體基頻ω,并計算單元阻尼矩陣,由各單元的變阻尼矩陣ce組裝得到總體阻尼矩陣C。

步驟10根據求出的結點位移un+1計算各單元的動剪應變γn+1和動剪應力τn+1。

步驟11重復步驟 5～10,得到各單元在每個時段內的動剪應變γd時程。

步驟12求出各單元γd時程中的最大值γdmax，根據等效動剪應變γd eff=0.65γdmax,查Gd/Gdmax～γd和λc～γd曲線得到新的Gd和λc。

步驟13重復步驟4～12,直到前后兩次采用的Gd相對誤差小于10%。

步驟14重復步驟3～13,直到各個時段全部計算結束,即整個地震歷程結束。

步驟15輸出計算結果。

2.5 地震加速度的輸入

壩體的動力反應計算需考慮“正常蓄水位+地震”工況。首先進行靜力分析,并將水庫水位蓄至正常蓄水位,隨后施加地震荷載,忽略壩體與地基的相互動力作用,加速度直接施加于壩體與壩基接觸面上,進行地震反應分析。

根據水庫壩址工程場地地震安全性評價報告中的場地基巖水平加速度反應譜曲線(圖5),采用人工合成法合成50年超越概率為10%的設計地震加速度,其中地震的持續時間取30 s。地震波輸入方向為:x向沿原河流方向水平加速度輸入;y向沿高程方向豎直加速度輸入,依據NB 35047—2015《水電工程水工建筑物抗震設計規范》,將其峰值折減2/3;z向為沿壩軸方向橫向加速度輸入。圖6為人工合成的設計地震三向加速度曲線。計算中將整個地震歷程劃分為30個大時段,每個大時段又劃分為50個小時段,因此,積分計算的時間步長為0.02 s。

圖5 場地基巖水平加速度反應譜

圖6 地震加速度時程曲線

計算過程中,記錄了壩體0+158斷面12個結點(圖7)的加速度反應和動位移反應及8個單元(圖8)的最大、最小動主應力反應。

圖7 0+158斷面12個結點位置及編號

圖8 0+158斷面8個單元位置及編號

3 三維非線性動力計算結果與分析

3.1 設計地震反應

由于試驗條件不允許,混凝土連接壩段和均質壤土擋水壩段永久變形計算參數參考了積石峽混凝土面板堆石壩[21]數據和相關經驗。典型壩體斷面(0-32斷面)連接壩段最大加速度幅值和最大動位移等值線分布見圖9和圖10。

圖9 0-32斷面動位移等值線分布(單位:cm)

圖10 0-32斷面加速度等值線分布(單位:m/s2)

設計提供的模擬地震動加速度曲線歷時達30 s,因此,在整理成果時給出了30 s的時程曲線,設計地震工況基巖輸入加速度取50年超越概率為10%的動力有限元分析計算成果:①x向和y向最大加速度分別為5.0 m/s2和4.0 m/s2,分別出現在壩體最大斷面下游壩頂附近和0-81斷面上游壩頂附近,x向放大系數為7.70,y向放大系數為6.16;②x向和y向最大位移分別為9.0 cm和1.8 cm,分別出現在壩體最大斷面壩頂附近和0-81斷面上游壩頂附近;③土體第一和第三最大主應力都為200 kPa,分別出現在0-81斷面上游壩底和0-32斷面上游壩底;④x向和y向地震永久變形分別為0.7 cm和2.5 cm,x向出現在0-32和0-81斷面下游壩頂以及壩體最大斷面下游壩坡1/3壩高處附近，y向出現在壩體最大斷面壩頂附近,地震沉陷量為壩高的0.06%。

圖11 0-81斷面接觸面結點分布

由計算成果可知:①最大地震動加速度、動位移反應位于壩頂局部位置,壩頂存在明顯的鞭鞘效應,需要在壩頂進行適當的抗震加固;②計算得到的最大放大系數為7.70,這可能是由于大壩很長,約束又加在兩端,正如橡皮筋一樣,導致計算得到的反應加速度放大倍數偏大,動力響應偏大,建議在壩頂加土工格柵。

圖12 端頭接觸面結點分布

3.2 上下游接觸面及端頭接觸面典型結點動位移

圖11為0-81斷面接觸面結點分布,圖12是端頭接觸面結點分布,由于這些典型結點的x、y、z向動位移反應歷時曲線過多,這里不再附圖展示。經統計分析,0-81斷面接觸面各結點在30 s設計地震過程中三向動位移極值范圍見表2,端頭接觸面第1～4列各結點在30 s設計地震過程中的三向動位移范圍見表3。由表2和表3可見,三向動位移反應均比較小,并且對稱分布,因此動力工況下接觸面工作狀態基本正常。

表2 0-81斷面接觸面結點動位移統計

表3 端頭接觸面第1～4列結點動位移統計

3.3 端頭接觸面典型單元法向動應力

端頭接觸面地震歷程30 s末法向動應力分布如圖13所示,可見設計地震工況下,端頭接觸面法向動應力最大值為3.70 MPa,出現在接觸面底部。經統計,端頭接觸面法向動應力均為壓應力,因此混凝土連接壩段在設計地震工況下是穩定的。

圖13 端頭接觸面典型單元法向動應力分布(單位:MPa)

4 結 論

a. 該大壩的最大地震動加速度、動位移反應位于壩頂局部位置,壩頂存在明顯的鞭鞘效應,需要在壩頂進行適當的抗震加固。

b. 最大地震動應力反應位于混凝土連接壩段上游壩底附近位置;計算得到的最大放大系數為7.70。

c. 各分區的設計與填筑標準、壩體分層填筑方案合理,壩體抗震安全性較好。

[1] 吳中如.水工建筑物安全監控理論及其應用[M].北京:高等教育出版社,2003.

[2] 龐林祥,傅華.高土石壩的地震動力響應分析方法研究與應用[J].水利水電技術,2015(7):3-9.(PANG Linxiang,FU Hua.Study on method for analysis on seismic dynamic response of high earth-rockfill dam and application[J].Water Resources and Hydropower Engineering,2015(7):3-9.(in Chinese))

[3] 李陽,任亮,王瑞駿,等.基于ABAQUS的面板砂礫石壩三維動力有限元分析[J].水資源與水工程學報,2014(6):39-43.(LI Yang,REN Liang,WANG Ruijun,et al.Finite element analysis of 3D dynamics for concrete face rock fill dam based on ABAQUS[J].Journal of Water Resources & Water Engineering,2014(6):39-43.(in Chinese))

[4] 任亮.考慮流變特性的面板砂礫石壩靜動力有限元分析[D].西安:西安理工大學,2013.

[5] 李炎隆,徐嬌,李守義.瀝青混凝土心墻堆石壩動力有限元分析[J].應用力學學報,2016,33(2):247-253.(LI Yanlong,XU Jiao,LI Shouyi.Finite element analysis of dynamic on rockfill dam with asphaltic concrete core wall[J].Chinese Journal of Applied Mechanic,2014(6):39-43.(in Chinese))

[6] 沈振中,崔娟,任華春.土石壩抗震穩定性分析的極限平衡水平條分法[J].南水北調與水利科技,2011,9(3):145-149(SHEN Zhenzhong,CUI Juan,REN Huachun.Quasi-static limit equilibrium horizontal slice method forseismic stability analysis of earth-rock dam[J].South to North Water Diversion and Water Science & Technology,2011,9(3):145-149(in Chinese))

[7] 劉小生,王鐘寧,趙劍明,等.面板堆石壩振動模型試驗及動力分析研究[J].水利學報,2002,33(2):29-35.(LIU Xiaosheng,WANG Zhongning,ZHAO Jianming,et al.Advancement of technology on shaking table model test and dynamic analysis of CFRD[J].Journal of Hydraulic Engineering,2002,33(2):29-35.(in Chinese))

[8] 竇興旺,林義興,夏頌佑,等.深覆蓋層上高堆石壩振動臺試驗與動力數值分析驗證研究[J].水利水電科技進展,2000,20(6):19-22.(DOU Xingwang,LIN Yixing,XIA Songyou,et al.Shaking-table test and dynamic numerical analysis of high rookfill dams constructed on deep overburdens [J].Advances in Science and Technology of Water resources,2000,20(6):19-22.(in Chinese))

[9] 竇興旺,夏頌佑,許百立.人工邊界法在土石壩動力分析中的應用[J].河海大學學報(自然科學版),2000,28(5):72-75.(DOU Xingwang,XIA Songyou,XU Baili.Application of artificial boundary method to dynamic analysis of embankment dams[J].Journal of Hohai University(Natural Sciences),2000,28(5):72-75.(in Chinese))

[10] 楊秀竹,李花珍,王星華.地震力作用下軟基土石壩的動力響應分析[J].鐵道科學與工程學報,2008,5(3):51-53.(YANG Xiuzhu,LI Huazhen,WANG Xinghua.Dynamical response analysis of the earth-rock dam on soft soil during earthquake[J].Journal of Railway Science and Engineering,2008,5(3):51-53.(in Chinese))

[11] 丁陸軍.面板堆石壩三維動力反應分析[D].西安:西安理工大學,2005.

[12] 曹學興,何蘊龍,熊堃,等.考慮圍壓效應的高土石壩動力響應分析[J].天津大學學報(自然科學與工程技術版),2014,47(2):116-123.(CAO Xuexing,HE Yunlong,XIONG Kun,et al.Dynamic response of high rockfill dam considering confining pressure effect [J].Journal of Tianjin University(Science and Technology),2014,47(2):116-123.(in Chinese))

[13] 孔憲京,周揚,鄒德高,等.汶川地震紫坪鋪面板堆石壩地震波輸入研究[J].巖土力學,2012,33(7):2110-2116.(KONG Xianjing,ZHOU Yang,ZOU Degao,et al.Study of seismic wave input of Zipingpu concrete face rockfill dam during Wenchuan earthquake[J].Rock and Soil Mechanics,2012,33(7):2110-2116.(in Chinese))

[14] 孔憲京,鄒德高,鄧學晶.高土石壩綜合抗震措施及其效果的驗算[J].水利學報,2006,37(12):1489-1495.(KONG Xianjing,ZOU Degao,DENG Xuejing.Comprehensive earthquake resistant measure of high earth-rockfill dams and effectiveness verification[J].Journal of Hydraulic Engineering,2006,37(12):1489-1495.(in Chinese))

[15] PIAO R,RIPPE A,MYERS B.Earth dam liquefaction and deformation analysis using numerical modeling[J].Rock Mechanics and Geotechnical Engineering,2006,8(8):80-85.

[16] 趙劍明,溫彥鋒,劉小生,等.深厚覆蓋層上高土石壩極限抗震能力分析[J].巖土力學,2010,31(增刊 1):41-47.(ZHAO Jianming,WEN Yanfeng,LIU Xiaosheng,et al.Study of maximum aseismic capability of high earth-rock dam on deep riverbed alluviums[J].Rock and Soil Mechanics,2010,31(Sup1):41-47.(in Chinese))

[17] 楊昕光,遲世春.土石壩坡抗震極限分析上限法[J].巖土力學,2013,34(3):721-726.(YANG Xinguang,CHI Shichun.Upper bound method for seismic stability limit analysis of earth-rock dams[J].Rock and Soil Mechanics,2013,34(3):721-726.(in Chinese))

[18] 顧淦臣.土石壩地震工程[M].南京:河海大學出版社,1989.

[19] 李子陽,谷艷昌,張磊.基于MSC.Marc的高拱壩蓄水期變形分析[J].水利水電科技進展,2008,28(4):15-19.(LI Ziyang,GU Yanchang,ZHANG Lei.MSC.Marc-based deformation analysis of high arch dam in water storage period [J].Advances in Science and Technology of Water resources,2008,28(4):15-19.(in Chinese))

[20] 楊貴,何敦明,劉漢龍,等.土石壩動力離心模型試驗顆粒流數值模擬[J].河海大學學報(自然科學版),2011,39(3):296-301.(YANG Gui,HE Dunming,LIU Hanlong,et al.Advancement of technology on shaking table model test anddynamic analysis of CFRD [J].Journal of Hohai University(Natural Sciences),2011,39(3):296-301.(in Chinese))

[21] 趙劍明,常亞屏,陳寧.強震區高混凝土面板堆石壩地震殘余變形與動力穩定分析[J].巖石力學與工程學報,2004,23(增刊 1):4547-4552.(ZHAO Jianming,CHANG Yaping,CHEN Ning.Study on earthquake-induced permanent deformation and dynamic stability of high concrete-faced rockfill dam [J].Chinese Journal of Rock Mechanics and Engineering,2004,23(Sup1):4547-4553.(in Chinese))