翟 明,李 亮，*,王上上,劉 旭,陳 富，楊正權

(1.青島理工大學 土木工程學院，青島 266033；2.中國水利水電科學研究院，北京 100048)

在眾多巖土問題中，邊坡失穩事故尤為突出。作為巖土工程三大經典問題之一，邊坡穩定問題引起了國內外學者的廣泛關注。目前常用的邊坡穩定分析方法主要有三種，即極限平衡法(Limit Equilibrium Method, LEM)、強度折減法(Strength Reduction Method, SRM)和極限分析上限法(Upper Bound Limit Analysis Method)。極限平衡法[1-2]是最早用于邊坡穩定分析的方法，目前已趨于完善和成熟。極限平衡法主要包括兩個步驟：①計算給定滑動面的安全系數；②基于智能優化算法變換不同的滑動面并計算其安全系數，確定最小安全系數及其對應的臨界滑動面。強度折減法思想最早由ZIENKIEWICZ等[3]提出，后來被許多學者廣泛應用于邊坡穩定分析。強度折減法的基本思想是將土體初始抗剪強度參數(黏聚力c、內摩擦角φ)不斷進行折減，直至邊坡達到臨界狀態，此時的折減系數即為安全系數。極限分析上限法[4]的基本思想為：對于給定的失效模式，從變形協調出發，根據外力功和內能消散相平衡的原理確定相應的安全系數，然后應用最優化方法，確定相應于最小安全系數的臨界失效模式。上述三種方法在邊坡穩定分析中均有應用，對于復雜荷載工況、復雜邊坡剖面的邊坡而言，如何檢驗三種方法的適用性并比較其計算性能對于邊坡穩定分析具有較強的指導意義。

近年來，光滑粒子流體動力學(Smoothed Particle Hydrodynamics，SPH)方法[5-9]逐漸被用于巖土工程領域，BUI等[5]較早地將Druck-Prager本構模型引入SPH方法中，進行了一系列巖土工程問題的模擬，為SPH方法在巖土工程問題中的應用奠定基礎。黃雨等[6]采用Bingham流體模型描述土體的流動大變形，并進行了SPH數值模擬，結果與模型試驗實測數據基本吻合，表明SPH方法對于土體大變形分析的有效性。唐宇峰等[7]研究了SPH在邊坡穩定計算中的失穩判據，LI等[8]用SPH方法并結合極限平衡方法的安全系數模擬了邊坡臨界狀態時的位移云圖，在此基礎上評估了極限平衡方法的臨界滑動面。由此可見，SPH方法能夠模擬滑坡從發生、運動直至最終土體堆積的整個過程，其位移場云圖可以結合強度折減法塑性剪應變增量云圖，從而全面直觀地比較傳統分析方法的臨界滑動面。在簡要敘述邊坡穩定分析方法的基礎上，采用三個邊坡算例，分別比較了三種邊坡穩定分析方法的結果，得出了有意義的結論。

1 方法介紹

1.1 極限平衡法

極限平衡法(LEM)是最早應用于邊坡穩定性分析的方法。極限平衡法大多以垂直條分法為基礎，通過引入條間力函數假定使問題變得靜定可解，從而利用力和力矩平衡方程得出邊坡穩定安全系數。目前已有多種基于極限平衡理論的邊坡穩定分析方法，如瑞典法、Bishop法、Janbu法、Spencer法、Morgenstern-Price法等。由于傳統的極限平衡法計算簡便、概念清晰，現在仍然是工程師所采用的主要方法。本文極限平衡法采用的計算程序為GEOStudio中的SLOPE/W，計算方法選取常用的Spencer法。

1.2 強度折減法

(1)

式中：c，φ為土體的初始黏聚力和內摩擦角；c′，φ′為折減后土體的黏聚力和內摩擦角；F為折減系數。

本文采用強度折減法結合有限差分的數值計算程序FLAC軟件求解邊坡穩定安全系數，采用塑性剪應變增量云圖來查看滑動區域。

