基于ABAQUS二次開(kāi)發(fā)平臺(tái)的邊坡有限元強(qiáng)度折減法研究

候 玲， 薛海斌， 周澤華， 劉海偉， 潘 峰

(1.西安理工大學(xué) 理學(xué)院，陜西 西安710054；2.西安理工大學(xué) 巖土工程研究所，陜西 西安710048)

?

基于ABAQUS二次開(kāi)發(fā)平臺(tái)的邊坡有限元強(qiáng)度折減法研究

候 玲1， 薛海斌2， 周澤華2， 劉海偉2， 潘 峰2

(1.西安理工大學(xué) 理學(xué)院，陜西 西安710054；2.西安理工大學(xué) 巖土工程研究所，陜西 西安710048)

有限元強(qiáng)度折減法在邊坡穩(wěn)定性分析中發(fā)揮著舉足輕重的作用，但通常的有限元強(qiáng)度折減法采用二分法來(lái)獲取邊坡的安全系數(shù)，該方法需要多次試算，效率較低。本文借助ABAQUS軟件的二次開(kāi)發(fā)平臺(tái)UFIELD程序，將軟件自帶的場(chǎng)變量定義為折減系數(shù)，建立場(chǎng)變量分別與軟件求解時(shí)步和強(qiáng)度參數(shù)的關(guān)系，方便地實(shí)現(xiàn)了基于場(chǎng)變量的邊坡有限元強(qiáng)度折減法，此方法僅需一次提交計(jì)算分析便可高效地實(shí)現(xiàn)邊坡安全系數(shù)的求解；在基于場(chǎng)變量的邊坡有限元強(qiáng)度折減法基礎(chǔ)上，依據(jù)滑面位置處的等效塑性應(yīng)變最大這一基本原理，利用不同垂直剖面逐點(diǎn)搜索法，實(shí)現(xiàn)了基于最大等效塑性應(yīng)變的邊坡潛在滑動(dòng)面位置精細(xì)化尋找；最后借助兩個(gè)經(jīng)典算例，將此套思路確定的潛在滑動(dòng)面位置、安全系數(shù)分別與傳統(tǒng)極限平衡M-P法的對(duì)應(yīng)結(jié)果進(jìn)行對(duì)比分析，驗(yàn)證了本文思路和方法的有效性。

邊坡有限元強(qiáng)度折減法； 場(chǎng)變量； 最大等效塑性應(yīng)變； 潛在滑動(dòng)面

有限元強(qiáng)度折減法在確定邊坡安全系數(shù)的同時(shí)還可以方便地確定出邊坡的滑動(dòng)面位置，是邊坡穩(wěn)定性分析中的一種重要方法。

自從Zienkiewicz O.C.[1]等提出有限元強(qiáng)度折減法以來(lái)，唐芬[2]、Yuan W.[3]、薛海斌[4]等提出了不同的強(qiáng)度折減方式；楊光華[5]、薛雷[6]和陳國(guó)慶[7]等對(duì)邊坡強(qiáng)度折減的范圍進(jìn)行了研究；張魯渝、鄭穎人[8]等研究了不同屈服準(zhǔn)則對(duì)強(qiáng)度折減安全系數(shù)的影響，并試圖探討這些安全系數(shù)之間的關(guān)系；鄭宏[9]、張培文[10]、楊光華[11]等分別提出了強(qiáng)度折減過(guò)程中其它配套參數(shù)的折減方式；宋二祥[12]、Griffiths D.V.[13]和連鎮(zhèn)營(yíng)[14]分別提出了三種不同的邊坡失穩(wěn)判據(jù)?？傊?，國(guó)內(nèi)外學(xué)者在強(qiáng)度折減法研究方面取得了豐碩的成果，但是這些學(xué)者大都通過(guò)人為地對(duì)強(qiáng)度參數(shù)進(jìn)行折減試算直到邊坡失穩(wěn)來(lái)確定邊坡的安全系數(shù)，每次折減都必須提交一次計(jì)算分析，過(guò)程繁復(fù)。為了在保證計(jì)算精度的同時(shí)，簡(jiǎn)化此計(jì)算過(guò)程，曹先鋒[15]、李寧[16]等分別提出了邊坡穩(wěn)定分析的溫控參數(shù)折減有限元法和基于優(yōu)化理論中的二分法來(lái)確定邊坡的安全系數(shù)。

