王志坤，陳艷江，楊　璐,*

(1.沈陽工業(yè)大學 建筑與土木工程學院，沈陽 110870;2.黑龍江大學 水利電力學院，哈爾濱 150080)

?

基于附加質(zhì)量法的三維重力壩地震響應(yīng)研究

王志坤1，陳艷江2，楊璐1,*

(1.沈陽工業(yè)大學 建筑與土木工程學院，沈陽 110870;2.黑龍江大學 水利電力學院，哈爾濱 150080)

應(yīng)用ABAQUS有限元軟件二次開發(fā)功能，編寫三維附加質(zhì)量單元子程序UEL，實現(xiàn)了三維模型中節(jié)點附加質(zhì)量的快速施加，并模擬Koyna混凝土重力壩地震破壞過程。研究發(fā)現(xiàn)：壩體中上部特別是壩面突變處是抗震薄弱部位，壩頂以下40 m左右出現(xiàn)貫穿裂縫；下游折坡點位置易發(fā)生應(yīng)力集中，瞬時豎向拉應(yīng)力達到12.738 MPa，應(yīng)盡量采用圓弧等形式避免壩面突變。地震過程中壩體上部震動劇烈，特別是壩頂處最大震動幅度為41.42 mm。計算結(jié)果基本符合1967年Koyna大壩震害記錄數(shù)據(jù)。數(shù)值結(jié)果證明三維附加質(zhì)量法子程序的開發(fā)是正確可靠的，用該二次開發(fā)子程序方法解決此類三維水體—結(jié)構(gòu)相互作用問題是可行的。

附加質(zhì)量法；動水壓力；重力壩；時程分析

近年來，我國水利事業(yè)蓬勃發(fā)展，建設(shè)了一批大型水利工程，但我國是一個地震多發(fā)國家，不少已建、在建、擬建的高壩大庫都位于強震區(qū)：雙江口堆石壩設(shè)計地震加速度0.205g，小灣拱壩設(shè)計地震加速度0.31g，大崗山拱壩設(shè)計地震加速度0.557g[1-2]，這些高壩大庫的發(fā)展也引起了專家學者對水壩抗震性能的重視，而其中一個重要的方面便是對地震作用下動水壓力的研究。

目前工程界處理地震作用下水壩動水壓力的問題采用的是Westergaard提出的不考慮水體可壓縮性的附加質(zhì)量法[3-4]。但目前所有的大型有限元軟件中，均沒有相應(yīng)的附加質(zhì)量單元類型，ABAQUS中僅提供了二維問題的附加質(zhì)量單元快速施加程序，但二維模型應(yīng)用范圍有限[5]。鑒于此本文基于ABAQUS軟件用戶子程序接口，編寫三維附加質(zhì)量單元子程序UEL，實現(xiàn)了三維模型中節(jié)點附加質(zhì)量的快速施加，為今后模擬此類三維水體-結(jié)構(gòu)相互作用問題提供了簡便快捷的處理方法。

1　計算原理

在地震荷載作用下，加速度激勵是隨時間不斷變化大小和方向的，壩體產(chǎn)生與之對應(yīng)的往復(fù)加速晃動，水體與壩體之間產(chǎn)生大小和方向也不斷變化的相對慣性作用力以及相對滑動。水體對壩體的作用可分為兩部分壓力的疊加：①靜水壓力；②地震作用下的慣性力，即動水壓力[6]。

Westergaard假定壩面最大動水壓力沿水深呈拋物線分布，并根據(jù)實際動水壓力對壩踵力距與近似動水壓力圖形對壩踵的力距相等的條件導(dǎo)出了韋氏附加質(zhì)量公式[7]：

(1)

式中ρω為水的密度；hω為設(shè)計水深；z為到水面的距離。

由于該方法簡單實用、易于計算、偏于安全，美國、日本等國家的建筑物抗震設(shè)計規(guī)范中至今仍沿用此忽略壩體變形和水體可壓縮性的動水壓力公式[8-10]，我國《水工建筑物抗震設(shè)計規(guī)范》亦是采用附加質(zhì)量模型作為重力壩抗震設(shè)計中動水壓力模擬方法[11]。

