邊坡穩定性分析中的有限元極限平衡法＊

郭子儀 范振華 朱云升 張謝東 曾格華

（內蒙古交通設計研究院有限責任公司1） 呼和浩特 010010） （武漢理工大學交通學院2） 武漢 430063）

0 引 言

當前，邊坡穩定性分析方法主要有2類：一類為建立在剛體極限平衡理論上的極限平衡法；另一類為以有限元法為代表的數值分析方法．極限平衡法能給出物理意義明確的邊坡穩定安全系數以及可能的破壞面．以有限元為代表的數值分析方法不但能考慮坡體變形對其穩定性的影響，而且能通過邊坡內應力應變分布分析邊坡失穩過程等，但很難給出明確的穩定安全系數以及可能的破壞面．邊坡應力應變一般是按未破壞時的邊界條件計算出來，而實際上邊坡內任何超過抗剪強度或抗拉強度的應力狀態都是不能穩定的，一旦發生局域破壞，應力將會重新調整，邊坡穩定性安全系數也隨之而有所變動．因此，簡單采用傳統的極限平衡法或者有限元法都不能真實模擬邊坡失穩破壞過程中穩定安全系數的變化［1－2］．

本文將傳統極限平衡方法和有限元方法相結合，提出有限元極限平衡法用于邊坡的穩定性分析，采用有限元方法計算邊坡巖土體內的應力應變分布，求出邊坡潛在滑動面上的應力，通過極限平衡法的概念，求出該潛在滑動面對應的穩定安全系數，通過邊坡潛在滑動面上的下滑力和抗滑力在邊坡失穩破壞過程中的變化規律探尋淺變質巖邊坡變形破壞機理．

1 邊坡失穩破壞有限元模擬原理

邊坡發生失穩破壞并不是瞬間發生整體滑動或劇滑失穩破壞，而是一個由局部破壞逐漸擴展以至貫通形成滑面的漸進蠕滑過程．當邊坡呈現局域性的集中應力超過材料的強度值時會發生局部破壞，而后發生應力釋放、應力轉移和應力重新調整．在破壞區域的鄰近區域所受到的影響最大，該鄰域可能由原先沒有超過強度值轉變為超過強度值而發生破壞，并進行應力釋放，又把多余的載荷轉加到其他區域．在不斷地發生應力釋放、轉移和調整過程中，邊坡破壞面也不斷地延伸，最終將有2種可能性：（1）是破壞面完全貫穿，滑體將在滑床上開始作加速運動；（2）是破壞面沒有貫穿，在伸展到某一區域后就停止擴展，其前方區域的應力應變均未超過強度值．前一種情況，顯然安全系數小于1，而在后一種情況安全系數大于1．由于應力調整的結果，按最終應力狀態計算出的安全系數必然與按初始邊界條件進行應力分析計算所得的安全系數不一致，按最終應力狀態計算出的安全系數更能反映邊坡穩定性的實際狀況［3－5］．

模擬邊坡失穩破壞過程在數學處理上不易實現，采用邊坡滑動面上的法向和切向彈簧模擬滑坡體和滑床之間的接觸摩擦問題，通過彈簧剛度變化模擬滑動面上拉應力開裂和巖土體剪應力超過其抗剪強度導致滑動面上應力重分配．將邊坡滑動面上切向彈簧力超過巖土體抗剪強度的區域定義為剪切塑性區，法向彈簧為拉力的區域定義為拉裂區，因邊坡潛在滑動面上的應力重分配會導致拉裂區和塑性區在邊坡失穩破壞過程中逐步擴展，通過邊坡潛在滑動面上塑性區的變化發展趨勢模擬邊坡失穩破壞過程．

2 有限元極限平衡法分析邊坡穩定性的實現

