基于有限元重度增加法的涉水岸坡穩定性分析研究

劉中才

(西藏眾興工程設計咨詢有限公司，西藏 山南 856100)

當前，已建和在建的水利工程數量數以萬計，我國水庫數量逐漸增加[1]。水庫建成蓄水后水位逐漸上升，造成山體長期浸泡在水中，水位反復升降直接削弱了山巖土體的各類力學強度，非常容易造成滑坡、泥石流等地質災害，這直接影響到水利樞紐工程的運行安全[2-4]。

因此，評估不同水位時涉水岸坡巖體的各項非線性特征顯得尤為重要。極限平衡法包含瑞典圓弧法、費倫紐斯(Fellenius)法和簡布(Janbu)法等[5-7]。而數值分析法包含了有限元法和離散元法等[8-9]。剛體極限平衡法由于僅能計算出穩定安全系數和滑弧，難以分析邊坡內部的應力應變特征而受到一定的限制。而數值分析法，按相同比例折減巖體內摩擦角和黏聚力[10]，理論上存在缺陷。本文結合了數值有限元分析法，提出采用重度增加法計算分析某涉水岸坡的穩定性，旨在為該邊坡工程的安全評估提供參考和借鑒。

1 研究方法

1.1 滲流計算原理

有限元法滲流計算基本機理：根據滲流作用力概念，將作用在滑動面上和劃分土塊表面的水壓力轉換為等價的體積力。把各節點水頭值換算成各單元滲透力，不需要考慮各單元體接觸面邊界上的孔隙水壓力和邊坡外水壓力。滑坡體中的地下水位隨水庫、河谷水位升降而變化，因而在滑坡體中就會形成飽和區和非飽和區。非飽和區土壤水與飽和區地下水運動密切相關，因此，由于水流運動導致邊坡巖土體狀態改變即為飽和與非飽和滲流問題，其滲透各向異性的飽和與非飽和滲流控制方程如式(1)、式(2)所示[11-13]：

(1)

(2)

式中：x、y表示位置坐標；Kx、Ky分別為巖土體的水平和垂直方向的滲透系數；ρw為水的密度；mw為釋水系數；Q為邊界流量；θw為體積含水率，ua為孔隙氣壓力，uw為孔隙水壓力，h為壓力水頭，t為時間。

1.2 穩定計算原理

計算工具采用巖土有限元軟件PLAXIS，對該涉水岸坡的應力應變和穩定性進行定量計算和分析。該軟件包含豐富的巖土體本構模型，可以應用于各類計算工況[14]。

鑒于當前邊坡穩定性方法均存在一定的理論缺陷問題，本文中邊坡穩定性分析時結合有限元數值分析法，采用有限元重度增加法[15-16]來分析邊坡的穩定性。由于有限元重度增加法僅需要梯度增加重度大小數值，然后不斷循環迭代計算，直到有限元計算不收斂或者節點位移突變作為判別條件，此時重度數值增加的倍數即定義為邊坡的安全系數，相比于強度折減法更為可靠，重度增加法計算式如式(3)所示：

(3)

式中：Fs為重度增加法定義的安全系數，γ為自然條件下邊坡巖土體的重度，γ′為根據梯級倍數增加后達到臨界狀態時邊坡巖土體的重度。

2 結果與分析

2.1 計算模型前處理

本次研究的邊坡對象位于西藏江達縣境內，由于含泥量較高，屬于強風化山體，山體的結構主要包含土、碎石和塊石。該邊坡頂部高程約3650 m，坡腳高程約2930 m，邊坡坡度平均33°。岸坡山體表層為碎石土，最大厚度超過50 m，邊坡基礎部分為弱風化巖體。該涉水邊坡剖面及網格劃分情況見圖1。參照巖土試驗結果，獲得了該岸坡碎石土和弱風化基巖的有限元計算力學參數，見下表1所示。

表1 涉水岸坡模型有限元計算力學參數

采用有限元重度增加法計算前，設置了不同水位的計算工況，分析了三種不同庫水位，高程分別為2950 m、2930 m、2910 m。本次計算主要考慮穩定滲流工況，不計算非穩定滲流，即認為巖土體內部水位隨著外面水位變化，與之保持一致，不存在滯后的現象。采用重度增加法計算邊坡穩定性時同時涉及表層碎石土和弱風化基巖的重度參數改變。巖土體材料模型選用摩爾庫倫模型，從而可以考慮涉水岸坡的非線性特征。

2.2 岸坡應力應變分析