除邊坡安全系數(shù)的確定之外，邊坡穩(wěn)定性分析的另一個(gè)關(guān)鍵問(wèn)題就是滑面的確定。事實(shí)上強(qiáng)度折減法中給出的是滑動(dòng)面位置所處的范圍而未能精確到一條曲線。在強(qiáng)度折減法計(jì)算結(jié)果基礎(chǔ)上，孫冠華等提出了基于等效塑性應(yīng)變的邊坡滑動(dòng)面搜索方法[17]；李劍在此基礎(chǔ)上借助FLAC3D軟件內(nèi)置的FISH語(yǔ)言實(shí)現(xiàn)了基于最大剪應(yīng)變?cè)隽康倪吰聺撛诨瑒?dòng)面搜索方法[18]；這些成果給出了確定滑動(dòng)面位置的方法，但是很少在非均質(zhì)邊坡工程中應(yīng)用。

本文擬借助ABAQUS的二次開(kāi)發(fā)平臺(tái)UFIELD程序，將軟件自帶的場(chǎng)變量定義為折減系數(shù)，建立場(chǎng)變量分別與軟件求解時(shí)步和強(qiáng)度參數(shù)之間的關(guān)系，來(lái)實(shí)現(xiàn)基于場(chǎng)變量的有限元強(qiáng)度折減法，此方法僅需一次提交計(jì)算分析便可確定邊坡的安全系數(shù)，較其它方法大大地提高了有限元強(qiáng)度折減法求解邊坡安全系數(shù)的效率；以此為基礎(chǔ)，依據(jù)滑面位置處的等效塑性應(yīng)變最大這一基本原理，利用不同垂直剖面逐點(diǎn)搜索法，基于等效塑性應(yīng)變對(duì)臨界狀態(tài)滑動(dòng)面位置進(jìn)行精細(xì)化尋找，給出了邊坡滑動(dòng)面的準(zhǔn)確位置；最后借助算例驗(yàn)證了本文方法的有效性。

1 基于ABAQUS二次開(kāi)發(fā)平臺(tái)的邊坡有限元強(qiáng)度折減法的基本原理

1.1 邊坡有限元強(qiáng)度折減法的基本原理

邊坡有限元強(qiáng)度折減法的基本原理就是通過(guò)不斷減小材料強(qiáng)度參數(shù)直到邊坡失穩(wěn)，其臨界狀態(tài)對(duì)應(yīng)的折減系數(shù)被定義為安全系數(shù)：

(1)

(2)

式中：SRF為折減系數(shù)；cini和c分別為折減前、后的黏聚力；φini和φ分別為折減前、后的內(nèi)摩擦角。

1.2 基于場(chǎng)變量的有限元強(qiáng)度折減法的基本原理

為了后續(xù)確定邊坡安全系數(shù)的方便，將場(chǎng)變量定義為折減系數(shù):

f=SRF

(3)

上述確定了場(chǎng)變量與折減系數(shù)之間的關(guān)系，于是材料參數(shù)與場(chǎng)變量的關(guān)系如下：

(4)

(5)

定義場(chǎng)變量與當(dāng)前分析步的時(shí)步之間的關(guān)系為：

f=a+bt

(6)

那么，由式(3)可得：

SRF=a+bt

(7)

通常情況下:t取為當(dāng)前分析步的時(shí)步，于是0≤t≤1；a參數(shù)決定著邊坡的初始折減系數(shù)，取為1；a+b參數(shù)決定著折減系數(shù)的范圍，為了計(jì)算的方便，b取為1；這樣將折減系數(shù)的范圍限定在[1,2]的閉區(qū)間內(nèi)，但如若想增大折減系數(shù)的范圍則可以通過(guò)增大b值來(lái)實(shí)現(xiàn)，具體a、b的取值視具體情況而定。

