Fellenius 法的一種新的解法

邱 超，王平安，劉建元

邊坡具有數(shù)量多、危害性大的特點，一直是巖土工程中研究的熱點和難點問題［1－4］。Fellenius法作為邊坡穩(wěn)定性分析的經(jīng)典方法，雖然該法不能完全滿足靜力平衡條件，但是其安全系數(shù)計算結(jié)果偏于安全［5］，被廣泛應用到邊坡穩(wěn)定性分析當中。

自Fellenius法被提出以來，不斷有學者對其進行研究與改進。Cao J G等［6］采用積分代替條分的方法，將安全系數(shù)表示為坐標和半徑的表達式，采用Powell法對問題進行了優(yōu)化，但是該法存在著得到局部最優(yōu)解的缺陷，蔣斌松等［7－8］采用解析計算的方法，得到了安全系數(shù)的解析表達式，然后利用函數(shù)求極值的方法得到安全系數(shù)的最優(yōu)解，但是其計算形式太過于復雜，不利于工程應用。李闖等［9］運用解析計算的方法得到了Fellenius法的安全系數(shù)關(guān)于邊坡參數(shù)的顯示表達式，該法無需試算與解方程，但是其計算太過于復雜，且在精度上不能保證。

針對以上問題，本文提出了Fellenius法的一種新的解法，采用解析方法得到安全系數(shù)的解析表達式，將豎向中心角作為優(yōu)化變量，將優(yōu)化問題轉(zhuǎn)化為有約束優(yōu)化問題，利用MATLAB中Fmincon函數(shù)對問題進行優(yōu)化，與傳統(tǒng)無約束優(yōu)化問題相比，該優(yōu)化方法可以更加容易且快速地得到優(yōu)化結(jié)果。提出的解析方法可以很好的得到安全系數(shù)的全局最優(yōu)解，其計算形式相對簡單，并不需要進行反復的迭代計算，有效地避免了傳統(tǒng)方法迭代繁瑣和得到局部最優(yōu)解的問題。

1 基本假定

如圖1所示，均質(zhì)土坡，潛在滑動面為圓弧滑動面，滑動面圓心 O，半徑 r0，邊坡高度 H，坡角θ，土體重度γ，內(nèi)摩擦角φ，黏聚力 c。x0、y0分別為邊坡滑動面圓心的橫、縱坐標。D、F分別為過C、A兩點豎直線與潛在滑動面的交點。μ為G點的橫坐標，α1、α2、η、β分別為G、F、E、B 對應于O 點的豎向中心角，φ為邊坡滑動面上任意一點對應于O點的豎向中心角，順時針為正。kv為豎向地震加速度系數(shù)，kh為水平地震加速度系數(shù)，Ti、Ni分別為土條上作用的切向力和正向力，Xi、Xi－1為條間作用的合力，Wi為土條重量。在分析問題之前做如下幾點基本假定：

（1）邊坡安全系數(shù)等于抗滑力與下滑力之比，且邊坡潛在滑動面上安全系數(shù)相同；

（2）地震作用下水平地震力為 khWi，豎向地震力為 kvWi；

（3）邊坡土體滿足摩爾 －庫侖屈服準則，即

式中：τf、σ分別為邊坡破壞時，滑動面上的切應力和正應力。

（4）土條條間合力滿足下列關(guān)系式：

圖1 滑動面參量示意圖

2 安全系數(shù)解析計算

對邊坡土體進行豎直條分，對任意豎直土條DF進行靜力平衡分析，有：

反解出 Ti，Ni有：

安全系數(shù)定義為邊坡潛在滑動面上抗滑力與下滑力之比，可以表示為：

式中：li為任意土條底面長度，n為劃分的土條數(shù)。為了更加精確的表達安全系數(shù)，采用積分的方法來代替條分法。

由圖1的幾何關(guān)系，邊坡滑動面上任意一點的土體上方高度h可以表示為：

將圓弧滑動面的圓心坐標x0，y0和半徑r0用豎向中心角表示，即：

聯(lián)立式（5）、式（6）、式（7），安全系數(shù)可以表示為：

式中：

