均質(zhì)土坡漸進(jìn)性破壞分析及變分法驗(yàn)證①

劉 蕓，李 林

(1.同濟(jì)大學(xué)土木工程學(xué)院地下建筑與工程系，上海200092;2.巖土及地下工程教育部重點(diǎn)實(shí)驗(yàn)室(同濟(jì)大學(xué))，上海200092)

0 引 言

隨著我國基礎(chǔ)建設(shè)投資力度的增大，邊坡工程的穩(wěn)定問題也日漸突出.每年我國發(fā)生的邊坡失穩(wěn)災(zāi)害不計(jì)其數(shù)，嚴(yán)重威脅人民的生命財(cái)產(chǎn)安全［1］.因此，為加強(qiáng)對邊坡及其作用形式的認(rèn)識，正確理解邊坡漸進(jìn)性失穩(wěn)的機(jī)理具有重大的現(xiàn)實(shí)指導(dǎo)意義.

就目前而言，分析邊坡穩(wěn)定性的主要方法有兩種，分別為傳統(tǒng)的基于極限平衡法和有限元數(shù)值模擬法［2～4］.但前者的分析方法中，需要做出許多近似假設(shè)［5～6］，因此結(jié)果中往往得不到滑移體內(nèi)土體的變形特性以及應(yīng)力遷移規(guī)律，也得不到支護(hù)結(jié)構(gòu)對土坡變形及穩(wěn)定性的影響.近年來，計(jì)算機(jī)技術(shù)不斷發(fā)展，基于有限元的強(qiáng)度折減系數(shù)法也開始受到重視.它可以靈活地結(jié)合強(qiáng)度折減法和三維彈塑性有限元分析法，在已知評判指標(biāo)的前提下，通過不斷地調(diào)整穩(wěn)定性折減系數(shù)，求出土坡發(fā)生破壞前提下的最小穩(wěn)定安全系數(shù)，進(jìn)而分析判定邊坡穩(wěn)定性.而邊坡穩(wěn)定性分析的變分法以變分原理和極限平衡法為基礎(chǔ)研究邊坡的穩(wěn)定性問題［7］.上述兩種方法在模擬邊坡漸進(jìn)破壞過程以及確定邊坡滑裂面位置方面具有不可比擬的優(yōu)勢.

本文以上海地區(qū)某均質(zhì)土坡工程為背景，構(gòu)建二維有限元數(shù)值模型，模擬土坡的破裂區(qū)域的漸進(jìn)性破壞發(fā)展過程，從而根據(jù)塑性區(qū)范圍搜索出潛在危險(xiǎn)的破裂滑動面.以變分法為基礎(chǔ)，研究平面應(yīng)變狀態(tài)下邊坡穩(wěn)定的平衡方程，求解出最危險(xiǎn)滑動面的函數(shù)表達(dá)式，并與數(shù)值模擬結(jié)果進(jìn)行了對比驗(yàn)證.

1 強(qiáng)度折減法

1.1 強(qiáng)度折減法基本原理

從Zienkiewicz［8］(1975)提出強(qiáng)度折減概念至今，絕大多數(shù)研究者都是沿著他的最初的思路進(jìn)行探討:

土的實(shí)際強(qiáng)度指標(biāo)為c 與φ，折減系數(shù)為Fr，則強(qiáng)度折減為:

令

則

參數(shù)φr與cr即為折減的強(qiáng)度指標(biāo)，將它們代入有限元計(jì)算，若Fr＞1，則有限元計(jì)算的位移就比實(shí)際大，應(yīng)力水平也大，破壞單元要多.令Fr從1.0 開始逐步增大，則算得的位移逐步增大，破壞單元逐步增多，乃至最后達(dá)到整體失穩(wěn).

邊坡失穩(wěn)的判據(jù)直接影響強(qiáng)度折減法計(jì)算的準(zhǔn)確性，現(xiàn)有的失穩(wěn)判據(jù)大致可分為3 種［9～11］:①折減后的強(qiáng)度參數(shù)使得計(jì)算不能收斂;②坡體內(nèi)塑性區(qū)從坡腳到坡頂貫通;③坡體中的特征點(diǎn)位移或應(yīng)變發(fā)生突變且無限發(fā)展.

1.2 邊坡漸進(jìn)性破壞計(jì)算思路

根據(jù)有限元強(qiáng)度折減法和土坡漸進(jìn)性破壞規(guī)律，動態(tài)模擬邊坡漸進(jìn)破壞過程.具體建模流程為:

