邊坡穩定分析的顆粒流方法研究

張小雪，王濱生，遲玉鵬，呂志強

(1.哈爾濱工程大學航天與建筑工程學院，黑龍江哈爾濱150001;2.黑龍江省六建建筑工程有限責任公司，黑龍江哈爾濱150076)

滑坡是一種較頻發的自然災害，滑坡運動問題屬于動力學問題，目前，研究人員對滑坡進行穩定性分析時大多采用極限平衡法［1］，這種方法沒有考慮土體中顆粒之間的相對運動及變形，只考慮靜力平衡條件和土的摩爾－庫侖準則，從而可知，這種方法不能很準確地表述滑坡運動的特點［2］。數值分析方法考慮了土體的應力－應變關系，其中有限單元法基于連續性假設，是數值模擬在邊坡穩定分析中應用最早的，但有限元法不能求解滑坡的大變形和大位移。離散元法其突出功能是在反應土顆粒之間接觸面的滑移、分離及側翻等大位移的同時，又能計算土顆粒內部的變形與應力分布。顆粒流(particle flow code，PFC)方法是離散元法中一種特殊的方法，在顆粒流中，把單元都看作是剛性的圓形或球形的顆粒，這種單元構成較符合散體介質的特性，所以對散體介質用顆粒流方法分析較為準確［3-4］。進行顆粒流分析時，在土樣應力－應變關系曲線已知的情況下，一個關鍵的問題是如何獲取能反映材料宏觀特性的細觀參數，目前應用最多的是通過雙軸壓縮試驗不斷調整細觀參數與宏觀特性，最終匹配得到合理的細觀參數［5］。近年來，徐小敏等［6-8］對粘性材料的抗剪強度參數與顆粒流細觀參數聯系起來，真實反映粘性材料剪切特性，為正確有效地選取顆粒細觀參數提供了依據。

1 接觸本構關系

在PFC2D(particle flow code in 2 dimensions)中，材料的本構特性是通過接觸本構模型來模擬的。接觸模型有3種，即接觸－剛度模型、滑動和分離模型和粘結模型。其中接觸－剛度模型有線性模型和Hertz-Mindlin模型2種。線性接觸模型通過法向和切向剛度定義，用假定接觸實體的接觸剛度是串聯的方法來計算聯合剛度。Hertz-Mindlin模型是非線性剛度模型，適用于模擬無粘結、應變小的情況。滑動和分離模型是2個接觸實體的內在特性，采用限制剪切力的方式，無法向抗拉強度，并允許滑動。粘結模型有2種，即接觸粘結模型和平行粘結模型。接觸粘結是球與球之間通過一點發生接觸，對于兩個相互接觸的顆粒，無法控制其相對轉動。所以，接觸粘結只可以傳遞力;平行粘結模擬顆粒接觸后還有其他粘結作用的情況，該作用粘結的有效剛度與顆粒接觸點剛度采用并聯連接，平行粘結可以傳遞力和力矩來阻止顆粒轉動［9-10］。接觸粘結和平行粘結可以同時設置，這樣可以使顆粒間有最小的相對運動。

2 數值試樣試驗

一般工程計算中，材料的特性是通過宏觀參數來描述的，這些宏觀參數往往是由實驗得到，或由實驗參數派生的，而對離散元理論，需要輸入的是土體的細觀參數，由細觀參數來反映宏觀參數，為能正確反映宏觀參數，首先必須建立正確的宏－細觀參數的關系。采用PFC2D程序進行雙軸試驗的顆粒流模擬，構造四周邊界約束的長方形試樣，頂部和底部約束模擬加載，左右兩邊的約束模擬試樣邊界施加的限制。試驗中，顆粒流程序自動記錄墻體和顆粒的運動過程和位置及顆粒之間的相互作用情況，得到顆粒和墻體運動的宏觀數據［11］。

2.1 試樣顆粒生成

先創建幾何尺寸為h×b=70 mm×35 mm的矩形試樣，該試樣由4道沒有摩擦的墻圍成(見圖1)，在PFC模型中，主要的控制參數見表1。