由式（8）可以看出，安全系數(shù)只與 c／（γh）、tanφ這兩個參數(shù)有關(guān)，對于不同的邊坡，如果這兩個參數(shù)相同，則安全系數(shù)與潛在滑動面均相同；如果參數(shù)c／（γhtanφ）相等，則潛在滑動面一致。對于一定的邊坡，安全系數(shù)只與α1、α2、η、β這4個變量有關(guān)，將這4個變量作為優(yōu)化變量對安全系數(shù)進行優(yōu)化，便可得到其全局最優(yōu)解。

3 MATLAB程序優(yōu)化

考慮到本文優(yōu)化問題為多維約束非線性優(yōu)化問題，本文采用MATLAB中Fmincon函數(shù)求解安全系數(shù)的最小值和其相對應的潛在滑動面，F(xiàn)mincon數(shù)學模型［10］為：

式中：Ax≤ b，線性不等式約束；Aeq．x ＝ beq，非線性等式約束；C（x）≤ 0，非線性不等式約束；Ceq（x）＝0，非線性等式約束；lb≤ x≤ ub，有界約束。

對于小規(guī)模優(yōu)化問題，F(xiàn)mincon采用序列二次規(guī)劃方法，該算法通過將原問題轉(zhuǎn)化為一系列二次規(guī)劃子問題的求解來獲得原問題的優(yōu)化解，對拉格朗日函數(shù)取二次近似，從而提高二次規(guī)劃子問題的精確程度，對非線性較強的問題也能優(yōu)化計算。序列二次規(guī)劃方法的基本原理為：在某個近似解處將原非線性規(guī)劃問題簡化為處理一個二次規(guī)劃問題，求取最優(yōu)解，如果有，則認為是原有非線性問題的最優(yōu)解，否則，用近似解代替構(gòu)成一個新的二次規(guī)劃問題，繼續(xù)迭代。

為了得到邊坡安全系數(shù)的最小值，以及其對應的邊坡潛在滑動面，應用MATLAB軟件，編制了相關(guān)的優(yōu)化程序進行求解，式（8）的安全系數(shù)為中心角α1、α2、η、β的函數(shù)，該問題為非線性約束問題，約束條件為：

在進行優(yōu)化的過程中，需要對初始值進行設(shè)置，初始值的選取對收斂速度或者是否收斂起著關(guān)鍵性的影響，初始值選取的不恰當，優(yōu)化問題的解可能陷入局部最優(yōu)解。在計算的過程中，α1，α2取值接近0，β值的取值范圍可以為1．2～1．5，η可以取α2和β的中間值，發(fā)現(xiàn)經(jīng)過最多3～4次嘗試便可以得到邊坡穩(wěn)定性分析的全局最優(yōu)解。

4 對比計算與分析

4．1 靜力條件下對比分析

為了驗證本文方法和優(yōu)化程序的正確性，對已有典型算例進行對比計算。均質(zhì)土坡［11］，土體重度γ＝18．5 kN／m3，邊坡高度 H ＝8 m，坡度分別為 1∶1和2∶1。為了分析強度參數(shù)對邊坡穩(wěn)定性的影響，土體內(nèi)摩擦角分別取 10°、15°、20°，黏聚力分別取 5 kPa、10 kPa、15 kPa、20 kPa、25 kPa、30 kPa。

本文計算方法得到的結(jié)果與Fellenius條分法計算得到的結(jié)果見表1。

表1 安全系數(shù)對比表

從表1可以看出，兩種方法計算的安全系數(shù)吻合得很好，驗證了提出的解析方法的正確性，同時也表明解析方法是條分法的下限解。利用本文編制的相關(guān)程序，可以很快的得到邊坡穩(wěn)定性分析的結(jié)果，同時也能保證得到全局最優(yōu)解，得到安全系數(shù)取最小值時的α1、α2、η、β值，帶入式（7），便可以得到臨界圓弧滑動面的圓心和半徑，從而確定邊坡潛在滑動面。

4．2 動力條件下對比分析

水平地震加速度系數(shù)與豎向地震加速度系數(shù)之間關(guān)系可以表示為：

式中：ξ表示kv相對于 kh的比例系數(shù)。研究表明［12－14］：一般取ξ ＝ 0．5，kh取 0．0 ～ 0．3 能夠滿足工程的要求。

