孫　平，陳　璽，王玉杰

（1.中國(guó)水利水電科學(xué)研究院 巖土工程研究所，北京　100048；2.西安理工大學(xué) 水利水電學(xué)院，陜西 西安　710048）

1　研究背景

極限平衡法與極限分析上限法是工程中解決邊坡穩(wěn)定問題的兩種常用方法。極限平衡法只考慮力與力矩平衡條件，不考慮變形協(xié)調(diào)條件，本質(zhì)上是一種近似方法，需要引入假定才能使問題變?yōu)殪o定。此外，極限平衡法獲得的解不受塑性力學(xué)上限或下限定理的支持，既不是一個(gè)上限解，也不是一個(gè)下限解［1］。極限分析上限法基于塑性力學(xué)的上限定理，通過在滑坡體內(nèi)部構(gòu)筑一個(gè)機(jī)動(dòng)許可的速度場(chǎng)，利用功能平衡方程求解邊坡安全系數(shù)，理論基礎(chǔ)嚴(yán)格。Donald等［2］提出了對(duì)土條進(jìn)行斜條分的極限分析上限方法，該方法將滑坡體離散為一系列具有傾斜界面的條塊，在假定滑坡體的底滑面與條塊界面同時(shí)達(dá)到極限狀態(tài)的條件下，應(yīng)用Mohr-Coulomb相關(guān)聯(lián)流動(dòng)法則建立一個(gè)機(jī)動(dòng)許可的速度場(chǎng)，再利用虛功原理求解安全系數(shù)。陳祖煜等［3］在此基礎(chǔ)上開發(fā)了巖質(zhì)邊坡穩(wěn)定分析上限法程序EMU，極大地推動(dòng)了這一方法在水利水電工程中的應(yīng)用［2-4］。

臨界滑動(dòng)模式的搜索是邊坡穩(wěn)定分析中十分關(guān)鍵的一步。在二維極限平衡法領(lǐng)域，眾多學(xué)者在應(yīng)用最優(yōu)化方法搜索臨界滑裂面方面開展了深入的研究工作，取得了豐富的成果［5-12］。應(yīng)該說，無論是圓弧滑裂面還是任意形狀滑裂面，搜索安全系數(shù)的整體極小值問題已經(jīng)得到較好的解決［5］。在二維極限分析斜條分上限法領(lǐng)域，由于待優(yōu)化變量中包括了滑裂面位置與條塊界面傾角，導(dǎo)致自由度與非線性程度大大增加，尋找安全系數(shù)的整體極小值變得十分困難。Donald［2］以滑面控制點(diǎn)坐標(biāo)與條塊界面傾角作為待優(yōu)化變量，應(yīng)用單純形法與隨機(jī)搜索法相結(jié)合搜索臨界滑動(dòng)模式。吳超等［13］將改進(jìn)遺傳算法應(yīng)用于地基承載力臨界滑動(dòng)模式的求解，但該方法不僅需要人為指定自變量搜索范圍，而且搜索過程中還需要剔除大量不合理的滑動(dòng)模式，計(jì)算效率較低。Leshchinsky［14］提出將離散組合優(yōu)化方法（Discontinuity layout optimization，即DLO）與極限分析上限法相結(jié)合（DLO-LA），尋找復(fù)雜邊坡最小安全系數(shù)對(duì)應(yīng)的臨界滑動(dòng)模式。該方法的基本思想是，在邊坡內(nèi)部與表面預(yù)先布置一系列均勻分布的網(wǎng)格點(diǎn)，試算滑裂面與條塊界面均由網(wǎng)格點(diǎn)的連線組成，并采用動(dòng)態(tài)規(guī)劃法進(jìn)行臨界滑動(dòng)模式的搜索。該處理方式試圖將連續(xù)的優(yōu)化問題轉(zhuǎn)化為離散優(yōu)化問題，使計(jì)算精度與結(jié)果嚴(yán)重依賴于網(wǎng)格點(diǎn)的布置。當(dāng)前，應(yīng)用最優(yōu)化方法搜索斜條分上限法滑動(dòng)模式的研究成果相對(duì)較少，且無法保證在任何情況下都收斂到全局最優(yōu)解［3］。

