基于應變強度分布準則的邊坡穩(wěn)定性分析

陳正垚，李世海，馮 春

(1.石家莊鐵道大學，河北石家莊 050043;2.中國科學院力學研究所，北京 100190)

0 引言

現(xiàn)有邊坡穩(wěn)定分析的理論與方法主要有以下幾種［1］:一是工程地質分析法;二是極限平衡分析法;三是極限分析法;四是數(shù)值分析方法;五是可靠性分析方法。目前較為常見的是極限平衡法、極限分析法和有限單元法。極限平衡法［2］是當前巖土工程界廣泛應用的方法，該方法基于剛塑性假設，且邊坡滑動面被簡化為圓弧形;對一般性邊坡，由于地質分層以及材料參數(shù)的多樣性，通過假定或簡化模型進行分析計算，其結果是不理想的。有限元強度折減法［3-5］無需事先假定滑動面的形狀和位置，只需對巖土體的強度參數(shù)不斷的折減，使邊坡巖土體因抗剪強度不能抵抗剪切應力而發(fā)生破壞，從而得到滑動面及其相應的安全系數(shù)。這使得有限元強度折減法與傳統(tǒng)的極限平衡法相比具有很大的優(yōu)勢，之所以沒有得到廣泛的應用，究其根源，有限元計算的精度問題是重要的原因之一。影響有限元計算精度的因素有很多，比如，本構模型、巖土力學參數(shù)、網格的疏密、邊坡失穩(wěn)的評判準則、迭代次數(shù)等。

本構模型作為影響有限元計算精度的重要因素之一，一直是專家們研究的重要課題。吳順川等［6］提出了遍布節(jié)理模型，認為巖體經強度折減后潛在破壞可能首先出現(xiàn)在巖體中或沿節(jié)理面或二者同時破壞，可同時考慮巖塊和節(jié)理屬性。錢建固等［7］在形分叉理論中引入非共軸彈塑性模型，構建了平面應變狀態(tài)下的土體軟化模型。楊光華等［8］基于Duncan-Chang非線性本構模型，提出在彈性階段對彈性模量也進行折減的變模量彈塑性模型強度折減法。藍航等［9］根據幾何損傷理論中的裂隙張量和FLAC-3D莫爾-庫侖模型的塑性流動格式建立考慮初始損傷的節(jié)理巖體彈塑性本構模型。姚仰平等［10］提出了土的統(tǒng)一硬化模型，該模型是基于劍橋模型，以變換應力方法和統(tǒng)一硬化參數(shù)為基本要素建立起來的彈塑性本構模型體系，它通過采用基于SMP準則、Lade準則或廣義非線性強度理論的變換應力方法，實現(xiàn)了模型與強度準則的有機結合和模型的三維化。

本文基于有限元和離散元有機結合的計算軟件CDEM［11-13］，通過引入能真實刻畫巖土體漸進破壞過程的應變強度分布模型［14］，對均質邊坡的穩(wěn)定性進行了分析，提出了在數(shù)值分析中能體現(xiàn)歷史沉積作用及土體應變硬化特征的土體強度隨應力水平的增加而提高的強化公式。研究了主要參數(shù)對邊坡穩(wěn)定性、滑動面位置、裂縫發(fā)展范圍以及破裂比的影響?？偨Y得出了一種用破裂比評判邊坡穩(wěn)定性的新方法。

1 應變強度分布準則的基本原理

李世海、周東等［14］提出的應變強度分布模型，從細觀層面上描述了巖土材料的損傷破裂機制。該模型假設巖土材料的破裂是其應變達到一定程度造成的，且材料中代表性體積單元的應變強度服從某種概率分布。以一個代表性體積單元為例，對單元施加剪切應變后，單元中低應變強度的部分被剪斷，服從庫倫摩擦定律，高應變強度的部分保持彈性，單元可以分為斷裂及連續(xù)兩個區(qū)域;隨著剪切應變的增加，不同應變強度的部分依次發(fā)生斷裂，單元的斷裂域不斷增大，當單元產生貫通的剪切面，即表明代表性體積單元已經完全剪斷。

CDEM中可通過單元邊界的斷裂實現(xiàn)材料漸進破壞過程的模擬。因此，文章在CDEM的單元邊界上引入上述應變強度分布模型式(1)、(2)。

其中，F(xiàn)n、Fs為界面上的法向、切向彈簧力，α、β為拉伸、剪切完整度因子，Kn=EA/L、Ks=GA/L為界面上的法向、切向彈簧剛度，Δun(大于0表示拉伸位移，小于0表示壓縮位移)、Δus為法向、切向彈簧位移，E、G為彈性模量及剪切模量，L為界面的特征厚度，A為界面的特征面積，φ為內摩擦角，F(xiàn)'s=-KsΔ us為不考慮損傷情況下的切向彈簧力。

