呂愛鐘 賈曉陽

(華北電力大學(xué)水利與水電工程學(xué)院， 北京 102206， 中國(guó))

0 引 言

邊坡穩(wěn)定問題是巖土力學(xué)領(lǐng)域的一大經(jīng)典問題，近百年來，學(xué)者們提出了一系列的穩(wěn)定性分析方法。

最早提出的瑞典圓弧法(Fellenius， 1936)假設(shè)邊坡的滑動(dòng)面為圓弧面，在發(fā)生滑動(dòng)破壞時(shí)滑動(dòng)面上達(dá)到Mohr-Coulomb 強(qiáng)度準(zhǔn)則的極限狀態(tài)，在進(jìn)行計(jì)算時(shí)，需要將滑動(dòng)體垂直分成若干條塊，并將條塊視為剛體，通過假定條塊間的內(nèi)力，列出靜力平衡方程對(duì)邊坡穩(wěn)定安全系數(shù)進(jìn)行求解。這也就是經(jīng)典極限平衡法的理論基礎(chǔ)，基于這一原理和不同的條塊間內(nèi)力的假定，Bishop(1955)、Morgenstern et al. (1965)、Janbu(1973)、Sarma(1973)等研究者們提出了一系列的極限平衡穩(wěn)定性分析方法。幾十年來，極限平衡法因其計(jì)算簡(jiǎn)單并且適用范圍廣，迅速衍生出了一系列改進(jìn)的分析方法。Maslov(1990)、Zhu et al. (2003)、Zheng(2009)基于滑動(dòng)面正應(yīng)力分布的假設(shè)，先后提出了滿足平衡條件而且不需要將滑動(dòng)體進(jìn)行豎向分割的分析方法。陳祖煜等(2001)、朱大勇等(2007)和鄭宏(2007)先后提出了邊坡的三維極限平衡計(jì)算方法。

以往的極限平衡條分法在被廣泛應(yīng)用的同時(shí)，也暴露出其理論上的缺陷，例如：在分析中認(rèn)為邊坡體是不可變形的剛體，這顯然是不符合實(shí)際的； 對(duì)條塊間內(nèi)力或者對(duì)滑動(dòng)面上正應(yīng)力分布的假定，以及對(duì)滑動(dòng)面形狀假定的是否合理？計(jì)算中只求解了平衡方程，所得結(jié)果是否滿足變形協(xié)調(diào)方程？這是經(jīng)典極限平衡法無法問答的問題(Lu et al.，2012)。

隨著計(jì)算機(jī)技術(shù)的發(fā)展，有限元法等數(shù)值方法在力學(xué)計(jì)算領(lǐng)域得到了廣泛應(yīng)用，這種將物體微元化的求解方法不僅滿足平衡方程而且還兼顧了變形協(xié)調(diào)方程，研究者們也將這種適用性很廣的方法引入到了邊坡穩(wěn)定性分析問題中，但是傳統(tǒng)的有限元法不能直接求得邊坡的穩(wěn)定安全系數(shù)。Giam et al. (1988)提出了一種根據(jù)有限元法計(jì)算應(yīng)力場(chǎng)確定臨界滑動(dòng)面及最小穩(wěn)定安全系數(shù)的模式搜索方法，Kim(1997)也提出了基于相似理論的方法。Ugai(1989)，鄭穎人等(2002)，韓偉歌等(2019)，高馮等(2020)，劉康琦等(2020)等學(xué)者提出了強(qiáng)度折減法，通過不斷增大折減系數(shù)，減小邊坡材料的強(qiáng)度參數(shù)，直至邊坡發(fā)生失穩(wěn)破壞，將此時(shí)的折減系數(shù)定義為邊坡的穩(wěn)定安全系數(shù)。此類方法計(jì)算工作量較大，定義的穩(wěn)定安全系數(shù)物理含義也不明確。