通過(guò)上述材料參數(shù)、場(chǎng)變量、計(jì)算時(shí)步及折減系數(shù)之間的相互關(guān)系，將基于場(chǎng)變量的有限元強(qiáng)度折減法基本原理歸納為：通過(guò)邊坡彈塑性分析中計(jì)算時(shí)步的增加，借助計(jì)算時(shí)步與場(chǎng)變量的關(guān)系，得到場(chǎng)變量，借助場(chǎng)變量與材料參數(shù)之間的關(guān)系來(lái)更新材料參數(shù)，通過(guò)場(chǎng)變量和折減系數(shù)之間的關(guān)系確定當(dāng)前的折減系數(shù)，那么邊坡失穩(wěn)時(shí)刻的折減系數(shù)即邊坡的最終安全系數(shù)。對(duì)于基于場(chǎng)變量的有限元強(qiáng)度折減法，本文采用ABAQUS二次開(kāi)發(fā)平臺(tái)上的UFIELD程序來(lái)實(shí)現(xiàn)。具體程序接口和變量說(shuō)明如下：

SUBROUTINE UFIELD(FIELD,KFIELD, NSECPT,KSTEP,KINC,TIME,NODE,

1 COORDS,TEMP,DTEMP,NFIELD)

C

INCLUDE ‘ABA_PARAM.INC’

C

DIMENSION FIELD(NSECPT,NFIELD), TIME(2), COORDS(3), KSTEP, KINC

1 TEMP(NSECPT), DTEMP(NSECPT)

C

user coding to define FIELD

RETURN

END

其中：FIELD表示場(chǎng)變量；KFIELD表示用戶自定義場(chǎng)變量數(shù)；NFIELD表示更新的場(chǎng)變量數(shù)； KSTEP 表示分析步數(shù)；KINC 表示增量步數(shù)；TIME(1)表示當(dāng)前時(shí)間步的時(shí)間；TIME(2) 表示總時(shí)間步的時(shí)間；NODE 表示節(jié)點(diǎn)數(shù)；COORDS 表示存放節(jié)點(diǎn)坐標(biāo)的數(shù)組；TEMP(NSECPT) 表示所有節(jié)點(diǎn)的當(dāng)前時(shí)刻溫度；DTEMP(NSECPT)表示節(jié)點(diǎn)處的溫度增量。

常規(guī)有限元強(qiáng)度折減法確定邊坡安全系數(shù)的流程見(jiàn)圖1(a)，從流程圖中可以發(fā)現(xiàn)，首先給邊坡賦予初始的抗剪強(qiáng)度參數(shù)，對(duì)邊坡進(jìn)行彈塑性分析，然后判斷邊坡是否失穩(wěn)，如果邊坡失穩(wěn)則減小折減系數(shù)，重復(fù)上一步的計(jì)算；如果邊坡未失穩(wěn)則增大折減系數(shù)，重復(fù)上一步計(jì)算；如此反復(fù)從而求得邊坡的安全系數(shù)。這樣必須提交n次計(jì)算才能獲得邊坡的安全系數(shù)。而基于場(chǎng)變量的有限元強(qiáng)折減法計(jì)算邊坡安全系數(shù)的流程如圖1(b)所示，在ti時(shí)刻各節(jié)點(diǎn)達(dá)到平衡后時(shí)間步向前推進(jìn)，隨著時(shí)間步的增加，強(qiáng)度參數(shù)不斷地折減，直到邊坡的塑性區(qū)貫通和特征點(diǎn)水平位移隨時(shí)間步的變化曲線出現(xiàn)拐點(diǎn)[19]，最終計(jì)算中斷。這樣僅需要一次提交計(jì)算便可以求取邊坡的強(qiáng)度折減安全系數(shù)。

圖1 常規(guī)有限元強(qiáng)度折減法和基于場(chǎng)變量的有限元強(qiáng)度折減法確定邊坡安全系數(shù)的對(duì)比分析流程圖Fig.1 Flow chart for calculating safety factor of slope using traditional strength reduction FEM and strength reduction FEM based on field variable

2 基于最大等效塑性應(yīng)變的邊坡潛在滑動(dòng)面精細(xì)化尋找方法