邊坡穩(wěn)定極限分析斜條分上限法中臨界滑動(dòng)模式的搜索，本質(zhì)上是一個(gè)工程的極小值問題，通過構(gòu)造合理的優(yōu)化模型，將這一工程上的極小值問題轉(zhuǎn)化為數(shù)學(xué)上的極值問題，是十分關(guān)鍵的一步。本文通過對(duì)滑裂面控制點(diǎn)坐標(biāo)與條塊界面傾角引入一系列的約束條件，避免在隨機(jī)搜索過程中生成不合理滑動(dòng)模式，進(jìn)而提出滑裂面通過與不通過軟弱夾層兩種情況下斜條分上限法滑動(dòng)模式的優(yōu)化模型，將該模型與全局優(yōu)化方法結(jié)合，尋找邊坡的臨界滑動(dòng)模式，并通過一系列經(jīng)典算例的對(duì)比分析，對(duì)本文提出的方法的可行性與有效性進(jìn)行驗(yàn)證。

2　邊坡穩(wěn)定極限分析斜條分上限法

圖1　雙塊體滑動(dòng)模式的速度場(chǎng)

極限分析斜條分上限法假定，當(dāng)滑坡體處于極限狀態(tài)時(shí)，底滑面與分界面同時(shí)達(dá)到極限狀態(tài)。本文以圖1所示的兩個(gè)塊體組合的平面滑動(dòng)問題為例進(jìn)行說明。

首先引入破壞面上“組合摩擦力”概念。破壞面上的抗剪力可分為兩部分：一部分為“凝聚力”，其值為ceA，A為破壞面面積；另一部分為法向力N與由N確定的摩擦阻力N tanφe的合力，稱之為“組合摩擦力”。變量中下標(biāo)‘e’表示破壞面上的抗剪強(qiáng)度指標(biāo)tanφ與c經(jīng)安全系數(shù)F折減后的值，即：tanφe=tanφ/F，ce=c/F。

對(duì)左、右塊體分別建立靜力平衡方程，有：

式中：Wl、Wr分別為左、右塊體的體積力；Pl,e、Pr,e、Pj,e分別為左、右塊體底滑面及條塊界面AB上的組合摩擦力；Cl,e、Cr,e、Cj,e分別為左、右塊體底滑面及條塊界面AB上的凝聚力。

極限分析定理假定材料服從相關(guān)聯(lián)流動(dòng)法則［1］?；谶@一假定，已經(jīng)證明Mohr-Coulomb材料發(fā)生剪切破壞時(shí)，破壞面上的塑性速度V與破壞面的夾角為內(nèi)摩擦角φ［1，3］。左、右塊體的塑性速度分別用Vl和Vr表示，與底滑面夾角分別用φl,e和φr,e表示，在條塊界面AB上，左塊體相對(duì)右塊體的塑性速度Vj為：

同理，Vj和AB的夾角為條塊界面的摩擦角φj，e。

可以證明，破壞面上的組合摩擦力與塑性速度垂直。本文結(jié)合圖2進(jìn)行說明。

如圖2所示，設(shè)線段AB為邊坡內(nèi)任意一個(gè)破壞面（即底滑面或條塊界面）。當(dāng)滑坡體達(dá)到極限狀態(tài)時(shí)，作用于該面上的法向力N與抗剪力T滿足Mohr-Cou?lomb屈服準(zhǔn)則，即：

圖2　破壞面上的組合摩擦力P與塑性速度V

式中：φe、ce分別為AB面上經(jīng)安全系數(shù)F折減后的摩擦角與凝聚力；A為AB面的面積。

如前所述，AB面上的組合摩擦力P為法向力N與摩擦阻力N tanφe的合力，則P與破壞面法向方向的夾角為φe。同時(shí)，因塑性速度V與AB面之間夾角為內(nèi)摩擦角φe，故P與V垂直。

令左、右塊體所受的力分別沿塑性速度Vl和Vr做功，由于作用在底滑面和條塊界面上的組合摩擦力、和分別與 Vl、Vr和 Vj垂直，故這些內(nèi)力沿該位移作的功均為零。

根據(jù)虛功原理，令：式（1）×Vl+式（2）×Vr，有

式（5）的標(biāo)量表達(dá)式為：