圖1 球面沖擊波作用于剛性平板示意圖Fig.1 The spherical shock wave acting on a rigid plate

為了使PFC模型更好的模擬土試樣的不均勻特性，顆粒單元半徑R采用從最小值Rmin到最大值Rmax的正態分布，取Rmin=0.25 mm，Rmax=0.5 mm。模型生成后，利用兩側向、頂部、底部墻體的運動，對試樣施加等向固結圍壓σ0，即σx=σy=σ0。σx和σy按下式計算:

式中:fx為左/右側墻體與顆粒單元在X方向接觸力;fy為頂/底部墻體與顆粒單元在Y方向接觸力;r為所有顆粒單元的平均半徑:

用下列關系計算X、Y方向上的應變［15］:

表1 粘性土PFC模型的輸入參數Table 1 PFC model input parameters of clay

2.2 土體宏、細觀力學參數關系分析

圖2為在圍壓σ3=100 kPa下，顆粒摩擦系數μ分別取0.5、0.9時，不同顆粒接觸強度C的土樣的應力－應變曲線。

圖2 圍壓100 kPa下的應力－應變曲線Fig.2 Stress-strain curves with 100 kPa confining

從圖2(a)可以看出在施加相同圍壓的情況下，具有相同摩擦系數，不同顆粒接觸強度的試樣，其峰值強度隨接觸強度的增大而增大，可見接觸粘結強度對峰值強度有較大的影響，說明粘結強度是影響粘土的主要因素之一，同時，隨著粘結強度的增大，曲線表現出較大的軟化特性。比較圖2(a)和(b)可以看出:具有相同顆粒接觸強度的試樣，隨著摩擦系數的增大峰值強度增大，可見摩擦系數主要影響土試樣的峰值強度。

圖3 試樣雙軸壓縮力學響應的應力－應變曲線Fig.3 Stress-strain behavior of biaxial testing

圖 3 為具有相同摩擦系數 μ=0.5、0.9，接觸強度C分別為50、100、200 kN，在不同的圍壓下得到的土試樣的應力－應變曲線。

圖4 摩爾－庫侖破壞包絡線Fig.4 Failure envelopes of Mohr-Coulomb

比較圖3中在不同摩擦系數下，其他參數相同的應力－應變曲線，可以明顯的看到，各條試驗曲線的切線彈性模量是一樣的，在該區域內，粘土的強度是由接觸粘結決定的，不論摩擦系數多少，都不會對粘土強度有影響，當接觸粘結破壞后，峰值強度隨著摩擦系數的增大而增大。

根據圖3的加載結果，得到了峰值強度，以此作為繪制摩爾應力圓的基本參數。圖4為顆粒摩擦系數相同，接觸強度不同的兩組摩爾－庫侖包絡線，從圖可以看出，每組圖的包絡線接近平衡，也就是說，試樣的摩擦角不變，這說明顆粒摩擦系數是影響試樣摩擦角的唯一因素;在對比兩組圖來看，試樣的粘聚力與顆粒摩擦系數和接觸強度有關。在摩擦系數相同的條件下，粘聚力隨接觸強度的增大而增大，而在接觸強度相同的條件下，粘聚力隨摩擦系數的增大而減小。

由此可見，在進行模擬分析時，可以根據得到的材料細觀參數與宏觀特性的關系調整細觀參數的數值，使用PFC方法模擬得到的應力－應變曲線與室內實驗得到的應力－應變曲線相符。

3 邊坡穩定顆粒流模擬

3.1 模型的建立

根據土體宏、細觀力學參數關系，按表2的參數取值，對粘性土坡進行數值模擬分析。

表2 邊坡的細觀參數Table 2 Micro-parameter of slope