在基于場(chǎng)變量的邊坡有限元強(qiáng)度折減方法計(jì)算所得等效塑性應(yīng)變的基礎(chǔ)上，依據(jù)滑面位置處的等效塑性應(yīng)變最大這一基本原理，可以實(shí)現(xiàn)一次自動(dòng)搜索出邊坡的潛在滑動(dòng)面位置。依據(jù)最大等效塑性應(yīng)變尋找邊坡潛在滑動(dòng)面的方法如下，首先在水平方向布置一系列的垂直線(見(jiàn)圖2中點(diǎn)劃線)，按照一定的空間步長(zhǎng)分別在每一條垂直線上尋找等效塑性應(yīng)變的最大值，將每條垂直線上等效塑性應(yīng)變最大值點(diǎn)相對(duì)應(yīng)的坐標(biāo)點(diǎn)記錄(見(jiàn)圖2紅色實(shí)心點(diǎn))，并依次將這一系列點(diǎn)連線從而形成邊坡臨界滑動(dòng)面的雛形，接著利用6次多項(xiàng)式采用最小二乘法將臨界滑動(dòng)面進(jìn)行光滑化處理，便可以得到最終的邊坡潛在滑動(dòng)面。特別注意上述在等效塑性應(yīng)變插值過(guò)程中采用等參元形函數(shù)進(jìn)行插值，而且水平向與垂直向布點(diǎn)的個(gè)數(shù)決定著邊坡滑動(dòng)面位置的精確程度。實(shí)際計(jì)算中，可以先擬定出滑動(dòng)面的大體位置，然后對(duì)初擬的滑動(dòng)面位置附近的搜索步長(zhǎng)進(jìn)行加細(xì)，可以起到既節(jié)省時(shí)間又提高精度的效果。

圖2 基于最大等效塑性應(yīng)變的邊坡潛在滑動(dòng)面精細(xì)化尋找方法示意圖Fig.2 Potential slip surface of the slope is determined by the maximum equivalent plastic strain

3 有限元強(qiáng)度折減法二次開(kāi)發(fā)程序的算例驗(yàn)證及滑動(dòng)面精細(xì)化尋找方法的驗(yàn)證

本文基于ABAQUS二次開(kāi)發(fā)平臺(tái)UFIELD程序開(kāi)發(fā)了基于場(chǎng)變量的有限元強(qiáng)度折減法，為了驗(yàn)證該方法的合理性及其基于最大等效塑性應(yīng)變的滑動(dòng)面精細(xì)化尋找方法的可靠性，分別選取均質(zhì)、非均質(zhì)兩個(gè)邊坡作為研究對(duì)象，將極限平衡M-P法計(jì)算結(jié)果作為衡量計(jì)算結(jié)果可靠性的標(biāo)準(zhǔn)。極限平衡M-P法[20]借助Geostudio軟件的SLOPE/W模塊來(lái)實(shí)現(xiàn)。

3.1 均質(zhì)邊坡算例分析

圖3所示均質(zhì)邊坡坡高20m，坡比為1∶1，模型總高度為40m，總寬度為105m；有限元計(jì)算模型共離散為5773個(gè)單元，5938個(gè)節(jié)點(diǎn)；計(jì)算模型兩側(cè)采用水平向約束邊界，模型底部采用全約束邊界[21]。計(jì)算中采用考慮Mohr-Coulomb強(qiáng)度準(zhǔn)則的理想彈塑性模型,具體計(jì)算過(guò)程中所采用的材料參數(shù)見(jiàn)表1。

圖3 計(jì)算網(wǎng)格圖Fig.3 Discrete grid of the calculation model

從計(jì)算結(jié)果表2可見(jiàn)，利用基于ABAQUS二次開(kāi)發(fā)平臺(tái)實(shí)現(xiàn)的基于場(chǎng)變量的有限元強(qiáng)度折減法，采用塑性區(qū)貫通與坡腳點(diǎn)水平位移出現(xiàn)拐點(diǎn)作為有限元強(qiáng)度折減法判別邊坡失穩(wěn)的判據(jù)，綜合確定的邊坡安全系數(shù)為1.210，此安全系數(shù)與傳統(tǒng)的有限元強(qiáng)度折減法計(jì)算結(jié)果完全一致，并且對(duì)應(yīng)安全系數(shù)位于極限平衡M-P法和有限元極限平衡法之間，同時(shí)三者安全系數(shù)之間的絕對(duì)誤差最大值為0.6%，這說(shuō)明本文基于ABAQUS二次開(kāi)發(fā)平臺(tái)實(shí)現(xiàn)的基于場(chǎng)變量的有限元強(qiáng)度折減法對(duì)于均質(zhì)邊坡是有效的。

表2 安全系數(shù)匯總表

