有限差分強度折減法求解邊坡安全系數

夏世友, 張電吉, 夏 亮, 曾夢秋, 王立光

(1.武漢工程大學環境與城市建設學院，湖北 武漢 430074；2.武漢市城鄉建設委員會，湖北 武漢 430015)

0 引 言

崩塌、滑坡和泥石流是山體在眾多因素作用下變形、失穩而給人類造成生命財產損失的地質現象.滑坡地質災害雖然不像地震災害那樣集中和強烈，但因其分布面比較廣、發生次數頻繁、影響時間長，所以災害的累計損失也比較大.當前面臨著嚴峻的地質災害形勢以及大量災害治理資金的投入，使得邊坡穩定性分析成為地質災害治理過程中的緊迫問題.許多學者在邊坡穩定性分析這個問題上，都提出了有效的方法，比如基于應力場和位移場的容重增加法、強度儲備法、能量法及各種搜索滑面的邊坡穩定性評價方法[1-3]，但這些方法大多數都比較復雜，不易被工程人員所掌握.而本文中所提到的強度折減法較為簡單，并且計算結果精度比較高，可以作為工程人員在處理邊坡穩定性問題時的首選之法.

1 強度折減法的理論背景

三維快速拉格朗日法(Fast Lagrangian Analysis of Continua)是由美國Itasca咨詢公司研發推出的連續性介質力學分析軟件，簡稱FLAC/FLAC3D.其理論基礎是數值分析方法中的顯式有限差分法，在模擬巖土或者其他材料的三維力學特性時運用比較多.

顯式有限差分方法在數值計算中不需要像有限單元法那樣將每個單元組裝成龐大的總剛度矩陣，其突出特點是計算速度快，存儲量小，不存在計算收斂性問題.目前，在連續性介質力學問題的范疇內，顯式有限差分法已經被證明具有與有限單元法同樣準確的計算結果[4].而且這個方法能夠很好地模擬巖土材料在達到強度極限或屈服極限時發生塑性流動的力學行為.FLAC和ANSYS等多場耦合有限元軟件相比較，在分析漸進破壞、失穩和模擬大變形等問題上有明顯的優勢，原因在于該軟件提供了各種可用于巖土體和支護體系的本構模型和結構單元，這些更加突出了FLAC的“專業”特性，因此在國際巖體工程界該軟件享有盛譽[5].

強度折減法理論是由Duncan提出的.該理論已經有比較長的歷史，在20世紀80年代國外就已經開始應用，但是因為力學尚處于起步階段和缺乏相應的失穩判據，一直以來在巖土工程界并沒有得到廣泛應用.隨著計算機性能的逐漸改善提高以及其他成熟商用軟件的推出，強度折減法受制于計算程序和計算精度的局面得到了根本性改變.20世紀90年代末重新成為巖土工程數值模擬研究的一個熱點.

2 強度折減法的基本原理

強度折減法中的邊坡穩定的安全系數定義[6]為：使邊坡剛好達到臨界破壞狀態時，對巖、土體的抗剪強度進行折減的程度，即定義安全系數為巖土體的實際抗剪強度與臨界破壞時的折減后的剪切強度的比值.強度折減法的要點是利用公式(1)和(2)來調整巖土體的強度指標C和Φ(式中，Cf為折減后的粘接力，Φf為折減后的摩擦力，Ftrial為折減系數)，然后對邊坡穩定性系數進行數值分析，不斷地增大折減系數，反復計算，直至達到臨界破壞狀態，此時所得的折減系數即為安全系數FS.

(1)

(2)

3 強度折減法分析邊坡穩定性的步驟

根據邊坡穩定性安全系數的定義，運用彈塑性理論及FLAC3D這個軟件求解邊坡穩定性安全系數的迭代步驟如圖1.

圖1 邊坡安全系數的迭代步驟圖Fig.1 Slope safety factor of iterative steps figure

在采用圖1迭代求解邊坡穩定性安全系數時，首先，確定好本構模型、屈服準則、邊界條件和材料參數等；其次，利用FISH語言來編寫計算模型的命令流程并建模；再次，設定邊坡穩定性安全系數的計算精度和上下限值；最后，利用軟件進行計算處理從而獲得邊坡穩定性安全系數和有關變量的等值陰影圖.

4 邊坡的失穩判據

目前，判斷邊坡是否失穩還沒有一個統一的標準，即安全系數求解過程的終止條件不統一.我們通過歸納總結得出的邊坡失穩判據主要有：

a.以數值分析計算不收斂性作為標準[7-9].目前，國內外學者經常采用數值分析計算結果的收斂性來判定邊坡是否發生破壞.但由于每個人所編寫的有限元程序的不同、迭代次數不同以及設置的收斂容許值大小不同，往往會導致計算的邊坡穩定性安全系數有所不一樣.此外，即使得到的數值計算分析結果收斂，這也不代表邊坡一定處于安全狀態[10].

b.以邊坡內某特征點位移的突變為標準[11].通過工程實際經驗發現任何邊坡的失穩破壞都伴隨著一定的位移變化，當邊坡內的某特征點位移發生突變時，我們就可以認為該邊坡處于臨界失穩狀態.但是，該準則受特征點所選位置以及位移容許限值的影響比較大，所以運用起來不方便.

