基于網(wǎng)格化應(yīng)力場的邊坡滑移面搜索模型及應(yīng)用

李 忠, 張 曦, 汪紅星

(1. 上海工程技術(shù)大學 城市軌道交通學院, 上海 201620; 2. 蘭州理工大學 土木工程學院, 甘肅 蘭州 730050)

隨著人口不斷增多，國內(nèi)外科技不斷發(fā)展，人們對交通運輸、旅游、住宅等方面的需求持續(xù)增加，使得自然邊坡、人工邊坡被大量使用.邊坡在地震、降水的作用下或人為施工過程中極易發(fā)生崩塌、滑坡、泥石流等自然災(zāi)害.因此，邊坡災(zāi)害防治問題越來越受到人們的關(guān)注，邊坡穩(wěn)定性分析變得尤為重要.正是因為大量學者持續(xù)地重視，邊坡穩(wěn)定性研究工作得到了不斷地發(fā)展.

邊坡穩(wěn)定性分析通常以極限平衡法和有限元法為主.極限平衡法是將邊坡穩(wěn)定問題當作剛體平衡問題來研究，即將組成滑坡體的巖土體視為在復雜狀態(tài)下的剛體.許多學者假定了巖土塊之間接觸面上作用力的方向，分布情況及作用位置等，進而衍生出了多種不同類型的極限平衡法[1].如Fellenius提出的瑞典條分法、Bishop法.之后，Janbu等[2]基于Bishop安全系數(shù)的定義又對條塊受力模式及滑面形狀進行過不同假設(shè)，提出了其他極限平衡方法.近年來，周豪等[3]提出了基于滑面正應(yīng)力修正的極限平衡法.王旭等[4]考慮土體多種條件參數(shù)的隨機性與變異性，在極限平衡法分析邊坡穩(wěn)定性時加入了可靠度理論.除此以外，陳祖煜等[5]提出了一種將二維Spencer法擴展到三維的邊坡穩(wěn)定性分析方法，該法在保證了滑坡體三個方向滿足靜力平衡的同時,還增加了一個整體力矩平衡條件,不再將行界面的條間力假定為水平,也不再將條底的剪力方向假定為平行于主滑平面方向.符貴軍等[6]發(fā)現(xiàn)現(xiàn)有的Sarma在進行邊坡分析時未考慮水平地震作用力的角度，進而提出一種考慮水平地震加速度方位角的改進Sarma法.鄧東平等[7]通過對一般情況下的條柱進行受力分析,從而建立了一種能同時滿足3個力和3個力矩平衡方程的準嚴格三維極限平衡計算公式.但由于這些方法都有著眾多的假定，尤其是假設(shè)了滑移面的形狀，極大限制了極限平衡法的應(yīng)用范圍.近年來隨著計算機的發(fā)展，國內(nèi)外在進行邊坡穩(wěn)定性分析時逐漸開始使用限元法與計算機相結(jié)合的方法，能更好地對滑移面進行確定，進而對邊坡穩(wěn)定性有了更加準確地分析.一些學者已經(jīng)做出了一些研究，如在用有限元分析建立應(yīng)力場的基礎(chǔ)上，采用螞蟻算法搜索出更接近實際情況的滑移面從而對邊坡穩(wěn)定性進行分析[8]，但這種算法有著收斂慢、容易陷入局部最優(yōu)的缺點.除此以外，在已知應(yīng)力場的情況下，使用多種群遺傳算法構(gòu)建適應(yīng)度函數(shù)搜索最危險滑移面，達到邊坡穩(wěn)定性分析的目的[9]；在二維的基礎(chǔ)上又提出了相應(yīng)的基于MPGA的三維邊坡穩(wěn)定性分析方法[10]；在邊坡穩(wěn)定性分析中應(yīng)用一種簡單的遺傳算法搜索臨界非圓弧破壞面，并用Mangnestern-Malk方法求解安全系數(shù)，分析了自然斜坡上既有建筑物和建筑物對地震和附加荷載的影響[11].但遺傳算法有著編程過程較為復雜、交叉率和變異率等參數(shù)通常需要靠經(jīng)驗選擇的缺點.