2　模型建立

2.1工程概況

迄今為止，全球范圍內(nèi)的大壩遭受強震震害的實例稀少，遭受Ⅷ度以上強震的百米級重力壩僅有4例，即中國的新豐江大壩、寶珠寺大壩，伊朗的西菲羅(Sefid Rud)大頭壩，印度的柯依那(Koyna)重力壩[1]。其中，Koyna大壩作為少數(shù)幾個在強震中破壞且有比較完整記錄的重力壩，一直是重力壩抗震分析中的經(jīng)典研究對象，已有大量學者[12-15]對該壩的地震破壞過程進行了模擬分析。

1967年12月11日Koyna壩址區(qū)域遭受了一次6.5級強烈地震，地震發(fā)生時的壩前水位91.75 m[12]。地震作用后，壩頂以下40 m左右產(chǎn)生了多條裂縫，下游面出現(xiàn)了嚴重的漏水現(xiàn)象。

2.2計算模型

Koyna混凝土重力壩位于印度的Koyna河上，壩長850 m，壩高103 m，壩底寬70 m，壩頂寬度14.8 m，地震發(fā)生時壩前水位為91.75 m，具體尺寸見圖1。

混凝土重力壩是大體積結(jié)構(gòu)，為防止開裂、滿足施工要求往往分縫分塊澆注，需沿壩軸線方向按約20 m設(shè)置橫縫，橫縫需在壩體冷卻至穩(wěn)定溫度時經(jīng)灌漿后形成整體，但灌漿的漿體僅能起傳遞壓應(yīng)力的填充作用，抗拉強度極低以至可以忽略。故而重力壩抗震分析一般取單個壩段進行計算，壩段間橫縫的影響可以不計。

依據(jù)上述Koyna混凝土重力壩尺寸選取單個壩段進行模擬計算，壩段寬度取為20 m，分別自上、下游壩踵、壩趾位置延伸2倍壩高，同時向下延伸1.5倍壩高，作為壩基，建立有限元模型見圖2，采用C3D8R單元，即8節(jié)點實體減縮單元，上游壩體與水體接觸表面采用附加質(zhì)量單元，模型中共33 792個節(jié)點，29 850個單元，其中29 200個C3D8R單元，650個附加質(zhì)量單元。

圖1　壩體剖面(單位m)Fig.1　Dam section profile(company:m)

圖2　計算模型Fig.2　Computational model

2.3計算參數(shù)

在模型中壩體材料采用混凝土損傷塑性本構(gòu)模型，壩基部分采用彈性模型。根據(jù)文獻[11]的規(guī)定，在抗震強度計算中，混凝土動態(tài)強度和動態(tài)彈性模量的標準值可較其靜態(tài)標準值提高30%，混凝土動態(tài)抗拉強度的標準值可取為動態(tài)抗壓強度標準值的10%[11]。依據(jù)相關(guān)文獻及Koyna大壩工程簡介，壩體混凝土各項力學參數(shù)見表1。

表1　壩體力學參數(shù)

另外，由于混凝土抗壓強度一般較大，故而沒有考慮壓縮引起的損傷[13]，混凝土拉伸屈服應(yīng)力、拉伸損傷因子和開裂位移之間的關(guān)系見圖3。

圖3　混凝土塑性指標Fig.3　Concrete plastic properties

2.4地震波

在進行壩體結(jié)構(gòu)動力分析中，同時計入壩體和地基的剛度，但只計入壩體的質(zhì)量，地基單元只考慮彈性不考慮質(zhì)量，以消除波的傳播效應(yīng)，避免人為的放大作用，這樣的處理方式較為接近實際工程情況。

地震波在空間的傳播方向十分復(fù)雜，不同方向的地震波對重力壩影響差異較大，模型中在壩基底面邊界施加水平和垂直方向地震加速度，選取1967年Koyna大壩實測的柯依那波進行計算，時間間隔0.01 s，地震總時程10 s。地震波時程曲線見圖4。

