劉 蕓,李 林
(1.同濟大學土木工程學院地下建筑與工程系,上海200092;2.巖土及地下工程教育部重點實驗室(同濟大學),上海200092)
隨著我國基礎建設投資力度的增大,邊坡工程的穩(wěn)定問題也日漸突出.每年我國發(fā)生的邊坡失穩(wěn)災害不計其數(shù),嚴重威脅人民的生命財產安全[1].因此,為加強對邊坡及其作用形式的認識,正確理解邊坡漸進性失穩(wěn)的機理具有重大的現(xiàn)實指導意義.
就目前而言,分析邊坡穩(wěn)定性的主要方法有兩種,分別為傳統(tǒng)的基于極限平衡法和有限元數(shù)值模擬法[2~4].但前者的分析方法中,需要做出許多近似假設[5~6],因此結果中往往得不到滑移體內土體的變形特性以及應力遷移規(guī)律,也得不到支護結構對土坡變形及穩(wěn)定性的影響.近年來,計算機技術不斷發(fā)展,基于有限元的強度折減系數(shù)法也開始受到重視.它可以靈活地結合強度折減法和三維彈塑性有限元分析法,在已知評判指標的前提下,通過不斷地調整穩(wěn)定性折減系數(shù),求出土坡發(fā)生破壞前提下的最小穩(wěn)定安全系數(shù),進而分析判定邊坡穩(wěn)定性.而邊坡穩(wěn)定性分析的變分法以變分原理和極限平衡法為基礎研究邊坡的穩(wěn)定性問題[7].上述兩種方法在模擬邊坡漸進破壞過程以及確定邊坡滑裂面位置方面具有不可比擬的優(yōu)勢.
本文以上海地區(qū)某均質土坡工程為背景,構建二維有限元數(shù)值模型,模擬土坡的破裂區(qū)域的漸進性破壞發(fā)展過程,從而根據(jù)塑性區(qū)范圍搜索出潛在危險的破裂滑動面.以變分法為基礎,研究平面應變狀態(tài)下邊坡穩(wěn)定的平衡方程,求解出最危險滑動面的函數(shù)表達式,并與數(shù)值模擬結果進行了對比驗證.
從Zienkiewicz[8](1975)提出強度折減概念至今,絕大多數(shù)研究者都是沿著他的最初的思路進行探討:
土的實際強度指標為c 與φ,折減系數(shù)為Fr,則強度折減為:
令
則
參數(shù)φr與cr即為折減的強度指標,將它們代入有限元計算,若Fr>1,則有限元計算的位移就比實際大,應力水平也大,破壞單元要多.令Fr從1.0 開始逐步增大,則算得的位移逐步增大,破壞單元逐步增多,乃至最后達到整體失穩(wěn).
邊坡失穩(wěn)的判據(jù)直接影響強度折減法計算的準確性,現(xiàn)有的失穩(wěn)判據(jù)大致可分為3 種[9~11]:①折減后的強度參數(shù)使得計算不能收斂;②坡體內塑性區(qū)從坡腳到坡頂貫通;③坡體中的特征點位移或應變發(fā)生突變且無限發(fā)展.
根據(jù)有限元強度折減法和土坡漸進性破壞規(guī)律,動態(tài)模擬邊坡漸進破壞過程.具體建模流程為:
①確定邊坡的地質概況、土體物理力學參數(shù)、初始應力場,根據(jù)上述信息構建相對應的數(shù)值計算模型.
②在折減系數(shù)Fr=1 時,進行邊坡的彈塑性力學計算,利用屈服接近度指標判斷坡體單元是否發(fā)生破損,如果沒有破損區(qū)域的發(fā)生,則相應酌情增加折減系數(shù)Fr的取值,再重新進行彈塑性力學分析,直至土坡發(fā)生局部破損為止停止.
③確定土坡漸進性破壞過程中的塑性破損區(qū),由折減系數(shù)Fr折減局部破損區(qū)內的強度參數(shù),然后將新的折減系數(shù)代入前次的折減有限元計算程序中,重新進行彈塑性力學模擬計算.
④參照步驟①~③,隨著折減次數(shù)的不斷增多以及折減系數(shù)Fr的逐漸增大,土坡塑性破損區(qū)不斷發(fā)生擴展,當土坡塑性破損區(qū)完全貫通時,土質邊坡即達到極限失穩(wěn)狀態(tài),至此土坡的漸進性破壞結束.
上海地區(qū)某土坡高H=10.0m,土體容重為18 kN·m-3,粘聚力c=28kPa,內摩擦角5°.采用四節(jié)點平面應變單元,邊坡左右兩側設置為水平約束,底面設置為水平和豎向約束,建立二維有限元模型.模型如圖1 所示.
圖1 計算模型示意圖
按照邊坡漸進破壞計算步驟,首先取折減系數(shù)Fr=1,此時模型所有單元都未出現(xiàn)破損,如圖2 所示.
圖2 未折減時土坡變形分布
調用第1 次折減計算完成的程序進行第2 次折減計算,此時折減系數(shù)Fr=1.02(圖3(a)),針對根據(jù)該折減模型,再次進行彈塑性力學計算.
按照上節(jié)所述的計算流程,此時可以對折減系數(shù)Fr進行反復自動搜索,折減計算過程中土坡破損區(qū)演化如圖3 所示.在第7 次折減計算時(此時折減系數(shù)Fr=1.15),有限元模擬計算結果正好顯示出土坡完全出現(xiàn)貫通塑性破損區(qū),此時土坡處于極限失穩(wěn)狀態(tài).由圖可見潛在破裂滑動面在起始時期變化發(fā)展比較緩慢,而在后期階段滑動面破裂貫通區(qū)快速發(fā)展,這一演化過程也符合邊坡滑坡的發(fā)展特點.
圖3 折減計算過程中破損區(qū)演化
邊坡穩(wěn)定性分析的變分法計算,其關鍵在于利用最危險滑動面函數(shù)和最危險滑動面上的正應力分布函數(shù)建立土條的靜力平衡方程,通過構造相應泛函并求解其極值,得到最危險滑動面函數(shù)和最危險滑動面上的正應力分布函數(shù)的表達式,進而研究邊坡的穩(wěn)定性問題.
為驗證基于強度折減法的均質土坡最危險滑裂面位置,以簡單邊坡坡角A 為原點建立直角坐標系,如圖4 所示.
圖4 直角坐標系邊坡計算模型
邊坡坡比1:n,高度H,均質邊坡土體黏聚力c,內摩擦角φ,容重γ,坡面函數(shù)ym(x)可寫為兩分段函數(shù),假設滑動面函數(shù)為y(x).取微元體土條進行受力分析,有
則:
土體重力:
式中,σ,τ 分別為沿滑動面函數(shù)y(x)分布的正應力和剪應力,α 為τ 與水平面的夾角,l 為滑動面長度.
本文考慮均質邊坡在土體自重、滑動面上的正應力和剪應力作用下達到平衡,則平衡方程為:
滑動面貫穿邊坡,端點A(x1,y1),C(x2,y2),土體屈服準則選用Mohr-Coulomb 準則,即
代入式(7),有
參考變分原理[12]和文獻[13],對長度和應力無量綱化,記
則有:
邊坡分析模型轉化為圖5 所示.
圖5 邊坡計算模型
由平衡方程知,最危險滑動面函數(shù)Y(X)及沿該滑動面分布的正應力函數(shù)σ(X)應滿足式(10).記泛函
式中λ1,λ2為待定常數(shù).邊坡最危險滑動面問題即轉化為泛函的極值問題:尋找Y(X)及σ(X)使G 取極值.由變分原理,上述問題的解Y(X)、σ(X)應滿足歐拉方程.
由式(14)可以確定滑動面及其正應力分布的控制微分方程,選用Mohr-Coulomb 屈服準則時,邊坡最危險滑動面與最危險滑動面上的正應力分布無關.求解式(14)得:
由式(15)即可確定滑動面函數(shù)Y(X)和沿該滑動面分布的正應力函數(shù)σ(X),分2 種情況討論極值函數(shù).
(1)當λ2=0 時,由式(15)得
當μ=tanφ 為常數(shù),即邊坡為均質情況時,式(16)確定的最危險滑動面函數(shù)Y(X)為直線,表達式為
為便于計算,引入關于未知極點O(X0,Y0)的極坐標系,如圖5 所示.令
式中,
代之入式(15)中
即
且由式(18)和式(19)得:
式中C1,C2為積分常數(shù).
式(23)和式(24)即為極坐標系下最危險滑動面及與之相對應的正應力表達式,為方便編程計算,將式(10)按照式(18)、(19)進行坐標變換,得
且由邊界條件,得
邊坡坡面函數(shù)表達式為:
邊坡B 點坐標為(nH,H),(n,1),(Rm(θ),arctan((λ1-nλ2)/(1+λ2))).假設邊坡最危險滑動面一端點坐標A 為(0,0),即最危險滑動面經過坡角[14],則
得
聯(lián)立式(25)和式(26),得
式中:
式(29)待求量共有5 個:待定常數(shù)λ1,λ2,積分常數(shù)C1,C2及未知坐標θ2,給定均質邊坡坡比1:n、坡高H、土體黏聚力c、內摩擦角φ 和容重γ 參數(shù)后,按照式(29)5 個方程即可求解.
輸入均質土坡已知參數(shù):坡比1:n,坡高H,土體黏 聚 力 c,內 摩 擦 角 φ 和 容 重 γ,利 用Mathematicas 軟件編程求解式(29),結果見表1.
表1 土坡參數(shù)及計算結果
利用坐標變換式(18)和(19)將邊坡最危險滑動面轉換到直角坐標系中.由此,即利用變分法得到了均質邊坡的最危險滑動面,極坐標系下邊坡最危險滑動面的表達式(23)說明其為對數(shù)螺旋面.并將理論結果與數(shù)值模擬結果進行對比,如圖6 所示.
圖6 最危險滑裂面對比圖
從圖6 中可以看出由強度折減法確定的土坡最危險滑裂面與采用變分法計算的解析結果趨勢大致相同,均為圓弧形滑裂面.但需要指出的是,運用變分法計算的最危險滑裂面位置要更靠近坡面,并且圓弧形滑裂面的曲率較強度折減法計算所得曲率要稍大.對比分析結果還表明,采用強度折減法模擬均質土坡漸進破壞過程及確定其最危險破裂面是準確、可行的.
本文以上海地區(qū)某均質土坡工程為背景,構建了二維有限元數(shù)值模型,動態(tài)模擬了土坡的塑性區(qū)發(fā)展和漸進破壞過程.以變分法為基礎,研究平面應變狀態(tài)下邊坡穩(wěn)定的平衡方程,求解出最危險滑動面的函數(shù)表達式,并與數(shù)值模擬結果進行了對比驗證.得出以下主要結論:
(1)隨著有限元模型折減次數(shù)的不斷增加,穩(wěn)定性折減系數(shù)Fr 也在逐漸增大,土質邊坡塑性破損區(qū)不斷動態(tài)發(fā)展破壞,邊坡的破裂面不斷發(fā)展形成,土體的穩(wěn)定性逐漸下降,當土坡滑裂面貫通時,土坡達到極限失穩(wěn)狀態(tài),漸進破壞完成.
(2)邊坡穩(wěn)定性分析的變分法計算,其關鍵在于利用最危險滑動面函數(shù)和最危險滑動面上的正應力分布函數(shù)建立土條的靜力平衡方程,通過構造相應泛函并求解其極值,得到最危險滑動面函數(shù)和最危險滑動面上的正應力分布函數(shù)的表達式.
(3)由強度折減法確定的土坡最危險滑裂面與采用變分法計算的解析結果趨勢大致相同,均為圓弧形滑裂面.運用變分法計算的最危險滑裂面位置要更靠近坡面,且圓弧形滑裂面的曲率較強度折減法計算所得曲率要稍大.
[1] 鄭穎人,趙尚毅,時為民.邊坡穩(wěn)定性分析的一些新進展[J].地下空間,2001,21(4):263-271.
[2] 丁樺,張均鋒,鄭哲敏.關于邊坡穩(wěn)定分析的通用條分法的探討[J].巖石力學與工程學報,2004,23(21):3684-3688.
[3] 張士兵,王練柱,張建.邊坡穩(wěn)定性彈塑性大變形有限元強度折減分析[J].巖石力學與工程學報,2004,23(1):4463-4467.
[4] 于玉貞,林鴻州,李榮建,等.非穩(wěn)定滲流條件下非飽和土邊坡穩(wěn)定分析[J].巖土力學,2008,29(11):2892-2898.
[5] 許建聰,尚岳全,陳侃福,等.順層滑彈塑性接觸有限元穩(wěn)定性分析[J].巖石力學與工程學報,2005,24(13):2231-2236.
[6] DUNCAN J M.Fellow State of Art:Limit Equilibrium and Finite Element Analysis of Slopes[J].Journal of Geotechnical Engineering,ASCE,1996,122(7):577-596.
[7] YAMAGAMI T,UETA Y.Search for Noncircular Slip Surfaces by the Morgenstern Preice Method[J].Proc.of 6th Int.Conf.on Numerical Methods in Geomechanics,1988,17(1):1335-1340.
[8] ZIENKIEWICZ O,HUMPHESON C,LEWIS R.Associated and Non-associated Visco-plasticity in Soil Mechanice[J].Geotechnique,1975,5(4):671 ~689.
[9] GRIFFITHS D V,LANE P A.Slope Stability Analysis by Finite Elements[J].Geotechnique,1999,49(3):387-403.
[10] 陳力華,靳曉光.有限元強度折減法中邊坡三種失效判據(jù)的適用性研究[J].土木工程學報,2012,45(9):136-146.
[11] 梁慶國,李德武.有限元強度折減法中邊坡三種失效判據(jù)的適用性研究[J].巖土力學,2008,29(11):3053-3 058.
[12] Dov Leshchinsky,San K C.Pseudostatic Seismic Stability of Slopes:Design Charts[J].Journal of Geotechnical Engineering,1994,120(9):1514-1532.
[13] Shen W Y,Teh C I.Practical Solution for Group Stiffness Analysis of Piles[J].Journal of Geotechnical and Geoenvironmental Engineering,2002,128(8):692-698.
[14] 柳厚祥,廖雪,李寧.公路邊坡穩(wěn)定性分析的二維變分方法[J].中國公路學報,2007,20(4):7-11.