基于以上分析，使用有限元建立應(yīng)力場，再結(jié)合計算機編程語言對應(yīng)力場進行分析，建立一種基于網(wǎng)格化應(yīng)力場的邊坡滑移面搜索模型能有效避免其他算法存在的不足.該方法不再局限于極限平衡法的種種假設(shè)限制,尤其是滑移面的假設(shè)，能適用于任意應(yīng)力狀態(tài)(極限狀態(tài)或非極限狀態(tài))下邊坡穩(wěn)定性分析，更加接近邊坡的真實受力情況.

1 邊坡穩(wěn)定性分析模型

1.1 建立模型

將數(shù)值應(yīng)力場與MATLAB編程相結(jié)合，在MATLAB軟件中利用坐標在平面內(nèi)將邊坡重新劃分一種單元大小均等的新網(wǎng)格，從而建立一種邊坡穩(wěn)定性分析通用模型.通過編程對數(shù)值應(yīng)力場進行分析，得到邊坡最小安全系數(shù)，從而達到邊坡穩(wěn)定性分析的目的.這種模型的關(guān)鍵在于，利用新網(wǎng)格通過搜索條件為安全系數(shù)的計算提供滑移面.應(yīng)力場與程序兩者相輔相成，緊密結(jié)合，最終得到最危險滑移面及其所對應(yīng)的最小安全系數(shù)，并進行邊坡穩(wěn)定性分析.

1.2 模型實現(xiàn)關(guān)鍵步驟

1) 根據(jù)土層參數(shù)使用MIDAS-GTS有限元分析軟件建立邊坡二維模型，導出模型相關(guān)的節(jié)點及單元坐標等信息.采用非線性靜力分析計算邊坡，計算后通過對結(jié)果進行后處理可以得到邊坡的數(shù)值應(yīng)力場情況.

2) 將保存的模型節(jié)點單元及計算得到的數(shù)值應(yīng)力場相關(guān)數(shù)據(jù)信息以數(shù)值矩陣形式導入MATLAB.在MATLAB中使用編程語言重新建立二維直角坐標系，通過坐標信息將模型重新劃分為新的網(wǎng)格式模型.

3) 建立原單元與新劃分網(wǎng)格模型單元的應(yīng)力對應(yīng)關(guān)系，將新的網(wǎng)格單元以每一豎條分為一組.利用土體單元評判準則在每一豎條里找到最危險的單元(最可能達到屈服條件的單元)進行標記，最終找到并繪制所有最危險單元組成的邊坡最危險滑移面，也就是最小安全系數(shù)所在的滑移面.

4) 根據(jù)最危險滑移面上所有單元的應(yīng)力情況，采用矢量和法計算邊坡最危險滑移面所對應(yīng)的安全系數(shù)，從而達到對二維邊坡進行最危險滑移面搜索及計算安全系數(shù)進行穩(wěn)定性評價的目的，從而對邊坡進行穩(wěn)定性分析.

1.3 模型實現(xiàn)關(guān)鍵問題

1.3.1新網(wǎng)格的劃分

根據(jù)有限元計算模型的尺寸大小重新建立二維直角坐標系，將邊坡寬度和高度所圍成的矩形范圍劃分為大小相等的新單元網(wǎng)格，如圖1所示.在x-y坐標內(nèi)，沿x方向?qū)⒛Ｐ蛣澐譃閍個網(wǎng)格，使得模型被劃分為a個豎條，再將每一豎條沿y方向劃分為b個網(wǎng)格，從而形成新的網(wǎng)格模型.為了使得計算結(jié)果相對準確，新劃分的網(wǎng)格單元尺寸應(yīng)略大于有限元計算時的網(wǎng)格，盡量保證每一個新單元內(nèi)存在一個完整的原有限元單元.

