(武漢工程大學(xué) 土木工程與建筑學(xué)院，湖北 武漢 430074)

1 概述

現(xiàn)有商業(yè)軟件[1]在進(jìn)行邊坡穩(wěn)定性分析時(shí)，一般采用滑面窮盡搜索算法，以獲取全局最小安全系數(shù)及其臨界滑面。在滑面選擇上，有圓弧滑面、多段折線滑面及多次樣條曲線滑面[2]。通常，這些商業(yè)軟件能較好地處理邊坡穩(wěn)定性計(jì)算問題。但是，在面對(duì)具有軟弱下臥層的邊坡時(shí)，圓弧滑面假設(shè)在搜索時(shí)有可能無(wú)法找到全局臨界滑面；多次樣條曲線滑面一般是在圓弧極限滑面的基礎(chǔ)上進(jìn)行局部微調(diào)得到的；此外，除非完全指定折線段的空間組成，商業(yè)軟件通常只能采用三段線搜索臨界滑面，且折線滑面需要工程師進(jìn)行事先的初步預(yù)估，大致選定臨界滑面折點(diǎn)的可能位置及首尾折線傾角范圍，整個(gè)計(jì)算過程存在明顯的人為因素影響，效率不高。

Malkawi于2001年提出了多段折線臨界滑面的Monte Carlo隨機(jī)搜索技術(shù)[3]，其采用逐步逼近策略，在坡體地質(zhì)較為簡(jiǎn)單時(shí)具備全自動(dòng)全局搜尋功能。張魯渝[4-5]將之改良后進(jìn)行了各種邊坡的隨機(jī)轉(zhuǎn)角及常值轉(zhuǎn)角的穩(wěn)定性分析計(jì)算，效果較好。白桃[6]也采用Monte Carlo隨機(jī)搜索技術(shù)進(jìn)行了邊坡等線段及等角度穩(wěn)定性分析計(jì)算，結(jié)果表明此方法在不太復(fù)雜土質(zhì)邊坡中應(yīng)用效果良好。而一旦在復(fù)雜土質(zhì)邊坡中使用時(shí)，基于其從頭至尾的掃描特性，掃描結(jié)果會(huì)依賴于初始滑面位置，易導(dǎo)致計(jì)算結(jié)果離散，使得滑面落入局部極值中。

為解決上述問題，筆者引入模擬退火法(SAA)這一全局極值解算方法，結(jié)合改進(jìn)后的Monte Carlo隨機(jī)搜索技術(shù)，采用Morgenstern-Price方法進(jìn)行含軟弱土層邊坡滑面的全局自動(dòng)搜索。

2 邊坡安全系數(shù)計(jì)算方法

一般而言，在進(jìn)行邊坡穩(wěn)定性的極限平衡法分析時(shí)，滑面條塊都會(huì)在納入邊坡拐點(diǎn)坐標(biāo)的情況下進(jìn)行水平方向上的等分。但在進(jìn)行多層邊坡穩(wěn)定性分析或者土體參數(shù)存在空間變異性時(shí)，條塊等分在計(jì)算上反而存在不便的地方。白桃[7]進(jìn)行了邊坡截面四面形劃分的穩(wěn)定性計(jì)算分析，但是對(duì)于含不規(guī)則土層截面的邊坡無(wú)法進(jìn)行網(wǎng)格劃分。三角形具備良好的網(wǎng)格適應(yīng)能力，本文以三角形網(wǎng)格劃分為例，進(jìn)行邊坡條塊劃分說明。圖1為澳大利亞計(jì)算機(jī)應(yīng)用協(xié)會(huì)(ACADS)發(fā)布的標(biāo)準(zhǔn)考題截面之一[8]。圖1(b)為圖1(a)在邊坡坡腳的局部放大圖。

提取多段折線滑面G(x)與邊坡截面三角網(wǎng)格相交后的所有交點(diǎn)，其與多段折線拐點(diǎn)合并后按x坐標(biāo)排序，就可以獲得多段折線滑面條塊劃分點(diǎn)位。如圖1(b)所示，G(x)上任意相鄰交點(diǎn)(含拐點(diǎn))i與j，其與頂面線S(x)之間的部分就構(gòu)成一個(gè)單獨(dú)的計(jì)算條塊。i、j的中點(diǎn)O與S(x)之間的垂直距離設(shè)定為條塊高度。對(duì)于分層土體，每個(gè)條塊高度范圍內(nèi)的土體重量進(jìn)行疊加計(jì)算。O點(diǎn)所在三角形網(wǎng)格的強(qiáng)度參數(shù)取定為條塊底面的強(qiáng)度參數(shù)，用于滑塊的穩(wěn)定性計(jì)算。

