邊坡失穩(wěn)破壞的顆粒流模擬

萬 江， 溫 泉， 昝世明

(1.十九冶成都建設(shè)有限公司，四川成都 610091； 2.東北大學(xué)資源與土木工程學(xué)院，遼寧沈陽 110819)

邊坡失穩(wěn)破壞的顆粒流模擬

萬 江1， 溫 泉1， 昝世明2

(1.十九冶成都建設(shè)有限公司，四川成都 610091； 2.東北大學(xué)資源與土木工程學(xué)院，遼寧沈陽 110819)

文章采用顆粒流離散元數(shù)值方法對典型邊坡算例失穩(wěn)變形破壞過程進(jìn)行了計(jì)算模擬，模擬中利用數(shù)值試驗(yàn)方法建立起宏觀與微觀力學(xué)參數(shù)關(guān)系，進(jìn)而確定計(jì)算用微觀參數(shù)。計(jì)算結(jié)果與典型考核算例基本吻合，從而說明了計(jì)算的有效性。同時(shí)計(jì)算表明利用PFC方法可以自動(dòng)搜索確定失穩(wěn)滑移面，表現(xiàn)破壞發(fā)展過程，在破壞過程中粘結(jié)力鏈不斷斷裂，但由于滑動(dòng)摩擦和顆粒間接觸作用仍能維持部分粘結(jié)力破壞區(qū)域巖土體的穩(wěn)定，因而粘結(jié)力斷裂區(qū)域范圍大于失穩(wěn)范圍，從而提高了對邊坡失穩(wěn)的認(rèn)識。

邊坡失穩(wěn)； 顆粒流( PFC)； 數(shù)值模擬； 平行粘結(jié)模型

邊坡廣泛存在于自然條件和生產(chǎn)條件下，將邊坡在遭遇地震、暴雨等觸發(fā)條件下失穩(wěn)將會產(chǎn)生滑坡，其會對人們的生命和財(cái)產(chǎn)安全帶來嚴(yán)重威脅。如“5·12”汶川地震誘發(fā)了大量的滑坡，給災(zāi)區(qū)帶來了重大的災(zāi)難[1]； 1979年12月,四川攀鋼蘭尖鐵礦滑坡摧毀運(yùn)輸主平硐,礦山停產(chǎn)半年,造成重大經(jīng)濟(jì)損；1975年6月云南省的一個(gè)錫礦滑坡,直接經(jīng)濟(jì)損失1 000萬元以上[2]。

邊坡穩(wěn)定計(jì)算問題是土力學(xué)研究中很重要的一個(gè)課題[3-4]。目前，邊坡穩(wěn)定性分析方法很多，其中傳統(tǒng)的方法有工程類比法、極限平衡分析法等[5-6]。其中，工程類比法需要具有豐富的實(shí)踐經(jīng)驗(yàn)，極限平衡分析法模型簡單、計(jì)算簡捷，被廣泛的應(yīng)用于工程實(shí)踐中，但該類方法并不能在復(fù)雜條件下深入探討滑坡的破壞過程和機(jī)理。而隨著計(jì)算技術(shù)的快速發(fā)展，以有限元法、無界元法、離散元法為代表的現(xiàn)代數(shù)值計(jì)算方法在邊坡穩(wěn)定性計(jì)算中得到了極大地發(fā)展。有限元法在邊坡穩(wěn)定性分析中應(yīng)用較多，其考慮了巖體的不連續(xù)性和非均質(zhì)性，能進(jìn)行線性分析和非線性分析，但是有限元法在解決大變形、應(yīng)力集中等問題時(shí)效果不如人意。1977 年，P Bettess提出了無界元法，是有限元方法的推廣，該法有效地解決了有限元方法的人為確定邊界的缺點(diǎn)，常常與有限元法聯(lián)合使用，但是不適合模擬非線性、非均質(zhì)邊坡穩(wěn)定性問題方面。1970 年由P A Cundall 提出并應(yīng)用于巖土體穩(wěn)定性分析的離散元法，該法可利用顯示時(shí)間差求解動(dòng)力平衡方程，適合于模擬非線性大位移和非連續(xù)介質(zhì)大變形問題。目前，國外內(nèi)學(xué)者利用離散元進(jìn)行了大量的研究[7-10]，離散元法得到很好的發(fā)展。其中塊體離散元計(jì)算程序適用于巖質(zhì)邊坡，而理想顆粒離散元方法則具有更廣的適用性，能夠模擬巖土體破壞演化過程，進(jìn)而研究其破壞的細(xì)觀機(jī)制。如學(xué)者 Bardet 和 Proubet[11-12]應(yīng)用理想二維顆粒集合體模擬了粒狀材料中剪切帶的結(jié)構(gòu)，對剪切帶的厚度、帶內(nèi)位移、孔隙比、體應(yīng)變及顆粒旋轉(zhuǎn)進(jìn)行了研究。周健等采用顆粒流法模擬了土的工程力學(xué)性質(zhì)，與室內(nèi)試驗(yàn)結(jié)果吻合較好，并用于土坡穩(wěn)定性分析中[13]。