近年來，巖土力學(xué)研究者嘗試了更多新方法來解決邊坡穩(wěn)定性分析問題。Huang et al. (1993)、楊明成等(2002)等學(xué)者提出了在有限元強(qiáng)度折減法或者條分法力學(xué)模型的基礎(chǔ)上，逐段尋找最危險(xiǎn)滑動(dòng)面的局部最小安全系數(shù)法，這類分析方法在計(jì)算和收斂速度上有優(yōu)越性，但無法回避強(qiáng)度折減法或條分法的局限性。王斌等(2017)，張巍等(2017)，劉鑫等(2019)等學(xué)者將物質(zhì)點(diǎn)法引入邊坡穩(wěn)定分析中，這種方法能夠很好地描述結(jié)構(gòu)初始失穩(wěn)破壞后的形態(tài)發(fā)展，但其計(jì)算十分復(fù)雜，無法得出意義明確的穩(wěn)定安全系數(shù)。蘇永華等(2019)假設(shè)滑動(dòng)面為對(duì)數(shù)螺旋面，基于塑性極限分析的上限分析法構(gòu)建邊坡可靠度計(jì)算模型，引入響應(yīng)面模型，求解得到邊坡可靠度的一個(gè)上限解，其結(jié)果有一定的工程借鑒意義，但也只是陳祖煜提出的塑性邊坡分析上下限原理的其一，不是準(zhǔn)確結(jié)果。

尋找一種物理意義明確，假設(shè)條件少，計(jì)算操作簡(jiǎn)單的邊坡穩(wěn)定分析方法，這正是本文要研究的問題。

圖 1 曲面上一點(diǎn)的法向正應(yīng)力σn和切向剪應(yīng)力τnFig. 1 The normal stress σn and the tangential shear stress τn at a point on the surface

本文提出的方法與經(jīng)典的極限平衡法不同，計(jì)算分析過程中不需要將滑體劃分成垂直條塊，不需要假設(shè)滑面上的法向應(yīng)力分布形式，而可以直接利用已有的坡體彈性應(yīng)力解。算例表明：本文方法所得到的與簡(jiǎn)化Bishop法和有限元強(qiáng)度折減法獲得的結(jié)果非常接近，且基本介于有限元法和簡(jiǎn)化Bishop法之間。

1 確定滑動(dòng)面的原理

1.1 滑動(dòng)面上的應(yīng)力分布

當(dāng)邊坡只受到重力作用時(shí)，坡體應(yīng)力場(chǎng)可直接通過彈性楔形體的彈性解給出(谷拴成等， 1994)：

(1)

式中：σx，σy分別為x、y向的正應(yīng)力，本文規(guī)定以壓為正；τxy為剪應(yīng)力；γ為巖體的容重；θ為邊坡的坡角。

當(dāng)為平面應(yīng)變問題時(shí)，AB曲面也可以理解為xOy平面內(nèi)的曲線。當(dāng)坡高為H時(shí)，可以求出通過坡角A的曲線x=f(y)上任一點(diǎn)的法向正應(yīng)力σn和切向剪應(yīng)力τn(圖 1)：

(2)

式中：

(3)

1.2 極限平衡原理和安全系數(shù)

τn是引起坡體滑動(dòng)的根源，而根據(jù)坡體材料滿足的Mohr-Coulomb準(zhǔn)則，σn所形成的內(nèi)摩擦力和材料固有的黏聚力可以阻止坡體的滑動(dòng)，則穩(wěn)定安全系數(shù)FS的表達(dá)式為：

(4)

2 滑動(dòng)面的參數(shù)表示

圖 2 滑動(dòng)面的分段表示Fig. 2 The piecewise description of sliding surface

如圖 2所示，對(duì)曲面AB沿豎直方向進(jìn)行等距分段，設(shè)每段曲面在豎直方向的高度都為Δh，且假定每段滑動(dòng)面的形狀都為一個(gè)二次多項(xiàng)式，當(dāng)分段間隔Δh足夠小時(shí)，這樣描述的曲面可以反映足夠多的曲面形狀。需要指出的是，這里進(jìn)行分段僅是為了描述滑動(dòng)面形狀，與經(jīng)典極限平衡法對(duì)滑動(dòng)體分條有本質(zhì)區(qū)別。若曲面AB共被分成n段，則滑動(dòng)面x=f(y)可以看作是如下的分段函數(shù)：