建立模型時先構建四面墻體，然后用擴大粒徑的方式生成顆粒，達到室內試驗所用式樣的孔隙率后，讓生成的顆粒在自重下運動，最終達到平衡，平衡后根據所模擬的邊坡形狀，將多余的顆粒刪除，最后形成符合宏觀參數的邊坡模型。為便于觀察滑裂面的發生發展及滑坡破壞的全過程，用不同的顏色標記顆粒［12］。此外，因為刪除多余的顆粒，使模型底部失去約束，應力釋放，故需要在進行一次迭代計算，達到平衡。在迭代計算過程中，體系的最大不平衡接觸力隨迭代次數的增加逐漸收斂于0，即體系達到了力學平衡狀態，運算結束［13］。

3.2 模擬結果及分析

采用PFC2D程序對滑坡的發生、發展的全過程進行仿真模擬。選擇時間為30 000步、70 000步、100 000步、140 000步的4張圖片作為典型圖片來分析，見圖5，從t=30 000步的圖可見，滑坡發生時先從坡頂出現微小的開裂，隨著微裂縫的發生，坡面上裂縫自上向下擴大，坡腳在上部塊體的重壓下滑動加劇，裂縫加大，直至出現剪切破壞，這時，邊坡頂部塊體開始向下快速滑落，邊坡底部土體因受重壓而變形，剪切帶自上而下貫通，最終形成圓弧形滑裂面。

圖5 粘性土坡破壞圖Fig.5 Destruction of cohesive soil slope

圖6是粘性土坡顆粒的位移矢量圖，從顆粒位移的變化及滑裂面形成的位置來看，用PFC2D程序對滑坡進行模擬符合滑坡演化規律，且與室內試驗的模擬結果吻合度較好。

圖6 粘性土坡顆粒位移矢量圖Fig.6 Displacement vectograph of cohesive soil slop

4 結論

1)利用PFC顆粒流軟件模擬雙軸試驗，以此繪制摩爾應力圓，對土試樣的宏觀強度指標與其對應的細觀參數之間的關系進行分析，結果表明，模型中各顆粒之間的摩擦系數，對試樣的摩擦角和試樣的粘聚力都有影響，而顆粒間的接觸強度只影響試樣的粘聚力，對試樣的摩擦角沒有影響。在摩擦系數相同的條件下，粘聚力隨接觸強度的增大而增大，而在接觸強度相同的條件下，粘聚力隨摩擦系數的增大而減小。

2)利用本文所得到的宏－細觀參數之間的關系，用PFC2D程序對宏觀參數已知的土試樣進行模擬，通過與土試樣的室內實驗結果比較，可以得出，用PFC2D程序進行雙軸試驗模擬，所建立的土試樣宏－細觀參數之間的關系是可行的。

3)采用PFC2D程序模擬滑坡變形破壞時，可以不需要先假設滑裂面的位置，而是顆粒根據受到的接觸力自行運動，在抗剪強度最弱的地方，顆粒分開，形成滑裂面。土顆粒在自重下形成滑坡，是土坡因為自身強度過小而形成的，首先在坡頂處形成小的開裂面，隨著時間的變化，開裂面逐漸擴大，最終形成貫穿上下的滑裂面，其受力特點是:滑坡上部受拉，中部受剪，下部受壓，與邊坡的實際受力情況相吻合，充分體現了采用離散元法進行邊坡穩定性分析的優越性。

［1］賈東遠，陰可，李艷華.巖石邊坡穩定性分析方法［J］.地下空間，2004，24(2):250-255.JIA Dongyuan，YIN Ke，LI Yanhua.Rock slope stability analysis method［J］.Underground Space，2004，24(2):250-255.

［2］杜永彬.太大公路徐家寨滑坡運動的數值模擬［J］.重慶交通大學學報，2008，27(6):1099-1102.DU Yongbin.Numerical simulation of the land slide movement in Xujiazhai in the high way from Taiyuan to Dalian［J］.Journal of Chongqing Jiaotong University，2008，27(6):1099-1102.

［3］廖雄華，周健，吳世明，等.粘性土平面應變試驗顆粒流模擬可行性研究［J］.佛山科學技術學院學報，2002，20(2):33-38.LIAO Xionghua，ZHOU Jian，WU Shiming，et al.Feasibility study of the simulation of plane strain test by particle flow code on clay［J］.Journal of Foshan University，2002，20(2):33-38.