本文將采用二維顆粒流程序( Par-ticle Follow Code PFC2D) 數(shù)值模擬軟件來模擬土體邊坡變形失穩(wěn)過程，并就邊坡失穩(wěn)破壞過程中的若干問題進(jìn)行探討。

1 PFC模型

PFC2D程序是基于離散單元法來模擬圓形顆粒介質(zhì)的運(yùn)動(dòng)及其相互作用的。其將巖土體看作由單粒、集?；蚰龎K等骨架單元共同構(gòu)成的結(jié)構(gòu)體系。PFC 顆粒流離散元是位移分析法，顆粒流理論在整個(gè)計(jì)算循環(huán)過程中，交替應(yīng)用力-位移定律和牛頓運(yùn)動(dòng)定律，通過力-位移定律更新接觸部分的接觸力，通過運(yùn)動(dòng)定律，更新顆粒與墻邊界的位置，構(gòu)成顆粒之間的新接觸。

在PFC2D中，材料的本構(gòu)特性是通過接觸本構(gòu)模型來模擬的。顆粒的接觸本構(gòu)模型有：(1)接觸剛度模型;(2)滑動(dòng)模型;(3)粘結(jié)模型。接觸剛度模型提供了接觸力和相對位移的彈性關(guān)系，滑動(dòng)模型則強(qiáng)調(diào)切向和法向接觸力使得接觸顆?？梢园l(fā)生相對移動(dòng)，而粘結(jié)模型是限制總的切向和法向力使得在粘結(jié)強(qiáng)度范圍內(nèi)發(fā)生接觸。PFC2D提供了 2種粘結(jié)模型，即接觸粘結(jié)模型(Contact-bond model)和平行粘結(jié)模型( Parallel-bond model) 。接觸粘結(jié)認(rèn)為粘結(jié)只發(fā)生在接觸點(diǎn)很小范圍內(nèi)，而平行粘結(jié)發(fā)生在接觸顆粒間圓形或方形有限范圍內(nèi)。接觸粘結(jié)只能傳遞力，而平行粘結(jié)同時(shí)能傳遞力矩。

我們選擇平行粘結(jié)模型來模擬土，離散的顆粒僅在接觸部分受力，當(dāng)接觸點(diǎn)受到的作用力大于接觸強(qiáng)度時(shí)，這些顆粒相互分離，此時(shí)平行粘結(jié)模型消失，接觸剛度模型和滑動(dòng)模型持續(xù)生效，顆粒發(fā)生變形或位移[7-9]。平行粘結(jié)模型的模型細(xì)觀參數(shù)有法向粘結(jié)強(qiáng)度、切向粘結(jié)強(qiáng)度、法向粘結(jié)剛度、切向粘結(jié)剛度及粘結(jié)半徑。

2 算例與分析

2.1 典型算例

計(jì)算采用澳大利亞計(jì)算機(jī)應(yīng)用協(xié)會(ACADS)調(diào)查所得的一道邊坡穩(wěn)定分析程序的考題。該考題于1989年4月在澳大利亞巖土工程協(xié)會會議上公布， “裁判程序”(Referee Program)答案邀請了多位國際著名專家參與調(diào)查，所得成果比較可靠，并在近年來較多應(yīng)用于裁定邊坡穩(wěn)定分析程序的有效性?？碱}尺寸及最危險(xiǎn)滑動(dòng)面如圖1所示。

圖1 考核算例邊坡形狀及最危險(xiǎn)滑動(dòng)面位置(單位:m)

2.2 計(jì)算參數(shù)選取