圖1 應變強度分布準則的剪應力－剪應變曲線(出自引文［14］)Fig.1 Shear stress-strain curve corresponding to strain distribution criterion

拉伸、剪切完整度因子α、β與界面上特征拉伸應變ε、特征剪切應變γ之間存在函數(shù)關系，如式3、4。其中 ε = Δun/L，γ = Δus/L，εmin、γmin為出現(xiàn)拉伸、剪切損傷時的最小應變值，εmax、γmax為完整度因子為零(完全破壞)時的拉伸、剪切應變值。F(ε)、F(γ)為與應變相關的完整度概率，介于0與1之間。

F(ε)、F(γ)的概率分布可以選擇不同的模式，如均勻分布、正態(tài)分布、韋伯分布等。以均勻分布為例，F(xiàn)(ε)、F(γ)可表述為5 式。

將理論公式給定具體參數(shù)得到如圖1所示的剪應力-剪應變曲線。

應變強度分布準則的主控參數(shù)為內摩擦角φ和最大剪應變強度γmax以及最小剪應變強度γmin，該準則可以很好的表述材料的屈服和軟化現(xiàn)象，當應變強度區(qū)間變化時可以得到脆斷模型、軟化模型以及理想彈塑性模型。在該準則下，細觀的材料屬性沒有變化，只是不同狀態(tài)下材料各部分作用方式發(fā)生變化，便可以得到復雜的宏觀材料力學行為。

2 邊坡穩(wěn)定性分析

2.1 模型尺寸和邊界條件

假定邊坡沿縱向有足夠的長度，認為是平面應變問題。采用Gmsh［15］進行建模和網格劃分，采用基于有限元和離散元有機結合的計算軟件CDEM進行數(shù)值計算，選用三角形平面應變單元建立邊坡有限元模型，通過單元邊界的斷裂，實現(xiàn)邊坡失穩(wěn)破壞過程的模擬。模型的邊界條件為，約束左側和右側邊界的X方向位移，約束底部的X和Y方向的位移。主要研究的參數(shù)是內摩擦角φ和最大剪應變強度γmax以及最小剪應變強度 γmin。此處將 γmax和 γmin做關聯(lián)處理，令 γ0=γmax=γmin×102，研究 φ 和 γ0對邊坡穩(wěn)定性的影響(圖2，表1)。

圖2 邊坡模型尺寸Fig.2 Model size

表1 巖土和模型參數(shù)Table 1 Parameters of geotechnical and model

2.2 應變分布準則強化公式

第一階段，令塊體和界面都是線彈性模型，算至穩(wěn)定。第二階段將界面改為應變強度分布準則模型，繼續(xù)算穩(wěn)定。

計算的結果顯示:第一階段邊坡內部會產生較大的應變，導致第二階段計算一開始邊坡內部較深部位就會有很多界面的應變超過γ0從而破壞，進而使內部產生大量的裂縫(圖3)，這些裂縫逐漸發(fā)展，最終導致邊坡失穩(wěn)，且失穩(wěn)時幾乎所有的界面都破壞。

圖3 邊坡內部裂縫圖Fig.3 Cracks of internal slope

運用這種方法，由于第一階段計算產生的應變是隨應力水平的增大而增大的，因此邊坡深部的應變總是大于邊坡表面的應變，在第二階段應用應變強度分布準則計算時，總是深部先達到應變極限值γ0從而破壞，深部完全破壞后，裂縫才會向表面發(fā)展。這樣顯然是不合理的。

在這里，引入一種強化機制，假設隨著應力水平的增大，土體的強度也會增大。其強化公式為:

對于具體的 γmax和 γmin，由于 γ0= γmax= γmin×102，所以6式可以寫成:

運用該強化公式，第二階段得到的邊坡臨界狀態(tài)的裂縫發(fā)展情況(圖4)?？梢钥吹剑谶吰率Х€(wěn)過程中，裂縫主要集中在滑動帶附近。當滑動帶中裂縫貫通并發(fā)展到一定程度后，邊坡失穩(wěn)破壞。

圖4 邊坡臨界狀態(tài)裂縫圖Fig.4 Slope cracks of critical state

2.3 破裂比的定義

