極限平衡有限元在邊坡穩定性分析中的應用

楊宏宇，白偉明，鄒洪海

（青島海洋地質研究所，山東 青島 266071）

0 引 言

邊坡失穩是工程安全的熱點問題[1-2]。實際工程中，邊坡穩定性分析方法主要采用剛體極限平衡分析方法[3]。由于天然巖土體應力狀態受巖土介質變形與強度特性的共同控制，采用剛體極限平衡分析方法評價邊坡穩定性時又需對滑體條間應力進行假設，計算結果往往不盡人意。同樣，采用基于經典數學物理方程與計算機技術的數值分析方法，雖然可得到邊坡應力場，且滿足應力變形協調條件，但卻又無法給出邊坡的安全系數。因此，若能將兩者有效聯系，則既可克服剛體極限平衡分析不滿足應力變形協調原理的缺陷[4-5]，又可得到更為合理的安全系數。

邊坡穩定性計算的關鍵是確定出滑面并獲得滑面法向、切向力[6]。當采用剛體極限平衡分析時，不同的條分法計算結果間往往存在一定偏差[7]，且僅當條分底面法向和切向力與真實值較為接近時，計算結果才比較理想。為此，國內外許多學者嘗試通過結合優化理論或基于有限元分析方法來提高計算的有效性，如動態規劃法[8]、遺傳算法等，均有效克服了安全系數局部最優的不足[9]。

本文以有限元數值計算所得的邊坡應力場為基礎，通過指定一個搜索區域，將其格柵化后，以格柵節點控制圓弧形滑移面圓心，不斷迭代試算，確定邊坡最危險滑移面位置及安全系數，然后將其與剛體極限平衡分析結果對比，以此驗證采用數值極限平衡有限元法分析邊坡穩定的適用性。

1 邊坡滑面搜索流程

邊坡失穩是邊坡體內部應力和強度變化的結果，故邊坡安全與其應力狀態密切相關。以有限元數值計算獲取邊坡應力場進行極限平衡分析，可彌補剛體極限平衡分析法不滿足應力變形協調原理的缺陷，使所得的邊坡最危險滑面及安全系數更為合理、準確。本文采用滑面圓心格柵化控制法確定最危險滑面，具體計算流程如下：

（1）預先指定滑面圓心搜索范圍，并對搜索范圍進行格柵化。

（2）依次選取格柵點作為滑面圓心，以坡腳可能發生剪出的范圍內任意一點作為滑面上一點，確定滑面及滑面方程。

（3）將步驟（2）中確定出的當前滑面離散化 （即按滑面穿過的單元大小將整個滑面劃分為小的弧段），以弧段中心為控制點搜尋每個弧段所穿過的單元位置，并根據單元節點應力計算出穿過該單元的弧段上的法向力和沿弧段的切向力。

（4）累計所有滑弧段的抗滑力與滑動力獲得滑面抗滑力與滑動力，從而計算出當前滑動面的安全系數，若當前計算結果小于先前計算值，更新全局安全系數及滑面位置。

（5）遍歷所有網格格柵點，即可確定出邊坡安全系數及最危險滑面位置。

2 工程實例

2.1 模型及參數

邊坡模型尺寸見圖1，斜坡坡腳為45°。為簡便模型構建，邊坡材料視為均質、各向同性的單一材料，其相應的材料參數見表1。有限元網格見圖2。模型由2 551個單元，5 338個節點構成。

圖1 邊坡幾何尺寸（單位：m）

表1 材料參數

圖2 有限元網格模型

2.2 計算結果

邊界約束采用位移約束，即固定模型左右及下邊界，然后應用經典的More-Columb屈服準則進行應力場的有限元數值模擬，所得主應力、應力分量見圖3。塑性區見圖4。

圖3 主應力、應力分量（單位：MPa）

計算結果表明，最大主應力－0.178 MPa，最小主應力－0.080 MPa，均表現為壓應力，坡腳位置應力變化較其他部位相對復雜；水平向 （X方向）和豎直向 （Y方向）應力變化表現出相同的規律；坡腳XY方向應力變化劇烈，出現剪切塑性區，故可斷定最危險滑面剪出口位置應該位于坡腳處，該結果同時可作為極限平衡穩定性分析滑面剪出口位置選擇。

圖4 剪切塑性區

2.3 穩定性分析