采用傳統極限平衡方法和有限元分析相結合的方法分析邊坡穩定性最關鍵的問題包括3個方面：（1）淺變質巖邊坡潛在滑動面的確定；（2）模擬潛淺變質巖邊坡蠕滑過程中潛在滑動面上應力重分配；（3）如何根據潛在滑動面上的應力計算邊坡穩定安全系數［6－7］．

2．1 計算邊坡的選擇

某高速公路JK104＋820～SK104＋840段最大路塹邊坡高度為37m，該路塹左側邊坡完成邊坡削方施工后，左側中線105m附近出現山體沿路線前進方向開裂，裂縫左側方向沿山間凹地近似垂直路線方向延伸，右側沿與路線近似37°方向向下發展，滑坡界限清楚，鋼筋混凝土框架局部斷裂，框架整體發生下滑，觀測顯示裂縫不斷向兩側及山頂發展．

根據深層側向位移監測成果，該滑坡后緣拉裂縫與中線水平距離最長192m，滑體平均厚度21．4m，最厚處達到39．5m，平均寬度為85m，滑體面積約為3萬m2，滑體體積約為51萬m3，屬于大型滑坡．滑動面發育在強風化與中風化的接觸面附近，主要物質成分為褐黃色－黃色強風化硅化絹云母板巖，裂隙極發育，其間充填粘土，鉆探巖芯隨含泥量的不同而呈泥柱狀和砂狀，泥柱狀天然狀態呈可塑－軟塑狀．

2．2 邊坡潛在滑動面的確定

邊坡潛在滑動面的確定可根據強度折減法結合現場地質勘探資料和深部位移測斜資料擬合分析得出［8－14］．邊坡巖土體材料強度系數折減后，經過非線性有限元分析，此時邊坡內將出現一塑性區，見圖1，塑性應變等值線示意圖見圖2．

圖1 主滑動面塑性區分布圖

圖2 最危險滑動面示意圖

滑坡塑性區主要出現在強風化層和弱風化層交界面上，屬于深層滑坡．以計算坡面貫通的塑性區作為潛在滑動面，結合邊坡鉆探地質資料，以強風化層、軟弱夾層為控制點采用擬合法對計算潛在滑動面進行修正，得出合理的潛在滑動面．

2．3 滑動面上應力重分配及穩定安全系數

本文在有限元分析過程中對淺變質巖邊坡潛在滑動面上施加切向和法向彈簧模擬滑坡體和滑床之間的接觸摩擦問題．假定滑動面上法向彈簧單元出現拉應力即認為該點出現拉裂縫，此時切向彈簧的應力即為該點抗剪強度，計算過程中通過改變彈簧剛度將超過土體抗剪強度的那部分剩余應力平均分配到其它彈簧上，由此模擬滑坡體在滑動過程中潛在滑動面上應力變化．在此過程中，大于最大切向力的彈簧力得到逐步釋放，剩余力由其他彈簧承擔，通過反復迭代計算即可實現淺變質巖邊坡在漸進破壞過程中潛在滑動面上的應力重分配．將邊坡滑動面上的塑性區基本貫通時各節點的法向彈簧力和切向彈簧力代入傳統極限平衡法邊坡穩定安全系數計算公式中計算邊坡穩定安全系數，見式（1）．這樣，一方面可以考慮邊坡漸進破壞過程中內部應力應變狀態對邊坡穩定性的影響，另一方面可以利用極限平衡理論邊坡穩定安全系數明確定義計算邊坡穩定安全系數．

式中：Fni和Fτi分別為有限元計算得到的破壞面上第i單元的法向力和切向力；Li為第i單元的長度．

2．4 邊坡主滑動面有限元模型建立