式中，ψl、ψr分別為左、右塊體的體積力與Vl和Vr的夾角。

式（6）中包含4個(gè)未知量，即Vl、Vr、Vj以及隱含于下標(biāo)‘e’中的安全系數(shù)F。由于Vr、Vj均可表達(dá)成Vl的線性函數(shù)［3］，故等式兩側(cè)的Vl、Vr、Vj均可消去，這樣式（6）中僅包含唯一的未知量F，可采用迭代法進(jìn)行求解。

上述雙塊體滑動(dòng)模式的求解方法可以很方便的推廣到多塊體滑動(dòng)模式。

3　斜條分上限法滑動(dòng)模式的優(yōu)化模型

圖3　邊坡幾何模型

3.1邊坡幾何模型定義與約束條件在如圖3所示的坐標(biāo)系中，一般將邊坡剖面簡(jiǎn)化為由若干線段組成的圖形，定義其幾何模型如下：y=o（x）為坡面線；y=w（x）為軟弱夾層線；y=b（x）為底邊界線；y=s（x）為滑裂面曲線。

在二維邊坡穩(wěn)定分析中，通常用m個(gè)控制點(diǎn)A1，A2，…，Am以直線或光滑曲線連接來模擬任意形狀滑裂面，其坐標(biāo)分別用（x1，y1），（x2，y2），…，（xm，ym）表示，令x1＜x2＜…＜xm。各控制點(diǎn)的條塊界面與坡面線的交點(diǎn)用Bi（i=2，3，…，m-1）表示，規(guī)定條塊界面傾角θ為AiBi與y軸正向的夾角，由正y方向轉(zhuǎn)向正x方向?yàn)檎?/p>

此外，軟弱夾層線y=w（x）被簡(jiǎn)化為由n個(gè)點(diǎn)M1，M2，…，Mn組成的多段線，其坐標(biāo)分別用（wx1，wy1），（wx2，wy2），…，（wxn，wyn）表示，令wx1＜wx2＜…＜wxn。當(dāng)軟弱夾層的厚度與邊坡高度相比可以忽略不計(jì)時(shí)，按無厚度處理［3］。規(guī)定在優(yōu)化的過程中，滑裂面控制點(diǎn)A2與Am-1分別在線段M1M2與Mn-1Mn上移動(dòng)。

邊坡抗滑穩(wěn)定安全系數(shù)F可表示為:

滑裂面的相鄰控制點(diǎn)之間采用等分方式進(jìn)一步細(xì)分條塊，并約定各細(xì)分點(diǎn)的條塊界面為各細(xì)分點(diǎn)與相鄰兩控制點(diǎn)條塊界面交點(diǎn)的連線，如圖3所示。

優(yōu)化過程中為避免構(gòu)造不合理滑動(dòng)模式，引入以下約束條件：

（1）滑裂面為下凸形。

（2）除剪入與剪出段外，滑裂面不能與坡面線相交。

（3）相鄰控制點(diǎn)所在的直線與x軸的夾角αi應(yīng)控制在一個(gè)合理的范圍內(nèi)［5］，即：-45°≤ αi≤ 80°。

（4）相鄰控制點(diǎn)的條塊界面不能在坡體內(nèi)相交。

（5）條塊界面傾角θi應(yīng)在一個(gè)合理的范圍內(nèi)［14］，即：-90°＜θi＜90°。

3.2滑裂面不通過軟弱夾層時(shí)滑動(dòng)模式的構(gòu)造步驟

3.2.1確定滑裂面剪入點(diǎn)與剪出點(diǎn)x坐標(biāo)的變化范圍如圖4所示，根據(jù)邊坡的幾何形狀，滑裂面的剪入點(diǎn)A1與剪出點(diǎn)Am在水平方向的變化范圍分別用［Lmin，Lmax］與［Rmin，Rmax］表示。用無量綱的標(biāo)準(zhǔn)變量表示，有：

當(dāng)點(diǎn)A1與Am的x坐標(biāo)確定后，其y坐標(biāo)可由函數(shù)y=o（x）唯一確定。