1.3.2單元內(nèi)離散點應(yīng)力計算

1) 確定新單元與原有限元單元的應(yīng)力關(guān)系

如圖2所示，取原有限元模型單元四個角點坐標的最大最小值，即單元的大小范圍，如圖2中矩形ABCD所示，得到搜索矩陣

A=[xmin,xmax,ymin,ymax]

找到新網(wǎng)格式模型單元如圖2中虛線框所示中心點i的坐標同時滿足以下條件的所有單元:

xmin

再根據(jù)SABCD=SA+SB+SC+SD找到滿足條件的唯一對應(yīng)的單元.

2) 計算滑面上點i的主應(yīng)力

在四邊形ABCD內(nèi)任意點i的主應(yīng)力可由四個角點的應(yīng)力通過形函數(shù)按比例等參差值的方法計算得到，可由下式表達：

(1)

上式中形函數(shù)Ln如下：

(2)

1.3.3土體單元判別準則及構(gòu)建判別函數(shù)

邊坡發(fā)生破壞時會在坡體內(nèi)形成一條相對軟弱的帶，稱之為滑移面.坡體沿滑移面土體發(fā)生剪切破壞，從而使邊坡發(fā)生失穩(wěn).滑移面是由無數(shù)發(fā)生破壞的土體單元組成的，所以可通過搜索出最有可能發(fā)生破壞的土體單元，從而進一步確定邊坡的最危險滑移面.通常采用摩爾-庫倫準則評判土體單元破壞問題[12]，即土體沿某一平面的剪切應(yīng)力與黏聚力的差值和正應(yīng)力的比值達到土體內(nèi)摩擦角的正切值時，表示土體單元在這個平面處于臨界破壞狀態(tài)，具體情況見下式：

(3)

在式(3)的基礎(chǔ)上構(gòu)建判別函數(shù)f如下式所示:

(4)

以該式作為靶單元的搜索標準，計算每一豎條所有單元的f值，則有：當f<1時土體單元未發(fā)生破壞；當f≥1時土體單元發(fā)生破壞.考慮到同一豎條中可能存在多個單元發(fā)生破壞，此時取f值最大的單元為該豎條最易發(fā)生破壞的單元(最可能發(fā)生屈服單元)進行標記，進而找到所有豎條的最危險單元所組成的邊坡最危險滑移面.

1.3.4安全系數(shù)的計算

在平面問題中，平面內(nèi)任意一點的應(yīng)力狀態(tài)，可用3個主應(yīng)力σ1、σ2、σ3來表示[12].故而邊坡滑移面上任一點的正應(yīng)力和切應(yīng)力可用下式計算：

(5)

式中：σi、τi為滑面上i點處的正應(yīng)力和切應(yīng)力；φ表示土體內(nèi)摩擦角.

矢量和法指邊坡在外部荷載的作用下，潛在滑移面上的極限抗滑力的總和與滑動力的總和之比[13].矢量和法由葛修潤院士提出，現(xiàn)在被越來越多的學者應(yīng)用于二維以及三維邊坡的穩(wěn)定性分析中，其中二維的矢量和法安全系數(shù)計算模型如圖3所示，安全系數(shù)由式(6)可得[14]：

圖3 矢量和法安全系數(shù)計算模型Fig.3 Vector sum method safety factor model

τfi=c+σitanφ

(6)

式中：τfi為i點的抗剪強度，其抗剪強度符合摩爾-庫倫定律，可由式(7)表示[15]：

(7)

式中:Δli表示滑面上i微段的長度；αi表示i微段與水平方向的夾角；n表示滑面的微段總數(shù)；θ表示投影方向與水平方向的夾角，可由式(8)表示：

(8)

由于單元劃分較密，所以滑移面上兩點之間連線可以近似的看作一條直線.計算安全系數(shù)時，如圖2所示，可以分別求得滑移面上i點與相鄰兩點i+1和i-1的中點ii和ii-1的坐標，在ii和ii-1所連成的直線上求得i點的應(yīng)力積分，根據(jù)式(6～8)可以求得邊坡的安全系數(shù).

