边坡失稳破坏的颗粒流模拟

万 江， 温 泉， 昝世明

(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