圖1 折線滑面條分法示意圖Figure 1 Diagram of Slices Method using n-segment Polyline

對(duì)于圓弧/非圓弧滑面的安全系數(shù)計(jì)算，數(shù)學(xué)方程較為完備的是Morgenstern-Price方法。為利于程序編制，本文采用朱大勇提出的簡(jiǎn)明Morgenstern-Price法[9]。

3 邊坡臨界滑面的掃描算法

3.1局部搜索掃描算法

圖2 局部臨界滑面搜索掃描算法[6]Figure 2 Searching Strategy for Local Critical Slip Surface

圖2為白桃[6]進(jìn)行多段折線滑面的局部臨界滑面掃描算法，其核心思想就是采用Monte Carlo隨機(jī)方法進(jìn)行等角/等邊逐段搜索最優(yōu)解的方法進(jìn)行邊坡穩(wěn)定的分析計(jì)算。從步驟1至步驟8為一輪，經(jīng)多輪掃描一般可獲取均值邊坡的極限滑面。在步驟2和步驟8中，采用等角度規(guī)避邊坡地面線拐點(diǎn)，其余旋轉(zhuǎn)步驟均采用等長(zhǎng)線段旋轉(zhuǎn)方式。這種方法存在明顯的短板，即面對(duì)含多層軟弱層土質(zhì)邊坡時(shí)，容易陷入局部極值的問題。文獻(xiàn)[6]中有專門的論述。

3.2全局搜索掃描算法

3.2.1多段折線模擬退火法

為解決上述掃描中的局部極值問題，本文通過引入模擬退火法[10]。模擬退火法用于描述系統(tǒng)在溫度t下從能量E(i)狀態(tài)i進(jìn)入具有能量E(j)狀態(tài)j的機(jī)制。t隨算法進(jìn)程逐步遞減，類似于固體退火過程中的溫度角色。計(jì)算持續(xù)進(jìn)行Metropolis算法：“產(chǎn)生新解—判斷—接受/舍棄”的迭代過程，對(duì)應(yīng)著固體在某一恒定溫度下趨于熱平衡的過程。模擬退火算法從某個(gè)初始解出發(fā)，經(jīng)大量解的變換最終求得最優(yōu)化問題的整體最優(yōu)解。本文對(duì)每次局部掃描獲得的極限滑面施加一定大小的抖動(dòng)，以概率接受的方式進(jìn)行全局極限滑面的位置篩選，計(jì)算流程見圖3。

模擬退火法是一種全局近似算法，在溫度下降過程中，系統(tǒng)接受惡化解的概率不斷降低，直到退火溫度后達(dá)到收斂解。在模擬退火的過程中，需要注意以下3個(gè)方面的技術(shù)問題：

a.初始溫度選擇。

退火初始溫度T需能保證平穩(wěn)分布中每個(gè)狀態(tài)的概率趨于相等，如式(1)，本文選取初始溫度為0.5 ℃。

exp[(F1-F2)/T]≈1

(1)

b.溫度下降方法。

研究表明，退火溫度下降方法的選取直接關(guān)系到計(jì)算量的大小和計(jì)算結(jié)果的準(zhǔn)確性。本文選用線性下降方法，降溫系數(shù)α=0.9，見式(2)。

Ti+1=αTi

(2)

c.算法終止條件。

模擬退火算法從初始溫度開始，通過在每一溫度的迭代和溫度的下降，最后達(dá)到終止原則而停止。本文設(shè)定的終止溫度為1×10-5℃，當(dāng)退火溫度達(dá)到終止溫度時(shí)，認(rèn)定掃描結(jié)束，此時(shí)的安全系數(shù)為全局最小安全系數(shù)。

圖3 模擬退火全局掃描計(jì)算流程圖Figure 3 Global searching flowchart using SAA

3.2.2圓弧窮盡掃描算法

為對(duì)比說明，對(duì)窮盡掃描算法[7]進(jìn)行移植，分析本文采用的三角網(wǎng)格邊坡截面，搜索圓弧臨界滑面，掃描示意見圖4。

圖4 三角網(wǎng)格圓弧掃描示意圖Figure 4 Diagram of circular slip surface searching in triangular mesh slope section