為了定量的研究邊坡失穩(wěn)過程中裂縫的發(fā)展情況，定義數(shù)值分析中破壞的界面數(shù)為破裂數(shù)d，破裂比，其中d1是當前階段的破裂數(shù)，d2為標準的破裂數(shù)。在本節(jié)中取d2為模型總的界面數(shù)。對某一給定的工程邊坡，當模型建成網格劃分好后，可以通過軟件統(tǒng)計出該模型總的界面數(shù)，即為標準破裂數(shù)。

2.4 內摩擦角對邊坡失穩(wěn)的影響

令γ0=210-3，研究φ對邊坡失穩(wěn)的影響。將邊坡破壞時，位移超過坡高2%(本例中為0.2m)的部分標記出來，得邊坡滑動范圍圖圖5?？梢缘贸觯瑑饶Σ两菍瑒用娴奈恢闷鹬饕刂谱饔?，內摩擦角越大，滑動范圍越小。統(tǒng)計不同內摩擦角時滑動的土方量，得出滑動方量與總土方量的比值——方量比，得到tanφ與方量比的關系曲線(圖6)。圖6說明方量比隨著內摩擦角的增大而減小，當內摩擦角趨近于0時，方量比趨近于1，當內摩擦角很大時，方量比趨近于0。

統(tǒng)計不同內摩擦角時的破裂比，得到tanφ與破裂比的關系曲線(圖7)。由圖7可知，隨著內摩擦角的增大，破裂數(shù)逐漸減少，破裂比逐漸減小。當內摩擦角趨近于0時，破裂比趨近于0.8，當內摩擦角很大時，破裂比趨近于0。并且在 φ～(5°～25°)時，tanφ與破裂比呈線性關系。

圖6 tanφ與方量比關系曲線Fig.6 Curves of tanφ and square ratio

綜上可知，內摩擦角主要控制滑動面位置，也即邊坡的滑動范圍。內摩擦角越小，邊坡的滑動范圍越大。

2.5 γ0對邊坡失穩(wěn)的影響

令φ=15°，研究γ0對邊坡失穩(wěn)的影響。將邊坡失穩(wěn)時，位移超過坡高2%(本例中為0.2m)的部分標記出來，得如下邊坡滑動范圍圖(圖8)并統(tǒng)計方量比(圖9)。圖8中，對于不同γ0滑動面的位置幾乎不變，可以得出γ0對滑動面的位置影響很小，圖9說明，γ0對邊坡失穩(wěn)時的滑動方量幾乎沒有影響。令φ=23°得出不同γ0對應的邊坡失穩(wěn)時的裂縫發(fā)展圖(圖10)。由圖10得出γ0主要控制的是裂縫的發(fā)展程度，γ0值越小，裂縫發(fā)展程度越大。由圖11(γ0與破裂比的關系曲線)可知，隨著γ0的增大，邊坡破裂比減小，兩者關系大體呈線性。綜上可知，在邊坡失穩(wěn)過程中γ0不影響滑動面位置，主要控制的是破裂比的大小也即裂縫的發(fā)展程度。

圖7 tanφ與破裂比關系曲線Fig.7 Curves of tanφ and rupture ratio

圖8 不同γ0時邊坡滑動范圍圖Fig.8 Slope sliding range of different γ0

2.6 內摩擦角和γ0值對邊坡穩(wěn)定性和破裂比的影響

改變內摩擦角和γ0，得到邊坡在不同tanφ值情況下對應的臨界γ0值，從而得出對特定邊坡由tanφ和γ0控制的邊坡臨界狀態(tài)曲線(圖12)。當已知某邊坡的tanφ和γ0時，可以通過該曲線得知邊坡是否穩(wěn)定。位于曲線之上，邊坡穩(wěn)定;位于曲線之下，不穩(wěn)定;在曲線上，則邊坡處于臨界狀態(tài)。

下面通過改變內摩擦角和γ0值研究不同tanφγ0組合情況下，邊坡的破裂情況。巖土參數(shù)同表1，通過數(shù)值分析得不同tanφ時，破裂比隨γ0變化圖(圖13)，不同γ0時，破裂比隨tanφ變化圖(圖14)。擬合得出tanφ-γ0-破裂比曲面圖(圖15)。

圖9 γ0與方量比關系曲線Fig.9 Curves of γ0and square ratio

圖10 不同γ0時邊坡裂縫發(fā)展圖Fig.10 Slope cracks of different γ0

圖11 γ0與破裂比關系曲線Fig.11 Curves of γ0and rupture ratio

圖12 邊坡臨界狀態(tài)曲線Fig.12 Slope curve of critical state

圖13 不同tanφ時，破裂比隨γ0變化圖Fig.13 Curves of γ0and rupture ratio in different tanφ

圖14 不同γ0時，破裂比隨tanφ變化圖Fig.14 Curves of tanφ and rupture ratio in different γ0