圖4　地震加速度時程曲線Fig.4　Seismic acceleration time history curve

3　模型分析

模態(tài)分析是各種動力學分析類型中的基礎(chǔ)內(nèi)容，在進行其他動力學分析之前首先要進行模態(tài)分析，采用ABAQUS有限元軟件中的Lanczos法提取大壩前4階頻率。表2、表3分別為Koyna重力壩空庫和滿庫狀態(tài)下前4階頻率，通過與文獻[14]對比可知誤差較小，證明建立的模型及開發(fā)的三維附加質(zhì)量單元用戶子程序是正確的。

表2　Koyna重力壩(空庫)前4階頻率

表3　Koyna重力壩(滿庫)前4階頻率

4　時程分析

運用ABAQUS有限元軟件對三維Koyna壩體進行動力時程運算并分析計算結(jié)果，分別提取壩體上關(guān)鍵位置點，并研究其位移、應(yīng)力等隨時間的變化情況，從而揭示混凝土重力壩在地震作用下的一般反應(yīng)規(guī)律。關(guān)鍵位置點布置見圖1，其中A、B分別為為上、下游壩頂，C為壩踵，D為壩趾，E為下游折坡點。

地震持續(xù)過程中壩體上各特征點在各個方向上應(yīng)力和位移的峰值及相應(yīng)峰值出現(xiàn)的時刻見表4，其中水平位移負值為上游方向，豎直位移負值為垂直向下方向，負應(yīng)力值代表拉應(yīng)力。由表4可見，壩頂處應(yīng)力很小，震后也未出現(xiàn)塑性破壞，只是在地震過程中晃動比較嚴重，出現(xiàn)較大位移。壩踵和壩趾處在地震中的晃動幾乎可以忽略，水平方向的應(yīng)力也很小，在地震波最劇烈的4～5 s內(nèi)壩踵處的豎直方向應(yīng)力較大，產(chǎn)生了塑性應(yīng)變，這是由于剛性地基沒有考慮地基的影響造成的，實際震害中不存在這種狀況，故不加分析。

相關(guān)震害記錄及研究資料表明，重力壩上部尤其是斷面突變處，即折坡點位置，是抗震薄弱部位。表4及圖5也表明地震對重力壩的破壞作用主要集中在斷面突變處，故著重分析折坡點E處應(yīng)力、位移等在地震中的變化規(guī)律。

表4　各特征點峰值應(yīng)力、位移

4.1應(yīng)力時程曲線

地震過程中特征點E的應(yīng)力時程曲線見圖5。由圖5可見，折坡點處水平和豎直方向應(yīng)力變化趨勢基本相同，只是峰值不同，水平應(yīng)力峰值為4.522 MPa，豎直應(yīng)力峰值為12.738 MPa，均超過了混凝土動態(tài)抗拉強度，出現(xiàn)了較嚴重的損傷破壞。兩個方向應(yīng)力峰值均出現(xiàn)在4.8 s左右，這也是壩體裂縫擴展最快的時刻，略滯后于地震波最大加速度出現(xiàn)時間。

在地震持續(xù)過程中，前期由于地震加速度小，壩體能夠抵御由震動引起的破壞作用，兩個方向的應(yīng)力均低于混凝土的抗拉強度，壩體只發(fā)生晃動，并未出現(xiàn)裂縫等不可逆轉(zhuǎn)的破壞現(xiàn)象。隨著地震波加強，地震的破壞作用超過了壩體自身抵抗能力，壩面突變處產(chǎn)生應(yīng)力集中，應(yīng)力值隨著地震的持續(xù)急劇增加，尤其是豎向應(yīng)力值，最大達到抗拉強度的4倍以上。地震作用末期，加速度減弱，拉應(yīng)力也相應(yīng)減小，但豎向拉應(yīng)力依然較大，直到地震結(jié)束。

圖5　折坡點E應(yīng)力時程曲線 Fig.5　Stress time history curve of slope point E

