考慮滲流體力的某大壩地震響應(yīng)仿真分析

韓 彰,蘇懷智,李 慧

(1.河海大學(xué)水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室,江蘇 南京 210098;2.河海大學(xué)水利水電學(xué)院,江蘇 南京 210098)

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

1 工程概況

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

2 數(shù)值模擬方法

采用動力計算本構(gòu)模型對該水庫均質(zhì)壤土擋水壩段及其混凝土連接壩段非線性有限元動力特性進(jìn)行數(shù)值模擬。

2.1 動力計算本構(gòu)模型

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

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

c. 接觸面的動力計算模型。接觸面單元的動力參數(shù)采用河海大學(xué)的試驗(yàn)成果[18]。接觸面單元的動力模型剪切勁度Kc與動剪應(yīng)變γd、阻尼比λc的關(guān)系分別為

(1)

(2)

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

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

2.2 壩體結(jié)構(gòu)網(wǎng)格、邊界條件與荷載

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

圖2 樁號位置示意圖

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

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

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

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

2.3 動力計算參數(shù)

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

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

表1 不同動剪應(yīng)變下壩料的動剪切模量比和動阻尼比

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

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

步驟2求出土體單元的初始動剪切模量Gdmax,0,土體單元的初始阻尼比根據(jù)經(jīng)驗(yàn)取為5%。

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

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

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

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

步驟10根據(jù)求出的結(jié)點(diǎn)位移un+1計算各單元的動剪應(yīng)變γn+1和動剪應(yīng)力τn+1。

步驟11重復(fù)步驟 5～10,得到各單元在每個時段內(nèi)的動剪應(yīng)變γd時程。

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

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

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

步驟15輸出計算結(jié)果。

2.5 地震加速度的輸入

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

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

圖5 場地基巖水平加速度反應(yīng)譜

圖6 地震加速度時程曲線

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

圖7 0+158斷面12個結(jié)點(diǎn)位置及編號

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

3 三維非線性動力計算結(jié)果與分析

3.1 設(shè)計地震反應(yīng)

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

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

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

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

圖11 0-81斷面接觸面結(jié)點(diǎn)分布

由計算成果可知:①最大地震動加速度、動位移反應(yīng)位于壩頂局部位置,壩頂存在明顯的鞭鞘效應(yīng),需要在壩頂進(jìn)行適當(dāng)?shù)目拐鸺庸?②計算得到的最大放大系數(shù)為7.70,這可能是由于大壩很長,約束又加在兩端,正如橡皮筋一樣,導(dǎo)致計算得到的反應(yīng)加速度放大倍數(shù)偏大,動力響應(yīng)偏大,建議在壩頂加土工格柵。

圖12 端頭接觸面結(jié)點(diǎn)分布

3.2 上下游接觸面及端頭接觸面典型結(jié)點(diǎn)動位移

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

表2 0-81斷面接觸面結(jié)點(diǎn)動位移統(tǒng)計

表3 端頭接觸面第1～4列結(jié)點(diǎn)動位移統(tǒng)計

3.3 端頭接觸面典型單元法向動應(yīng)力

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

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

4 結(jié) 論

a. 該大壩的最大地震動加速度、動位移反應(yīng)位于壩頂局部位置,壩頂存在明顯的鞭鞘效應(yīng),需要在壩頂進(jìn)行適當(dāng)?shù)目拐鸺庸獭?/p>

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

c. 各分區(qū)的設(shè)計與填筑標(biāo)準(zhǔn)、壩體分層填筑方案合理,壩體抗震安全性較好。

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

[2] 龐林祥,傅華.高土石壩的地震動力響應(yīng)分析方法研究與應(yīng)用[J].水利水電技術(shù),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].水資源與水工程學(xué)報,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].西安:西安理工大學(xué),2013.

[5] 李炎隆,徐嬌,李守義.瀝青混凝土心墻堆石壩動力有限元分析[J].應(yīng)用力學(xué)學(xué)報,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] 沈振中,崔娟,任華春.土石壩抗震穩(wěn)定性分析的極限平衡水平條分法[J].南水北調(diào)與水利科技,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] 劉小生,王鐘寧,趙劍明,等.面板堆石壩振動模型試驗(yàn)及動力分析研究[J].水利學(xué)報,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] 竇興旺,林義興,夏頌佑,等.深覆蓋層上高堆石壩振動臺試驗(yàn)與動力數(shù)值分析驗(yàn)證研究[J].水利水電科技進(jìn)展,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] 竇興旺,夏頌佑,許百立.人工邊界法在土石壩動力分析中的應(yīng)用[J].河海大學(xué)學(xué)報(自然科學(xué)版),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] 楊秀竹,李花珍,王星華.地震力作用下軟基土石壩的動力響應(yīng)分析[J].鐵道科學(xué)與工程學(xué)報,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] 丁陸軍.面板堆石壩三維動力反應(yīng)分析[D].西安:西安理工大學(xué),2005.

[12] 曹學(xué)興,何蘊(yùn)龍,熊堃,等.考慮圍壓效應(yīng)的高土石壩動力響應(yīng)分析[J].天津大學(xué)學(xué)報(自然科學(xué)與工程技術(shù)版),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] 孔憲京,周揚(yáng),鄒德高,等.汶川地震紫坪鋪面板堆石壩地震波輸入研究[J].巖土力學(xué),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] 孔憲京,鄒德高,鄧學(xué)晶.高土石壩綜合抗震措施及其效果的驗(yàn)算[J].水利學(xué)報,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].巖土力學(xué),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].巖土力學(xué),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].南京:河海大學(xué)出版社,1989.

[19] 李子陽,谷艷昌,張磊.基于MSC.Marc的高拱壩蓄水期變形分析[J].水利水電科技進(jìn)展,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] 楊貴,何敦明,劉漢龍,等.土石壩動力離心模型試驗(yàn)顆粒流數(shù)值模擬[J].河海大學(xué)學(xué)報(自然科學(xué)版),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] 趙劍明,常亞屏,陳寧.強(qiáng)震區(qū)高混凝土面板堆石壩地震殘余變形與動力穩(wěn)定分析[J].巖石力學(xué)與工程學(xué)報,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))