薛 逵，付英杰，張永學

（1.普蘭店市水土保持辦公室，遼寧普蘭店116200；2.沈陽經(jīng)濟技術開發(fā)區(qū)農(nóng)業(yè)發(fā)展局，遼寧沈陽110016）

地震荷載下土石壩變形分析

薛 逵1，付英杰2，張永學2

（1.普蘭店市水土保持辦公室，遼寧普蘭店116200；2.沈陽經(jīng)濟技術開發(fā)區(qū)農(nóng)業(yè)發(fā)展局，遼寧沈陽110016）

土石壩在地震荷載下的內(nèi)部應力應變變化越來越受到重視，針對此問題以遼寧省撫順市清原縣小孤家水庫土石壩為例，使用ABAQUS軟件，模擬其在地震荷載下的應力應變變化情況，分析計算結(jié)果，最后得出該土石壩在此荷載條件下x方向位移和永久變形分布圖及加速度變化曲線圖，其結(jié)果與理論相符，利用ABAQUS軟件可以模擬土石壩在實際工程中的受力情況，對于今后土石壩內(nèi)部應力研究有一定的意義。

土石壩；地震荷載；應力應變；變形分析

0　引言

土石壩是目前世界上建設數(shù)量最多的壩型，據(jù)國內(nèi)外建壩數(shù)量統(tǒng)計，土石壩占已建壩總數(shù)的90%以上，我國的土石壩建設主要集中在20世紀50到70年代。隨著近年地震多發(fā)，各地的專家開始對于土石壩在地震荷載下的受力及形變情況更加重視。因為在發(fā)生地震的時段，含有大片非粘性土料（如砂）的土石壩，在飽和狀態(tài)下受地震震動時，其強度會降低，從而導致過大的有害位移。受地震破壞的主要原因是土體中的孔隙壓力增大，導致土的強度降低而發(fā)生液化，是導致垮壩或嚴重破壞的主要原因。

本文針對以上現(xiàn)象，對大壩觀測和現(xiàn)場試驗所得到的及時準確的外部數(shù)據(jù)進行整理和分析，使用ABAQUS軟件，應用有限元方法計算出大壩內(nèi)部的應力應變情況，荷載分布狀態(tài)。根據(jù)有限元法的計算結(jié)果、水庫水位和壩體土壤層結(jié)構模擬地震時土石壩內(nèi)部孔隙水壓力的變化以及壩體內(nèi)個單元的變形。

1　工程概況

遼寧省撫順市清原縣小孤家水庫是一座以城鎮(zhèn)供水為主、兼顧灌溉、防洪、水產(chǎn)養(yǎng)殖等中型水利樞紐工程。水庫大壩為粘土心墻壩，壩長347.0 m，最大壩高29.5 m，正常蓄水位481.5 m。壩體主要由中粗砂組成，壩殼材料為砂礫料，壩址區(qū)基巖為前震旦紀鞍山群變質(zhì)巖系片麻巖，巖性為角閃斜長片麻巖：暗灰色或灰褐色，中細粒變晶結(jié)構，片麻狀構造，巖石整體性差，但風化不嚴重，新鮮巖石較堅硬。第四系沖積層廣泛分布于河道兩側(cè)、河漫灘、二級階地之上。巖性自上而下由中細砂和粗砂及部分礫石、卵石組成。主壩基礎開挖深3～4 m，底寬5.0 m，建基面為弱風化巖，副壩基開挖深達29 m才見基巖的探坑；壩端處理即壩肩基槽最底蓋板高程297 m，順巖坡到305 m高程，心墻插入巖石的深度一般是1.2～1.4 m，底寬1.5 m。

根據(jù)以上地質(zhì)資料將小孤家水庫主壩及壩址區(qū)土體和巖石屬性分為5種：1）壩前粗砂土；2）壩后中粗砂土；3）心墻粘土；4）壩下強風化巖層，厚度8 m；5）壩下中風化巖層，厚度4 m。

2　計算模型及參數(shù)

2.1計算模型

本文應用土的靜力有限元計算模型、土的滲流有限元計算模型和土的地震有限元計算模型分別得出土石壩在靜力、蓄水和地震3種狀態(tài)下的受力和形變情況。計算方程如下：