4.2位移時程曲線

下游折坡點E的時間位移時程曲線見圖6，其中水平方向主要是隨著地震過程的持續(xù)壩體沿上下游方向晃動，晃動幅度基本相同，均隨著地震加速度的加強而增大，但略滯后于加速度，在4.48 s時達到最大位移11.43 mm。隨后地震作用減弱，壩體晃動幅度減小，至地震作用結(jié)束時壩體基本回到原來位置，僅向上游方向偏移1.45 mm。

在地震作用前期，壩體在豎直方向僅發(fā)生上下晃動，晃動幅度在2 mm以內(nèi)，但在地震加速度峰值過后，下游壩面突變處即E點處出現(xiàn)塑性損傷，4 s后E點的豎直方向出現(xiàn)2 mm左右的位移，并一直持續(xù)到地震結(jié)束，表明在E處出現(xiàn)了寬度為2 mm左右的裂縫，與1967年的Koyna大壩震害記錄吻合。

圖6　折坡點E位移時程曲線Fig.6　Displacement time curve of slope point E

4.3損傷分析

Koyna大壩是國外遭受Ⅷ度以上強震的兩個百米級重力壩之一，震災(zāi)情況為12～18號、24～30號壩段壩頂以下40 m左右位置產(chǎn)生了多條水平裂縫，下游面出現(xiàn)了嚴重的漏水現(xiàn)象，即裂縫已經(jīng)貫通了上下游。說明重力壩上部，尤其是斷面突變處，是其抗震薄弱部位，在強震作用下，上部壩體易開裂[16]。震后Koyna壩體保持了整體穩(wěn)定性，壩基面未出現(xiàn)拉裂和剪切損壞[17]。

模擬的地震過程中壩體斷面拉伸損傷云紋圖見圖7。由圖7可見，壩體破壞主要集中在下游折坡點處。地震前期地震動加速度較小，未對壩體造成損傷破壞，圖7(a)為2 s時的壩體斷面拉伸損傷云紋圖，損傷值為零。2 s后地震加速度變大，震動劇烈，壩體局部區(qū)域主要是下游折坡點位置及壩頂以下40 m左右處，4 s左右壩體開始出現(xiàn)明顯的破壞，損傷區(qū)域由壩面處向壩體內(nèi)部擴展，逐漸形成貫穿上下游的水平裂縫。5 s后地震作用減弱，破壞面積沒有進一步擴大，主要是已經(jīng)形成的貫穿裂縫繼續(xù)發(fā)展，嚴重的貫穿裂縫導(dǎo)致下游面出現(xiàn)漏水現(xiàn)象，符合1967年的震害記錄。

壩踵處出現(xiàn)了輕微的損傷破壞，原因是剛性地基沒有考慮地基的影響[18]，但震后壩體保持了整體穩(wěn)定性，亦符合當時的實際情況。說明三維附加質(zhì)量法的運用是合理可行的，能相對準確地模擬地震作用下水體對壩體的相互作用情況，為今后工程界處理此類問題提供了簡便可靠的解決方法。

圖7　壩體斷面損傷云紋圖Fig.7　Damage moire of dam section

5　結(jié)　論

依據(jù)1967年印度Koyna混凝土重力壩地震相關(guān)資料數(shù)據(jù)，建立壩體三維計算模型，采用Westergaard提出的附加質(zhì)量法考慮地震作用下動水壓力的影響，利用ABAQUS軟件用戶子程序接口UEL，開發(fā)編寫三維附加質(zhì)量單元子程序。進行地震時程計算，并分析壩體損傷、應(yīng)力、位移等在地震過程中的響應(yīng)規(guī)律，得出以下結(jié)論：

1) 計算結(jié)果與1967年震害記錄基本吻合，說明三維附加質(zhì)量子程序的開發(fā)是合理可行的，為工程界處理類似三維水體-結(jié)構(gòu)動水相互作用問題提供了新的簡單快捷的方法。