從圖4中基于最大等效塑性應(yīng)變搜索得到的滑動(dòng)面位置與邊坡失穩(wěn)時(shí)刻的塑性區(qū)對(duì)比分析圖可以看出，基于最大等效塑性應(yīng)變搜索得到的滑動(dòng)面位置位于等效塑性應(yīng)變等值線密集的滑帶范圍內(nèi)，并且可以清晰地看出塑性應(yīng)變等值線的特征點(diǎn)剛好落在此滑動(dòng)面上。

圖4 基于最大等效塑性應(yīng)變的邊坡潛在滑動(dòng)面確定Fig.4 Potential slip surface of the slope is determined by the maximum equivalent plastic strain

從圖5中不同方案確定的邊坡潛在滑動(dòng)面位置可以看出，離散數(shù)據(jù)點(diǎn)組成的滑動(dòng)面經(jīng)過(guò)光滑處理后與滑動(dòng)面雛形相近。極限平衡M-P法確定的邊坡滑動(dòng)面位置與基于最大等效塑性應(yīng)變尋找方法確定的滑動(dòng)面位置在下半段非常接近，但是在上半段本文方法相對(duì)于極限平衡法要深緩，但整體偏差較小。綜上所述可說(shuō)明基于最大等效塑性應(yīng)變的搜索對(duì)均質(zhì)邊坡臨界滑動(dòng)面的精細(xì)化尋找是有效的。

圖5 不同方案確定的邊坡潛在滑動(dòng)面位置Fig.5 Positions of the potential slip surface of the slope are determined using the different schemes

3.2 非均質(zhì)邊坡算例分析

如圖6所示,非均質(zhì)邊坡取自澳大利亞計(jì)算機(jī)應(yīng)用協(xié)會(huì)(ACADS)的考核題[22]，為了在有限元計(jì)算中減小邊界效應(yīng)對(duì)計(jì)算結(jié)果的影響，將左側(cè)邊界延伸10m，向下延伸10m。具體模型尺寸為：坡高10m，坡比1 ∶2，模型總高度為25m，總寬度為60m；有限元計(jì)算模型共離散為16 556個(gè)單元，16 824個(gè)節(jié)點(diǎn)；模型兩側(cè)采用水平向約束邊界，模型底部采用全約束邊界。計(jì)算采用考慮Mohr-Coulomb強(qiáng)度準(zhǔn)則的理想彈塑性模型,具體計(jì)算過(guò)程中所采用的材料參數(shù)見(jiàn)表3。

圖6 計(jì)算網(wǎng)格圖Fig.6 Discrete grid of the calculation model

材料編號(hào)土體密度/(kg/m3)彈性模量/MPa泊松比內(nèi)摩擦角/()黏聚力/kPa側(cè)壓力系數(shù)①195010025380065②1950100252353065③1950100252072065

從安全系數(shù)計(jì)算結(jié)果匯總表4可以看出，利用本文方法采用塑性區(qū)貫通與坡腳點(diǎn)水平位移出現(xiàn)拐點(diǎn)作為有限元強(qiáng)度折減法判別邊坡失穩(wěn)的判據(jù)，綜合確定的邊坡安全系數(shù)為1.364，此安全系數(shù)位于前人采用傳統(tǒng)的有限元強(qiáng)度折減法確定的安全系數(shù)之間[21,23]，并且多種方法計(jì)算所得安全系數(shù)之間的彼此差別不大，這說(shuō)明本文基于ABAQUS二次開(kāi)發(fā)平臺(tái)實(shí)現(xiàn)的基于場(chǎng)變量的有限元強(qiáng)度折減法在確定非均質(zhì)邊坡對(duì)應(yīng)安全系數(shù)時(shí)也是適用的。

表4 安全系數(shù)匯總表

從圖7中基于最大等效塑性應(yīng)變搜索得到的滑動(dòng)面位置與邊坡失穩(wěn)時(shí)刻的塑性區(qū)對(duì)比分析圖可以看出，基于最大等效塑性應(yīng)變搜索得到的滑動(dòng)面位置位于邊坡失穩(wěn)時(shí)刻的塑性區(qū)內(nèi)部，并且基于最大等效塑性應(yīng)變搜索得到的滑動(dòng)面基本上通過(guò)了等效塑性應(yīng)變等值線的特征點(diǎn)。

圖7 基于最大等效塑性應(yīng)變的邊坡潛在滑動(dòng)面確定Fig.7 Potential slip surface of the slope is determined by the maximum equivalent plastic strain