(5)

式中：ai，bi，ci(i=1， 2，..，n)為待定系數(shù)，式(5)中共有3n個(gè)未知系數(shù)，這3n個(gè)未知系數(shù)并不是相互獨(dú)立的，根據(jù)對(duì)曲面AB的連續(xù)性和光滑性要求，可以證明只有n+1個(gè)是獨(dú)立的。

曲面AB的連續(xù)性和光滑性是指分段函數(shù)在相鄰兩分段交界點(diǎn)處函數(shù)值及一階導(dǎo)數(shù)值是相等的，由此通過式(5)可列出在n-1個(gè)交點(diǎn)處，所必須滿足的2(n-1)個(gè)方程為：

(6)

(7)

式(6)和式(7)共2n-1個(gè)方程，其中的yi(i=0， 1，..，n-1)是已知量。在這些方程中含有3n個(gè)未知系數(shù)ai，bi，ci(i=1， 2，..，n)，所以共有n+1個(gè)獨(dú)立的待求未知量，本文將a1，b1，a2，a3…，an-1，an作為待求的未知量時(shí)，則由式(6)和式(7)可以求出其他未知量為：

(8)

式中：1≤i≤n-1，這樣將滑動(dòng)面分成n段共需要n+1個(gè)獨(dú)立的未知量就可以描述整個(gè)滑動(dòng)面了，變量個(gè)數(shù)的減少將更有利于最危險(xiǎn)滑動(dòng)面的優(yōu)化計(jì)算。

3 最危險(xiǎn)滑動(dòng)面的搜索

3.1 計(jì)算前的理論準(zhǔn)備

將式(3)代入式(2)可得：

(9)

式中的dx/dy=f′(y)可通過式(5)獲得。式(4)可寫為：

(10)

將式(1)代入式(9)，再代入式(10)便可得到計(jì)算任一曲面AB對(duì)應(yīng)的FS的公式：

(11)

使用數(shù)值積分方法計(jì)算式(11)中相關(guān)積分值，將曲面AB上的每個(gè)二次曲面段沿豎直方向等距劃分為m個(gè)微段，沿整個(gè)曲面共有m×n個(gè)微段，m×n+1個(gè)端點(diǎn)。當(dāng)每個(gè)微段的長(zhǎng)度足夠小時(shí)，在每個(gè)微段上可以認(rèn)為應(yīng)力呈直線分布，式(11)中相關(guān)積分值可以使用梯形積分公式計(jì)算。

為了保證搜索過程中滑動(dòng)面的合理性，還需要增加一些約束條件，本文僅對(duì)數(shù)值積分的分段點(diǎn)施加約束條件，第i個(gè)分段點(diǎn)的y坐標(biāo)值為hi，則有

(12)

-f′(hi)>0 1≤i≤m

(13)

式(12)中，-hi/tanθ為坡面線AO上y坐標(biāo)值為hi的點(diǎn)的x坐標(biāo)值，f(hi)代表滑動(dòng)面AB上y坐標(biāo)值為hi的點(diǎn)的x坐標(biāo)值，此約束條件的幾何意義為限制坡體內(nèi)任意深度處滑動(dòng)面AB均在坡面線AO的右側(cè)，保證滑動(dòng)面上的點(diǎn)在邊坡體的內(nèi)部。式(13)中的約束條件的幾何意義為保證滑動(dòng)面曲線方程有單調(diào)性，防止曲線產(chǎn)生向坡肩方向的彎曲。

3.2 最危險(xiǎn)滑動(dòng)面的優(yōu)化計(jì)算

