基于有限元極限平衡法的三維邊坡穩(wěn)定性

蘇振寧，邵龍?zhí)?/p>

大連理工大學工業(yè)裝備結構分析國家重點實驗室，大連 116024

邊坡穩(wěn)定分析是巖土工程中的重要課題.通常三維空間中的邊坡選取典型斷面并假設斷面處于平面應變條件，從而進行二維穩(wěn)定分析[1].二維分析導致結果偏于保守，且不能反映邊坡破壞的三維特性.許多學者開展了三維邊坡穩(wěn)定分析方法的研究，有滿足三個方向上的力和力矩平衡的嚴格三維極限平衡法[2-6]，也有用于三維邊坡穩(wěn)定分析的強度折減法[7-10].

葛修潤[11]指出極限平衡法和強度折減法都建立在強度折減的概念上，存在許多不合理之處，故在“真實”應力狀態(tài)的基礎上提出了矢量和法并進行了一系列研究[12-13].“真實”應力狀態(tài)是指采用數(shù)值方法根據(jù)未折減的土體強度參數(shù)計算得到的應力場.也有學者在真實應力狀態(tài)下提出了安全系數(shù)為剪應力比形式的三維穩(wěn)定分析方法[14-16]，該形式的安全系數(shù)，理論基礎不牢固，僅在滑動面是球形或橢球形時，安全系數(shù)才具有力矩比的物理意義[17].

在二維邊坡穩(wěn)定分析中，邵龍?zhí)逗屠罴t軍[18]將一點極限平衡條件擴展到沿滑動面的整體極限平衡條件，根據(jù)該條件定義了安全系數(shù)，明確了其物理意義，建立了有限元極限平衡法，并將其用于各類土工結構的穩(wěn)定分析[19].該方法理論體系嚴密，應力場真實，能適用于任意形狀的滑動面.

本文將有限元極限平衡法擴展到三維，提出了三維滑動面上一點在滑動方向上的極限平衡條件，證明了沿滑動面的整體極限平衡等價條件.給出了主滑方向的確定方法以及基于主滑方向的滑動方向計算方法.基于極限平衡條件定義了局部和整體安全系數(shù).最后通過算例驗證了本文方法的合理性、有效性以及對任意形狀滑動面的適應性.

1 極限平衡條件

1.1 一點滑動方向上的極限平衡條件

極限平衡狀態(tài)是研究對象即將失去而未失去平衡的狀態(tài)，極限平衡法和強度折減法都通過折減強度參數(shù)的方法，使邊坡達到極限平衡狀態(tài).極限平衡法中的極限平衡狀態(tài)為滑動面上每個點在滑動切平面上的剪應力都等于抗剪強度，強度折減法的極限平衡狀態(tài)判斷則根據(jù)失穩(wěn)判據(jù)主要分成3 種[20]：數(shù)值解非收斂、塑性區(qū)貫通和位移突變.采用有限元極限平衡法對邊坡進行穩(wěn)定分析，因為應力場真實（非極限平衡狀態(tài)），所以需要對極限平衡狀態(tài)進行預估.

在二維平面應變問題中，滑動面切平面決定剪應力方向和抗剪強度方向.剪應力方向沿滑動面切向指向滑動方向.無論是非極限平衡狀態(tài)還是極限平衡狀態(tài)，剪應力只有大小變化沒有方向變化，如圖1(a).

圖1 滑動面應力分析.(a)二維；(b)三維Fig.1 Stress on a slip surface: (a) two-dimensional;(b) threedimensional

在三維中，滑動面上一點的法向正應力和剪應力可根據(jù)該點處應力張量和滑動面切平面法向量計算：

式中，σ為該點處應力張量，n為滑動面切平面的法向單位矢量，T為應力矢量，σn為法向正應力矢量，τ為剪應力矢量.當滑動面確定后，由非極限平衡狀態(tài)向極限平衡狀態(tài)的變化過程中，法向正應力矢量方向不變，但剪應力矢量方向在滑動面的切平面上變化.根據(jù)Chen 的研究[21]，土體發(fā)生塑性應變時剪應變與剪應力可假定為同向.所以假設極限平衡狀態(tài)時剪應力方向與位移方向一致，即為滑動方向.需注意的是，各點的滑動方向是在滑動面切平面上的方向，每個點切平面的法方向不同，滑動方向也不同.若確定了滑動方向，非極限平衡狀態(tài)下的滑動面應力如圖1(b)所示.其中τf是極限平衡狀態(tài)的抗剪強度矢量.在極限平衡狀態(tài)下，考慮在滑動方向上有：

式中，d是一點滑動方向矢量，因為 σn垂直于滑動面，d在滑動面上，所以 σn·d=0，保持左右兩側一致，并考慮到τf與d方向相反，式（2）變?yōu)?/p>