分析計算了庫水位2950 m、2930 m、2910 m時該涉水邊坡的應力應變特征。其中，應力為剪切應力，應變為剪切總應變。由于坡體強度主要由抗剪強度進行反應和體現，因此，邊坡應力應變分析以剪應力和總應變為主。圖2為不同水位時邊坡的剪應力分布情況，分析可知，三種情況下山體巖體內部的剪應力分布情況比較一致，表層碎石土整體出現應力松弛現象，剪應力非常小，介于200～-200 kPa之間，正值受拉，負值受壓。當水位為2950 m和2930 m時，在高程3430 m附近出現了較大的正值剪應力，大約1900 kPa，此處巖層分界坡度較陡，極易出現斷層錯位。當水位為2910 m時，在高程3430 m處的剪應力相對減小，僅為1800 kPa。三種水位下剪應力的最大值均位于模型底部居中位置，越往山體內部走剪應力數值越小，這對該岸坡的表層滑動影響很小，基本可以忽略不計。以上分析說明，在不同蓄水位時涉水岸坡的應力狀態差別很小，且表層碎石土應力松弛現象較為明顯，當在外界自然條件作用下，在坡度較陡的位置容易發生局部垮塌，應引起關注。

圖2 不同水位時涉水岸坡剪應力分布

圖3為不同水位時邊坡的總應變分布情況，陰影部位范圍大小反映了總應變大小。分析可知，三種水位情況下山體巖體內部的總應變分布情況差別較大。庫水位為2950 m時總應變分布范圍最大，主要分布在巖層分界線兩側，剪應變集中高程在3380～3480 m，3380 m以下往坡腳位置總應變范圍逐漸減小。當庫水位為2930 m時此處的總應變分布量相比2950 m時大大減小，但坡體中部基本沒有發生應變。當庫水位為2910 m時，在高程3380～3480 m的總應變分布量與2930 m時相差不大，但坡腳位置附近的剪應變基本消失，中部同樣也沒有發生剪切應變。分析說明，在不同蓄水位時涉水坡體內部的剪應變分布范圍區別相差較大，主要體現在涉水邊坡的中上部位、巖層分界線坡度較陡的位置。

圖3 不同水位時涉水岸坡總應變分布

2.3 岸坡穩定性分析

基于數值有限元分析基本原理，采用重度增加法計算了不同水位時該岸坡的穩定安全系數，見表2。庫水位在2950 m時，岸坡的安全系數僅為1.145，穩定性較差。庫水位2930 m時，邊坡穩定性系數為1.248，庫水位為2910 m時，邊坡的安全系數最大，為1.259。在穩定滲流條件下，蓄水位越高，涉水邊坡的穩定性反而越差，因為蓄水位越高，水體浸泡的邊坡巖體也越多，導致巖體抗剪強度降低。同時，根據重度增加法的計算結果，導出了三種水位時涉水岸坡的最不利滑裂面，發現不同水位時的邊坡滑裂面完全相同，滑裂面起于邊坡臨空面頂部，內部與巖層分界線相切，剪出口位于坡面中下部。圖4僅給出了蓄水水位2930 m時該涉水岸坡的最不利滑裂面。

表2 涉水岸坡穩定安全系數

圖4 岸坡水位2930 m時最不利滑裂面

3 結 論

本文針對西藏某涉水岸坡的穩定性問題，以有限元軟件PLAXIS為計算平臺，結合數值有限元分析法，提出采用有限元重度增加法分析邊坡在不同蓄水位時的巖土體應力應變特征，并且定量計算不同水位工況時的岸坡穩定安全系數，得到以下結論。

(1)不同蓄水位時涉水岸坡的剪應力狀態差別很小，僅在巖層分界線處存在局部差異。

(2)不同蓄水位時坡體內部剪應變分布區別相差較大，蓄水位越高，剪應變分布量越大，在最高水位2950 m時剪應變集中在高程3380～3480 m的巖層分界線處。

(3)當山體內部處于穩定滲流狀態時，蓄水位越高，邊坡的穩定性越差，最高水位和最低水位安全系數相差可達10%。

(4)三種水位時的邊坡滑裂面完全相同，起于邊坡臨空面頂部，內部與巖層分界線相切，剪出口位于坡面中下部。

證明了該岸坡在高水位時穩定性欠佳，遇到暴雨地震等特殊工況時存在滑動的可能，且滑動模式屬于高位滑坡，因此該涉水岸坡的治理應重點關注坡體中上部位。