以邊坡主滑動方向上的剖面作為主滑動剖面建立二維有限元模型，在滑動面每個節點上施加2個彈簧單元，其中1個彈簧單元與該節點處滑動面的切線方向相同，另1個彈簧單元與該節點處滑動面的法線方向相同，可以通過潛在滑動面上相鄰兩節點的坐標計算出這2節點中任意1個節點的切線方向和法線方向．在邊坡潛在滑動面上某節點施加的切向彈簧單元應力模擬該點滑動面上的剪應力大小，法向彈簧單元應力模擬該點潛在滑動面上正應力大小．滑動面每個節點上的2個彈簧單元一端與節點相連，另一端的位移自由度和轉動自由度全部施加約束．通過這兩個彈簧單元的位移大小可以較好地模擬潛在滑體和滑床之間的相對滑動趨勢．

2．5 材料參數和荷載的確定

本文計算材料抗剪強度參數主要根據室內試驗結果和參數反演法結果綜合取值確定，根據潛在滑動面以上坡體巖土層處于天然狀態時不同厚度進行加權平均得出潛在滑動面以上巖土體綜合容重、彈性模量和泊松比，邊坡荷載僅考慮潛在滑動面以上滑體的自重作用，根據滑體綜合容重來計算自重，計算參數見表1．

表1 計算參數表

3 計算結果分析

3．1 邊坡應力和變形分析

通過有限元計算可以得出主滑動剖面第一主應力云分布圖，合成位移圖以及剪應力云圖，分別見圖3～5．同時，可以得到主滑動面內拉應力區分布和關鍵點拉應力值，見表2，表中所列拉應力區是指從起點到終點連續出現拉應力的區域．

表2 主滑動坡面拉應力值分布及其高程

圖3 主滑動剖面第一主應力云

圖4 主滑動剖面合成位移云

圖5 主滑動剖面剪應力云

從上述表中邊坡在自重作用下的第一主應力計算結果來看，主滑動剖面的表面和潛在潛在滑動面上不同高程處均會出現連續拉應力區，拉應力基本上分布在邊坡中后緣部位，由于邊坡巖土體抗拉強度往往較低，在這些拉應力作用下，邊坡后緣往往會率先產生拉裂縫．潛在滑動面拉裂區平均拉應力均小于坡面拉裂區平均拉應力，表明邊坡拉裂區均是由坡面產生，隨著滑體滑動逐步向坡體內部擴展．

若遇降水，滑坡土體處于一定程度的飽水狀態，除土體自重增加以外，土體的抗剪強度也要下降，滑坡體的穩定性降低，上述拉應力區的任何部位都有出現拉裂縫的可能性．邊坡坡體內一旦產生拉裂縫，裂縫會伴隨著滑坡體滑動過程逐步擴展，加速滑坡體下滑過程，最終在坡體內形成連續剪切破壞面，導致邊坡失穩破壞．邊坡上部向下蠕滑過程中對坡腳部分巖土體起到推擠作用，使得坡腳位置承受較大的壓應力，因此，邊坡下部主要承受壓應力作用．

計算剖面最大合成位移出現在人工開挖臨空面頂部和坡體后緣之間，人工開挖臨空面上出現擠出變形．因此，人工削方臨空面直接影響滑體位移走向，臨空面過高過陡會導致過大位移產生，坡體除了沿潛在滑動面產生主滑動破壞外，還有可能由于臨空面位移過大產生淺表層次生滑動．

從計算剖面剪應力云圖分布來看，坡體內最大剪應力出現在潛在滑動面上，沿著潛在滑動面方向邊坡的剪切應力逐漸增大，尤其是坡腳部位剪應力較大，當潛在滑動面上剪切應力超過了該點的抗剪強度，則在該點產生剪切破壞，潛在滑動面上出現剪切破壞點后隨著滑坡體滑動過程的逐步發展，潛在滑動面上發生應力重分配，破壞點逐步增多形成連續破壞面．

3．2 滑動面受力特征及穩定安全系數分析

考慮邊坡潛在滑動面上部分區域剪應力超過巖土體抗剪強度而產生應力重分配的受力特點，進行主滑動面應力重分配的迭代計算，可以得出失穩邊坡潛在滑動面上下滑力和抗滑力的變化規律以及邊坡穩定安全系數變化規律．