式（3）為一點在滑動方向上處于極限平衡狀態(tài)的條件.

1.2 滑動面整體極限平衡條件

任意滑動面上土體整體達到極限平衡狀態(tài)與滑動面上土體各處在滑動方向上處于極限平衡狀態(tài)等價.滑動面整體達到極限平衡狀態(tài)定義為：

式（3）在滑動面上處處成立與式（4）的等價，證明如下：

如果滑動面上各處都在其自身的滑動方向上處于極限平衡狀態(tài)，即滿足式（3），則對滑動面上微元有

對式（5）兩側進行積分，所以式（4）成立.

反過來，如果式（4）成立，考慮到積分符號內為標量，并且矢量點乘有結合律，則有

因為對于穩(wěn)定或者處于極限平衡狀態(tài)的土體，每一點都必有

且 ds＞0，則要使式（6）成立，必須積分域上處處滿足：

所以滑動面上每一點都滿足式（3）.

需要注意，即使當滑動面上各點的剪應力矢量的大小都等于抗剪強度，整個滑動面也不是處于極限平衡狀態(tài).滑動方向意味著一個運動許可方向，只有當各點的剪應力矢量方向指向滑動方向且大小等于抗剪強度時，整體滑動面才處于極限平衡狀態(tài).

2 滑動方向與主滑方向

滑動面上不同點有不同的滑動方向d.滑動面上各點滑動方向的計算方法是本文方法的關鍵點之一.需要根據(jù)現(xiàn)有的有限元計算結果，得到合理的滑動方向.

滑動方向的計算需要用到主滑方向的概念.Kalatehjari 等[22]根據(jù)極限平衡法的基本假定之一，剛體滑動體假定，闡明了三維邊坡滑動存在唯一的主滑方向.主滑方向被定義為極限平衡狀態(tài)下，滑坡體滑動方向的水平投影方向.主滑方向不同于滑動方向，主滑方向是定義在XOY平面上，每一個滑動面只對應一個主滑方向.現(xiàn)對唯一主滑方向的合理性進行說明：

對滑動體進行垂直條分，根據(jù)假定每個條分柱體都是剛體.在極限平衡狀態(tài)下，假設每個柱體運動方向的水平投影有3 種情況：平行、聚攏和分離.剛體假設阻止了柱體聚攏，而柱體分離則意味著柱體間沒有力的相互作用，力和力矩平衡被打破，這與實際情況不符.所以互相平行的柱體運動方向水平投影是唯一的選項.

在三維極限平衡法中，很多方法對條柱底面的極限剪應力方向也有類似假定，比如要求剪應力平行于XOZ平面[23-24]或剪應力在XOY平面內的投影方向保持一致[5-6,22].極限平衡法是通過折減材料強度參數(shù)獲得滿足力和力矩平衡的極限狀態(tài)從而計算安全系數(shù)的方法，底面剪應力方向就是極限狀態(tài)時的滑動方向.在非極限狀態(tài)下，本文根據(jù)真實應力場計算滑動面上的剪應力沒有這種方向一致性，但通過主滑方向計算得到的滑動方向，正是對極限狀態(tài)時滑動方向的預估，也是對極限狀態(tài)時剪應力方向的預估.

滑動面上各點的滑動方向可以根據(jù)唯一的主滑方向計算.根據(jù)滑動方向處于滑動面上；滑動方向在XOY平面的投影是主滑方向；滑動方向是單位矢量，可以聯(lián)立方程計算滑動方向單位矢量：

式中，du是主滑方向矢量，k是Z方向的單位矢量.3 個方程可以求解滑動方向矢量d的3 個元素.

通過非極限平衡狀態(tài)的滑動面應力場對極限平衡狀態(tài)推測得到主滑方向，可以將滑動面剪應力矢量的和矢量的水平投影方向作為主滑方向，用公式表示為：

主滑方向的選取具有一定的主觀性，比如設計人員認為邊坡在未來會受到某一方向的荷載，從而產生沿某一特定方向的破壞，那么要考慮這一特定方向的安全系數(shù)，就可以將這一特定方向設為主滑方向進行計算.

3 局部和整體安全系數(shù)定義

安全系數(shù)是對研究對象距離極限平衡狀態(tài)的描述，安全系數(shù)等于1 則認為研究對象處于極限平衡狀態(tài).根據(jù)對一點在滑動方向上處于極限平衡狀態(tài)的定義，一點處的局部安全系數(shù)被定義為：

式（11）表示一點安全系數(shù)是抗剪強度與剪應力在滑動方向上的投影的比，其物理意義是一點在滑動方向與極限平衡狀態(tài)的距離，也是一點在滑動方向上強度的儲備.