3.2.2確定滑裂面中間控制點(diǎn)A2，A3，…，Am-1的x坐標(biāo)的變化范圍為使任意滑裂面構(gòu)造更為靈活，將（x1，xm）進(jìn)行m-2等分（圖5），各分點(diǎn)ai的計(jì)算式為：

滑裂面中間控制點(diǎn)Ai（i=2，…，m-1）x坐標(biāo)的變化范圍為［ai，ai+1］，用無量綱的標(biāo)準(zhǔn)變量表示，有

圖4　滑裂面剪入與剪出點(diǎn)坐標(biāo)的確定　

圖5　滑裂面中間控制點(diǎn)x坐標(biāo)變化范圍

圖6　滑裂面中間控制點(diǎn)y坐標(biāo)變化范圍

3.2.3確定滑裂面中間控制點(diǎn)Ai的y坐標(biāo)的變化范圍［yi，min，yi，max］如圖6所示，過滑裂面控制點(diǎn)Ai（i=2，…，m-1）的豎直線用Li表示。

下限yi，min的確定方法為：

（1）控制點(diǎn)位于底邊界y=b（x）的上方。直線Li與曲線y=b（x）求交，交點(diǎn)的y坐標(biāo)用y1表示；

（2）根據(jù)約束條件1，直線Ai-2Ai-1與Li求交，交點(diǎn)的y坐標(biāo)用y2表示；

（3）根據(jù)約束條件3，過點(diǎn)Ai-1且與水平線夾角為-45°的直線與直線Li求交，交點(diǎn)的y坐標(biāo)用y3表示。則 yi，min=max（y1，y2，y3）。

上限 yi，max的確定方法為：

（1）直線Li與曲線y=o（x）求交，交點(diǎn)的y坐標(biāo)用y4表示；

（2）根據(jù)約束條件1，直線Ai-1Am與直線Li求交，交點(diǎn)的y坐標(biāo)用y5表示；

（3）根據(jù)約束條件2，令P表示x坐標(biāo)在［xi-1，xi］之間的坡面線端點(diǎn)的集合，即：

過點(diǎn)Ai-1與集合P中所有點(diǎn)的直線與Li求交，交點(diǎn)y坐標(biāo)的最小值用y6表示；

（4）根據(jù)約束條件3，過點(diǎn)Ai-1且與水平線夾角為80°的直線與直線Li求交，交點(diǎn)的y坐標(biāo)用y7表示。則yi,max=min（y4，y5，y6，y7）。用無量綱的標(biāo)準(zhǔn)變量表示，有：

3.2.4確定條塊界面傾角θi的變化范圍［θi,min，θi,max］條塊界面傾角的變化范圍如圖7所示。

下限θi,min的確定方法為：

（1）線段AiAi-1與y軸正向的夾角，記為β1；

（2）根據(jù)約束條件5，令β2=-90°；

（3）根據(jù)約束條件4，記線段AiBi-1與y軸正向的夾角，記為β3。則θi,min=max（β1，β2，β3）。

上限θi,max的確定方法為：

（1）根據(jù)約束條件4，令β4=θi-1；

（2）根據(jù)約束條件5，令β5=90°；

（3）線段AiAm與y軸正向的夾角，記為β6。則θi，max=min（β4，β5，β6）。用無量綱的標(biāo)準(zhǔn)變量表示，有:

3.3滑裂面通過軟弱夾層時(shí)滑動(dòng)模式的構(gòu)造步驟

3.3.1確定位于軟弱夾層上的滑裂面如圖8所示，滑裂面控制點(diǎn)A2在線段M1M2上移動(dòng)，其x坐標(biāo)的變化范圍為［wx1，wx2］，用無量綱的標(biāo)準(zhǔn)變量表示，有:

圖7　條塊界面傾角的變化范圍

圖8　滑裂面剪入與剪出段傾角的變化范圍

同理，點(diǎn)Am-1在線段Mn-1Mn上移動(dòng)，其x坐標(biāo)的變化范圍為［wxn-1，wxn］，用無量綱的標(biāo)準(zhǔn)變量表示，有:

3.3.2確定滑裂面剪出段A2A1與x軸負(fù)方向的夾角α的變化范圍［αmin，αmax］滑裂面剪入與剪出段傾角的變化范圍，如圖8所示。

下限αmin的確定方法為：