計算剖面潛在滑動面經過每次迭代計算過程均發生了應力重分配，導致潛在滑動面上塑性區隨著迭代次數增加而逐步擴展，直至基本貫通，迭代計算過程中，邊坡各計算剖面上抗滑力和下滑力也隨之發生變化，邊坡穩定安全系數逐步減小，滑坡體產生完全失穩破壞．潛在滑動面上抗滑力、下滑力以及穩定安全系數變化見表3．

表3 計算剖面潛在滑動面受力表

通過潛在滑動面上的法向和切向彈簧模擬滑坡體和滑床間的接觸摩擦問題，采用本文有限元極限平衡法對邊坡進行計算分析可知：

1）計算剖面潛在滑動面上出現明顯分段，潛在滑動面上邊坡后緣頂部產生拉應力集中，出現“卸荷拉裂段”，邊坡前緣坡腳處或者人工開挖臨空面上產生剪應力集中或者邊坡坡頂后緣由于坡頂拉裂產生蠕滑而出現剪切塑性區，定義為“蠕變剪滑段”，邊坡中部應力接近而尚未超過該處巖體抗剪強度，僅有變形而未剪滑，為邊坡穩定區，該段對保持斜坡整體穩定起到鎖固作用，定義為“鎖固段”．隨著各計算剖面迭代計算次數增加，潛在滑動面上“卸荷拉裂段”、“蠕變剪滑段”和“鎖固段”的分布長度在逐步改變，“剪變蠕滑段”分布長度隨著潛在滑動面上應力重分配過程逐步增大，超過巖土體抗剪強度的那部分剪應力逐步向“鎖固段”轉移，當鎖固段巖體內的剪應力積聚到超過了該處巖體的抗剪強度時，“鎖固段”則逐步轉換成“剪變蠕滑段”，“鎖固段”分布長度則逐步減小，當“鎖固段”全部變成“剪變蠕滑段”后坡體剪切滑動面貫通，邊坡產生整體失穩．

2）邊坡初始穩定安全系數均等于1左右，處于臨界極限平衡狀態，潛在滑動面上均分布有長度不等的剪切塑性區．路塹邊坡開挖后，形成了較高的臨空面，改變了坡體受力狀態，由于邊坡處于斷層破碎影響帶內，強風化層巖石結構松散，風化裂隙發育，巖質較軟，坡體上水田內蓄積大量積水，地下水埋藏淺，致使巖土體抗剪強度降低．滑坡體裂縫已發展到寬2～20cm，長3～84m，滑坡周界裂縫仍在不斷發展，已基本連通．深層側向位移監測數據表明，滑坡體上不同部位的13個孔均有位移，滑坡體一直存在滑動速率．

3）從潛在滑動面上應力重分配的重復迭代計算結果可看出，各計算剖面的潛在滑動面上剪切塑性區均隨著迭代次數增加而逐步擴展，鎖固段長度逐漸減小，當鎖固段的剪切應變能賦存到與鎖固段的抗剪強度相等時，此時整個坡體處于穩定不平衡狀態，滑體沿著潛在滑動面作時而停止時而運動的緩慢移動．在這交替運動過程中，潛在滑動面巖體受到強烈的擾動，使得巖體的抗剪強度驟然下降，當鎖固段失穩斷裂后與坡體的后緣拉裂段和前緣的滑移段形成貫通面，滑坡體就形成了完整的滑動面，坡體緩慢地向邊坡前緣軟弱結構面尋找出滑口，前緣臨空面的巖體開始與母體脫離，滑向坡下，滑坡體產生完全失穩破壞．

4）隨著迭代計算次數的增加，潛在滑動面上下滑力基本保持穩定，抗滑力隨著潛在滑動面上剪切塑性區的擴展而逐步減小，穩定安全系數也隨之減小，滑坡體逐步產生蠕變滑動．

3．3 計算結果對比分析