經(jīng)過上述準(zhǔn)備之后，對(duì)最危險(xiǎn)滑動(dòng)面的搜索轉(zhuǎn)化成尋找一組使穩(wěn)定安全系數(shù)FS達(dá)到最小的一組系數(shù)a1，b1，a2，a3…，an-1，an，可以通過混合罰函數(shù)優(yōu)化方法來實(shí)現(xiàn)這個(gè)過程(Fletcher， 1981； 萬耀青， 1983)。

罰函數(shù)法是處理約束優(yōu)化問題的一個(gè)有效方法，其基本思想是：將原問題的目標(biāo)函數(shù)和約束條件按一定的方式構(gòu)造出帶參數(shù)的增廣目標(biāo)函數(shù)，把約束優(yōu)化問題轉(zhuǎn)化為一系列無約束優(yōu)化問題來求解。按照約束形式和構(gòu)造的函數(shù)及罰因子的遞推方法等不同，罰函數(shù)法又分為內(nèi)點(diǎn)罰函數(shù)法，外點(diǎn)罰函數(shù)法和混合罰函數(shù)法3種(Himmelblau， 1972)。內(nèi)點(diǎn)法容易處理不等式約束優(yōu)化設(shè)計(jì)問題，外點(diǎn)法容易處理等式約束優(yōu)化設(shè)計(jì)問題，把外點(diǎn)法和內(nèi)點(diǎn)法結(jié)合起來，則可同時(shí)處理既具有等式約束條件，又具有不等式約束條件的優(yōu)化問題，這種方法稱為混合懲罰函數(shù)法。

在本文的優(yōu)化計(jì)算中，設(shè)計(jì)變量為滑動(dòng)面的形狀參數(shù)即X=[a1，b1，a2，a3…an-1，an]，使用的優(yōu)化模型為：

(14)

式中：bi(X)>0表示不等式約束條件，可由式(11)，式(12)獲得，即：

(15)

在本文的優(yōu)化模型中，罰函數(shù)的形式為：

(16)

式中：r為罰因子，下標(biāo)集合I1，I2定義為：

(17)

式中：X(0)為給出的原問題的計(jì)算初始點(diǎn)。優(yōu)化過程可以通過編寫的計(jì)算機(jī)程序來完成。

4 算例分析

以坡體高度為8m的土坡為例(Huang， 1993； 楊明成等， 2002)，坡體材料容重為18.5kN·m-3，并取兩種不同坡率的邊坡進(jìn)行分析，一種是1︰1的邊坡，另一種是1︰2的邊坡。邊坡材料的黏聚力c從2.5kPa變化到30kPa，內(nèi)摩擦角φ從10°變化到20°，以考察強(qiáng)度參數(shù)變化對(duì)邊坡穩(wěn)定性的影響。對(duì)每一種情況，分別用本文方法，鄭穎人等(2002)中的有限元方法和邊坡設(shè)計(jì)中常用而且我國(guó)《公路路基設(shè)計(jì)規(guī)范》推薦的簡(jiǎn)化Bishop法確定臨界滑動(dòng)面并計(jì)算最小穩(wěn)定安全系數(shù)，計(jì)算結(jié)果見表 1。

表 1 不同分析方法所獲得的穩(wěn)定安全系數(shù)Table 1 Factor of safety obtained by different methods

表中，本文方法取m=100，n=10計(jì)算。簡(jiǎn)化Bishop法對(duì)應(yīng)的整體最小穩(wěn)定安全系數(shù)由文獻(xiàn)《土坡穩(wěn)定分析》(黃仰賢，1988)的穩(wěn)定圖表給出，這種方法假設(shè)滑動(dòng)面為圓弧面，認(rèn)為條塊間豎直方向的剪力為零，條塊間只存在水平方向的作用力，通過不同圓弧面的圓心和半徑組合試算尋找最危險(xiǎn)滑動(dòng)面。有限元法對(duì)應(yīng)的整體最小穩(wěn)定安全系數(shù)使用局部最小安全系數(shù)法(Huang et al.,1993)得到，該法使用基于準(zhǔn)則的有限元方法求解邊坡體的有效應(yīng)力狀態(tài)， 然后使用莫爾-庫侖準(zhǔn)則估算坡體內(nèi)一點(diǎn)處的剪切強(qiáng)度s，根據(jù)莫爾圓計(jì)算剪應(yīng)力τ，并定義局部安全系數(shù)為FSL=s/τ，通過強(qiáng)度折減，認(rèn)為局部穩(wěn)定安全系數(shù)小于或等于1的點(diǎn)為破壞點(diǎn)，將這些破壞點(diǎn)連接起來，就得到臨界滑動(dòng)面，而強(qiáng)度折減系數(shù)就是所要求的最小穩(wěn)定安全系數(shù)。

