周桂云
(1.金陵科技學(xué)院建筑工程學(xué)院,江蘇南京 211169;2.中國人民解放軍理工大學(xué)工程兵工程學(xué)院,江蘇南京 210007)
基于強度折減的邊坡穩(wěn)定安全系數(shù)有限元迭代解法
周桂云1,2
(1.金陵科技學(xué)院建筑工程學(xué)院,江蘇南京 211169;2.中國人民解放軍理工大學(xué)工程兵工程學(xué)院,江蘇南京 210007)
討論有限元強度折減法的失穩(wěn)判據(jù),結(jié)果表明:失穩(wěn)判據(jù)對安全系數(shù)影響較大,對臨界滑動面形狀和位置的確定影響不大。在此基礎(chǔ)上,提出基于強度折減的邊坡穩(wěn)定安全系數(shù)有限元迭代解法,先用強度折減法搜索滑動面,再按滑動面已知邊坡用迭代解法求解穩(wěn)定安全系數(shù)。該方法充分結(jié)合了強度折減法和迭代解法的優(yōu)點,求解邊坡穩(wěn)定安全系數(shù)時效率快、精度高,可為實際工程參考應(yīng)用。
強度折減法;迭代解法;邊坡穩(wěn)定安全系數(shù);有限元法
20世紀(jì)60年代以來,有限元法被廣泛應(yīng)用于邊坡穩(wěn)定分析,已成為近年來邊坡穩(wěn)定分析的新趨勢。有限元法克服了極限平衡法的不足,并有其獨特的優(yōu)點:①能夠充分考慮土體的非線性本構(gòu)關(guān)系、土體的非均質(zhì)性和各向異性等復(fù)雜特性;②能模擬施工過程,可考慮變形體中泥化夾層的滲流效應(yīng)、孔隙水壓力等;③能給出邊坡內(nèi)的應(yīng)力場、位移場,可以了解邊坡的逐步破壞機理和塑性區(qū)的開展情況;④能適應(yīng)復(fù)雜的邊界條件,比如復(fù)雜的幾何形狀、結(jié)構(gòu)物的作用、有支護的邊坡穩(wěn)定;⑤可以真實考慮水庫水位驟降、地震荷載等對邊坡穩(wěn)定性的影響。Duncan[1]對有限元法在邊坡和堤壩穩(wěn)定分析中的應(yīng)用研究進行了詳細(xì)的歸納總結(jié),按照所采用的本構(gòu)模型將它們分為4類,即:①線彈性模型;②折線彈性模型;③雙曲線彈性模型;④彈塑性或彈黏塑性模型。從這些總結(jié)可以看出,大部分研究工作是對邊坡進行有限元分析,對計算范圍內(nèi)各單元上積分點的應(yīng)力進行強度判別,凡應(yīng)力狀態(tài)達到剪切破壞或拉破壞判別標(biāo)準(zhǔn)的部位稱之為破壞區(qū),根據(jù)破壞區(qū)的分布位置和范圍的大小可以對邊坡的穩(wěn)定性進行評價,為邊坡的治理、施工方法提供依據(jù)。隨著計算機的發(fā)展,有限元法作為一種強有力的數(shù)值計算方法,在邊坡穩(wěn)定分析中已經(jīng)取得了大量的成果。目前主要的方法有2類:①有限元迭代解法;②有限元強度折減法。迭代解法和強度折減法在邊坡穩(wěn)定安全系數(shù)的定義上是相同的,不同的是迭代解法只適用于滑動面預(yù)先知道的情況下求解邊坡的穩(wěn)定安全系數(shù),而抗剪強度折減法主要用于解決滑動面未知的情況(當(dāng)然也可以用于滑動面已知的情況)。對滑動面未知情況可以利用破壞時刻的塑性區(qū)分布圖來確定滑動面的位置。2種有限元法各有優(yōu)缺點,本文在對有限元強度折減法的失穩(wěn)判據(jù)進行討論的基礎(chǔ)上,結(jié)合有限元強度折減法和迭代解法,提出基于強度折減的邊坡穩(wěn)定安全系數(shù)有限元迭代解法,使得求解邊坡穩(wěn)定安全系數(shù)時速度快且精度高。
有限元強度折減法將強度折減技術(shù)與彈塑性有限元方法相結(jié)合,首先進行強度折減系數(shù)下邊坡的彈塑性有限元分析,得到邊坡內(nèi)的應(yīng)力場、應(yīng)變場、位移場,然后以位移、廣義剪應(yīng)變等描述變形程度的某種物理量作為評判指標(biāo),定量地描繪邊坡的潛在塑性破壞區(qū)域及其破壞程度與發(fā)展趨勢,據(jù)此基于一定的經(jīng)驗評判準(zhǔn)則確定邊坡的極限平衡狀態(tài),并將由此所確定的相應(yīng)強度折減系數(shù)作為邊坡的穩(wěn)定安全系數(shù)[2-4]。
式中:K為安全系數(shù);c和φ為土體的實際抗剪強度參數(shù);ce和φe為土體達到極限狀態(tài)時的抗剪強度參數(shù)。
有限元強度折減法不需要事先確定滑動面的形狀和位置,并可以按一定的規(guī)則搜索出最危險臨界滑動面,因而能解決復(fù)雜的工程問題。然而,該方法的關(guān)鍵問題是如何根據(jù)有限元計算結(jié)果來判斷邊坡是否處于破壞狀態(tài),即如何定義邊坡的失穩(wěn)判據(jù),至今工程界還沒有統(tǒng)一的標(biāo)準(zhǔn)。目前,主要有3種失穩(wěn)判據(jù):
a.以有限元迭代求解過程的不收斂作為邊坡失穩(wěn)的判據(jù)[5-9]。在非線性有限元迭代求解時,給出迭代收斂標(biāo)準(zhǔn),若在給定的迭代次數(shù)和收斂準(zhǔn)則內(nèi)有限元計算仍未收斂,則認(rèn)為發(fā)生破壞。迭代收斂準(zhǔn)則有位移準(zhǔn)則、失衡力準(zhǔn)則和能量準(zhǔn)則3種。
b.以塑性區(qū)(或等效塑性應(yīng)變)從坡腳到坡頂貫通作為邊坡失穩(wěn)的判據(jù)[10-12]。根據(jù)強度折減彈塑性有限元分析結(jié)果,繪制坡體內(nèi)的塑性應(yīng)變(廣義剪應(yīng)變)分布圖,當(dāng)某一幅值的塑性應(yīng)變(廣義剪應(yīng)變)開展區(qū)貫通時則認(rèn)為邊坡失穩(wěn)破壞。
c.以坡體內(nèi)特征部位位移發(fā)生突變作為邊坡失穩(wěn)的判據(jù)[13-14]。邊坡的變形破壞總會具有一定的位移特性,位移結(jié)果是邊坡失穩(wěn)最直觀的表達形式。根據(jù)有限元計算結(jié)果建立特征點的位移和折減系數(shù)關(guān)系曲線,若位移發(fā)生突變且無限發(fā)展則認(rèn)為邊坡失穩(wěn)破壞。
通過研究分析,上述3種失穩(wěn)判據(jù)各有優(yōu)缺點:由于巖土體的復(fù)雜性,影響迭代計算不收斂的因素很多,如預(yù)先定義的收斂標(biāo)準(zhǔn)、本構(gòu)方程形式、荷載步長、初始地應(yīng)力、土體材料特性等,且收斂準(zhǔn)則的定義具有一定的人為任意性,所以以迭代過程的不收斂作為失穩(wěn)判據(jù)是不合理的,適用性差;以塑性區(qū)貫通作為失穩(wěn)判據(jù)具有較明確的物理意義,然而塑性應(yīng)變的大小受到材料類型的影響,而且以某一定量值作為標(biāo)準(zhǔn)缺乏廣泛的適用性,邊坡是在塑性應(yīng)變貫通前失穩(wěn)還是貫通后失穩(wěn),缺乏依據(jù);以坡體內(nèi)特征部位位移發(fā)生突變作為失穩(wěn)判據(jù)的物理意義明確且直觀,然而,其缺陷是究竟用哪個點以及哪個方向的位移沒有統(tǒng)一的認(rèn)識,如何從曲線上給出安全系數(shù)也沒有明確的方法??紤]到實用性與簡便性,文獻[15][16]建議聯(lián)合采用特征部位位移突變和塑性區(qū)貫通作為失穩(wěn)判據(jù)。
為了分析失穩(wěn)判據(jù)對穩(wěn)定安全系數(shù)和滑動面的影響,假定一均質(zhì)土坡如圖1所示,土的計算參數(shù)為:彈性模量E=2.0×105Pa,泊松比 ν=0.25,密度ρ=2000kg/m3,黏聚力c=1.0×104Pa,摩擦角 φ=20°。
圖1 均質(zhì)土坡計算簡圖
有限元網(wǎng)格剖分如圖2所示,由有限元計算可知,坡體是由下往上逐漸貫通破壞的,所以選擇處于塑性區(qū)范圍內(nèi)的A點為特征點。以有限元迭代求解過程的不收斂為失穩(wěn)判據(jù)時求得穩(wěn)定安全系數(shù)為1.1236,以塑性區(qū)貫通作為失穩(wěn)判據(jù)時求得穩(wěn)定安全系數(shù)為1.1494,以特征部位位移發(fā)生突變?yōu)槭Х€(wěn)判據(jù)求得穩(wěn)定安全系數(shù)為1.1628。圖3為3種失穩(wěn)判據(jù)下邊坡失穩(wěn)時的塑性區(qū)分布。
圖2 均質(zhì)土坡有限元模型
經(jīng)對比分析可知,有限元計算迭代是否收斂與邊坡失穩(wěn)沒有必然的關(guān)系,迭代不收斂時邊坡不一定發(fā)生失穩(wěn);塑性區(qū)貫通是破壞的必要條件,但不是充分條件,塑性區(qū)貫通并不一定意味著破壞;位移突變是破壞的充分條件,邊坡破壞總有一定的位移特征。所以建議結(jié)合塑性區(qū)貫通和位移突變來判別邊坡是否破壞。計算結(jié)果表明,不同失穩(wěn)判據(jù)求得的穩(wěn)定安全系數(shù)不同,且失穩(wěn)判據(jù)對穩(wěn)定安全系數(shù)的影響較大。但是塑性區(qū)分布相差不大,若以塑性區(qū)分布圖上塑性應(yīng)變較大點的連線作為臨界滑動面,則3種失穩(wěn)判據(jù)下的臨界滑動面很接近,說明失穩(wěn)判據(jù)對臨界滑動面形狀和位置的確定影響不大。
圖3 3種失穩(wěn)判據(jù)下的塑性區(qū)分布
針對有限元強度折減法的上述問題,為了避免由于失穩(wěn)判據(jù)對安全系數(shù)的影響,筆者提出改進方法,即基于強度折減的邊坡穩(wěn)定安全系數(shù)有限元迭代解法。強度折減法之所以方便實用是因為此法不需要事先假定滑動面的形式和位置,然而該方法的關(guān)鍵是失穩(wěn)判據(jù)的定義,且失穩(wěn)判據(jù)對穩(wěn)定安全系數(shù)影響較大,對臨界滑動面位置影響不大,而迭代法適用于已知滑動面的情況求解邊坡的穩(wěn)定安全系數(shù)。因此,基于強度折減的邊坡穩(wěn)定安全系數(shù)有限元迭代解法的思路為:對未知滑動面的邊坡先用有限元強度折減法分析,得到臨界滑動面,再用有限元迭代解法按已知滑動面情況求解穩(wěn)定安全系數(shù)。在已知滑動面的有限元迭代求解時,關(guān)于滑動面的單元形式在文獻[17]中有詳細(xì)的討論,提出了常規(guī)的矩形(平面)或立方體(空間)單元的統(tǒng)一模式。對厚度為零的滑動面采用較小厚度取代零厚度使得問題有解。計算分析結(jié)果表明,一般可近似地取結(jié)構(gòu)面的厚度t=1/100l,其中l(wèi)為與滑動面相連的實體單元的最小邊長。滑動面的材料參數(shù)取與滑動面相連的2種材料中較小尺寸相對應(yīng)的材料參數(shù)。滑動面上的應(yīng)力強度準(zhǔn)則一般采用摩爾-庫侖準(zhǔn)則,同時滿足相應(yīng)的抗拉條件。
圖4為一垂直坡覆蓋了一寬度為B的剛性條形基礎(chǔ),基礎(chǔ)中部施加垂直向下的集中荷載F,不考慮自重,右側(cè)和底部施加約束,其理論極限荷載值為F/(Bcu)=2.0,其中cu為土體的黏聚力。由強度折減法的抗滑穩(wěn)定安全系數(shù)定義有F/[B(cu/K*)]=2.0,則該坡體理論安全系數(shù)為K*=2.0Bcu/F。本次計算中取F=1500N,B=5m,cu=200Pa,摩擦角 φ=0°,彈性模量E=1.0×105Pa,泊松比ν=0.495。理論破壞面與水平面成45°角方向,理論安全系數(shù)K*=1.3333。
圖4 垂直剪問題模型尺寸
先用有限元強度折減法搜索滑動面,未知滑動面網(wǎng)格剖分如圖5所示,計算所得塑性區(qū)分布如圖6所示。根據(jù)折減法計算所得的塑性區(qū)分布確定滑動面,作為已知滑動面問題用有限元迭代解法求解穩(wěn)定安全系數(shù),已知滑動面網(wǎng)格剖分如圖7所示,取迭代初始穩(wěn)定安全系數(shù)K1=1.0,經(jīng)過4次迭代計算,結(jié)果為K3=K4=1.3309,與理論解僅相差1.8%,破壞時的位移矢量如圖8所示,破壞滑移變形如圖9所示。
圖5 未知滑動面網(wǎng)格剖分
圖6 塑性區(qū)分布
圖7 已知滑動面網(wǎng)格剖分
圖8 位移矢量
圖9 破壞變形
采用圖1所示的均質(zhì)土坡,土的計算參數(shù)如下:彈性模量E=2.0×105Pa,泊松比 ν=0.25,密度 ρ=2000kg/m3,黏聚力c=1.0×104Pa,摩擦角 φ=20°。
未知滑動面計算時網(wǎng)格剖分如圖10所示,這里網(wǎng)格剖分較密,可以比較準(zhǔn)確地確定滑動面,塑性區(qū)分布如圖11所示,根據(jù)塑性區(qū)分布圖確定滑動面后剖分網(wǎng)格如圖12所示。經(jīng)過3次迭代滿足誤差要求,穩(wěn)定安全系數(shù)K=1.1732,與文獻[18]中不同網(wǎng)格計算所得穩(wěn)定安全系數(shù)基本一致。破壞時的位移矢量如圖13所示,破壞變形如圖14所示。
圖10 未知滑動面網(wǎng)格剖分
圖11 塑性區(qū)分布
圖12 已知滑動面網(wǎng)格剖分
圖13 位移矢量
圖14 破壞變形
對有限元強度折減法的失穩(wěn)判據(jù)進行討論,建議結(jié)合使用塑性區(qū)貫通和位移突變作為邊坡失穩(wěn)判據(jù),并認(rèn)為失穩(wěn)判據(jù)對穩(wěn)定安全系數(shù)影響較大,但對搜索臨界滑動面的形狀和位置影響不大。在此基礎(chǔ)上,提出基于強度折減的邊坡穩(wěn)定安全系數(shù)有限元迭代解法。算例結(jié)果表明,該方法方便有效,充分利用了強度折減法能夠方便準(zhǔn)確地搜索滑動面的優(yōu)點和迭代法計算穩(wěn)定安全系數(shù)效率高、精度高的優(yōu)點,是一種求解邊坡穩(wěn)定安全系數(shù)的快速、簡便方法,可為實際工程應(yīng)用提供參考。
[1]DUNCAN J M.State of the art:limit equilibrium and finite element analysis of slopes[J].Journal of Geotechnical Engineering,ASCE,1996,122(7):477-596.
[2]ZIENKIEWICZ O C,HUMPHESON C,LEWIS R W.Associated and non-associated visco-plasticity and plasticity in soil mechanics[J].Geotechnique,1975,25(4):671-689.
[3]宋二祥.土工結(jié)構(gòu)安全系數(shù)的有限元計算[J].巖土工程學(xué)報,1997,19(2):1-7.
[4]張魯渝,鄭穎人,趙尚毅,等.有限元強度折減系數(shù)法計算土坡穩(wěn)定安全系數(shù)的精度研究[J].水利學(xué)報,2003,34(1):21-27.
[5]GRIFFITH D V,LANE P A.Slope stability analysis by finite element[J].Geotechnique,1999,49(3):387-403.
[6]DAWSON E M,ROTH W H,DRESCHER A.Slope stability analysis by strength reduction[J].Geotechnique,1999,49(6):835-840.
[7]趙尚毅,鄭穎人,時衛(wèi)明,等.用有限元強度折減法求邊坡穩(wěn)定安全系數(shù)[J].巖土工程學(xué)報,2002,24(3):343-346.
[8]趙尚毅.用有限元強度折減法及其在土坡與巖坡中的應(yīng)用[D].重慶:后勤工程學(xué)院,2004.
[9]鄭穎人,趙尚毅.用有限元強度折減法求滑(邊)坡支擋結(jié)構(gòu)的內(nèi)力[J].巖石力學(xué)與工程學(xué)報,2004,23(20):3552-3558.
[10]鄭宏,李春光,李焯芬,等.求解安全系數(shù)的有限元法[J].巖土工程學(xué)報,2002,24(5):626-628.
[11]連鎮(zhèn)營,韓國城,孔憲京.強度折減有限元法研究開挖邊坡的穩(wěn)定性[J].巖土工程學(xué)報,2001,23(4):406-411.
[12]欒茂田,武亞軍,年延凱.強度折減有限元法中邊坡失穩(wěn)的塑性區(qū)判據(jù)及其應(yīng)用[J].防災(zāi)減災(zāi)工程學(xué)報,2003,23(3):1-8.
[13]宋二祥.土工結(jié)構(gòu)安全系數(shù)的有限元計算[J].巖土工程學(xué)報,1997,19(2):1-7.
[14]周桂云,李同春.飽和-非飽和非穩(wěn)定滲流作用下巖質(zhì)邊坡穩(wěn)定性分析[J].水電能源科學(xué),2006,24(5):79-82.
[15]姚緯明,李同春,任旭華,等.帶軟弱結(jié)構(gòu)面的彈塑性有限元分析[J].河海大學(xué)學(xué)報:自然科學(xué)版,1999,27(3):34-38.
[16]劉金龍,欒茂田,趙少飛,等.關(guān)于強度折減有限元方法中邊坡失穩(wěn)判據(jù)的討論[J].巖土力學(xué),2005,26(8):1345-1348.
[17]李同春,盧智靈,姚緯明,等.邊坡抗滑穩(wěn)定安全系數(shù)的有限元迭代解法[J].巖石力學(xué)與工程學(xué)報,2003,22(3):446-450.
[18]張丹.邊坡穩(wěn)定及抗滑樁內(nèi)力分析方法研究[D].南京:河海大學(xué),2004.
Finite element iteration method for safety factor of slope stability based on strength reduction method
ZHOU Gui-yun1,2(1.Institute of Civil Engineering and Architecture of Jinling Institute of Technology,Nanjing211169,China;2.Engineering Institute of PLAUniversity of Science and Technology,Nanjing210007,China)
The instability criterion of strength reductionmethodwas discussed.The instability criterionhad a great influence onthe safety factor,but has a little influence on the shape and position of the critical sliding surface.On this basis,the finite element iteration method for the safety factor of slopesbased on the strength reductionmethod was proposed.The critical sliding surface was searched by use of the strength reductionmethod,and then the safety factor of given slopes of the sliding surface was solved by the finite element iteration method.A conclusion is drawn that the proposed method fully combines the advantages of the strength reduction method and the iteration method.It is fast,accurate and practical in solving the safety factor for the slope stability.
strength reduction method;iteration method;safety factor of slopes;finite element method
TU311.3
A
1006-7647(2010)03-0058-04
10.3880/j.issn.1006-7647.2010.03.015
中國博士后科學(xué)基金(20080431343);金陵科技學(xué)院科研基金(jit-n-2009-014)
周桂云(1980—),女,江蘇泰州人,講師,博士,從事土木工程教學(xué)及科研工作。E-mail:guiyunzhouhhu@yahoo.com.cn
book=94,ebook=275
2009-07-17 編輯:高建群)