其中：［B］是應變位移矩陣；［C］是本構矩陣；［a］是節(jié)點x與y位移增量的列向量；A是沿元素邊界的區(qū)域；v是單元素體積；b是單位體荷載強度；〈N〉是行向量插值函數(shù)；p是表面壓力增量；｛Fn｝是節(jié)點集中荷載增量；q為邊界節(jié)點單位流量；u為孔隙水壓力；γ為水的容重；［M］是質(zhì)量矩陣；［D］是阻尼矩陣；［K］是剛度矩陣；［F］是載體荷載；｛·a·｝是節(jié)點矢量加速度；｛·a｝是節(jié)點矢量速度；｛a｝是節(jié)點矢量位移。

方程（1）為土的靜力有限元計算模型，用于計算對于一個給定的時間增量內(nèi)通過荷載增量產(chǎn)生的位移增量，并將位移增量累積最終得出土石壩各單元的變形值。方程（2）為土的滲流有限元計算模型，它可以得到任意時間段內(nèi)土體內(nèi)孔隙水壓力變化情況及分布，將方程（1）與（2）耦合分析就能得出土石壩在蓄水狀態(tài)下各單元的變形值。方程（3）為土的地震有限元計算模型，方程為有限元分析得到的動態(tài)響應體系方程，與（1）和（2）耦合分析能夠得出土石壩各單元在地震荷載作用下的永久形變。

2.2模型參數(shù)

本文選擇摩爾-庫倫模型為土的本構模型，土石壩計算模型參數(shù)由現(xiàn)場取得土樣并在實驗室內(nèi)做三軸試驗獲得，根據(jù)固結(jié)排水試驗結(jié)果，計算模型所需參數(shù)見表1。

表1　各土層計算參數(shù)表

在計算地震荷載作用時，還要引入兩組地震加速度時程，兩組時程數(shù)值相同，但是方向分別是x軸方向和y軸方向。兩組時程基準期為100年，超越概率為2%，加速度最大值是116 gal，時程采樣間隔0.02 s，持續(xù)時間40 s。根據(jù)動態(tài)周期荷載作用下土體剪切模量隨著周期剪切應變的增長而增長的關系，周期剪應變振幅與剪切模量的變化關系可以由G-Reduction函數(shù)來表示，由此函數(shù)推得動態(tài)周期荷載條件下的剪切模量G。由于G隨著周期剪應變的增加而減少，區(qū)域內(nèi)應力應變滯后回線所包圍的面積增加，導致了阻尼比隨著動應變振幅的增加而增加。因此，在動荷載作用下的非線性的剪切模量和阻尼比可以近似地用等效線性分析來模擬，即常數(shù)G和阻尼比在動荷載分析中使用，新的模量G和阻尼比在動態(tài)分析中通過周期或者等效周期剪應變計算獲得。然后，用新模量和阻尼比重新進行動態(tài)荷載分析。

3　計算與分析

3.1幾何模型建立及其定解條件

本文對小孤家水庫土石壩進行二維有限元分析，幾何模型所選的是主壩最大斷面縱剖面，壩頂最大高程495.5 m，底部中風化層最低高程447.0 m。

此次計算為竣工期，大壩的正常蓄水位為481.5 m，在靜力與蓄水狀態(tài)下，壩下土層側(cè)向水平位移受限制，土層底部水平位移及垂直位移受限制，心墻進入土層部分不做主從面處理，在蓄水狀態(tài)時心墻迎水面做主從面處理以模擬不透水心墻在水壓力及孔隙水壓力作用下產(chǎn)生的形變。靜力分析時設置作用時間為1 s，時間步長0.025 s；蓄水時設置作用時間為10 s，時間步長0.000 1 s；地震時設置作用時間為40 s，時間步長0.02 s。在地震分析中，將二維有限元分析用一維有限元分析代替，即地震分為兩部分，分別導入一組x方向加速度時程A1和一組y方向加速度時程A2。

3.2計算結(jié)果