圖13說明:對同一tanφ，破裂比隨γ0的增大而減小，并且在臨界值時，會有明顯的突變，破裂比從較大的值突變到0。不同的tanφ決定了滑動范圍的大小，也就決定了破裂比的最大值，tanφ值越小，破裂比的最大值越大。圖14說明:對同一γ0值，破裂比隨tanφ的增大而減小。對同一 tanφ，隨著 γ0的增大，裂縫發(fā)展程度減小，其破裂比也減小。

圖 15 tanφ-γ0-破裂比曲面圖Fig.15 Surfaces of tanφ，γ0and damage ratio

3 運用破裂比評價邊坡穩(wěn)定性的可行性

在第二節(jié)中，研究了應變強度分布準則下，參數(shù)tanφ和γ0對邊坡滑動范圍、裂縫發(fā)展程度和破裂比的影響。據此提出以下評價邊坡穩(wěn)定性的方法。通過實地勘察，得到邊坡的φ和γ0值的大體范圍，以及地表裂縫的發(fā)展情況，對重要的坡體可通過鉆孔了解坡體內部裂縫的發(fā)展情況。在勘察得到的φ和γ0值的范圍內取值，進行數(shù)值分析，得出多組參數(shù)下，坡體裂縫的發(fā)展情況，與實地勘察得到的裂縫發(fā)展情況進行對比，找出比較接近的算例若干組。保持上述若干組中每組的φ值不變，折減γ0使坡體裂縫發(fā)展與實際最接近，此時的破裂數(shù)為當前破裂數(shù)d1，繼續(xù)折減γ0使坡體失穩(wěn)破壞，得到該滑動面所對應的最大破裂數(shù)d2。則破裂比(破裂比越接近1，表明邊坡越靠近臨界狀態(tài))。上述破裂比中的最大值，表征了坡體當前所處的最不利狀態(tài)，根據該最不利狀態(tài)下破裂比的具體數(shù)值，即可確定該邊坡所處的穩(wěn)定性級別。

該方法在坡體地質條件比較簡單的情況下，得出的結果比較準確;當?shù)刭|條件復雜，特別是坡體分層較多且層間巖土性質相差較大時，在用數(shù)值分析模擬坡體裂縫發(fā)展的階段，不同層的參數(shù)有各自的取值范圍，組合后工作量會增加很多，要得出裂縫發(fā)展情況與實際近似的算例，難度也會增加。對于大型坡體，其內部裂縫的發(fā)展情況無法全面的掌握，對于最后的結果也會有影響。

對于層理較多的坡體，可以主要研究薄弱層，其他的層取具有代表性的幾組參數(shù)。對大型坡體，可以重點勘探主要的裂縫，輔以坡體表面的裂縫與數(shù)值結果進行對比分析。

4 結論

本文基于有限元和離散元有機結合的計算軟件CDEM，通過引入能真實刻畫巖土體漸進破壞過程的應變強度分布模型，對均質邊坡的穩(wěn)定性進行了分析。在數(shù)值分析中，針對自然界中土體的應變硬化效應引入一種土體強度隨應力水平的增加而提高的強化機制，從而較好地模擬了邊坡失穩(wěn)過程中，裂縫的發(fā)生和發(fā)展過程。

研究了內摩擦角和γ0對邊坡穩(wěn)定性的影響，得出了可以判斷特定邊坡穩(wěn)定性的邊坡臨界狀態(tài)曲線。通過單因素研究，得出內摩擦角對滑動面的位置起主要控制作用，內摩擦角越大，滑動范圍越小。γ0對裂縫的發(fā)展起主要控制作用，γ0值越小，裂縫發(fā)展程度越大。得出了tanφ-γ0-破裂比曲面圖。

提出了運用破裂比評價邊坡穩(wěn)定性的方法，該方法通過邊坡表面及內部的裂縫信息反演獲得應變強度分布準則需要的參數(shù)(內摩擦角、應變強度等)，進而通過分析不同潛在滑動模式下的破裂比給出邊坡所處的最不利狀態(tài)，最后根據最不利狀態(tài)下邊坡破裂比的具體數(shù)值給出邊坡的穩(wěn)定性級別。

［1］楊天鴻，張鋒春，于慶磊，等.露天礦高陡邊坡穩(wěn)定性研究現(xiàn)狀及發(fā)展趨勢［J］.巖土力學，2011，32(5):1437-1472.YANG Tianhong，ZHANG Fengcun，YU Qinglei，et al.Research situation of open-pit mining high and steep slope stability and its developing trend［J］.Rock and Soil Mechanic，2011，32(5):1437-1472.