①確定邊坡的地質(zhì)概況、土體物理力學(xué)參數(shù)、初始應(yīng)力場，根據(jù)上述信息構(gòu)建相對應(yīng)的數(shù)值計(jì)算模型.

②在折減系數(shù)Fr=1 時(shí)，進(jìn)行邊坡的彈塑性力學(xué)計(jì)算，利用屈服接近度指標(biāo)判斷坡體單元是否發(fā)生破損，如果沒有破損區(qū)域的發(fā)生，則相應(yīng)酌情增加折減系數(shù)Fr的取值，再重新進(jìn)行彈塑性力學(xué)分析，直至土坡發(fā)生局部破損為止停止.

③確定土坡漸進(jìn)性破壞過程中的塑性破損區(qū)，由折減系數(shù)Fr折減局部破損區(qū)內(nèi)的強(qiáng)度參數(shù)，然后將新的折減系數(shù)代入前次的折減有限元計(jì)算程序中，重新進(jìn)行彈塑性力學(xué)模擬計(jì)算.

④參照步驟①～③，隨著折減次數(shù)的不斷增多以及折減系數(shù)Fr的逐漸增大，土坡塑性破損區(qū)不斷發(fā)生擴(kuò)展，當(dāng)土坡塑性破損區(qū)完全貫通時(shí)，土質(zhì)邊坡即達(dá)到極限失穩(wěn)狀態(tài)，至此土坡的漸進(jìn)性破壞結(jié)束.

2 算例分析

2.1 有限元模型構(gòu)建

上海地區(qū)某土坡高H=10.0m，土體容重為18 kN·m－3，粘聚力c=28kPa，內(nèi)摩擦角5°.采用四節(jié)點(diǎn)平面應(yīng)變單元，邊坡左右兩側(cè)設(shè)置為水平約束，底面設(shè)置為水平和豎向約束，建立二維有限元模型.模型如圖1 所示.

圖1 計(jì)算模型示意圖

2.2 結(jié)果分析

按照邊坡漸進(jìn)破壞計(jì)算步驟，首先取折減系數(shù)Fr=1，此時(shí)模型所有單元都未出現(xiàn)破損，如圖2 所示.

圖2 未折減時(shí)土坡變形分布

調(diào)用第1 次折減計(jì)算完成的程序進(jìn)行第2 次折減計(jì)算，此時(shí)折減系數(shù)Fr=1.02(圖3(a))，針對根據(jù)該折減模型，再次進(jìn)行彈塑性力學(xué)計(jì)算.

按照上節(jié)所述的計(jì)算流程，此時(shí)可以對折減系數(shù)Fr進(jìn)行反復(fù)自動搜索，折減計(jì)算過程中土坡破損區(qū)演化如圖3 所示.在第7 次折減計(jì)算時(shí)(此時(shí)折減系數(shù)Fr=1.15)，有限元模擬計(jì)算結(jié)果正好顯示出土坡完全出現(xiàn)貫通塑性破損區(qū)，此時(shí)土坡處于極限失穩(wěn)狀態(tài).由圖可見潛在破裂滑動面在起始時(shí)期變化發(fā)展比較緩慢，而在后期階段滑動面破裂貫通區(qū)快速發(fā)展，這一演化過程也符合邊坡滑坡的發(fā)展特點(diǎn).

圖3 折減計(jì)算過程中破損區(qū)演化

3 變分法驗(yàn)證

邊坡穩(wěn)定性分析的變分法計(jì)算，其關(guān)鍵在于利用最危險(xiǎn)滑動面函數(shù)和最危險(xiǎn)滑動面上的正應(yīng)力分布函數(shù)建立土條的靜力平衡方程，通過構(gòu)造相應(yīng)泛函并求解其極值，得到最危險(xiǎn)滑動面函數(shù)和最危險(xiǎn)滑動面上的正應(yīng)力分布函數(shù)的表達(dá)式，進(jìn)而研究邊坡的穩(wěn)定性問題.

3.1 土坡穩(wěn)定性分析模型

為驗(yàn)證基于強(qiáng)度折減法的均質(zhì)土坡最危險(xiǎn)滑裂面位置，以簡單邊坡坡角A 為原點(diǎn)建立直角坐標(biāo)系，如圖4 所示.

圖4 直角坐標(biāo)系邊坡計(jì)算模型