應用ABAQUS對小孤家水庫土石壩進行有限元分析，在靜力狀態(tài)下產(chǎn)生的x方向和y方向變形如圖1和圖2所示，在蓄水狀態(tài)下產(chǎn)生的x方向和y方向變形如圖3和圖4所示，在地震荷載作用下產(chǎn)生的x方向和y方向變形如圖5和圖6所示。

圖1　土石壩x方向位移等值線圖

圖2　土石壩沉降等值線圖

圖3　蓄水后土石壩x方向位移等值線圖

圖4　蓄水后土石壩沉降等值線圖

圖5　 A1條件下土石壩x方向永久變形等值線圖

圖6　 A2條件下土石壩y方向永久變形等值線圖

從圖1中反映出土石壩迎水坡和背水坡分別產(chǎn)生了向左的3.926 cm最大位移和向右的3.75 cm最大位移，這是由于上面土體下壓使得坡面土體向心墻移動符合受力規(guī)律。圖2中顯示土石壩的最大沉降為12.17 cm，位置是在土石壩最上端，出現(xiàn)這種情況的原因是由于小孤家水庫土石壩的心墻不是一直修筑至壩頂?shù)?，而是與土石壩頂部有0.5 m的空隙，這部分空隙由砂土來填充所以會出現(xiàn)沉降較大的現(xiàn)象。

從圖3中可以看到在水庫蓄水后土石壩的x方向最大位移為1.02 cm，發(fā)生在迎水坡內(nèi)部。由于孔隙水壓力的作用土石壩在x方向上的位移與無水條件下的數(shù)據(jù)相比有很大的差別，不僅是數(shù)值上的改變，在孔隙水壓力和浮托力作用下土石壩迎水坡和背水坡x方向位移方向由原來的指向心墻變?yōu)榕c心墻背離方向，并且背水坡最大位移中心由于心墻迎水面受到水壓力而向左移動了很多出現(xiàn)在心墻位置。圖4中，土石壩最大沉降發(fā)生在大壩頂部砂土與心墻的空隙部分，只是由于孔隙水壓力和浮托力的作用最大位移中心向左移動了5 m，最大沉降值12.96 cm比無水條件下略有增大。從以上結(jié)果可以看出x方向位移的改變最明顯，這是由于蓄水后水對土石壩的滲透和沖擊作用所導致的，心墻作為不透水面攔截了大部分滲流，同時也承受了較大的水壓力使其位移發(fā)生顯著變化。

在地震荷載作用下土石壩形變改變較大，圖5中顯示土石壩在A1作用下x方向的永久變形最大值為57.42 cm，位置在迎水坡上，從迎水坡到背水坡變形逐漸減小，最小減至19.14 cm。在圖6中，豎向變形主要發(fā)生在心墻兩側(cè)，最大沉降為34.48 cm在迎水坡面，心墻處沉降為8.41 cm。

在地震荷載計算過程中 40 s內(nèi)的數(shù)據(jù)ABAQUS都以表格形式儲存，利用此數(shù)據(jù)可以用來繪制加速度峰值沿心墻軸線隨深度變化圖。首先將數(shù)據(jù)中加速度最值選出，然后將分布在心墻軸線上的數(shù)據(jù)點連成曲線得圖7和圖8。

圖7　加速度A1峰值沿高度的分布情況

從加速度分分布可以看出，加速度峰值的最大值出現(xiàn)在壩頂部位，土石壩30 m以上加速度變化較大波動明顯，30 m以下加速度值程下降趨勢且浮動變小。這也說明了砂土相對于風化巖層的不穩(wěn)定性和易液化性。

圖8　加速度A2峰值沿高度的分布情況

4　結(jié)論與展望

4.1結(jié)論

本文以摩爾-庫倫模型為土體本構模型，采用有限插值法對小孤家水庫土石壩進行了有限元分析及計算。計算情況分為靜力、蓄水和地震3種：靜力部分是為了得到土石壩最初的形態(tài)和應力狀態(tài)，各應力位移的分布情況，為后面的計算分析提供依據(jù)和基礎數(shù)據(jù)；蓄水是模擬水庫正常運行期間孔隙水壓力和浮托力對迎水坡和心墻的影響；地震部分綜合了前兩步的數(shù)據(jù)和結(jié)果，在給定的加速度時程曲線下分析出土石壩永久變形和加速度分布情況。