c.以塑性區從坡腳至坡頂貫通為標準[12-13].通俗的說法是，當邊坡的塑性應變區從坡腳貫通至坡頂時,則認為邊坡發生了失穩破壞.但是,塑性區的貫通并不意味著該區域內的所有塑性破壞點都貫通了，往往還存在著局部塑性應變發展滯后區域，這些區域還可能進一步發展.文獻[14]認為：塑性區的貫通只是邊坡破壞的必要條件而不是充分條件，即在塑性區貫通的同時,還要參考其他的判據，比如，該邊坡是否產生了很大的并且無限發展的塑性變形或者是該邊坡內某特征點位移是否產生了突變.

綜合考慮上面的分析，可以發現這三種失穩判據都存在一定的局限性.判據a、b受人為因素影響較大，判據c無法單獨地判斷邊坡是否失穩，因此在下面的算例分析中，將聯合采用塑性區從坡腳貫通至坡頂和速度矢量圖來判斷邊坡是否失穩.因為當分析的結果中出現貫通的塑性區時，僅僅能表示該邊坡的某些區域已出現了塑性破壞，但是無法明確地判斷邊坡是否處于危險狀態，而速度矢量圖不僅可以明顯地看見貫通的塑性區，而且可以看見滑動面的外側區域各網格點的速度明顯大于其他區域，因此可以確定該邊坡已經發生了破壞.

5 運用強度折減法計算邊坡穩定算例

某均質巖石邊坡坡高12 m，巖體容重γ=22 kN/m3,黏聚力c=12 kPa,內摩擦角φ=20°,坡角β=45°,彈性模量E=6 MPa,泊松比μ=0.3.該勻質巖石邊坡斷面圖如圖2.

圖2 勻質巖石邊坡斷面圖Fig.2 Homogeneous rock slope section in figure

利用FLAC3D軟件計算，在計算中，a.建立的計算模型按下列要求進行：將所使用的坐標系的垂直方向定為Y向，水平方向定為X方向，本構模型選用Mohr-Coulomb本構模型；b.計算邊界條件設置如下：邊坡的左邊界和右邊界作為X方向的約束，即橫向位移；底部邊界作為Z方向的約束，即縱向位移；整體約束Y方向，即平面法線方向位移.材料參數值按算例所給值設置；邊坡穩定性安全系數的上、下限值各設置定為2和0，計算精確度設定為0.000 1；同時，在數值計算過程中只考慮土體自身重力作用，計算結果如圖3.

從圖3中可以明顯看到塑性貫通區域，即坡體的潛在滑動面；同時該速度矢量圖也有力地證明了該邊坡外側區域已明顯出現了滑動，即坡體發生了失穩破壞，此時的安全系數Fs=1.008 71，因為滑動面外側區域各網格點的速度明顯大于其他區域.

圖3 剪切應變增量云圖及速度矢量圖Fig.3 Shear strain increment imageries and velocity vector

6 結 語

綜上所述，筆者發現在分析邊坡穩定時，運用以強度折減法作為理論基礎的顯式拉格朗日有限差分法，具有以下優點：a.計算簡單方便、結果可靠；b.避免了極限平衡法需要事先假設滑動面以及人為假定等不足；c.將強度折減法原理和FLAC3D相結合，形象的描繪出了邊坡的變形過程及其滑動面形狀.

參考文獻：

[1] 徐衛亞，周家文，鄧俊曄，等.基于Dijkstra算法的邊坡極限平衡有限元分析[J].巖土工程學報，2007，29(8)：1159-1172.

[2] 鄭宏，李春光，李焯芬，等.求解安全系數的有限元法[J].巖土工程學報，2002，24(5)：626-628.

[3] 鄭宏，劉德富，羅先啟.基于變形分析的邊坡潛在滑面的確定[J].巖石力學與工程學報，2004，23(5)：092716.

[4] 陳小康.露天坑下殘礦回收安全控制技術研究[D].長沙：中南大學，2010.

[5] 黃貫虹，方剛.系統工程方法與應用[M].廣州：暨南大學出版社，2005.

[6] 陳育民，徐鼎平.FLAC/FLAC3D基礎與工程實例[M].北京：中國水利水電出版社，2009：260-261.

[7] 趙尚毅，鄭穎人，時衛民，等.用有限元強度折減法求邊坡穩定安全系數[J].巖土工程學報，2002，24(3)：343-346.

[8] GriffithsD V, Lane P A.Slope stability analysis by finite elements [J]. Geotechnique,1999,49(3);387-403.

[9] 林杭，曹平，趙延林，等.強度折減法在Hoek-Brown準則中的應用[J].中南大學學報:自然科學版，2007，38(6)：1219-1224.

[10] 董璞，劉金龍，李亮輝.強度折減有限元法分析邊坡穩定性的精度探討[J].四川建筑科學研究，2009，35(2)：146-150.

[11] 宋二祥.土工結構安全系數的有限元計算[J].巖土工程學報，1997，19(2)：1-7.

[12] 鄭宏，李春光，李焯芬.求解安全系數的有限元法[J].巖土工程學報， 2002，24(5)：626-628.

[13] 欒茂田，武亞軍，年廷凱.強度折減有限元法中邊坡失穩的塑性區判據及其應用[J].防災減災工程學報，2003，23(3)：1-8.

[14] 趙尚毅，鄭穎人，張玉芳.極限分析有限元法講座-II有限元強度折減法中邊坡失穩的判據探討[J]，巖土力學，2005,26(2)：332-336.