2) 重力壩在地震過程中的損傷破壞主要集中在壩體中上部，壩面突變處是其抗震薄弱部位，特別是下游折坡點位置易出現(xiàn)貫穿裂縫，應(yīng)引起足夠重視。

3) 下游折坡點處易產(chǎn)生應(yīng)力集中，應(yīng)力值遠大于混凝土的動態(tài)抗拉強度值，最大拉應(yīng)力達到12.738 MPa，應(yīng)避免出現(xiàn)壩面突變，折坡點處盡量采用圓弧等形式。

4) 地震中壩頂震動劇烈，主要沿上下游方向晃動，最大振幅達到41.42 mm，壩頂處構(gòu)造設(shè)施應(yīng)采取相應(yīng)措施避免發(fā)生破壞。

[1]潘家錚. 水電與中國[M]. 北京：中國水利水電出版社,2007.(PAN Jia-Zheng. Hydro-power and China[M]. Beijing: China Water Conservancy and Hydro-power Press,2007.(in Chinese))

[2]陳厚群. 混凝土高壩抗震研究[M]. 北京：高等教育出版社，2011.(CHEN Hou-Qun. Study on seismic behavior of concrete dam[M]. Beijing: Higher Education Press,2011.(in Chinese))

[3]Lee J, Fenves G L. A plastic-damage concrete model for earthquake analysis of dams[J].Earthquake Engineering and Structural Dynamics,1998,27(4)：937-956.

[4]Omidi V S. Seismic cracking of concrete gravity dams by plastic-damage model using different damping mechanisms[J].Finite Elements in Analysis and Desing,2013,63(6)：80-97.

[5]徐磊. 基于ABAQUS的四節(jié)點三維附加質(zhì)量單元開發(fā)[J]. 水電能源科學，2011,29(6)：100-102.(XU Lei. Development of four nodes 3-D added mass element in ABAQUS[J]. Water Resources and Power, 2011,29(6)：100-102.(in Chinese))

[6]宋波，劉浩鵬，張國明. 基于附加質(zhì)量法的橋墩地震動水壓力分析與實例研究[J]. 土木工程學報，2010,43(增刊)：102-107.(SONG Bo, LIU Hao-Peng, ZHANG Guo-Ming. Analysis and case study of earthquake-induced hydrodynamic pressure on deep water bridge based on added masses approach[J]. China Civil Engineering Journal,2010,43(s):102-107.(in Chinese))

[7]Westergaard H M. Water pressures on dams during earthquakes[J].Transactions of the American Society of Civil Engineers, 1933,98(3)：418-433.

[8]王銘明，陳健云，范書立. 重力壩地震動水壓力試驗研究[J]. 水電能源科學，2012,30(5)：51-53.(WANG Ming-Ming, CHEN Jian-Yun, FAN Shu-Li. Experimental study on seismic water pressure of gravity dam[J]. Water Resources and Power, 2012,30(5)：51-53.(in Chinese))

[9]Pailoplee, Santi. Earthquake hazard of dams along the Mekong mainstream[J]. Natural Hazards,2014,74(3)：1813-1827.

[10] Das, R., Cleary, P. W. A mesh-free approach for fracture modeling of gravity dams under earthquake[J].International Journal of Fracture,2013,179(2)：9-33.

[11] 中華人民共和國水利部.水工建筑物抗震設(shè)計規(guī)范：SL203-97[S].北京：中國水利水電出版社，1997. (Ministy of Water Resources of the People’s Republic of China.Code for seismic design of hydraulic structures：SL203-97[S]. Beijing:China Water Conservancy and Hydropower Press,1997.(in Chinese))

[12] 徐金英，李德玉，郭勝山. 基于ABAQUS的兩種庫水附加質(zhì)量模型下重力壩動力分析[J]. 中國水利水電科學研究院學報，2014,12(1)：98-103.(XU Jin-Ying, LI De-Yu, GUO Sheng-Shan. Dynamic analysis of gravity dam under two kinds of reservoir water quality model based on ABAQUS[J]. Journal of China Institute of Water Resources and Hydropower Research, 2014,12(1)：98-103.(in Chinese))