圖4中，截面上方的方格交點(diǎn)構(gòu)成搜索滑面的圓心集，L1和L2線段分別表示掃掠圓弧最下端的上限和下限。將L1至L2之間的距離L進(jìn)行n等分，對(duì)于每一個(gè)圓心O，則需要進(jìn)行n次掃描，每次掃描的半徑為：

(3)

式中:RL1為圓心到L1的垂直距離；i為針對(duì)某個(gè)圓心的第i次掃描。

提取圖4中圓弧與網(wǎng)格線的交點(diǎn)，每相鄰兩交點(diǎn)之間的土條作為極限平衡法的一個(gè)計(jì)算單元。圖4中ABCD就是一個(gè)計(jì)算單元，CD所在三角網(wǎng)格強(qiáng)度參數(shù)作為CD線上的強(qiáng)度參數(shù)。對(duì)每個(gè)圓弧滑面，采用朱大勇提出M-P簡(jiǎn)明算法進(jìn)行安全系數(shù)求解計(jì)算[9]。

4 算例

4.1算例1

圖5(a)為澳大利亞計(jì)算機(jī)應(yīng)用協(xié)會(huì)(ACADS)發(fā)布的標(biāo)準(zhǔn)考題截面之一，截面尺寸、各層土體重度γ、粘聚力c及內(nèi)摩擦角φ均在圖中標(biāo)示。首先采用本文所提出方法進(jìn)行邊坡的折線滑面穩(wěn)定性分析，經(jīng)試算可以確定當(dāng)采用8節(jié)點(diǎn)7段線滑面時(shí)安全系數(shù)收斂，此時(shí)安全系數(shù)FsL=1.345。圖5(b)中虛線為任意隨機(jī)生成的8節(jié)點(diǎn)7段線初始滑面線。

(a) 邊坡截面尺寸及物理參數(shù)

(b) 搜索得到的臨界滑面圖5 算例1參數(shù)Figure 5 The Parameters of exapnple 1

同時(shí)，對(duì)本文圓弧窮盡掃描算法分析本文采用的三角網(wǎng)格邊坡截面，搜索圓弧臨界滑面，計(jì)算所得安全系數(shù)Fso=1.370。GeoStudio 2007商業(yè)軟件計(jì)算圓弧滑面安全系數(shù)Fsg=1.389。8節(jié)點(diǎn)臨界滑面及圓弧滑面見圖5(b)。圖5(b)說明本文提出的截面網(wǎng)格劃分方法具有良好的適應(yīng)性，能夠針對(duì)含不規(guī)則形狀土層邊坡進(jìn)行穩(wěn)定性分析。雖然本例中多段折線滑面與圓弧滑面很接近，但也很明顯：多段折線滑面假設(shè)相對(duì)于圓弧滑面假設(shè)而言，能搜索得到更為危險(xiǎn)的臨界滑面，因而更具有使用優(yōu)勢(shì)。為清楚地顯示，后續(xù)示意圖不再畫出三角網(wǎng)格。

4.2算例2

本算例為兩土層中間夾雜一層軟弱斜坡層邊坡[11]，截面尺寸及物理強(qiáng)度參數(shù)見圖6(a)。各滑面計(jì)算安全系數(shù)見表1，各方法得到的臨界滑面見圖6(b)。

(a) 邊坡截面尺寸及物理參數(shù)

表1 各臨界滑面計(jì)算安全系數(shù)Table1 FactorofSafetyforCriticalSlipSurfaces滑面類型安全系數(shù)4節(jié)點(diǎn)0.422多段折線滑面5節(jié)點(diǎn)0.4136節(jié)點(diǎn)0.4117節(jié)點(diǎn)0.410圓弧滑面本文0.425GeoStudio20070.425

結(jié)果表明多段折線滑面在折點(diǎn)數(shù)量n=6時(shí)安全系數(shù)達(dá)到收斂。圖6(b)顯示了隨機(jī)生成的不同節(jié)點(diǎn)數(shù)量初始滑面掃描得到的臨界滑面。多次試算顯示，無(wú)論初始滑面位置如何，在給定節(jié)點(diǎn)數(shù)量的情況下，都可以搜尋得到相應(yīng)的臨界滑面，表明本文所提方法對(duì)初始滑面位置不敏感。只要生成初始滑面，就可以完全由程序自身完成臨界滑面的全自動(dòng)快速搜索。通過對(duì)比多段折線滑面和圓弧滑面，同樣發(fā)現(xiàn)多段折線滑面展現(xiàn)出較圓弧滑面更為優(yōu)異的搜索能力。