1.3 極限分析上限法

極限分析上限法有很多實現途徑，本文采用最新提出的不連續布局優化[10-11]方法(Discontinuity Layout Optimization，DLO)。不連續布局優化使用嚴格的線性優化技術識別分布在土體中節點之間的不連續discontinuity的關鍵布局，它可以識別塊體的滑動和平移，能夠利用數學優化技術識別關鍵的破壞滑移線。

DLO方法以節點來離散求解域，節點之間相互連接，形成多個不連續discontinuities(以下簡稱disc)，然后將土體參數賦予每個不連續disc，在荷載作用下，位移增大，外力做功也相應增大，當外力做功與系統的能量耗散相等時，達到極限狀態。采用數學優化技術規劃多條可能的路徑，對路徑中每一個不連續disc進行搜索，計算路徑的能量耗散，能量耗散最小的路徑即為破壞路徑。DLO方法確定破壞模式的過程如圖1所示，首先描述所要求解的問題(圖1(a))，然后以一定密度的節點離散所求解問題的幾何體(圖1(b))，節點之間相互連接形成不連續disc(圖1(c))，最后應用有效的數學優化技術得到破壞時的滑動面(圖1(d))。需要注意的是，DLO的計算過程中，所布設的節點、方式、密度、本構模型參數取值以及邊界條件的設定對結果都會有影響，本文所采用的計算程序為LimitState：GEO，采用其默認設置，網格密度取中等。

圖1 DLO方法確定破壞模式的過程

1.4 光滑粒子流體動力學方法(SPH)

SPH[5-9]是一種無網格的數值方法，與有限元等需要劃分網格的數值方法不同，SPH方法用有限個粒子來離散問題域，不需要網格劃分和小變形假設，這使得SPH可以不受網格畸變問題的影響，因此，SPH方法可以直接用于模擬巖土工程的大變形問題。SPH的概念圖如圖2所示，每個粒子被賦予狀態變量(含有速度、位移、應力、應變等)和材料屬性(含有重度、黏聚力、內摩擦角、彈性模量、泊松比等)。

圖2 SPH的概念

圖2中黑色實心圓圈所示的粒子i的信息通過核函數W與影響區域內臨近粒子k相互作用來獲得，與影響區域外的粒子無關，影響區域通過變量h來確定。SPH方法中粒子的運動信息由質量守恒方程(式(2))和動量守恒方程(式(3))確定。

(2)

(3)

式中：N為粒子總數；ρ為粒子中土體的密度；m為粒子質量；v為粒子的速度；t為時間；Wik為粒子i與粒子k間的光滑核函數；σ為單個粒子的總應力張量，其表達式由土體的本構模型確定；x為粒子位置；α，β為不同的坐標軸；fα為外力引起的加速度分量。

根據式(1)～(7)可得:第10個考核期Cpv為3.976萬元，Cav為4.704萬元，Ev為3.96萬元，Vs為-0.016萬元，Vc為-0.744萬元，Is為1.034，Ic為0.842。由計算結果可知:該項目在第10個考核期實際安全成本投入處于超支狀態，實際安全保障水卻沒有達到計劃水平，可以排除是由于安全成本節約導致了安全度的降低，很明顯是由于管理措施不到位而造成的，項目經理部應結合專家對當月安全評價體系的評分情況和項目實際安全投入情況進行深入分析，找出原因，加強安全管理力度的同時嚴格控制安全成本的支出，保證項目施工安全順利開展。

SPH方法模擬巖土工程問題時，算法參數取值對計算結果有重要影響。譬如，核函數W，影響半徑h以及邊界問題處理時所需的人工黏性比例因子等。本團隊基于SPH方法的基本原理與公式自編了Fortran程序，并進行了3年的調試，與前人的研究結果進行了對比，并發表了相關文章，可參閱文獻[8]。

2 SPH驗證