從表 1可見，邊坡的穩(wěn)定安全系數(shù)隨著坡體材料內(nèi)摩擦角φ和黏聚力c的增大而增大。3種方法獲得的穩(wěn)定系數(shù)相差很小。當(dāng)坡比為1︰1時(shí)，在12個(gè)算例中，基本上有這樣的規(guī)律：有限元法獲得的穩(wěn)定安全系數(shù)最大，簡(jiǎn)化Bishop法獲得的穩(wěn)定安全系數(shù)最小，而本文方法介于兩者之間； 而當(dāng)坡比為1︰2時(shí)，獲得的規(guī)律與坡比為1︰1時(shí)獲得的規(guī)律完全不同，在12個(gè)算例中，除了相等的穩(wěn)定安全系數(shù)之外，簡(jiǎn)化Bishop法獲得的穩(wěn)定安全系數(shù)有5個(gè)為最大，有限元法和本文方法獲得的穩(wěn)定安全系數(shù)各有2個(gè)為最大。綜上所述，本文方法獲得的穩(wěn)定安全系數(shù)基本介于有限元法和簡(jiǎn)化Bishop法之間。

圖 3和圖 4給出了1︰1的邊坡，兩種內(nèi)摩擦角和黏聚力情形下，3種方法確定的最可能滑動(dòng)面形狀。

圖 3 1︰1邊坡的最可能滑動(dòng)面(c=10kPa，φ=15°)Fig. 3 The most dangerous sliding surface in 1︰1 slope(c=10kPa，φ=15°)

圖 4 1︰1邊坡的最可能滑動(dòng)面(c=15kPa，φ=10°)Fig. 4 The most dangerous sliding surface in 1︰1 slope(c=15kPa，φ=10°)

5 結(jié) 論

對(duì)于巖質(zhì)或土質(zhì)邊坡，只要其材料組成大致均勻，就可以簡(jiǎn)化成均質(zhì)的楔形彈性體，若坡體只受到重力作用，則可利用彈性力學(xué)的逆解法，求出坡體的各個(gè)應(yīng)力分量。本文基于坡體為庫侖材料的假定，提出了一種任意滑動(dòng)面的均質(zhì)邊坡滑動(dòng)穩(wěn)定分析的極限平衡方法，該方法具有以下特點(diǎn)：

(1)物理意義比較明確，假設(shè)條件較少，在求解過程中也不存在數(shù)值計(jì)算上的困難。

(2)通過算例分析表明，本文方法所得到的邊坡穩(wěn)定安全系數(shù)與簡(jiǎn)化的Bishop法和有限元強(qiáng)度折減法的結(jié)果非常接近，且基本介于有限元法和簡(jiǎn)化Bishop法之間。

(3)計(jì)算得到的最危險(xiǎn)滑動(dòng)面可以反映上陡下緩的特點(diǎn)，滑動(dòng)面是向上凹的光滑曲線，且在坡頂?shù)男螤钆c實(shí)際情形更加接近，可以很好地解釋坡頂所出現(xiàn)的陡峭裂縫。

當(dāng)坡體的材料組成不均勻，甚至有多層土以及某些軟弱土層時(shí)，難以求出應(yīng)力的解析解，但可以利用有限元等數(shù)值方法求出坡體內(nèi)的應(yīng)力，再利用本文對(duì)滑動(dòng)面形狀的處理方式，便可解決這種復(fù)雜情形下的坡體穩(wěn)定分析問題。