按照上述傳遞系數法的基本原理對王家寨和南約溝滑坡各計算剖面進行穩定安全系數進行計算，結果見表4．

表4 邊坡穩定安全系數對比表

從表中可以看出，傳遞系數法由于沒有考慮滑坡體內部和應力狀態，計算的主滑動剖面穩定安全系數大于有限元極限平衡法計算結果。該邊坡本來屬于潛在不穩定體，處于極限平衡狀態，如果采用傳遞系數法計算的穩定安全系數進行穩定性評價，則可以認為是穩定體．采用有限元極限平衡法計算各剖面的穩定安全系數，考慮了邊坡潛在滑動面上的應力在邊坡蠕變滑動過程中的重分配問題，這種方法計算的穩定安全系數更符合邊坡的實際情況．

4 結 論

1）采用彈塑性有限元強度折減和淺變質巖邊坡現場地質勘探資料相結合的方法確定淺變質巖風化層邊坡潛在滑動面是合理可行的．

2）采用淺變質巖風化層邊坡潛在滑動面上的法向和切向彈簧單元模擬滑坡體和滑床之間的接觸摩擦問題能夠真實反映淺變質巖風化層邊坡漸進破壞的基本過程和特點．

3）淺變質巖風化層邊坡潛在滑動面在邊坡漸進破壞過程中發生應力重分配，坡破壞從邊坡后緣拉裂開始，隨著潛在滑動面應力重分配過程的發生，潛在滑動面上的剪切塑性區逐步擴展，直至貫通，淺變質巖風化層邊坡主要破壞模式為漸進蠕滑拉裂破壞．

4）通過具體工程實例分析可知，本文提出的有限元極限平衡法分析邊坡穩定性是一種有效可行的方法．

［1］黃昌乾，丁恩保．邊坡工程常用穩定性分析方法［J］．水電站設計，1999，15（1）：45－48．

［2］朱云升，李洪亮．道路滑坡穩定性分析有限元方法研究［J］．公路交通科技，2007，24（4）：55－59．

［3］戴自航，盧才金．邊坡失穩機理的力學解釋［J］．巖土工程學報，2006，28（10）：104－107．

［4］江學平，于遠忠．邊坡漸進失穩破壞的概率分析［J］．西南科技大學學報，2007，22（2）：78－82．

［5］LIU Y Q，LI H B，DAI H C，et al．Preliminary study on the progressive failure of a layered rock slope under explosions［J］．Key Engineering Materials，2006，306－308：1461－1466．

［6］朱典文，唐小兵，謝支鋼．邊坡穩定性分析的數值解研究［J］．武漢理工大學學報：交通科學與工程版，2003，22（4）：856－859．

［7］EBERHARDT E，STEAD D，COGGAN J S．Numerical analysis of initiation and progressive failure in natural rock slopes－the 1991Randa rocksfide［J］．International Journal of RockMechanics and Mining Sciences，2004，41：69－87．

［8］DAWSON E M，ROTH W H，DRESCHER A．Slope stability anaysis by strength reduction ［J］．Geotechnique，1999，49（6）：835－840．

［9］鄭穎人，趙尚毅．有限元強度折減法在土坡和巖坡中的應用［J］．巖石力學與工程學報，2004，23（19）：3381－3388．

［10］MATSUI T，SAN K C．Finite element slope stability analysis by shear strength reduction technique［J］．Soils and Foundation，1992，32（1）：59－70．

［11］靳付成．巖質邊坡失穩判據及臨界滑動面確定方法研究［D］．南京：河海大學，2007．

［12］張金保．推移式滑坡對施工中橋梁的影響以及處理對策［J］．青海交通科技，2013（4）：55－58．

［13］岳尚全．巖體穩定和區域穩定數值模擬邊界條件研究方法［J］．巖石力學與工程學季，1999（2）：262－266．

［14］張魯渝，時衛民，鄭穎人．平面應變條件下土坡穩定有限元分析［J］．巖土工程學報，2002（4）：487－490．