（1）滑裂面剪出點(diǎn)A1在坡面線上。記坡面線的左端點(diǎn)為E，線段A2E與x軸負(fù)方向的夾角記為α1；

（2）根據(jù)約束條件1，線段A3A2與x軸負(fù)方向的夾角記為α2。則αmin=max（α1，α2）。

上限αmax的確定方法為：

由約束條件3確定上限αmax=80°。用無量綱的標(biāo)準(zhǔn)變量表示，有：

3.3.3確定滑裂面剪入段Am-1Am與x軸正方向的夾角β的變化范圍［βmin，βmax］

下限βmin的確定方法為：

（1）滑裂面控制點(diǎn)Am在坡面線上。記坡面線的右端點(diǎn)為F，線段Am-1F與x軸正方向的夾角記為β1；

（2）根據(jù)約束條件1，線段Am-2Am-1與x軸正向的夾角記為β2。則βmin=max（β1，β2）。

上限βmax的確定方法為：

由約束條件3確定上限βmax=80°。用無量綱的標(biāo)準(zhǔn)變量表示，有：

條塊界面傾角θi的變化范圍的確定方法與前文相同，故不再贅述。

3.4最優(yōu)化數(shù)學(xué)模型與優(yōu)化求解綜上所述，對(duì)于滑裂面不通過軟弱夾層的情況，建立的優(yōu)化模型為：

對(duì)于滑裂面通過軟弱夾層的情況，其最優(yōu)化數(shù)學(xué)模型為：

可以看出，通過建立斜條分上限法滑動(dòng)模式的數(shù)學(xué)模型，滑動(dòng)模式的搜索問題轉(zhuǎn)化為一個(gè)多自由度的有界約束極小值問題。下文將遺傳算法（GA）和粒子群算法（PSO）這兩種在工程中應(yīng)用廣泛的全局優(yōu)化算法與本文方法相結(jié)合，開展臨界滑動(dòng)模式的搜索計(jì)算。

4　算例分析

圖9　邊坡剖面及不同方法獲得的臨界滑動(dòng)模式

4.1算例1：ACADS考核題1（c） 本算例為一個(gè)無軟弱夾層的非均質(zhì)土坡，如圖9所示，各土層的物理力學(xué)參數(shù)如表1所示。文獻(xiàn)［5］與文獻(xiàn)［3］分別采用STAB程序（Bishop簡(jiǎn)化法+單純形法）與EMU程序（斜條分上限法+單純形法）對(duì)該邊坡進(jìn)行了穩(wěn)定性計(jì)算。建立斜條分上限法的優(yōu)化模型時(shí)，滑裂面控制點(diǎn)數(shù)為5，總自由度個(gè)數(shù)為11，各控制點(diǎn)之間采用光滑曲線連接。圖9列出了不同方法獲得的臨界滑裂面位置，表2列出了不同方法的計(jì)算結(jié)果。

由計(jì)算結(jié)果可知，本文方法獲得的臨界滑動(dòng)模式與STAB程序解十分接近，安全系數(shù)也明顯小于EMU程序解。此外，本文方法得到的安全系數(shù)略大于STAB程序解，表明本文方法獲得的是一個(gè)合理的、略大于極限平衡解的上限解。

表1　ACADS考核題1（c）的物理力學(xué)參數(shù)

表2　不同計(jì)算方法的計(jì)算結(jié)果

4.2算例2：ACADS考核題EX3如圖10所示的非均質(zhì)邊坡，在邊坡底部發(fā)育有一層產(chǎn)狀水平、力學(xué)性質(zhì)較差的軟弱夾層，即土層2。各土層的物理力學(xué)性質(zhì)見表3。文獻(xiàn)［5］采用STAB程序（Spencer法+單純形法）對(duì)該邊坡進(jìn)行了穩(wěn)定性分析。由于本算例為典型的沿軟弱夾層滑動(dòng)的問題，本文將土層2按兩種方式處理：（1）土層2作為一種“普通”土層，滑裂面控制點(diǎn)數(shù)為4，自由度為8，滑裂面控制點(diǎn)之間采用直接連接；（2）土層2作為無厚度的軟弱夾層，假定滑裂面經(jīng)過其底邊界，滑裂面控制點(diǎn)為4，自由度為6。為進(jìn)行對(duì)比研究，采用EMU程序（上限解+單純形法）對(duì)該算例進(jìn)行了計(jì)算。圖10～圖11顯示了不同方法獲得的臨界滑裂面，表4列出了不同方法計(jì)算得到的安全系數(shù)。

