邊坡穩(wěn)定極限分析斜條分上限法的全局優(yōu)化方法

孫　平，陳　璽，王玉杰

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

1　研究背景

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

臨界滑動模式的搜索是邊坡穩(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］以滑面控制點坐標與條塊界面傾角作為待優(yōu)化變量，應(yīng)用單純形法與隨機搜索法相結(jié)合搜索臨界滑動模式。吳超等［13］將改進遺傳算法應(yīng)用于地基承載力臨界滑動模式的求解，但該方法不僅需要人為指定自變量搜索范圍，而且搜索過程中還需要剔除大量不合理的滑動模式，計算效率較低。Leshchinsky［14］提出將離散組合優(yōu)化方法（Discontinuity layout optimization，即DLO）與極限分析上限法相結(jié)合（DLO-LA），尋找復(fù)雜邊坡最小安全系數(shù)對應(yīng)的臨界滑動模式。該方法的基本思想是，在邊坡內(nèi)部與表面預(yù)先布置一系列均勻分布的網(wǎng)格點，試算滑裂面與條塊界面均由網(wǎng)格點的連線組成，并采用動態(tài)規(guī)劃法進行臨界滑動模式的搜索。該處理方式試圖將連續(xù)的優(yōu)化問題轉(zhuǎn)化為離散優(yōu)化問題，使計算精度與結(jié)果嚴重依賴于網(wǎng)格點的布置。當前，應(yīng)用最優(yōu)化方法搜索斜條分上限法滑動模式的研究成果相對較少，且無法保證在任何情況下都收斂到全局最優(yōu)解［3］。

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

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

圖1　雙塊體滑動模式的速度場

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

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

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

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

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

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

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

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

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

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

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

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

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

式（5）的標量表達式為：

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

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

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

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

圖3　邊坡幾何模型

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

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

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

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

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

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

（1）滑裂面為下凸形。

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

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

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

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

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

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

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

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

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

圖4　滑裂面剪入與剪出點坐標的確定　

圖5　滑裂面中間控制點x坐標變化范圍

圖6　滑裂面中間控制點y坐標變化范圍

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

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

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

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

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

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

（1）直線Li與曲線y=o（x）求交，交點的y坐標用y4表示；

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

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

過點Ai-1與集合P中所有點的直線與Li求交，交點y坐標的最小值用y6表示；

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

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）。用無量綱的標準變量表示，有:

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

3.3.1確定位于軟弱夾層上的滑裂面如圖8所示，滑裂面控制點A2在線段M1M2上移動，其x坐標的變化范圍為［wx1，wx2］，用無量綱的標準變量表示，有:

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

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

同理，點Am-1在線段Mn-1Mn上移動，其x坐標的變化范圍為［wxn-1，wxn］，用無量綱的標準變量表示，有:

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

下限αmin的確定方法為：

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

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

上限αmax的確定方法為：

由約束條件3確定上限αmax=80°。用無量綱的標準變量表示，有：

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

下限βmin的確定方法為：

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

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

上限βmax的確定方法為：

由約束條件3確定上限βmax=80°。用無量綱的標準變量表示，有：

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

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

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

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

4　算例分析

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

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

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

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

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

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

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

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

圖10　土層2作為普通土層時的臨界滑動模式

圖11　土層2作為軟弱夾層時的臨界滑動模式

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

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

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

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

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

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

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

5　結(jié)論

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

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