通過以上成果不難看出，由于砂土是由土顆粒構成的骨架及空隙中的水和空氣組成的，顆粒之間連接較弱，土顆粒骨架具有不穩(wěn)定性，在靜力作用下動力荷載小所以主要表現(xiàn)為彈性變形；在孔隙水壓力和循環(huán)荷載作用下動應變增大時，本來表現(xiàn)的彈性變形逐漸轉(zhuǎn)變?yōu)樗苄宰冃?，引起土結(jié)構的改變，從而使土的殘余變形和強度喪失，尤其對于飽和砂土而言極易發(fā)生振動液化。雖然本文研究對象是砂土中相對穩(wěn)定的中粗砂，但是在結(jié)果中可以看到其在靜力與動力過程中應力應變和位移的變化程度及波動是非常明顯的，甚至發(fā)生了位移反向移動的現(xiàn)象。其原因就是在土顆粒中應力是通過觸點來傳播的，但是在浸水后觸點力就要先傳遞給孔隙中的水，由水來繼續(xù)承擔。這樣在循環(huán)荷載的作用下，觸點力的增加都轉(zhuǎn)化成了孔隙水壓力從而引起了孔隙水壓力的增加，當孔隙水壓力超過靜水壓力時，就是會使土體的有效應力減少，造成動力抗剪強度降低或者完全喪失，使土發(fā)生液化。

綜上所述，砂土在靜力狀態(tài)下所表現(xiàn)強度可以滿足土石壩安全要求，在動荷載作用下所產(chǎn)生土體性質(zhì)變化和液化問題應該值得我們注意。在地震多發(fā)區(qū)建設水工建筑物時應謹慎考慮砂土及類似土型的選用。

4.2展望

在對比驗證方面，本文研究的是竣工期水庫大壩，所以在能夠得到有力驗證的沉降方面由于產(chǎn)生沉降值太小沒有實現(xiàn)與真實值的對比。希望以后能夠以施工期土石壩為研究對象進一步驗證模型分析的正確性。也希望能夠有更多的研究方法和數(shù)據(jù)來實現(xiàn)這一目的。

在分析方面，文采用二維方法建模，這種方法只能表現(xiàn)一個斷面的數(shù)據(jù)及形態(tài)，但是土石壩是三維立體式的，對于Z方向和XY平面的表現(xiàn)也很重要。在今后的工作中希望能夠使用三維建模對水工建筑進行分析計算。

在數(shù)據(jù)方面，本文中大部分數(shù)據(jù)來源于現(xiàn)場采樣在實驗室進行的三軸實驗和水庫資料數(shù)據(jù)，由于沒有做固結(jié)儀實驗所以在剪脹角方面是估計值。地震波數(shù)據(jù)也由于資源有限只采用一組加速度時程曲線來完成整個計算過程。隨著近期地震活動頻繁，各地接連發(fā)生強震并造成大量損失，人們對于地震對于建筑的災害越來越重視，希望能夠有更多數(shù)據(jù)和研究為今后建筑物抗震提供有效的依據(jù)。

［1］孔令偉，郭愛國，陳善雄，等.膨脹土的承載強度特征與機制［J］.水利學報，2004（11）：54－61.

［2］李廣信.高等土力學［M］.北京：清華大學出版社，2004.

［3］時衛(wèi)民，鄭穎人.庫水位下降情況下滑坡的穩(wěn)定性分析［J］.水利學報，2004（3）：67－71.

［4］王厚德.大壩安全監(jiān)測［C］.北京：中國水利水電出版社，2004.

［5］白冰.強夯荷載作用下飽和土層孔隙水壓力簡化計算方法［J］.巖土力學與工程學報，2003，22（9）：1469—1473.

［6］米占寬，高面板壩壩體流變性狀研究［D］.南京水利科學研究院，2001，6.

［7］汝乃華，牛運光.大壩事故與安全·土石壩［M］.北京：中國水利水電出版社，2001.

TV698.1+3

A

1002-0624（2016）01-0038-04

2015-02-03