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

大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024

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

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

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

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

1 極限平衡條件

1.1 一點(diǎn)滑動(dòng)方向上的極限平衡條件

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

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

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

在三維中，滑動(dòng)面上一點(diǎn)的法向正應(yīng)力和剪應(yīng)力可根據(jù)該點(diǎn)處應(yīng)力張量和滑動(dòng)面切平面法向量計(jì)算：

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

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

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

1.2 滑動(dòng)面整體極限平衡條件

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

式（3）在滑動(dòng)面上處處成立與式（4）的等價(jià)，證明如下：

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

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

反過(guò)來(lái)，如果式（4）成立，考慮到積分符號(hào)內(nèi)為標(biāo)量，并且矢量點(diǎn)乘有結(jié)合律，則有

因?yàn)閷?duì)于穩(wěn)定或者處于極限平衡狀態(tài)的土體，每一點(diǎn)都必有

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

所以滑動(dòng)面上每一點(diǎn)都滿足式（3）.

需要注意，即使當(dāng)滑動(dòng)面上各點(diǎn)的剪應(yīng)力矢量的大小都等于抗剪強(qiáng)度，整個(gè)滑動(dòng)面也不是處于極限平衡狀態(tài).滑動(dòng)方向意味著一個(gè)運(yùn)動(dòng)許可方向，只有當(dāng)各點(diǎn)的剪應(yīng)力矢量方向指向滑動(dòng)方向且大小等于抗剪強(qiáng)度時(shí)，整體滑動(dòng)面才處于極限平衡狀態(tài).

2 滑動(dòng)方向與主滑方向

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 計(jì)算方法

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

式中，σn為 法向正應(yīng)力，φ為摩擦角，c為黏聚力.

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

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

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

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

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

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

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

5 算例

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

5.1 算例1

算例1 來(lái)自Fredlund 和Krahn[25]對(duì)二維邊坡的研究，有學(xué)者將這一算例拓展到三維情況并進(jìn)行了研究.該算例包含2 種工況，一種是均質(zhì)邊坡，滑動(dòng)面為橢球形；一種是帶軟弱薄層的邊坡，滑動(dòng)面為橢球型和位于軟弱薄層的平面組合.橢圓的長(zhǎng)軸為133.8 m，有限元模型y軸方向的長(zhǎng)度為110 m，其他相關(guān)參數(shù)如圖4 所示，圖中γ為重度，E為彈性模量，v為泊松比，H為坡高，β為坡角，R為短軸半徑.橫向邊界條件為約束Y軸方向位移（平面應(yīng)變）.圖5 展示滑動(dòng)面網(wǎng)格和有限元網(wǎng)格尺寸對(duì)安全系數(shù)的影響.隨著滑動(dòng)面網(wǎng)格尺寸從4 減小到0.75 m，安全系數(shù)逐漸減小并穩(wěn)定，有限元網(wǎng)格尺寸從4 減小到1 m，安全系數(shù)略有減小.圖6 和圖7 分別展示工況1 的滑動(dòng)方向和剪應(yīng)力矢量方向，通過(guò)對(duì)比可見(jiàn)兩者方向基本一致，只在側(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ù)對(duì)安全系數(shù)的影響.(a)滑動(dòng)面；(b)有限元Fig.5 Influence of the number of units on the safety factor: (a) slip surface;(b) FEM

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

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

5.2 算例2

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

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

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

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

圖10 非對(duì)稱楔形體的剪應(yīng)力矢量方向與主滑方向.(a)等軸測(cè)圖；(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é)果可見(jiàn)，對(duì)非球形滑動(dòng)面采用剪應(yīng)力比形式的安全系數(shù)F1較嚴(yán)格極限平衡法偏差較大.因?yàn)镕1計(jì)算方法中默認(rèn)滑動(dòng)方向與非極限平衡狀態(tài)下的剪應(yīng)力矢量方向一致，這與嚴(yán)格極限平衡法的假定不同.

在本文中，雖然安全系數(shù)F2與嚴(yán)格極限平衡法的計(jì)算過(guò)程不同，但結(jié)果仍具有一致性.原因是：（1）本文方法計(jì)算中通過(guò)剛體假定和主滑方向的假設(shè)，預(yù)估了極限平衡狀態(tài)下的滑動(dòng)方向.因?yàn)椴捎孟嗤募僭O(shè)，所以滑動(dòng)方向與極限平衡法計(jì)算得到的剪應(yīng)力矢量方向一致.（2）本文方法中一點(diǎn)安全系數(shù)代表在滑動(dòng)方向上強(qiáng)度的儲(chǔ)備，整體安全系數(shù)是滑動(dòng)面上各點(diǎn)的安全系數(shù)的中值.而極限平衡法通過(guò)對(duì)強(qiáng)度進(jìn)行折減達(dá)到極限平衡狀態(tài)得到的安全系數(shù)也代表著強(qiáng)度儲(chǔ)備.

6 結(jié)語(yǔ)

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

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

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