應用本團隊自編SPH程序來模擬砂柱垮塌試驗，模擬結果與HUNGR[12]所做的模型試驗結果進行對比，來驗證自編SPH程序的可靠性。砂柱高度為0.2 m，寬度為0.4 m，密度ρ=1630 kg/m3，黏聚力c=0 kPa，內摩擦角φ=30.9°。在SPH參數設置中，兩個主要的參數為粒子的間距與影響半徑h，MAO等在文獻[9]中對SPH方法的參數做了細致的研究，本文參數取值參考文獻[9]研究結果進行。將粒子半徑取為0.005 m，影響半徑取為0.006 m，共產生3200個粒子，分析模型如圖3所示。

圖3 砂柱試驗SPH分析模型

SPH方法的計算過程包含兩個步驟，首先模擬荷載作用下土的固結過程，固結完成以后，再模擬土體在荷載作用下的失穩滑動過程。SPH方法模擬砂柱垮塌后的形態與文獻[12]的試驗結果對比如圖4所示。

圖4 砂柱試驗的SPH模擬與試驗對比

圖4中，紅色線為砂柱垮塌試驗完成后的坡面曲線。由圖4可知，SPH方法模擬的結果與實驗室砂柱垮塌試驗結果表現出了良好的一致性，所以，本團隊自編SPH程序可以用于模擬滑坡過程中的大變形問題，其位移場云圖可以與強度折減法塑性剪應變增量云圖一起來綜合衡量邊坡穩定性問題。

3 算例分析

3.1 算例1

考慮圖5所示的簡單均質邊坡，模型寬度為20 m，總高度為10 m，其中坡高為6 m，坡角45°。假設土體材料服從Mohr-Coulomb強度準則，黏聚力c=20 kPa，內摩擦角φ=17°，重度γ=19 kN/m3。

圖5 算例1模型

對于算例1，分別使用極限平衡法、強度折減法、不連續布局優化法進行分析，得到的安全系數分別為1.623，1.674和1.742。其臨界滑動面和塑性貫通區如圖6所示。

圖6 算例1臨界滑動面及塑性剪應變增量云圖

通過圖6中三種方法結果的對比可知，極限平衡法和強度折減法所得安全系數基本一致，DLO所得安全系數略高；極限平衡法和不連續布局優化法得到的臨界滑動面非常接近，均落在強度折減法的塑性破壞區內(由塑性剪應變增量云圖確定)。所以，對于簡單的均質邊坡，三種方法所得臨界滑動面幾乎一致，安全系數略有差別。

應用SPH方法分析該問題，將粒子半徑定義為0.1 m，影響半徑取為0.12 m共產生3965個粒子，SPH分析模型如圖7所示。利用原始參數進行邊坡自重應力的模擬，進而將原始參數仿照強度折減法公式(1)進行折減，折減系數取為LEM所得安全系數，即F=1.623，將折減后的參數輸入SPH程序，得到計算穩定時刻的位移場云圖與LEM，DLO方法的臨界滑動面比較如圖8所示。

圖7 算例1的SPH分析模型

圖8 算例1的臨界滑動面及SPH位移場云圖

由圖8可以看出，SPH得到的位移場云圖明顯分為藍色的不動區與青色/黃色以及紅色所示的滑動區，二者之間的分界線與LEM，DLO所得臨界滑動面基本吻合，從而間接驗證了強度折減法所得塑性剪應變增量云圖與SPH位移場云圖在識別邊坡滑動面方面的一致性。

3.2 算例2

邊坡的幾何尺寸與算例1保持一致，將土體黏聚力改為c=0 kPa，此時問題成為無黏性土坡的穩定性問題，此時的安全系數應為

(4)

式中：φ為土體的內摩擦角；β為邊坡坡角。

極限平衡法得到的安全系數為F=0.33，破壞模式為表層破壞，與理論分析的結果一致。用強度折減法計算得到的安全系數為Fs=0.30，塑性剪應變增量云圖如圖9所示，同樣為表層破壞模式。

圖9 算例2的SRM塑性區云圖

SPH方法得到的位移場云圖如圖10所示，同樣為沿表層的滑移破壞。