邊坡宏觀物理力學(xué)參數(shù)主要有楊氏模量、泊松比、黏聚力和內(nèi)摩擦角等，本算例的宏觀力學(xué)參數(shù)如表1所示。而程序需要輸入微觀力學(xué)參數(shù)，為了得到能夠表現(xiàn)宏觀力學(xué)性質(zhì)的微觀力學(xué)參數(shù)，需要進(jìn)行雙軸數(shù)值試驗(yàn)，通過一系列雙軸試驗(yàn)得到PFC模擬與巖土宏觀特征相應(yīng)的微觀力學(xué)參數(shù)。通過數(shù)值試驗(yàn)得到PFC數(shù)值模擬所用微觀力學(xué)參數(shù)如表2所示，其對應(yīng)的宏觀力學(xué)參數(shù)如表1所示。從表1中可以看到表2中微觀力學(xué)參數(shù)所表現(xiàn)的宏觀參數(shù)與算例基本相符。

表1 土的宏觀材料

表2 土的微觀參數(shù)

2.3 土質(zhì)邊坡計(jì)算過程與結(jié)果

按照表1、表2中列出的宏細(xì)觀參數(shù)建立黏性土坡模型，首先用wall命令建立算例中邊坡的輪廓，并利用no_shadow命令，通過控制空隙率在輪廓內(nèi)生成致密顆粒材料；然后設(shè)置重力、密度、摩擦系數(shù)并計(jì)算至平衡；第三步賦予材料其他的微觀參數(shù)使其具有形成具有粘結(jié)強(qiáng)度的土，圖2為邊坡的顆粒模型；第四步刪掉頂部束縛顆粒的3面墻，相當(dāng)于對邊坡進(jìn)行了開挖或卸載，所以需要計(jì)算達(dá)到新的平衡；最后將顆粒的速度、位移清零并計(jì)算至滑坡過程完畢，并記錄該狀態(tài)的顆粒的位移及變形情況。

圖2 邊坡模型

離散元是位移分析法，用臨界位移來確定破壞標(biāo)準(zhǔn)是比較合適的方法，但是PFC模擬常常以顆粒最大不平衡力的發(fā)展情況作為停止計(jì)算依據(jù)。綜合考慮，在有明顯位移時(shí)，以累計(jì)位移導(dǎo)致邊坡失穩(wěn)作為破壞標(biāo)準(zhǔn)；在無明顯位移時(shí)，以顆粒的平均不平衡力小于10-1N，并且最大不平衡力與平均不平衡力之比小于10作為計(jì)算結(jié)束的標(biāo)準(zhǔn)[4]。

圖3、圖4為顆粒流程序?qū)ι鲜鏊憷M(jìn)行模擬后的最終結(jié)果的位移等值線，其中圖3是由邊坡未失穩(wěn)時(shí)顆粒的位置與失穩(wěn)后的位移繪制成的，圖4是邊坡失穩(wěn)后顆粒的位置與位移繪制成的。滑動(dòng)面根據(jù)該算例的尺寸，所建立的模型中約有15 000，平均半徑為0.1 m，位移小于平均半徑的顆粒只發(fā)生了形變，因此以位移等于0.1 m為準(zhǔn)則確定滑動(dòng)面。圖3、圖4的中光滑曲面為算例中的滑動(dòng)面，兩個(gè)滑動(dòng)面大部分對應(yīng)的很好。由于建立模型時(shí)顆粒半徑不夠接近實(shí)際土顆粒以及整體模型不夠精細(xì)，坡頂位置的滑動(dòng)面對應(yīng)的不太好。對比兩張位移等值線圖可以發(fā)現(xiàn)，滑動(dòng)面上的球的位置基本不變，也可以證明該數(shù)值模擬的正確性。

圖3 位移等值線

圖4 位移等值線

邊坡失穩(wěn)時(shí)平行粘結(jié)模型破壞過程如圖5所示。但是由于有滑動(dòng)模型與接觸剛度模型依然存在，粘結(jié)強(qiáng)度只是減小而不是消失，所以粘結(jié)力鏈斷裂范圍大于失穩(wěn)范圍。

圖5 粘結(jié)力鏈破壞過程