邊坡穩定分析的關鍵問題是確定出合理的最危險滑動面位置及相應的安全系數。許多情況下，在采用剛體極限平衡分析法進行邊坡穩定性分析時發現，搜尋出的邊坡安全系數小于1.0，而實際邊坡仍然處于穩定狀態，計算結果與實際情況矛盾，此時的剛體極限平衡分析結果具有局限性。本文將有限元計算所得單元應力、應變、位移等結果引入到求解邊坡臨界滑動面安全系數中，建立了數值計算與滑坡穩定評價的聯系，考慮了滑體條間應力假設的不足，克服了剛體極限平衡分析法不滿足應力變形協調的缺陷，可有效地確定出邊坡最危險滑移面位置及安全系數。

以有限元數值計算所得的應力場為基礎，應力變形分析結果為參照，即以坡腳位置作為剪出口，邊坡左上側區域格柵節點控制圓弧形滑面圓心，進行邊坡最危險潛在滑面位置搜索及安全系數的迭代計算，最終確定出邊坡最危險滑面位置及安全系數。剛體極限平衡法 （包括Bishop法、Janbu修正、Spencer法和M-P法）與極限平衡有限元法最危險滑面搜尋結果見圖5。由于剛體極限平衡計算結果非常相近，搜尋滑面重合，故圖5中剛體極限平衡滑面僅1條。將其與極限平衡有限元法搜尋結果對比發現，極限平衡有限元所得滑面位置較剛體極限平衡結果更靠近邊坡內，主要是由于剛體極限平衡分析法將土條視為剛體，沒有考慮邊坡土體在天然狀態下產生塑性變形時應力應變場的變化，從而導致了分析結果間的差異。

表2為采用不同的方法計算所得的安全系數。從表2可知，剛體極限平衡分析計算的安全系數均值為1.469，極限平衡有限元計算所得安全系數為1.516，兩者結果相近，具有一致性，說明數值極限平衡有限元法是有效的，可合理確定出邊坡最危險滑面位置及安全系數，特別是對已經產生塑性變形的邊坡將更具適用性。而對于土釘、錨桿、錨索等支護邊坡的穩定性計算及優化設計，采用極限平衡有限元計算分析亦可很好地解決其在復雜應力狀態下的安全分析與評價問題。

圖5 不同方法計算的最危險滑面位置

表2 不同方法計算的安全系數

3 結 論

（1）邊坡極限平衡有限元分析考慮了滑體條間應力假設的不足，克服了剛體極限平衡分析法不滿足應力變形協調的缺陷，可確定出邊坡最危險滑面位置及安全系數。

（2）對已產生塑性變形的邊坡進行安全性評價時，極限平衡有限元法比剛體極限平衡分析法更具適用性。

［1］ 王玉平，曾志強，潘樹林.邊坡穩定性分析方法綜述［J］.西華大學學報：自然科學版,2012，31（2）:101-105.

［2］ SHI Chong,ZHANG Yulong,ZHANG Qiang.Safety factor calculation method of anchored slope based on pseudo cohesion theory ［J］.Applied Mechanics and Materials,2012,170-173:453-456.

［3］ 陳祖煜.土質邊坡穩定分析原理·方法·程序［M］.北京：中國水利水電出版社，2003.

［4］ TURNER J P,JENSEN W G.Landslide stabilization using soil nail and mechanically stabilized earth walls： Case study［J］.J.of Geotech.and Geoenvirn.Eng.,ASCE,2005,131（2）:141-150.

［5］ COTTON P E,DAVID M,LUARK P F.Seismic response and extended life analysis of the deepest top-down soil nail wall in the U.S［J］.Geotechnical Special Publication,2004,124（3）:723-740.

［6］ 徐衛亞，周家文，鄧俊曄，等.基于Dijkstra算法的邊坡極限平衡有限元分析［J］.巖土工程學報， 2007， 29（8）:1159-1172.

［7］ 張鵬.巖土邊坡剛體極限平衡法的誤差根源與范圍研究［D］.西安：西安理工大學，2003.

［8］ 范鵬賢，朱大勇，郭志昆，等.動態規劃法求解邊坡安全系數最小上限解［J］.巖土工程學報， 2007， 29（3）:467-470.

［9］ 陳云敏，李育超，凌道盛.蒙特卡洛法與有限元結合搜索邊坡臨界滑動面［J］.巖土力學， 2004， 25（S2）:75-80.