［2］張發(fā)明，陳祖煜，彌宏亮.三維極限平衡理論及其在塊體穩(wěn)定分析中的應用［J］.水文地質工程地質，2002，4:33-35.ZHANG Faming，CHEN Zuyu，MI Hongliang.3-D limit equilibrium theory and its application in block stability analysis［J］.Hydrogeology and Engineering Geology，2002，4:33-35.

［3］李紅衛(wèi)，馬惠民，張忠平.強度折減法在高含水滑坡穩(wěn)定性分析中的應用［J］.中國地質災害與防治學報，2009，20(3):27-30.LI Hongwei，MA Huimin，ZHANG Zhongping.Application of strength reduction method in stability analysis of a landslide with high water content［J］.The Chinese Journal of Geological Hazard and Control，2009，20(3):27-30.

［4］陳昌富，翁敬良.基于廣義Hoek-Brown準則邊坡穩(wěn)定性分析強度折減法［J］.中國地質災害與防治學報，2010，21(1):13-18.CHEN Changfu，WENG Jingliang.Strength reduction method based on generalized Hoek-Brown criterion for slope stability analysis［J］.The Chinese Journal of Geological Hazard and Control，2010，21(1):13-18.

［5］候連成，王笑二，王家偉，等.基于強度折減法的某變電所邊坡穩(wěn)定性分析［J］.水文地質工程地質，2012，39(1):72-74.HOU Liancheng，WANG Xiaoer，WANG Jiawei，et al.Slope stability analysisby strength reduction fora substation［J］.Hydrogeology and Engineering Geology，2012，39(1):72-74.

［6］吳順川，金愛兵，高永濤.基于遍布節(jié)理模型的邊坡穩(wěn)定性強度折減法分析［J］.巖土力學，2006，27(4):537-542.WU Shunchuan，JIN Aibing，GAO Yongtao.Slope stability analysis by strength reduction method based on ubiquitous-joint model［J］.Rock and Soil Mechanic，2006，27(4):537-542.

［7］錢建固，呂璽琳，黃茂松.平面應變狀態(tài)下土體的軟化特性與本構模擬［J］.巖土力學，2009，30(3):617-622.QIAN Jiangu，LV Xilin，HUANG Maosong.Softening characteristics of soils and constitutive modelingunder plane strain condition［J］.Rock and Soil Mechanic，2009，30(3):617-622.

［8］楊光華，張玉成，張有祥.變模量彈塑性強度折減法及其在邊坡穩(wěn)定分析中的應用［J］.巖石力學與工程學報，2009，28(7):1506-1512.YANG Guanghua，ZHANG Yucheng，ZHANG Youxiang.Variable modulus elastoplastic strength reduction method and its application to slope stability analysis［J］.Chinese Journal of Rock Mechanics and Engineering，2009，28(7):1506-1512.

［9］藍航，姚建國，張華興，等.基于 FLAC-3D的節(jié)理巖體采動損傷本構模型的開發(fā)及應用［J］.巖石力學與工程學報，2008，27(3):572-579.LAN Hang，YAO Jianguo，ZHANG Huaxing， et al.Development and application of constitutive model of jointed rock mass damage due to mining based on FLAC3D［J］.Chinese Journal of Rock Mechanics and Engineering，2008，27(3):572-579.

［10］姚仰平，侯偉，羅汀.土的統(tǒng)一硬化模型［J］.巖石力學與工程學報，2009，28(10):2135-2151.YAO Yangping，HOU wei，LUO ting.Unified hardening model foe soils［J］.Chinese Journal of Rock Mechanics and Engineering，2009，28(10):2135-2151.

［11］LI Shihai，ZHAO Manhong，WANG Yuannian，et al.A new numerical method for dem-block and particle model［J］.International Journal of Rock Mechanics and Mining Sciences，2004，41(3):436-436.

［12］Wei Zuoan，Li Shihai，J.G.Wang，et al.A dynamic comprehensive method for landslide control［J］.Engineering Geology，2006，84(1-2):1-11.

［13］Continuum-based Discrete ElementMethod and its Application.Li Shihai 2008.9，Beijing DEM’08.

［14］LI Shihai，ZHOU Dong.Progressive failure constitutive model of fracture plane in geomaterial based on strain strength distribution［J］.International Journal of Solids and Structures，2013，50:570-577.

［15］C.Geuzaine，J.-F.Remacle.A three-dimensional finite elementmesh generatorwith built-in pre-and postprocessing facilities［J］.International Journal for Numerical Methods in Engineering，2009，79(11):1309-1331.