均質(zhì)土坡，邊坡參數(shù)：邊坡高度 H＝30 m，c／（γH）＝ 0．15，φ ＝ 40°，ξ＝ 0．5。當邊坡安全系數(shù)Fs＝1時對應的水平向加速度系數(shù)即為水平向屈服加速度系數(shù)，記為 kh。對應不同的邊坡坡角，本文計算得到的水平向屈服加速度系數(shù)與已有計算結(jié)果見表2。從表2可以看出，本文計算的水平向屈服加速度系數(shù)與已有計算結(jié)果很接近，但是由于本文考慮了豎向地震作用系數(shù)的影響，二者的計算結(jié)果還是有一點不同，豎向地震加速度系數(shù)的影響不可忽視。同時，由表2可以得到，隨著邊坡坡角的增大，水平向屈服加速度系數(shù)不斷減小。

表2 水平向屈服加速度對比結(jié)果

通過以上算例分析表明，靜力條件下和動力條件下案例分析表明了本文計算結(jié)果和優(yōu)化程序的正確性。

5 結(jié) 論

通過以上分析，得到如下結(jié)論：

（1）提出了Fellenius法的一種新的解法，以豎向中心角為自變量，推導出了安全系數(shù)的解析表達式，并給出了相應的優(yōu)化程序。

（2）隨著邊坡坡角的增大，水平向屈服加速度系數(shù)不斷減小，同時豎向地震加速度系數(shù)的影響不可忽視。

（3）提出的解析法能夠有效避免得到局部最優(yōu)解的缺陷，具有更強的全局搜索能力和收斂能力。

［1］ 鄧東平，李 亮，趙煉恒．地震作用下邊坡穩(wěn)定性分析的擬靜力解法研究［J］．中南大學學報（自然科學版），2014，45（10）：3578-3588．

［2］ 趙 博，趙亞瓊．邊坡穩(wěn)定性能的判定方法研究及其工程應用［J］．水利與建筑工程學報，2017，15（1）：25-29．

［3］ 王華俊，卿翠貴，姚文杰．一種改進的遺傳算法在邊坡工程中的應用［J］．水利與建筑工程學報，2015，13（3）：195-199．

［4］ Baker R，Garber M．Theoretical analysis of the Stabi-lity of slopes［J］．Geotechnique，1978，28（4）：195-411．

［5］ 方玉樹．邊坡穩(wěn)定性分析條分法最小解研究［J］．巖土工程學報，2008，30（3）：331-335．

［6］ Gao J G，Zaman M M．Analytical metnod for analysis of slope stablity［J］．International Jounal for Numerical and Analytical Methods in Geomechanics，1999，23（5）：439-449．

［7］ 蔣斌松，蔡美峰，呂愛鐘．邊坡穩(wěn)定性的解析計算［J］．巖石力學與工程學報，2004，33（6）：2726-2729．

［8］ 蔣斌松，蔡美峰，呂愛鐘．純黏土邊坡穩(wěn)定性解析計算［J］．工程力學，2003，20（5）：204-208．

［9］ 李 闖，趙盛杰，董 曄，等．Fellenius法的解析算法［J］．水利與建筑工程學報，2016，14（5）：202-210．

［10］ 李 婭．MATLAB R2015b最優(yōu)化計算［M］．北京：清華大學出版社，2017．

［11］ Huang S L，Yamasaki K．Slope failure analysis using local minimum factor－ of－ safety approach［J］．Journal of Geotechnical engineering ASCE，1993，119（12）：1974-1989．

［12］ Ausilio E，Conte E，Dente G．Seismic stability analysis of reinforced slopes［J］．Soil Dynamics and Earthquake Engineering，2000，19（3）：159-172．

［13］ Michalowskir R L．Soil reinforcement for seismic design of geotechnical structures［J］．Computers and Geotechnics，1998，23（1）：1-17．

［14］ 中華人民共和國建設(shè)部．建筑抗震設(shè)計規(guī)范：GB 50011—2001［S］．北京：中國建筑工業(yè)出版社，2002．

［15］ Chang C J，Chen W F，Yao J T．Seismic displacements in slopes by limit analysis［J］．Journal of Geotechnical Engineering ASCE，1984，110（7）：860-874．

［16］ 羅 強，趙煉恒，李 亮，等．地震效應和坡頂超載對均質(zhì)土坡穩(wěn)定性影響的擬靜力分析［J］．巖土力學，2010，31（12）：3835-3841．