從圖8中不同方案確定的邊坡潛在滑動(dòng)面位置可以看出，用極限平衡M-P法確定的邊坡滑動(dòng)面位置與基于最大等效塑性應(yīng)變尋找方法確定的滑動(dòng)面位置在滑動(dòng)面下半段幾乎重合，但在滑動(dòng)面上半部分有一定的偏差。從本文方法的離散數(shù)據(jù)點(diǎn)明顯可以看出塑性區(qū)在材料邊界處的方向調(diào)整，然而基于極限平衡M-P法確定的邊坡滑動(dòng)面位置和本文方法擬合曲線均將此信息淹沒(méi)。

圖8 不同方案確定的邊坡潛在滑動(dòng)面位置Fig.8 Positions of the potential slip surface of the slope are determined using the different schemes

但綜合對(duì)比分析，兩種方案的計(jì)算成果，可以發(fā)現(xiàn)兩種方法確定的邊坡潛在滑動(dòng)面位置仍然在計(jì)算可接受的范圍內(nèi)，說(shuō)明基于最大等效塑性應(yīng)變的邊坡滑動(dòng)面精細(xì)化尋找方法對(duì)于非均質(zhì)邊坡是有效的。

4 結(jié) 論

1) 借助ABAQUS軟件的二次開(kāi)發(fā)平臺(tái)的UFIELD程序，將軟件自帶的場(chǎng)變量定義為折減系數(shù)，建立場(chǎng)變量分別與軟件求解時(shí)步和強(qiáng)度參數(shù)之間的關(guān)系，方便地實(shí)現(xiàn)基于場(chǎng)變量的邊坡有限元強(qiáng)度折減法。此方法相對(duì)于其它方法的優(yōu)勢(shì)在于：僅需一次提交計(jì)算分析便可高效地實(shí)現(xiàn)邊坡安全系數(shù)的求解。

2) 在基于場(chǎng)變量的邊坡有限元強(qiáng)度折減方法計(jì)算所得等效塑性應(yīng)變的基礎(chǔ)上，利用不同垂直剖面逐點(diǎn)搜索法，實(shí)現(xiàn)了基于最大等效塑性應(yīng)變的邊坡潛在滑動(dòng)面位置精細(xì)化尋找。

3) 通過(guò)兩個(gè)經(jīng)典算例，展示了此方法和思路進(jìn)行邊坡穩(wěn)定性分析的效果，同時(shí)將此思路確定的潛在滑動(dòng)面位置、安全系數(shù)分別與傳統(tǒng)極限平衡M-P法的對(duì)應(yīng)結(jié)果進(jìn)行對(duì)比分析，發(fā)現(xiàn)不同方法確定的安全系數(shù)相差不大，兩種方法確定的滑動(dòng)面位置下半段相近，上半段本文方法確定滑動(dòng)面相對(duì)于M-P法確定滑動(dòng)面要深緩，但是兩種方法的偏差較小，這說(shuō)明本文所提方法是有效的。

[1]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.

[2]唐芬, 鄭穎人, 趙尚毅. 土坡漸進(jìn)破壞的雙安全系數(shù)討論[J].巖石力學(xué)與工程學(xué)報(bào), 2007, 26(7):1402-1407.

TANG Fen, ZHENG Yingren, ZHAO Shangyi. Discussion on two safety factors for progressive failure of soil slope[J].Chinese Journal of Rock Mechanics and Engineering, 2007, 26(7):1402-1407.

[3]YUAN Wei, BAI Bing, LI Xiaochun, et al. A strength reduction method based on double reduction parameters and its application[J].Journal of Central South University, 2013, 20(9):2555-2562.

[4]薛海斌, 黨發(fā)寧, 尹小濤, 等. 邊坡強(qiáng)度參數(shù)非等比例相關(guān)聯(lián)折減法研究[J].巖石力學(xué)與工程學(xué)報(bào), 2015, 34(S2):4005-4012.

XUE Haibin, DANG Faning, YIN Xiaotao, et al. Research on method of slope strength parameters non-proportional associated reduction[J].Chinese Journal of Rock Mechanics and Engineering, 2015, 34(S2): 4005-4012.