通過觀察該土質(zhì)邊坡失穩(wěn)過程(圖6)，我們發(fā)現(xiàn)由于黏性較小，邊坡自身的彈性變形、塑性變形或流變變形較大、流變趨勢持續(xù)時(shí)間很長，整個(gè)坡體表現(xiàn)塑性流動(dòng)狀態(tài)，可發(fā)現(xiàn)無明顯的滑動(dòng)面和裂縫。但是，通過對速度和、位移的檢測，可以發(fā)現(xiàn)土質(zhì)邊坡的大部分是受剪破壞(圖7)。

圖6 邊坡失穩(wěn)過程

圖7 速度、位移矢量

3 結(jié) 論

(1)通過顆粒元平行粘結(jié)模型，采用數(shù)值試驗(yàn)的手段，可以很好地模擬了土體滑坡的漸進(jìn)性破壞過程，并且不用設(shè)定滑動(dòng)面的位置與形狀；可以很直觀地模擬邊坡失穩(wěn)的發(fā)展和運(yùn)動(dòng)以及失穩(wěn)面，對失穩(wěn)開始的位置的理解、失穩(wěn)的發(fā)展和最后的失穩(wěn)形態(tài)有很大的實(shí)際意義，如土坡加固等。

(2)對于黏性較小的土質(zhì)邊坡，變形邊坡自身的彈性變形、塑性變形或流變變形較大、流變趨勢持續(xù)時(shí)間過長，整個(gè)坡體表現(xiàn)塑性流動(dòng)狀態(tài)，沒有明顯的裂縫出現(xiàn)。

(3)為提高模擬準(zhǔn)確性，邊坡表平面顆粒半徑應(yīng)當(dāng)由大到小進(jìn)行逐漸細(xì)化，邊坡模型卸載、安全系數(shù)的折減應(yīng)當(dāng)逐步進(jìn)行。

[1] 李秀珍，孔紀(jì)名，鄧紅艷，等. “5·12”汶川地震滑坡特征及失穩(wěn)破壞模式分析[J]. 四川大學(xué)學(xué)報(bào)：工程科學(xué)版，2009，41(3):72-77.

[2] 阮善菊，張福生. 露天礦山的生態(tài)問題及解決對策探討[J]. 現(xiàn)代礦業(yè)，2010，26(8):81-84.

[3] 周健，池永. 土的工程力學(xué)性質(zhì)的顆粒流模擬[J]. 固體力學(xué)學(xué)報(bào), 2004, 25(4):377-382.

[4] 周健，王家全，曾遠(yuǎn)，等. 顆粒流強(qiáng)度折減法和重力增加法的邊坡安全系數(shù)研究[J]. 巖土力學(xué)，2009，30(6):1549-1554.

[5] SARMA S K. Stability analysis of embankments and slopes[J]. Journal of Geotechnical Engineering Division, ASCE, 1979, 105(12): 1511-1524.

[6] BISHOP A W. The use of the slip circle in the stability analysis of slopes[J]. Geotechnique, 1955, 5(1): 7-17.

[7] CUNDALL P E, HART R G. Numerical modeling of discontinua[J]. Engineering Computations, 1992, 9(2):101-113.

[8] Itasca Consulting Group. PFC2D user’s manual (version3.1)[M]. Minneapolis, Minnesota: Itasca Consulting Group, Inc, 2004.

[9] Itasca Consulting Group. PFC2d theory and back-ground[M]. Minnesota, Minneapolis: Itasca Consulting Group, 2004.

[10] 焦玉勇，葛修潤. 基于靜態(tài)松弛法求解的三維離散單元法[J]. 巖石力學(xué)與工程學(xué)報(bào)，2000，19(4): 451-456.

[11] BARDET J P, PROUBET J. Shear-band analysis in idealized granular material[J]. Journal of Engineering Mechanics, 1992, 118(8): 397-415.

[12] 賀續(xù)文，劉忠，廖彪，等. 基于離散元法的節(jié)理巖體邊坡穩(wěn)定性分析[J].巖土力學(xué)，2011,32(7):2199-2204.

[13] PFC2D (Particle Flow Code in 2 Dimensions), Version 1.1. Itasca Consulting Group, Inc. 1999a, Minneapolis, MN: ICG.

十九冶成都建設(shè)有限公司科研開發(fā)項(xiàng)目(SJY-44)

萬江(1988～)，男，工學(xué)碩士，助理工程師，主要從事邊坡穩(wěn)定性研究工作。

TU413.6+2

A

[定稿日期]2016-09-28