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

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

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

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

在二維邊坡穩(wěn)定分析中，邵龍?zhí)逗屠罴t軍[18]將一點極限平衡條件擴展到沿滑動面的整體極限平衡條件，根據(jù)該條件定義了安全系數(shù)，明確了其物理意義，建立了有限元極限平衡法，并將其用于各類土工結(jié)構(gòu)的穩(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ā)生塑性應變時剪應變與剪應力可假定為同向.所以假設(shè)極限平衡狀態(tài)時剪應力方向與位移方向一致，即為滑動方向.需注意的是，各點的滑動方向是在滑動面切平面上的方向，每個點切平面的法方向不同，滑動方向也不同.若確定了滑動方向，非極限平衡狀態(tài)下的滑動面應力如圖1(b)所示.其中τf是極限平衡狀態(tài)的抗剪強度矢量.在極限平衡狀態(tài)下，考慮在滑動方向上有：

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

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

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

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

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

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

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

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

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

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

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

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

2 滑動方向與主滑方向

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

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

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

在三維極限平衡法中，很多方法對條柱底面的極限剪應力方向也有類似假定，比如要求剪應力平行于XOZ平面[23-24]或剪應力在XOY平面內(nèi)的投影方向保持一致[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è)計人員認為邊坡在未來會受到某一方向的荷載，從而產(chǎn)生沿某一特定方向的破壞，那么要考慮這一特定方向的安全系數(shù)，就可以將這一特定方向設(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）右側(cè)變?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)相同的設(shè)定，對傾向滑動方向的滑動面，剪應力矢量方向變化不大，但對傾向方向偏離滑動方向的滑動面，剪應力矢量方向在非極限平衡狀態(tài)和極限平衡狀態(tài)差距很大，這會導致誤差，這也是采用式（16）定義安全系數(shù)的限制.

4 計算方法

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

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

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

圖2 滑動面構(gòu)造步驟.(a)三角形網(wǎng)格；(b)網(wǎng)格節(jié)點調(diào)整；(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)建三角形網(wǎng)格；

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

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

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

滑動面三角形網(wǎng)格中心點應力張量由有限元網(wǎng)格節(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 種工況，一種是均質(zhì)邊坡，滑動面為橢球形；一種是帶軟弱薄層的邊坡，滑動面為橢球型和位于軟弱薄層的平面組合.橢圓的長軸為133.8 m，有限元模型y軸方向的長度為110 m，其他相關(guān)參數(shù)如圖4 所示，圖中γ為重度，E為彈性模量，v為泊松比，H為坡高，β為坡角，R為短軸半徑.橫向邊界條件為約束Y軸方向位移（平面應變）.圖5 展示滑動面網(wǎng)格和有限元網(wǎng)格尺寸對安全系數(shù)的影響.隨著滑動面網(wǎng)格尺寸從4 減小到0.75 m，安全系數(shù)逐漸減小并穩(wěn)定，有限元網(wǎng)格尺寸從4 減小到1 m，安全系數(shù)略有減小.圖6 和圖7 分別展示工況1 的滑動方向和剪應力矢量方向，通過對比可見兩者方向基本一致，只在側(cè)面靠近坡面處略有差別.因此表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]來驗證本文方法的合理性.算例包括具有對稱結(jié)構(gòu)表面的楔形體和具有非對稱結(jié)構(gòu)表面的楔形體.算例的材料參數(shù)和幾何信息列于表2 和表3中.通過圖9 和圖10 的對比可以明顯看出，在非極限平衡狀態(tài)下，剪應力矢量方向與滑動方向不一致，所以導致表4 中F1結(jié)果小于其他方法安全系數(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ù)結(jié)果可見，對非球形滑動面采用剪應力比形式的安全系數(shù)F1較嚴格極限平衡法偏差較大.因為F1計算方法中默認滑動方向與非極限平衡狀態(tài)下的剪應力矢量方向一致，這與嚴格極限平衡法的假定不同.

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

6 結(jié)語

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

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

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