2 典型算例分析

2.1 算例介紹

以澳大利亞計算機協(xié)會(ACADS)所提供的考試題目為例[16].該邊坡為一個均質(zhì)邊坡，坡比為1∶2.該算例通常被國內(nèi)外學者們用來檢驗二維邊坡計算模型的可行性.算例邊坡的尺寸、邊界條件等具體情況如圖4所示.土層參數(shù):黏聚力c=3.0 kPa,內(nèi)摩擦角φ=19.6°,彈性模量E=10 MPa,容重γ=20 kN/m3,泊松比μ=0.3.

利用MADIS-GTS建立的有限元計算模型如圖4所示，寬度方向為40 m，高度方向為20 m.該二維模型荷載只考慮自重，其邊界約束條件為,上邊界自由無約束條件，左右邊界進行水平約束，下邊界進行水平及豎向約束.采用非線性靜力分析，計算得到邊坡的應(yīng)力場情況.該模型共劃分3 400個單元，3 526個節(jié)點.

圖4 邊坡有限元計算模型 (m)

2.2 結(jié)果對比分析

將該算例用極限平衡法計算得到最危險滑移面及最小安全系數(shù)，再將其結(jié)果與本文基于應(yīng)力場的網(wǎng)格式計算模型對該算例進行穩(wěn)定性分析所得的結(jié)果進行對比.對滑移面與安全系數(shù)分別進行比較.滑移面形狀及位置具體情況如圖5所示，兩種方法所得邊坡安全系數(shù)對比結(jié)果見表1.

圖5 滑移面示意圖Fig.5 Schematic diagram of slip surface

表1 安全系數(shù)對比結(jié)果

圖中分別繪出極限平衡法及本文方法所計算的結(jié)果以及邊坡的坡線，從圖5所示滑移面可知，基于應(yīng)力場的滑移面搜索模型所搜索到的滑移面和極限平衡法相比，前者得到的滑移面曲線始終在極限平衡法所得的滑移面曲線附近波動，表明了基于應(yīng)力場的網(wǎng)格式邊坡滑移面搜索模型進行滑移面搜索及邊坡穩(wěn)定性分析的可行性.

從表1可以看出基于應(yīng)力場的滑移面搜索方法計算出的安全系數(shù)和極限平衡法給出的答案相差不大.其中極限平衡法所計算的結(jié)果較小，這是由于極限平衡法的種種假設(shè)限制導致其結(jié)果較為保守，進一步說明了基于應(yīng)力場的網(wǎng)格式滑移面搜索模型進行滑移面搜索的合理性.

3 結(jié)論

1) 網(wǎng)格化應(yīng)力場邊坡滑移面搜索模型無需對滑移面的形狀與位置等進行假定，通過單元土體應(yīng)力分析的方式，在平面內(nèi)搜索出潛在破壞土體單元，構(gòu)建二維滑移面，為二維邊坡潛在最危險滑移面的自動搜索創(chuàng)造條件.該方法的優(yōu)點在于能適用于任意應(yīng)力狀態(tài)(極限狀態(tài)或非極限狀態(tài))下邊坡穩(wěn)定性分析，使得該方法能應(yīng)用于主動加固條件下邊坡穩(wěn)定狀態(tài)評價.

2) 通過對ACADS經(jīng)典算例的計算和分析，將兩種方法計算結(jié)果進行對比分析，驗證了該種方法可以適用于邊坡穩(wěn)定性分析，滑移面形式為不規(guī)則曲線，更符合邊坡的真實受力狀態(tài).

3) 如何將基于應(yīng)力場的二維邊坡滑移面搜索及邊坡穩(wěn)定性分析方法應(yīng)用于復雜應(yīng)力狀態(tài)下的邊坡，及將該方法擴展到三維邊坡的穩(wěn)定性分析，后續(xù)還需進一步深入的研究.