[13] 楊璐，金峰．水工混凝土結(jié)構(gòu)中的數(shù)值計算與實例[M]. 北京：科學出版社,2015.(YANG Lu, JIN Feng. Numerical calculation and examples of hydraulic concrete structures[M]. Beijing: Science Press,2015.(in Chinese))

[14] 趙光恒. 結(jié)構(gòu)動力學[M]. 北京：中國水利水電出版社，1996.(ZHAO Guang-Heng. Structural dynamics[M]. Beijing: China Water Conservancy and Hydro-power Press,1996.(in Chinese))

[15] 黃耀英，孫大偉，田斌. 兩種庫水附加質(zhì)量模型的重力壩動力響應(yīng)研究[J]. 人民長江，2009,40(7)：64-66.(HUANG Yao-Ying, SUN Da-Wei, TIAN Bin. Study on dynamic response of gravity dam with two kinds of reservoir water quality model[J]. Yangtze River, 2009,40(7)：64-66.(in Chinese))

[16] 方修君，金峰，王進廷. 基于擴展有限元法的Koyna重力壩地震開裂過程模擬[J]. 清華大學學報:自然科學版. 2008,48(12)：2065-2069.(FANG Xiu-Jun, JIN Feng, WANG Jin-Ting. Simulation of seismic cracking of Koyna gravity dam based on extended finite element method[J]. Journal of Tsinghua University(Science and Technology), 2008,48(12)：2065-2069.(in Chinese))

[17] 邵長江，錢永久. Koyna混凝土重力壩的塑性地震損傷響應(yīng)分析[J]. 振動與沖擊，2006,25(4)：129-131.(SHAO Chang-Jiang, QIAN Yong-Jiu. Analysis of plastic seismic damage response of Koyna concrete gravity dam[J]. Journal of Vibration and Shock, 2006,25(4)：129-131.(in Chinese))

[18] Kartal, Murat Emre,Bayraktar. Non-linear response of the rock-fill in a concrete-faced rock-fill dam under seismic excitation[J]. Mathematical and Computer Modelling of Dynamical Systems,2014,21(1)：77-101.

Research on seismic response of 3-D gravity dam based on added mass method

WANG Zhi-Kun1, CHEN Yan-Jiang2, YANG Lu1,*

(1.SchoolofArchitectureandCivilEngineering,ShenyangUniversityofTechnology,Shenyang,110870,China;2.SchoolofHydrauicandElectricPower,HeilongjiangUniversity,Harbin150080,China)

3-D added mass subroutine UEL based on the redevelop function of ABAQUS is developed，fast application of added mass in 3-D model is realized, and seismic damage process of Koyna concrete gravity dam is simulated. According to the result: middle upper part especially the mutation of dam body is weak part in Seismic, and through cracks appear below the dam crest of 40m around; The position of the downstream slope occurs stress concentration, the instantaneous vertical tensile stress reaches 12.738 MPa, and the circular arc and other forms should be used to avoid the mutation of the dam surface. In earthquake process, the upper shakes violently especially the maximum vibration amplitude of the dam crest reaches 41.42 mm. The calculation results are in accord with the damage record of Koyna dam in 1967, the development of 3-D added mass subroutine is correct and reliable, and the method to deal with the interaction problem of 3-D water-structure is feasible.

added mass method；hydrodynamic pressure；gravity dam；time history analysis

10.13524/j.2095-008x.2016.03.035

2016-05-06

國家自然科學基金資助項目(11102118)

王志坤(1990-)，男，河北昌黎人，碩士研究生，研究方向：水工結(jié)構(gòu)抗震及ABAQUS數(shù)值模擬， E-mail：xxzhikun@163.com;*通訊作者：楊璐(1973-),女，山東掖縣人，教授，博士，碩士研究生導(dǎo)師，研究方向：混凝土彈塑性損傷本構(gòu)和ABAQUS的數(shù)值模擬，E-mail:yanglu515@163.com。

TV312

A

2095-008X(2016)03-0012-08