邊坡坡比1:n，高度H，均質(zhì)邊坡土體黏聚力c，內(nèi)摩擦角φ，容重γ，坡面函數(shù)ym(x)可寫為兩分段函數(shù)，假設(shè)滑動面函數(shù)為y(x).取微元體土條進(jìn)行受力分析，有

則:

土體重力:

式中，σ，τ 分別為沿滑動面函數(shù)y(x)分布的正應(yīng)力和剪應(yīng)力，α 為τ 與水平面的夾角，l 為滑動面長度.

本文考慮均質(zhì)邊坡在土體自重、滑動面上的正應(yīng)力和剪應(yīng)力作用下達(dá)到平衡，則平衡方程為:

滑動面貫穿邊坡，端點(diǎn)A(x1，y1)，C(x2，y2)，土體屈服準(zhǔn)則選用Mohr－Coulomb 準(zhǔn)則，即

代入式(7)，有

參考變分原理［12］和文獻(xiàn)［13］，對長度和應(yīng)力無量綱化，記

則有:

邊坡分析模型轉(zhuǎn)化為圖5 所示.

圖5 邊坡計(jì)算模型

3.2 最危險(xiǎn)滑動面變分法求解

由平衡方程知，最危險(xiǎn)滑動面函數(shù)Y(X)及沿該滑動面分布的正應(yīng)力函數(shù)σ(X)應(yīng)滿足式(10).記泛函

式中λ1，λ2為待定常數(shù).邊坡最危險(xiǎn)滑動面問題即轉(zhuǎn)化為泛函的極值問題:尋找Y(X)及σ(X)使G 取極值.由變分原理，上述問題的解Y(X)、σ(X)應(yīng)滿足歐拉方程.

由式(14)可以確定滑動面及其正應(yīng)力分布的控制微分方程，選用Mohr－Coulomb 屈服準(zhǔn)則時(shí)，邊坡最危險(xiǎn)滑動面與最危險(xiǎn)滑動面上的正應(yīng)力分布無關(guān).求解式(14)得:

由式(15)即可確定滑動面函數(shù)Y(X)和沿該滑動面分布的正應(yīng)力函數(shù)σ(X)，分2 種情況討論極值函數(shù).

(1)當(dāng)λ2=0 時(shí)，由式(15)得

當(dāng)μ=tanφ 為常數(shù)，即邊坡為均質(zhì)情況時(shí)，式(16)確定的最危險(xiǎn)滑動面函數(shù)Y(X)為直線，表達(dá)式為

為便于計(jì)算，引入關(guān)于未知極點(diǎn)O(X0，Y0)的極坐標(biāo)系，如圖5 所示.令

式中，

代之入式(15)中

即

且由式(18)和式(19)得:

式中C1，C2為積分常數(shù).

式(23)和式(24)即為極坐標(biāo)系下最危險(xiǎn)滑動面及與之相對應(yīng)的正應(yīng)力表達(dá)式，為方便編程計(jì)算，將式(10)按照式(18)、(19)進(jìn)行坐標(biāo)變換，得

且由邊界條件，得

邊坡坡面函數(shù)表達(dá)式為:

邊坡B 點(diǎn)坐標(biāo)為(nH，H)，(n，1)，(Rm(θ)，arctan((λ1－nλ2)/(1+λ2))).假設(shè)邊坡最危險(xiǎn)滑動面一端點(diǎn)坐標(biāo)A 為(0，0)，即最危險(xiǎn)滑動面經(jīng)過坡角［14］，則

得

聯(lián)立式(25)和式(26)，得

式中:

式(29)待求量共有5 個(gè):待定常數(shù)λ1，λ2，積分常數(shù)C1，C2及未知坐標(biāo)θ2，給定均質(zhì)邊坡坡比1:n、坡高H、土體黏聚力c、內(nèi)摩擦角φ 和容重γ 參數(shù)后，按照式(29)5 個(gè)方程即可求解.

3.3 結(jié)果分析

輸入均質(zhì)土坡已知參數(shù):坡比1:n，坡高H，土體黏 聚 力 c，內(nèi) 摩 擦 角 φ 和 容 重 γ，利 用Mathematicas 軟件編程求解式(29)，結(jié)果見表1.

表1 土坡參數(shù)及計(jì)算結(jié)果

利用坐標(biāo)變換式(18)和(19)將邊坡最危險(xiǎn)滑動面轉(zhuǎn)換到直角坐標(biāo)系中.由此，即利用變分法得到了均質(zhì)邊坡的最危險(xiǎn)滑動面，極坐標(biāo)系下邊坡最危險(xiǎn)滑動面的表達(dá)式(23)說明其為對數(shù)螺旋面.并將理論結(jié)果與數(shù)值模擬結(jié)果進(jìn)行對比，如圖6 所示.