圖10 算例2的SPH方法滑動位移

所以，極限平衡法、強度折減法以及SPH方法適用于無黏性土邊坡的穩定性分析。而應用DLO方法計算該問題時，折減強度后輸出的結果為“unknown”，折減荷載后的結果為“unstable”，均沒能給出一個確定的安全系數。將黏聚力調整為c=1 kPa時，結果仍然不能確定。這說明DLO方法不能很好地處理c接近于0的幾乎無黏性、純摩擦情況的邊坡穩定問題。在實際工程中遇到此類情況時應慎用DLO。

3.3 算例3

采用文獻[11]中的算例，算例幾何模型如圖11所示，材料參數見表1。高度為8 m的大壩(土層A)修建于成層土地基上，其中第一層土為強度較高的砂土層(土層B)，第二層為軟弱層(土層C)，最下層為基巖。

圖11 算例3的模型

表1 算例3土層參數

分別應用LEM，SRM，DLO三種方法分析該問題，得到的臨界滑動面、安全系數和塑性剪應變增量云圖如圖12所示。需要指出的是：由于該邊坡為失穩狀態，在SPH分析時，需要人為提高土層的強度參數以便獲取邊坡的自重應力，然后再利用原始參數進行滑坡過程的模擬。

圖12 算例3的LEM，DLO臨界滑動面及SRM塑性區云圖

由圖12可知，LEM，SRM和DLO得到的安全系數分別為0.875，0.86和1.106，LEM和SRM所得安全系數基本一致。DLO方法的安全系數與LEM比較，計算誤差為26.4%；與SRM比較，計算誤差為28.6%，誤差過大。

再分析該問題的臨界滑動面，DLO方法得到的滑動面在土層交界處均出現一個比較大的凸(凹)角，這與SPH位移場云圖吻合較好，而LEM的臨界滑動面左右兩部分的趨勢變化不太明顯，這種較為常規的滑動面在軟弱交替土層中出現似乎是不合理的。另外，滑動面在軟弱帶中的部分，LEM的滑動面沒有穿過軟弱帶到達基巖上部，而DLO方法的滑動面以及SRM，SPH方法的塑性區均到達基巖表面(圖12、圖13)，這說明LEM尋找復雜邊坡問題臨界滑動面的能力有待提高。這是因為LEM通常需要人為設置滑動面形狀和滑入、滑出范圍，對于這種復雜的情況，潛在的滑入、滑出范圍不明顯，人為設置過于主觀，存在一定誤差。所以，對于確定復雜條件下邊坡失穩的臨界滑動面位置，DLO與SPH方法比LEM更加具有優勢。

圖13 算例3的LEM，DLO臨界滑動面及SPH塑性區云圖

對于幾何模型復雜、存在軟弱夾層等復雜地質條件的邊坡穩定性問題，建議應用多種理論或軟件進行分析，綜合考量以得出正確的結論來指導工程設計。

4 結論

1) 對于簡單的均質邊坡，LEM，SRM和DLO方法均能得到令人滿意的結果，可以得出合理的安全系數以及臨界滑動面，SPH模擬得到的位移場云圖與強度折減法所得塑性剪應變增量云圖在識別臨界滑動面方面具有一致性。

2) DLO方法對于接近無黏性、純摩擦土體的情況會出現不易收斂，甚至計算失敗的現象，遇到此類情況應慎用DLO方法。

3) 對于幾何模型復雜、存在軟弱夾層等復雜地質條件的邊坡穩定性問題，在計算安全系數上，LEM與SRM所得安全系數基本一致，DLO方法所得安全系數偏高，誤差較大。在確定臨界滑動面上，DLO方法的結果與SPH方法的位移場云圖十分一致，這是因為兩者均不需對滑動面的位置做過多假設，對于確定復雜邊坡的臨界滑動面表現出了明顯的優勢。所以，對于復雜條件的邊坡穩定問題，建議應用多種理論進行分析，綜合考量后得出合理結論。