[5]楊光華, 鐘志輝, 張玉成, 等. 用局部強(qiáng)度折減法進(jìn)行邊坡穩(wěn)定性分析[J].巖土力學(xué), 2010, 31(S2):53-58.

YANG Guanghua, ZHONG Zhihui, ZHANGYucheng, et al. Slope stability analysis by local strength reduction method[J].Rock and Soil Mechanics, 2010, 31(S2):53-58.

[6]薛雷, 孫強(qiáng), 秦四清, 等. 非均質(zhì)邊坡強(qiáng)度折減法折減范圍研究[J].巖土工程學(xué)報(bào), 2011, 33(2):275-280.

XUE Lei, SUN Qiang, QIN Siqing, et al. Scope of strength reduction for inhomogeneous slopes[J].Chinese Journal of Geotechnical Engineering, 2011, 33(2):275-280.

[7]陳國(guó)慶, 黃潤(rùn)秋, 周輝, 等. 邊坡漸進(jìn)破壞的動(dòng)態(tài)強(qiáng)度折減法研究[J].巖土力學(xué), 2013, 34(4):1140-1146.

CHEN Guoqing, HUANG Runqiu, ZHOU Hui, et al. Research on progressive failure for slope using dynamic strength reduction method[J].Rock and Soil Mechanics, 2013, 34(4):1140-1146.

[8]張魯渝, 鄭穎人, 趙尚毅, 等. 有限元強(qiáng)度折減系數(shù)法計(jì)算土坡穩(wěn)定安全系數(shù)的精度研究[J].水利學(xué)報(bào), 2003, (1):21-27.

ZHANG Luyu, ZHENG Yingren, ZHAO Shangyi, et al. The feasibility study of strength reduction method with FEM for calculating safety factors of soil slope stability[J]. Journal of Hydraulic Engineering, 2003, (1):21-27.

[9]鄭宏, 李春光, 李焯芬, 等. 求解安全系數(shù)的有限元法[J]. 巖土工程學(xué)報(bào), 2002, 24(5):626-628.

ZHENG Hong, LI Chunguang, LEE C F, et al. Finite element method for solving the factor of safety[J].Chinese Journal of Geotechnical Engineering, 2002, 24(5):626-628.

[10]張培文, 陳祖煜. 剪脹角對(duì)求解邊坡穩(wěn)定的安全系數(shù)的影響[J].巖土力學(xué), 2004, 25(11):1757-1760.

ZHANG Peiwen，CHEN Zuyu. Finite element method for solving safety factor of slope stability[J].Rock and Soil Mechanics, 2004, 25(11):1757-1760.

[11]楊光華, 張玉成, 張有祥. 變模量彈塑性強(qiáng)度折減法及其在邊坡穩(wěn)定分析中的應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報(bào), 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.

[12]宋二祥. 土工結(jié)構(gòu)安全系數(shù)的有限元計(jì)算[J].巖土工程學(xué)報(bào), 1997, 19(2):1-7.

SONG Erxiang. Finite element analysis of safety factor for soil structures[J].Chinese Journal of Geotechnical Engineering, 1997, 19(2):1-7.

[13]GRIFFITHS D V, LANE P A. Slope stability analysis by finite element[J].Geotechnique, 1999, 49(3):387-403.

[14]連鎮(zhèn)營(yíng), 韓國(guó)城, 孔憲京. 強(qiáng)度折減有限元法研究開(kāi)挖邊坡的穩(wěn)定性[J].巖土工程學(xué)報(bào), 2001, 23(4):407-411.

LIAN Zhenying，HAN Guocheng，KONG Xianjing. Stability analysis of excavation by strength reduction FEM [J].Chinese Journal of Geotechnical Engineering, 2001, 23(4):407-411.

[15]曹先鋒, 徐千軍. 邊坡穩(wěn)定分析的溫控參數(shù)折減有限元法[J]. 巖土工程學(xué)報(bào), 2006, 28(11):2039-2042.

CAO Xianfeng, XU Qianjun. Temperature driving strength reduction method for slope stability analysis [J]. Chinese Journal of Geotechnical Engineering, 2006, 28(11):2039-2042.

[16]李寧, 許建聰. 基于場(chǎng)變量的邊坡穩(wěn)定分析有限元強(qiáng)度折減法[J]. 巖土力學(xué), 2012, 33(1):314-318.