F(s)是 在滑動面s上的局部安全系數(shù)函數(shù)，在各點符合式（11）的定義，可使滑動面上土體微元均達到極限平衡狀態(tài)，那么土體在滑動面s上整體達到極限平衡就變?yōu)?/p>

應用積分中值定理，式（12）右側變?yōu)椋?/p>

FG是整體安全系數(shù)，表征滑動面上各點局部安全系數(shù)的中值.在推導整體安全系數(shù)的定義中，對滑動面形狀和對稱性未做要求，所以方法可以用于任意形狀滑動面的安全系數(shù)計算.

二維情況下，因為 τ和d方向相同，τf和d方向相反，式（9）退化為

這是在基于應力的二維邊坡穩(wěn)定分析方法中被廣泛使用的全局安全系數(shù)定義式，也是有限元極限平衡法的安全系數(shù)定義[18].

三維情況下，剪應力比形式的安全系數(shù)定義為：

對比安全系數(shù)定義式（16）和式（14），如果假定 τ和d方 向相同，τf和d方向相反，式（14）可以退化為式（16）.這種假定隱含剪應力矢量方向在非極限平衡狀態(tài)和極限平衡狀態(tài)相同的設定，對傾向滑動方向的滑動面，剪應力矢量方向變化不大，但對傾向方向偏離滑動方向的滑動面，剪應力矢量方向在非極限平衡狀態(tài)和極限平衡狀態(tài)差距很大，這會導致誤差，這也是采用式（16）定義安全系數(shù)的限制.

4 計算方法

本文方法需要計算邊坡區(qū)域的應力場，并基于應力計算邊坡的安全系數(shù).有限元法可以根據(jù)本構關系計算應力場.在本文中，有限元本構使用線性彈性理想塑性模型，計算軟件為商業(yè)有限元軟件Abaqus 6.14.理想的塑性屈服準則與Mohr-Coulomb 破壞準則一致.抗剪強度采用Mohr-Coulomb 計算：

式中，σn為 法向正應力，φ為摩擦角，c為黏聚力.

在計算過程中需要對滑動面上的變量進行積分.數(shù)值積分需要對滑動面進行離散，對滑動面進行網格劃分.在這里選擇將滑動面劃分為三角形網格，原因是三角形網格比四邊形網格更容易計算滑動面的法線方向并更準確地描述滑動面，在邊界處也更容易分割.通過3 個步驟構造滑動面，如圖2 所示.

圖2 滑動面構造步驟.(a)三角形網格；(b)網格節(jié)點調整；(c)滑動面與坡面相交Fig.2 Construction steps of the slip surface: (a) triangle mesh;(b) mesh node adjustment;(c) slip surface intersects the slope surface

第一步：在XOY平面中創(chuàng)建三角形網格；

第二步：通過插值將三角形網格節(jié)點的z坐標從XOY平面調整到構造表面的位置；

第三步：獲得坡面與三角形網格的交點，三角形網格被分割，模型內部的部分被保留，作為滑動表面.

另一種方法是，如果滑動面在有限元模型中已經被劃分（比如滑動面是巖體結構面或確定的軟弱層），則可以直接將滑動面位置的四面體單元表面的三角形網格提取出來，作為滑動面網格使用；如果有限元采用的是六面體單元，則可以將四邊形表面從中剖分，得到兩個三角形網格單元.

滑動面三角形網格中心點應力張量由有限元網格節(jié)點處應力張量插值計算，插值方法選用三維線性插值.本文方法安全系數(shù)計算流程如圖3所示.

圖3 安全系數(shù)計算流程Fig.3 Flowchart of the safety factor calculation

5 算例

算例中將采用式（16）剪應力比形式定義的安全系數(shù)記作F1，將采用式（14）定義的安全系數(shù)記作F2.

5.1 算例1

算例1 來自Fredlund 和Krahn[25]對二維邊坡的研究，有學者將這一算例拓展到三維情況并進行了研究.該算例包含2 種工況，一種是均質邊坡，滑動面為橢球形；一種是帶軟弱薄層的邊坡，滑動面為橢球型和位于軟弱薄層的平面組合.橢圓的長軸為133.8 m，有限元模型y軸方向的長度為110 m，其他相關參數(shù)如圖4 所示，圖中γ為重度，E為彈性模量，v為泊松比，H為坡高，β為坡角，R為短軸半徑.橫向邊界條件為約束Y軸方向位移（平面應變）.圖5 展示滑動面網格和有限元網格尺寸對安全系數(shù)的影響.隨著滑動面網格尺寸從4 減小到0.75 m，安全系數(shù)逐漸減小并穩(wěn)定，有限元網格尺寸從4 減小到1 m，安全系數(shù)略有減小.圖6 和圖7 分別展示工況1 的滑動方向和剪應力矢量方向，通過對比可見兩者方向基本一致，只在側面靠近坡面處略有差別.因此表1 中F1與F2在工況1 下基本一致.

