邊坡失穩破壞的顆粒流模擬

萬 江， 溫 泉， 昝世明

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

邊坡失穩破壞的顆粒流模擬

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

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

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

邊坡失穩； 顆粒流( PFC)； 數值模擬； 平行粘結模型

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

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

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

1 PFC模型

PFC2D程序是基于離散單元法來模擬圓形顆粒介質的運動及其相互作用的。其將巖土體看作由單粒、集粒或凝塊等骨架單元共同構成的結構體系。PFC 顆粒流離散元是位移分析法，顆粒流理論在整個計算循環過程中，交替應用力-位移定律和牛頓運動定律，通過力-位移定律更新接觸部分的接觸力，通過運動定律，更新顆粒與墻邊界的位置，構成顆粒之間的新接觸。

在PFC2D中，材料的本構特性是通過接觸本構模型來模擬的。顆粒的接觸本構模型有：(1)接觸剛度模型;(2)滑動模型;(3)粘結模型。接觸剛度模型提供了接觸力和相對位移的彈性關系，滑動模型則強調切向和法向接觸力使得接觸顆粒可以發生相對移動，而粘結模型是限制總的切向和法向力使得在粘結強度范圍內發生接觸。PFC2D提供了 2種粘結模型，即接觸粘結模型(Contact-bond model)和平行粘結模型( Parallel-bond model) 。接觸粘結認為粘結只發生在接觸點很小范圍內，而平行粘結發生在接觸顆粒間圓形或方形有限范圍內。接觸粘結只能傳遞力，而平行粘結同時能傳遞力矩。

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

2 算例與分析

2.1 典型算例

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

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

2.2 計算參數選取

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

表1 土的宏觀材料

表2 土的微觀參數

2.3 土質邊坡計算過程與結果

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

圖2 邊坡模型

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

圖3、圖4為顆粒流程序對上述算例進行模擬后的最終結果的位移等值線，其中圖3是由邊坡未失穩時顆粒的位置與失穩后的位移繪制成的，圖4是邊坡失穩后顆粒的位置與位移繪制成的。滑動面根據該算例的尺寸，所建立的模型中約有15 000，平均半徑為0.1 m，位移小于平均半徑的顆粒只發生了形變，因此以位移等于0.1 m為準則確定滑動面。圖3、圖4的中光滑曲面為算例中的滑動面，兩個滑動面大部分對應的很好。由于建立模型時顆粒半徑不夠接近實際土顆粒以及整體模型不夠精細，坡頂位置的滑動面對應的不太好。對比兩張位移等值線圖可以發現，滑動面上的球的位置基本不變，也可以證明該數值模擬的正確性。

圖3 位移等值線

圖4 位移等值線

邊坡失穩時平行粘結模型破壞過程如圖5所示。但是由于有滑動模型與接觸剛度模型依然存在，粘結強度只是減小而不是消失，所以粘結力鏈斷裂范圍大于失穩范圍。

圖5 粘結力鏈破壞過程

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

圖6 邊坡失穩過程

圖7 速度、位移矢量

3 結 論

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

(2)對于黏性較小的土質邊坡，變形邊坡自身的彈性變形、塑性變形或流變變形較大、流變趨勢持續時間過長，整個坡體表現塑性流動狀態，沒有明顯的裂縫出現。

(3)為提高模擬準確性，邊坡表平面顆粒半徑應當由大到小進行逐漸細化，邊坡模型卸載、安全系數的折減應當逐步進行。

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

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

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

[4] 周健，王家全，曾遠，等. 顆粒流強度折減法和重力增加法的邊坡安全系數研究[J]. 巖土力學，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] 焦玉勇，葛修潤. 基于靜態松弛法求解的三維離散單元法[J]. 巖石力學與工程學報，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] 賀續文，劉忠，廖彪，等. 基于離散元法的節理巖體邊坡穩定性分析[J].巖土力學，2011,32(7):2199-2204.

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

十九冶成都建設有限公司科研開發項目(SJY-44)

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

TU413.6+2

A

[定稿日期]2016-09-28