LI Ning, XU Jiancong. Strength reduction FEM for slope stability analysis based on field variable [J]. Rock and Soil Mechanics, 2012, 33(1):314-318.

[17]孫冠華, 鄭宏, 李春光. 基于等效塑性應(yīng)變的邊坡滑面搜索[J]. 巖土力學(xué), 2008, 29(5):1159-1163.

SUN Guanhua, ZHENG Hong, LI Chunguang. Searching critical slip surface of slopes based on equivalent plastic strain [J]. Rock and Soil Mechanics, 2008, 29(5):1159-1163.

[18]李劍, 陳善雄, 余飛. 基于最大剪應(yīng)變?cè)隽康倪吰聺撛诨瑒?dòng)面搜索[J]. 巖土力學(xué), 2013, 34(supp.1):371-378.

LI Jian, CHEN Shanxiong, YU Fei. A method for searching potential failure surface of slope based on maximum shear strain increment [J]. Rock and Soil Mechanics, 2013, 34(supp.1):371-378.

[19]張媛, 李榮建, 李錦, 等. 超載作用下土質(zhì)邊坡失穩(wěn)與破壞模式分析[J]. 西安理工大學(xué)學(xué)報(bào), 2015, 31(3): 347-352.

ZHANG Yuan, LI Rongjian, LI Jin, et al. An analysis of instability and failure mode of soil slope under overloading action [J]. Journal of Xi’an University of Technology, 2015, 31(3): 347-352.

[20]MORGENSTERN N R, PRICE V E. The analysis of the stability of general slip surfaces [J]. Geotechnique, 1965, 15(1):79-93.

[21]邵龍?zhí)? 李紅軍. 土工結(jié)構(gòu)穩(wěn)定分析—有限元極限平衡法及其應(yīng)用[M]. 北京: 科學(xué)出版社, 2011.

[22]陳祖煜. 土質(zhì)邊坡穩(wěn)定性分析—原理.方法.程序[M]. 北京: 中國(guó)水利水電出版社, 2003.

[23]楊有成. 強(qiáng)度折減法在斜坡穩(wěn)定性分析中的適用性分析[D]. 武漢: 中國(guó)地質(zhì)大學(xué), 2008.

YANG Youcheng. Research on applicability of shear strength reduction method in slope stability analysis [D]. Wuhan: China University of Geosciences, 2008.

(責(zé)任編輯 楊小麗)

Research on the strength reduction FEM of the slope by applying the secondary development platform of ABAQUS

HOU Ling1， XUE Haibin2， ZHOU Zehua2， LIU Haiwei2， PAN Feng2

(1.School of Sciences，Xi’an University of Technology，Xi’an 710054，China；2.Institute of Geotechnical Engineering，Xi’an University of Technology，Xi’an 710048，China)

The strength reduction FEM plays an important role in the slope stability analysis. But the strength reduction FEM often gets the safety factor of slope by the dichotomy, which is less efficient because of its need for trial. With the secondary development platform UFIELD program of ABAQUS, the strength reduction FEM based on the field variable is realized by defining the field variable of the software as the reduction factor and establishing the relationship between the time step and strength parameters and field variable. This method only submits an analysis to calculate the safety factor of the slope effectively. Based on the principle that the equivalent plastic strain of the sliding zone is maximum, the searching for the location of the potential slip surface of the slope is realized using point by point searching method of different vertical profiles. Finally, taking the two typical slope as examples, the comparative analysis of the potential slip surface and safety factor is made by the traditional limit equilibrium method of M-P and the method proposed in this paper. Thus, the validity of the method proposed in this paper is verified.

strength reduction FEM of the slope；field variable；maximum equivalent plastic strain；potential slip surface

10.19322/j.cnki.issn.1006-4710.2016.04.013

2016-02-16

國(guó)家自然科學(xué)基金資助項(xiàng)目(51679199)；水利部公益性行業(yè)科研專項(xiàng)基金資助項(xiàng)目(201501034-04);陜西省科技統(tǒng)籌創(chuàng)新工程重點(diǎn)實(shí)驗(yàn)室資助項(xiàng)目(2014SZS15-Z01)

候玲，女，副教授，研究方向?yàn)橛?jì)算力學(xué)。E-mail：houling@xaut.edu.cn

TU43

A

1006-4710(2016)04-0449-06