由計(jì)算結(jié)果可知，針對(duì)土層2的兩種處理方式，基于本文提出的優(yōu)化模型，采用兩種優(yōu)化方法得到的臨界滑動(dòng)模式十分接近，安全系數(shù)略大于STAB程序解，且小于EMU程序解。此外，土層2按第（1）種處理方式得到的安全系數(shù)略小于第（2）種處理方式的計(jì)算結(jié)果，表明土層2的厚度對(duì)邊坡穩(wěn)定安全系數(shù)有一定的影響，但影響不大?？偟目磥?，采用本文方法獲得的是一個(gè)相對(duì)較優(yōu)的、合理的上限解。

表3　ACADS考核題EX3的物理力學(xué)參數(shù)

圖10　土層2作為普通土層時(shí)的臨界滑動(dòng)模式

圖11　土層2作為軟弱夾層時(shí)的臨界滑動(dòng)模式

表4　不同方法獲得的最小安全系數(shù)計(jì)算結(jié)果

4.3算例3：無重力地基極限承載力算例圖12所示為無重量承載垂直表面荷載的地基極限承載力問題。1960年代，Skolovvskii［15］采用滑移線理論提出了這一問題的解析解。地基承載力是土力學(xué)的經(jīng)典穩(wěn)定問題之一，常作為極限分析方法的標(biāo)準(zhǔn)考題而被廣泛引用［3］。在地基承載力領(lǐng)域，通常用加載系數(shù)代替?zhèn)鹘y(tǒng)的安全系數(shù)來表征其安全儲(chǔ)備能力。若邊坡表面作用有荷載q0，通過不斷增加這個(gè)荷載，直至邊坡達(dá)到極限狀態(tài)的荷載q，則定義加載系數(shù)η為：

圖12　無重量地基承載力問題

表5　不同方法獲得的計(jì)算結(jié)果對(duì)比

本算例計(jì)算參數(shù)：基礎(chǔ)寬度B=17m，地基土體內(nèi)摩擦角φ=0°，c=30 kPa，容重γ=0 kN/m3，Prandtl解獲得的地基土極限承載力qu=154.25 kPa。

令q0=qu=154.25 kPa，以加載系數(shù)η作為目標(biāo)函數(shù)，滑裂面控制點(diǎn)個(gè)數(shù)為5，用光滑曲線連接，最右側(cè)滑裂面控制點(diǎn)坐標(biāo)在優(yōu)化過程中保持不動(dòng)，自由度個(gè)數(shù)為10。文獻(xiàn)［3］采用上限解+單純形法對(duì)本算例進(jìn)行了計(jì)算。不同方法得到的計(jì)算結(jié)果如圖12與表5所示。

從計(jì)算結(jié)果可以看出，采用本文方法得到的臨界滑動(dòng)模式中，條塊界面均收斂到一點(diǎn)，加載系數(shù)η小于0.04，與理論解非常接近，且明顯優(yōu)于采用上限解+單純形法的解。

5　結(jié)論

本文提出了一種邊坡穩(wěn)定斜條分上限法滑動(dòng)模式的優(yōu)化模型，該模型針對(duì)滑裂面通過與不通過軟弱夾層兩種情況，對(duì)滑裂面的控制點(diǎn)坐標(biāo)與條塊界面傾角引入一系列的約束條件，使滑動(dòng)模式的優(yōu)化問題轉(zhuǎn)化為一個(gè)具有多自由度的有界約束的極小值問題。

多個(gè)典型算例的驗(yàn)證與對(duì)比分析表明，本文提出的數(shù)學(xué)模型與遺傳算法、粒子群算法等全局優(yōu)化方法相結(jié)合，具有較好的全局收斂性。特別地，對(duì)于破壞機(jī)構(gòu)復(fù)雜的地基承載力問題，本文方法能夠給出一個(gè)合理的、與理論解十分接近的上限解。