表1 算例1 的安全系數(shù)比較Table 1 Comparison of the safety factors computed for example 1

圖4 算例1.(a)材料參數(shù)和斷面；(b)工況1；(c)工況2Fig.4 Example 1: (a) material parameters and cross-section;(b) case 1;(c) case 2

圖5 單元數(shù)對安全系數(shù)的影響.(a)滑動面；(b)有限元Fig.5 Influence of the number of units on the safety factor: (a) slip surface;(b) FEM

圖6 工況1 的滑動方向與主滑方向.(a)等軸測圖；(b)俯視圖Fig.6 Sliding direction and the main sliding direction of case 1: (a)isometric view;(b) top view

圖7 工況1 的剪應力矢量方向與主滑方向.(a)等軸測圖；(b)俯視圖Fig.7 Shear stress vector direction and the main sliding direction of case 1: (a) isometric view;(b) top view

5.2 算例2

楔形體破壞是三維巖石邊坡破壞的常見情況.如圖8 所示，這里使用兩個典型的楔形破壞算例[26]來驗證本文方法的合理性.算例包括具有對稱結構表面的楔形體和具有非對稱結構表面的楔形體.算例的材料參數(shù)和幾何信息列于表2 和表3中.通過圖9 和圖10 的對比可以明顯看出，在非極限平衡狀態(tài)下，剪應力矢量方向與滑動方向不一致，所以導致表4 中F1結果小于其他方法安全系數(shù).Hoek-Bray 解是一個靜力平衡解，一般被認為是該算例的理論解，本文方法和嚴格極限平衡法都與理論解一致，但F1小于理論解，有較大誤差.

表4 算例2 的安全系數(shù)比較Table 4 Comparison of the safety factors for example 2

圖8 算例2 楔形體破壞.(a)對稱楔形體；(b)非對稱楔形體Fig.8 Example 2: wedge failure: (a) symmetric wedge;(b) asymmetric wedge

圖9 非對稱楔形體的滑動方向與主滑方向.(a)等軸測圖；(b)俯視圖Fig.9 Sliding direction and the main sliding direction of an asymmetrical wedge: (a) isometric view;(b) top view

圖10 非對稱楔形體的剪應力矢量方向與主滑方向.(a)等軸測圖；(b)俯視圖Fig.10 Shear stress vector direction and the main sliding direction of an asymmetrical wedge: (a) isometric view;(b) top view

表2 算例2 的材料參數(shù)Table 2 Mechanical parameters used in example 2

表3 算例2 的幾何信息Table 3 Geometric information used in example 2 (°)

5.3 討論

從算例1 和算例2 的安全系數(shù)結果可見，對非球形滑動面采用剪應力比形式的安全系數(shù)F1較嚴格極限平衡法偏差較大.因為F1計算方法中默認滑動方向與非極限平衡狀態(tài)下的剪應力矢量方向一致，這與嚴格極限平衡法的假定不同.

在本文中，雖然安全系數(shù)F2與嚴格極限平衡法的計算過程不同，但結果仍具有一致性.原因是：（1）本文方法計算中通過剛體假定和主滑方向的假設，預估了極限平衡狀態(tài)下的滑動方向.因為采用相同的假設，所以滑動方向與極限平衡法計算得到的剪應力矢量方向一致.（2）本文方法中一點安全系數(shù)代表在滑動方向上強度的儲備，整體安全系數(shù)是滑動面上各點的安全系數(shù)的中值.而極限平衡法通過對強度進行折減達到極限平衡狀態(tài)得到的安全系數(shù)也代表著強度儲備.

6 結語

基于滑動面唯一主滑方向計算滑動方向，通過對三維滑動面極限平衡狀態(tài)的分析，將滑動面上局部安全系數(shù)定義為抗剪強度與剪應力在滑動方向上投影的比，并通過證明滑動面整體極限平衡條件，定義滑動面整體安全系數(shù).通過與經典算例對比，驗證了本文方法的可行性，并獲得如下結論：

（1）本文方法安全系數(shù)定義合理，物理意義明確，計算過程簡單，可以適用于任意形狀滑動面，能應用于三維邊坡穩(wěn)定性分析中.

（2）采用剪應力比形式的三維安全系數(shù)定義沒有考慮滑動方向，會錯誤的估計極限平衡狀態(tài)的剪應力矢量方向，導致安全系數(shù)誤差.本文方法考慮了滑動方向上的極限平衡狀態(tài)，且滑動方向假設合理，計算結果與嚴格極限平衡法一致.