4.3算例3

圖7(a)為文獻(xiàn)[6]中未能直接進(jìn)行全局極限滑面搜索的復(fù)雜土層邊坡。在改良前，隨機(jī)生成的初始滑面在經(jīng)歷多輪循環(huán)后極易落入局部極值，即圖7(b)中的局部臨界滑面S′，對(duì)應(yīng)安全系數(shù)Fs′=1.82。因此，為完成此類復(fù)雜邊坡的極限滑面搜索，需要通過引入工程師經(jīng)驗(yàn)，預(yù)判滑面的可能位置，然后比較各可能位置計(jì)算之后的安全系數(shù)值大小，解算過程比較繁瑣。

(a) 邊坡截面尺寸及物理參數(shù)

(b) 搜索得到的臨界滑面圖7 算例3參數(shù)Figure 7 The Parameters of exapnple 3

引入模擬退火法之后，再次對(duì)這一復(fù)雜邊坡進(jìn)行了多次試算。為方便說明，圖7(b)中的初始隨機(jī)滑面選為①②③，分別代表初始滑面位置跨越不同土層區(qū)間。計(jì)算結(jié)果表明，只需要隨機(jī)生成初始滑面的任意位置，而無(wú)論初始滑面落入哪個(gè)區(qū)間，本文算法總能得到收斂解，快速且自動(dòng)搜索得到全局臨界滑面S，對(duì)應(yīng)安全系數(shù)Fs=1.55。

5 結(jié)論

本文采用Monte Carlo等角/等邊逐段掃描與模脫空時(shí)梁段的撓度變化呈縱向分布，有脫空發(fā)生時(shí)梁段的撓度變化呈橫向分布，無(wú)滑道脫空時(shí)梁段縱向兩邊的撓度基本為0，有些許下?lián)系内厔?shì)，中間梁段的撓度有向上翹起的趨勢(shì)；滑道脫空后，從該滑道至臨近區(qū)域段的梁體呈下?lián)馅厔?shì)，其余梁體出現(xiàn)上翹現(xiàn)象。單滑道脫空時(shí)邊滑道脫空引起梁體下?lián)陷^大，多滑道脫空中同側(cè)邊中滑道脫空引起梁體下?lián)厦黠@大于其他脫空工況，達(dá)到6 mm以上。有限元分析結(jié)果表明：箱梁頂推時(shí)，支撐反力及梁體受力對(duì)滑道高程差非常敏感，毫米級(jí)的誤差都將導(dǎo)致支撐反力發(fā)生劇烈變化，梁體局部有應(yīng)力集中的風(fēng)險(xiǎn)，同時(shí)存在箱梁截面發(fā)生翹曲變形的可能。在實(shí)際頂推過程中，需對(duì)滑道高程偏差進(jìn)行實(shí)時(shí)監(jiān)控，通過動(dòng)態(tài)調(diào)整保持各滑道橫向始終處于同一高程，以減輕在頂推過程中由于滑道高程差導(dǎo)致的支撐反力差過大造成梁體局部受力集中的問題。

擬退火法相結(jié)合進(jìn)行含軟弱夾層邊坡的多段折線全局臨界滑面搜索。同時(shí)，編制了圓弧滑面窮盡掃描算法程序作為對(duì)比。以商業(yè)軟件Fortran和Matlab為計(jì)算手段進(jìn)行聯(lián)合編程實(shí)現(xiàn)了Morgenstern-Price極限平衡搜索算法。計(jì)算結(jié)果表明：

a.截面三角網(wǎng)格劃分方法具有良好的適應(yīng)性，能較好地對(duì)不規(guī)則層位幾何形狀進(jìn)行離散，進(jìn)而方便后續(xù)計(jì)算分析。

b.多段折線滑面具有較圓弧滑面更為優(yōu)秀的極限滑面定位功能，其安全系數(shù)隨著滑面段數(shù)的增加會(huì)逐步達(dá)到收斂。

c.不再依賴于工程師經(jīng)驗(yàn)進(jìn)行初始滑面大致范圍選擇，避免了人為因素影響，本文所提算法能快速地、全自動(dòng)地進(jìn)行由弱到強(qiáng)不同復(fù)雜土質(zhì)情況的邊坡穩(wěn)定性計(jì)算分析，并能找到全局臨界滑面與安全系數(shù)。