圖6 最危險(xiǎn)滑裂面對比圖

從圖6 中可以看出由強(qiáng)度折減法確定的土坡最危險(xiǎn)滑裂面與采用變分法計(jì)算的解析結(jié)果趨勢大致相同，均為圓弧形滑裂面.但需要指出的是，運(yùn)用變分法計(jì)算的最危險(xiǎn)滑裂面位置要更靠近坡面，并且圓弧形滑裂面的曲率較強(qiáng)度折減法計(jì)算所得曲率要稍大.對比分析結(jié)果還表明，采用強(qiáng)度折減法模擬均質(zhì)土坡漸進(jìn)破壞過程及確定其最危險(xiǎn)破裂面是準(zhǔn)確、可行的.

4 結(jié) 論

本文以上海地區(qū)某均質(zhì)土坡工程為背景，構(gòu)建了二維有限元數(shù)值模型，動態(tài)模擬了土坡的塑性區(qū)發(fā)展和漸進(jìn)破壞過程.以變分法為基礎(chǔ)，研究平面應(yīng)變狀態(tài)下邊坡穩(wěn)定的平衡方程，求解出最危險(xiǎn)滑動面的函數(shù)表達(dá)式，并與數(shù)值模擬結(jié)果進(jìn)行了對比驗(yàn)證.得出以下主要結(jié)論:

(1)隨著有限元模型折減次數(shù)的不斷增加，穩(wěn)定性折減系數(shù)Fr 也在逐漸增大，土質(zhì)邊坡塑性破損區(qū)不斷動態(tài)發(fā)展破壞，邊坡的破裂面不斷發(fā)展形成，土體的穩(wěn)定性逐漸下降，當(dāng)土坡滑裂面貫通時(shí)，土坡達(dá)到極限失穩(wěn)狀態(tài)，漸進(jìn)破壞完成.

(2)邊坡穩(wěn)定性分析的變分法計(jì)算，其關(guān)鍵在于利用最危險(xiǎn)滑動面函數(shù)和最危險(xiǎn)滑動面上的正應(yīng)力分布函數(shù)建立土條的靜力平衡方程，通過構(gòu)造相應(yīng)泛函并求解其極值，得到最危險(xiǎn)滑動面函數(shù)和最危險(xiǎn)滑動面上的正應(yīng)力分布函數(shù)的表達(dá)式.

(3)由強(qiáng)度折減法確定的土坡最危險(xiǎn)滑裂面與采用變分法計(jì)算的解析結(jié)果趨勢大致相同，均為圓弧形滑裂面.運(yùn)用變分法計(jì)算的最危險(xiǎn)滑裂面位置要更靠近坡面，且圓弧形滑裂面的曲率較強(qiáng)度折減法計(jì)算所得曲率要稍大.

［1］ 鄭穎人，趙尚毅，時(shí)為民.邊坡穩(wěn)定性分析的一些新進(jìn)展［J］.地下空間，2001，21(4):263－271.

［2］ 丁樺，張均鋒，鄭哲敏.關(guān)于邊坡穩(wěn)定分析的通用條分法的探討［J］.巖石力學(xué)與工程學(xué)報(bào)，2004，23(21):3684－3688.

［3］ 張士兵，王練柱，張建.邊坡穩(wěn)定性彈塑性大變形有限元強(qiáng)度折減分析［J］.巖石力學(xué)與工程學(xué)報(bào)，2004，23(1):4463－4467.

［4］ 于玉貞，林鴻州，李榮建，等.非穩(wěn)定滲流條件下非飽和土邊坡穩(wěn)定分析［J］.巖土力學(xué)，2008，29(11):2892－2898.

［5］ 許建聰，尚岳全，陳侃福，等.順層滑彈塑性接觸有限元穩(wěn)定性分析［J］.巖石力學(xué)與工程學(xué)報(bào)，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］ 陳力華，靳曉光.有限元強(qiáng)度折減法中邊坡三種失效判據(jù)的適用性研究［J］.土木工程學(xué)報(bào)，2012，45(9):136－146.

［11］ 梁慶國，李德武.有限元強(qiáng)度折減法中邊坡三種失效判據(jù)的適用性研究［J］.巖土力學(xué)，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］.中國公路學(xué)報(bào)，2007，20(4):7－11.