［4］劉洋，吳順川，周健.單調荷載下砂土變形過程數值模擬及細觀機制研究［J］.巖土力學，2008，29(12):3199-3216.LIU Yang，WU Shunchuan，ZHOU Jian.Numerical simulation of sand deformation under monotonic loading and mesoscopic mechanism analysis［J］.Rock and Soil Mechanics，2008，29(12):3199-3216.

［5］王培濤，楊天鴻，于慶磊，等.節理邊坡巖體參數獲取與PFC2D應用研究［J］.采礦與安全工程學報，2013，30(4):560-565.WANG Peitao，YANG Tianhong，YU Qinglei，et al.On obtaining jointed rock slope geo-parameters and the application of PFC2D［J］.Journal of Mining and Safety Engineering，2013，30(4):560-565.

［6］徐小敏，凌道盛，陳云敏，等.基于線性接觸模型的顆粒材料細-宏觀彈性常數相關關系研究［J］.巖土工程學報，2010，32(7):991-999.XU Xiaomin，LING Daosheng，CHEN Yunmin，et al.Correlation of microscopic and macroscopic elastic constants of granular materials based on linear contact model［J］.Chinese Journal of Geotechnical Engineering，2010，32(7):991-999.

［7］周博，汪華斌，趙文峰，等.粘性材料細觀與宏觀力學參數相關性研究［J］.巖土力學，2012，33(10):3171-3178.ZHOU Bo，WANG Huabin，ZHAO Wenfeng，et al.Analysis of relationship between particle mesoscopic and macroscopic mechanical parameters of cohesive materials［J］.Rock and Soil Mechanics，2012，33(10):3171-3178.

［8］尹成薇，梁冰，姜利國.基于顆粒流方法的砂土宏-細觀參數關系分析［J］.煤炭學報，2011，36(2):264-267.YIN Chengwei，LIANG Bing，JIANG Liguo.Analysis of relationship between macro-micro-parameter of sandy soil based on particle flow theory［J］.Journal of China Coal Society，2011，36(2):264-267.

［9］張曉平，吳順川，張志增，等.含軟弱夾層土樣變形破壞過程的細觀數值模擬分析［J］.巖土力學，2008，29(5):1200-1204.ZHANG Xiaoping，WU Shunchuan，ZHANG Zhizeng，et al.Numerical simulation and anaiysis of failure process of soil with weak intercalated layer［J］.Rock and Soil Mechanics，2008，29(5):1200-1204.

［10］王宇，李曉，王聲星，等.滑坡漸進破壞運動過程的顆粒流仿真模擬［J］.長江科學院院報，2012，29(12):46-52.WANG Yu，LI Xiao，WANG Shengxing，et al.The particle flow simulation of dynamic process of landslide failure［J］.Journal of Yangtze River Scientific Research Institute，2012，29(12):46-52.

［11］周健，池毓蔚，池永，等.砂土雙軸試驗的顆粒流模擬［J］.巖土工程學報，2000，22(6):701-704.ZHOU Jian，CHI Yuwei，CHI Yong，et al.Simulation of biaxial test on sand by particle flow code［J］.Chinese Journal of Geotechnical Engineering，2000，22(6):701-704.

［12］周健，廖雄華，池永，等.土的室內平面應變試驗的顆粒流模擬［J］.同濟大學學報，2002，30(9):1044-1050.ZHOU Jian，LIAO Xionghua，CHI Yong，et al.Simulating plane strain test of soils by particle flow code［J］.Journal of Tongji University，2002，30(9):1044-1050.

［13］周健，王家全，曾遠，等.土坡穩定分析的顆粒流模擬［J］.巖土力學，2009，30(1):86-90.＇ZHOU Jian，WANG Jiaquan，ZENG Yuan，et al.Simulation of slope stability analysis by particle flow code［J］.Rock and Soil Mechanics，2009，30(1):86-90.