引 言
边坡失稳问题广泛存在于岩土工程、环境科学、交通运输等诸多领域．边坡受到雨水冲刷、地震以及冲击作用时会发生失稳造成滑坡，对人们造成巨大的生命威胁，对社会造成巨大的经济损失，而对边坡堆积体的运动状态和力学状态的研究能更清楚地对边坡失稳问题进行了解．因此，对土颗粒边坡堆积体模型失稳问题进行研究具有非常重要的实际意义．
对于边坡的稳定性分析，常用极限平衡法与有限元法．但极限平衡法在对各刚性块体进行下滑力和抗滑力计算时忽略了其内部的应力-应变关系，有限元法计算中对边坡不同部位上的力学性质和应力状态等考虑不足，两种方法很难进一步地去描述边坡离散型破坏和运动破坏，难以反映边坡的渐进失稳过程[1-3]．而离散元法计算直接从细观参数出发分析边坡的失稳模式，将边坡中的土颗粒以离散的方式表现出来，更能直观地了解到裂变区域位置、宽度、演化过程、能量转化等．周健等[4]运用颗粒离散元法对砂性土坡和黏性土坡分别进行数值模拟，得到砂性土坡整体破坏时为塑性流动状态，边坡上无明显裂缝，而黏性土坡为脆性破坏，变形超过一定量就发生失稳现象．Fakhimi等[5-8]通过模拟单轴压缩试验，对岩体的变形破坏特性、宏观性质和微观特征进行了研究．杨冰等[9-13]基于PFC2D和EDEM软件，对边坡失稳过程的宏观响应进行了深入研究和探讨．
复杂网络能够有效表征颗粒的介观结构且有利于颗粒的尺寸效应研究，因此复杂网络理论被广泛应用于颗粒填充场的研究中．对于颗粒边坡堆积体的研究，需要从细观上定义颗粒之间的接触关系，而复杂网络能够解决这一问题，对边坡堆积体的位移变形、裂隙发育等，复杂网络能从微观接触研究边坡堆积体模型失稳现象，进行量化分析．近年来，许多学者基于复杂网络理论对颗粒填充场进行了研究，并取得了较为满意的成果[14-16]．但现阶段复杂网络的分析方法多用于分析密实颗粒填充体试样压缩过程，而从复杂网络视角对颗粒边坡堆积体失稳过程中接触力网络的演化分析鲜有研究，这在边坡失稳研究中是一种新的数学分析方法．
本文采用颗粒离散元法，对边坡堆积体稳定性进行了研究，通过金属板的下压对边坡堆积体上部施加一个静荷载以模拟边坡堆积体在上部存在荷载时的破坏失稳过程，并分析了不同时间段边坡堆积体坡顶的力学行为，将模拟结果与相关实验做了对比，最后引进复杂网络理论对黏性土颗粒边坡堆积体坡顶的力链网络进行分析，结合强度折减法，验证了复杂网络方法用于研究边坡堆积体的准确性与可行性．
1 数值仿真方法和网络分析方法
1.1 颗粒离散元算法
本文采用Hertz-Mindlin with bonding黏结模型，模型通过一定尺寸的“黏结键”将颗粒黏结在一起．这种“黏结键”可以承受切向和法向的应力，当达到最大法向应力或切向应力时，该黏结键断裂．黏结键断裂形成开裂面，导致边坡堆积体坡体出现滑移，边坡堆积体失稳破坏．
颗粒间“黏结键”形成之前，颗粒的接触计算采用标准的Hertz-Mindlin接触模型．当颗粒间“黏结键”在时刻tBOND形成之后，颗粒的法向力和力矩的比值(Fn,t/Tn,t)归零，并且应用下列等式在每个计算时间步长对法向力和力矩的值进行更新：
δFn=-vnSnAδt，
(1)
δFt=-vtStAδt，
(2)
δMn=-ωnStJδt，
(3)
(4)
(5)
(6)
式中A表示接触区域面积，Sn和St分别为法向和切向刚度，RB是“黏结键”的半径；δt是时间步长；vn，vt分别是颗粒的法向和切向速度；ωn，ωt分别是法向和切向角速度．当法向和切向应力超过预先设定的阈值时，黏结断裂．因此定义法向应力和切向应力的最大值为
(7)
(8)
1.2 复杂网络理论
复杂网络是指具有自组织、自相似、吸引子、小世界、无标度中部分或全部性质的网络．复杂网络的主要参数包括：度及度分布，集聚程度及其分布特征，最短距离及其分布特征，以及连通集团的规模分布．
1.2.1 节点的度
节点υi的度ki定义为与该节点连接的边数．从直观上看， 一个节点的度越大就意味着这个节点在某种意义上越重要．网络中所有节点υi的度ki的平均值称为网络的平均度， 记为〈k〉，即
(9)
式中，N为网络节点数．
1.2.2 聚类系数
网络的聚类系数度量一个复杂网络的集团化程度，网络的聚类系数越大，说明网络的集团化程度越高，节点之间的关系越紧密．一般地，假设网络G中一个节点i有ki条边和其他节点相连，这ki个节点就称为节点i的“邻居”．显然，这ki个节点间最多可能有ki(ki-1)/2条边，ki个节点之间实际存在的边数Ei与总的可能边数的比值定义为节点i的聚类系数Ci，网络中所有节点i的聚类系数的平均值就是网络的聚类系数，用C表示，即
(10)
(11)
1.2.3 平均路径长度
节点的距离dij定义为连接两个节点i和j的最短路径上的边数．网络中的平均路径长度L定义为任意两个节点间距离的平均值：
(12)
式中，N为网络节点数．网络的平均路径长度也称为网络的特征路径长度，用来衡量网络节点间的信息传递效率．
2 数值仿真模拟
2.1 计算模型
本文采用离散元软件对边坡堆积体失稳过程进行数值模拟计算，离散元模型如图1所示．边坡堆积体模型由1 637个直径为3 mm的均匀颗粒构成，对边坡堆积体左侧、右侧及底面施加固定端约束，在边坡堆积体坡顶不同位置A，B，C，D(x=-52.5，-62.5，-72.5，-82.5 mm)处沿y轴以一定的速度对土体边坡进行加载(原点为模型的中心，即(200/2,120/2)处)，施加作用力的金属板的宽度为10 mm,坡顶测点E为监测点．
图1 边坡堆积体模型
Fig. 1 The slope pile model
本文中的颗粒材料选取硬黏土，土颗粒材料的基本特征参数及黏结模型参数通过压缩试验确定[17]，如表1、2所示．
表1 颗粒材料参数
Table 1 Granular material parameters
density ρ/(kg/m3)shear modulus G/MPastiffness ratio K/(N/m)Poisson’s ratio νstatic friction coefficient μrolling friction coefficient σ2 0004691.50.350.50.001
表2 黏结颗粒模型参数
Table 2 BPM parameters
normal stiffness Sn/(GN/m3)tangential stiffness St/(GN/m3)normal critical pressure σmax/MPatangential critical pressure τmax/MPaglue radius RB/mm1 67066736240.8
2.2 失稳过程
边坡堆积体模型在上部荷载作用下逐渐发生失稳破坏，图2是当金属板作用在位置B(x=-62.5 mm)时，边坡堆积体土颗粒速度矢量随加载时间变化的速度矢量图．在加载开始时，金属板两侧的土体由于挤压作用开始往左右两边运动，左侧土体由于左侧存在的固定约束，坡顶土体出现隆起的趋势，而右侧土体由于靠近自由端，在自重和下压板的作用下出现向右下方运动的趋势．随着金属板的持续向下作用，土体颗粒沿边坡方向出现裂纹和剪切滑移破坏．当停止加载后，边坡堆积体大部分土颗粒逐渐恢复稳定状态，只有少量的土颗粒在重力作用下，脱离边坡堆积体向下滑动．
(a) t=0.2 s (b) t=0.6 s
(c) t=0.8 s (d) t=1 s
图2 边坡堆积体颗粒速度矢量图
Fig. 2 Particle velocity vector diagrams for the slope pile
图3为颗粒堆积体边坡的两组离散元实验结果[18]，将其与图2进行对比．从模拟得出的边坡堆积体的变化趋势及滑移开裂面的形成位置可以看出，失稳路径主要为陡倾的断层和向坡外的缓倾角裂隙组合成的阶梯状形式，失稳模式主要是局部的剪切-拉张失稳模式,与实验结果基本相似，验证了本离散元模型模拟符合滑坡演化规律．边坡强度一定的条件下，边坡破坏的基本形式是滑移面向坡面方向移动，结果与Katz等[19]采用DEM得到的边坡滑移规律相同．
图3 实验结果[18]
Fig. 3 The results of experiment[18]
2.3 荷载作用位置的影响
为探究荷载作用位置对边坡堆积体失稳的影响，对A，B，C，D不同位置的边坡堆积体失稳过程进行离散元模拟，结果如图4所示，边坡堆积体滑移开裂面的角度随着金属板作用位置的变化而变化．当对位置A进行加载时，边坡堆积体滑移开裂面的角度为63°，略大于土体边坡堆积的角度；当对位置B进行加载时，边坡堆积体滑移开裂面的角度为40°；当对位置C进行加载时，边坡堆积体滑移开裂面的角度为31°；当对位置D进行加载时，边坡堆积体仅仅在加载板正下方受到挤压，加载板两侧土体保持稳定状态，相对运动的趋势较小，没有出现滑移开裂面．
根据边坡模型在坡顶局部载荷作用下的土体变化规律实验[20]，土体边坡发生破坏时，水平向张拉变形最大值发生在坡体上部，上层土体在局部载荷作用下更容易发生变形．再由本文研究的黏性土颗粒堆积体施加的约束方向以及颗粒特性，所以选择坡顶E点作为监测点，结合之前四组边坡堆积体失稳的模拟效果，对坡顶测点E的y方向平均速度进行实时监测．
图5是不同作用位置时边坡堆积体坡顶测点E的y方向平均速度-时间曲线图．对比A，B，C三个不同作用位置时边坡堆积体坡顶测点E的平均速度-时间曲线可知，随着金属板作用位置离边坡堆积体坡顶测点E越远，测点E向上运动速度的峰值出现得越晚，这是因为加载点越远，金属下压板对土体的挤压效应传递时间越长；且随着金属板作用位置离边坡堆积体坡顶测点E越远，测点E向下运动速度的峰值也出现得越晚，这是因为金属板加载点越远，边坡堆积体坡顶的土体出现滑移开裂面的时间越晚．由作用位置为D时的边坡堆积体坡顶测点E的平均速度-时间曲线可知，当金属板作用位置离边坡堆积体坡顶测点足够远时，金属板对土体的挤压效应虽然可以传递到边坡堆积体坡顶，但是不会出现滑移开裂面．
(a) A(x=-52.5 mm) (b) B(x=-62.5 mm)
(c) C(x=-72.5 mm) (d) D(x=-82.5 mm)
图4 不同作用位置边坡堆积体土颗粒速度云图
Fig. 4 Velocity contours of soil particles in slope piles for different loading positions
注 为了解释图中的颜色，读者可以参考本文的电子网页版本．
图5 不同作用位置时测点E的y方向平均速度-时间曲线
Fig. 5 The average velocity-time curves in the y-direction at measure point E of the slope pile for different loading positions
3 基于复杂网络边坡堆积体失稳分析
本文使用复杂网络理论的统计工具来分析土颗粒填充结构的网络结构，量化边坡堆积体土体颗粒在外荷载作用下介观结构的演化规律，从而能更好地表征边坡堆积体从静态堆积、土体开裂再到滑移失稳的力学特性，以及边坡堆积体失稳过程的特征．利用PAJEK软件计算边坡堆积体土体颗粒复杂网络拓扑特性的相关参数．
3.1 基于复杂网络的失稳研究
图6是边坡堆积体土体颗粒自然堆积下的力链网络．为了找到滑移面出现的位置，本文将接触力网络划分为Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ共5个区域，并分别对网络参数平均度、平均最短路径、聚类系数进行计算和分析．
图7～9分别是金属板作用在位置B时不同区域的平均度、平均最短路径、聚类系数随时间的演化过程．从图7可知在加载开始时，Ⅱ区颗粒之间的接触开始出现少量的断裂，导致颗粒接触力链网络的平均度出现较为平缓的下降．随着荷载的持续作用，接触力出现较大的缺失，边坡堆积体坡顶土颗粒接触力链网络的平均度下降速度加快．而Ⅰ、Ⅲ、Ⅳ、Ⅴ区由于固定端约束，挤压的变形较小，颗粒接触力网络变化也较小，网络的平均度基本没有发生变化．从图8可知加载刚开始时，Ⅱ区的颗粒接触力网络的平均最短路径相对稳定，在加载一段时间后，边坡堆积体坡顶土颗粒接触力网络的平均最短路径出现跳跃式骤增，最后网络的平均最短路径又保持相对稳定．从图9可以看出，在加载的整个过程中，Ⅱ区的接触力网络聚类系数随时间的推移逐渐增大，而Ⅰ、Ⅲ、Ⅳ、Ⅴ区的颗粒接触力网络变化较小，网络的聚类系数基本没有发生变化．
图6 边坡堆积体土体颗粒自然堆积下的力链网络 图7 不同区域的接触力网络平均度-时间曲线
Fig. 6 The force chain network under natural accumulation Fig. 7 Average degree-time curves of the contact of soil particles in the slope pile force network in different areas
图8 不同区域的接触力网络平均最短路径-时间曲线 图9 不同区域的接触力网络聚类系数-时间曲线
Fig. 8 Average shortest path-time curves of the Fig. 9 Clustering coefficient-time curves of the contact force network in different areas contact force network in different areas
分析图7～9网络参数的演化趋势，发现Ⅱ区的网络变化最大且最能反映失稳过程，所以在之后只选取Ⅱ区的网络进行研究．
图10是加载过程中边坡堆积体坡顶土颗粒接触力网络平均度-时间曲线．由图10可知，当金属板作用位置为A，B，C(x=-52.5，-62.5，-72.5 mm)时，加载开始后，由于金属板对土体的挤压，颗粒之间的接触出现少量的断裂，边坡堆积体坡顶土颗粒接触力链网络的平均度出现较为平缓的下降；在加载一段时间后，边坡堆积体出现裂纹，延展开裂，接触力出现较大的缺失，边坡堆积体坡顶土颗粒接触力链网络的平均度下降速度加快；在形成滑移开裂面后，边坡堆积体坡顶形成相对稳定的整体，接触网络变化较小，平均度变化也较小．当金属板作用位置为D(x=-82.5 mm)时，由于金属板加载位置离边坡堆积体坡顶位置较远，传递到边坡堆积体坡顶位置的土体挤压效应比较微弱，颗粒接触力断裂的数量较少，因此，边坡堆积体坡顶土颗粒接触网络变化较小，平均度变化也较小．对比A，B，C三条边坡坡顶颗粒接触力网络平均度-时间曲线，可以看出，随着金属板作用位置离边坡堆积体坡顶的顶点越远，当出现滑移开裂面时，最终接触力网络的平均度越大．
图10 平均度-时间曲线 图11 聚类系数-时间曲线
Fig. 10 Average degree-time curves Fig. 11 Clustering coefficient-time curves
图11是加载过程中边坡堆积体坡顶土颗粒接触力网络聚类系数-时间曲线．由图可知，在整个加载过程中，聚类系数随时间的推移而增大，而且随着金属下压板作用位置与边坡堆积体坡顶顶点距离的靠近，增长的速率也在增大，最终网络聚类系数的值也越大．加载过程中，对边坡堆积体坡体的土颗粒有一个局部的挤压作用，导致坡顶接触力网络会发生局部团聚，且金属板离边坡堆积体坡顶的距离越近，坡体变形越大挤压作用越明显，网络局部团聚越紧密．
图12是加载过程中边坡堆积体坡顶土颗粒接触力网络平均最短路径-时间曲线．由图可知，当金属板作用位置为A,B,C(x=-52.5,-62.5,-72.5 mm)时，加载刚开始时，传递到边坡堆积体坡顶位置的挤压效应较小，边坡堆积体坡顶土颗粒接触力网络的平均最短路径相对稳定;在加载一段时间后，传递到边坡堆积体坡顶位置的挤压效应增大，边坡堆积体出现裂纹，延展开裂，边坡堆积体的力传递效率降低，边坡堆积体坡顶土颗粒接触力网络的平均最短路径出现断崖式突变;在形成滑移开裂面后，边坡堆积体坡顶形成相对稳定的整体，网络的平均最短路径保持相对稳定．当金属板作用位置为D(x=-82.5 mm)时，由于金属板加载位置离边坡堆积体坡顶位置较远，传递到边坡堆积体坡顶位置的土体挤压效应比较微弱，因此， 边坡堆积体坡顶土颗粒接触网络变化较小， 平均最短路径变化也较小．对比A， B， C(x=-52.5，-62.5，-72.5 mm)三条边坡堆积体坡顶颗粒接触力网络平均最短路径-时间曲线，可以看出，金属板作用位置离边坡堆积体坡顶顶点越远，平均最短路径的值出现断崖式突变的时间越晚，变化量也越小．
图12 平均最短路径-时间曲线 图13 加载位置-时间散点图
Fig. 12 Average shortest path-time curves Fig. 13 Loading position-time scatter plots
图13是不同加载位置，坡顶颗粒沿y方向向下运动速度出现峰值的时间散图和平均最短路径出现骤增的时间散图．由图可知，坡顶颗粒运动速度出现峰值的时间和平均最短路径出现骤增的时间几乎是同时的，因此我们可以用平均最短路径出现骤增的时间来预测坡顶颗粒向下运动速度出现最大值的时间．
3.2 基于强度折减法的安全系数
强度折减法即强度参数降低法[21]，基本原理就是不断降低材料的强度参数直至破坏．因模拟使用的是离散元法，在考虑细观颗粒时，切向黏结强度与法向黏结强度劣化作用是不同的，切向黏结强度的劣化作用远大于法向黏结强度[22]．因此折减时将法向黏结系数定义为1，对切向黏结系数进行折减，在颗粒流离散元中，摩擦因数的劣化作用也远远小于切向黏结强度，但为了提高精度方便其他方法验证也纳入折减计算中．
(13)
(14)
式中，c，φ分别为黏性土颗粒真实黏聚力、内摩擦角；k为折减系数；ck，φk分别为折减后的黏性土颗粒的黏聚力、内摩擦角．边坡的安全系数即边坡临界状态的折减系数：
Fs=kcritical．
(15)
对于堆积体边坡失稳如何判断有很多种判别方法，其中坡顶竖直方向的位移是否突变可作为边坡失稳的判据[23]．对于不同加载位置A，B，C(由以上分析可知，D点加载位置太远，影响较小，故不做分析)，通过关键点E的速度时程曲线的趋势判断边坡的安全系数[23]．三个不同加载位置对应不同加载时间计算所得折减系数为1.12，而当系数为1.13时，三个加载位置对应速度时程在加载时段末不归于零，而是保持一固定状态基本不变．由此可以看出，当折减系数为1.12时，边坡处于临界状态．而不同加载位置A，B，C加载时E的竖直和水平方向速度时程曲线分别在0.2，0.3，0.6 s时趋于零，这与平均最短路径出现骤增的时间几乎相同，误差小于0.01 s．验证了运用平均最短路径骤增起点来预测边坡堆积体失稳的准确性以及可行性．
4 结 论
1) 本文将一种新的分析方法——复杂网络，引入边坡堆积体失稳宏观响应与边坡堆积体土颗粒接触力网络的拓扑特性研究，结合实验模型和强度折减法，验证了本文方法的准确性与可行性．采用离散元法和复杂网络理论结合研究边坡堆积体失稳问题，既克服了有限元法仅适用于小变形的缺陷，也解决了传统离散元法计算岩土粒子连通性的高运算成本问题，相较于极限平衡法，边坡变形的演化过程得到了呈现．
2) 在金属板作用下的边坡堆积体坡顶颗粒接触力链网络演化过程中，接触力网络平均度随时间的变化趋势大致一致，但是随着金属下压板作用位置离边坡堆积体坡顶顶点越近，接触力网络的平均度下降趋势越明显，最终的接触力网络的平均度越大．在一定范围内，金属板作用位置离边坡堆积体坡顶顶点越远，平均最短路径的值出现跳跃式骤增的时间越晚，变化量也越小．聚类系数随着加载的持续，呈线性增长．加载位置离边坡堆积体坡顶顶点越近，聚类系数的增长速度越快，当加载位置远离顶点时，边坡堆积体相对稳定，聚类系数变化较小．
3) 边坡堆积体坡顶测点土颗粒运动速度出现峰值的时间和平均最短路径出现骤增的时间几乎是同时的，且与得出边坡安全系数的速度时程趋于零时间一致．因此，可以用平均最短路径出现骤增的时间来预测坡顶颗粒向下运动速度出现最大值的时间，从而预测滑移面出现的时间．
4) 复杂网络参数中的平均最短路径是以上三个网络参数中，在描述边坡堆积体失稳过程时最有效的．平均最短路径的骤增可以对边坡堆积体失稳进行预警，当边坡颗粒接触力网络的平均最短路径出现急剧增长时，边坡堆积体的传力效率也会发生急剧的下降，边坡堆积体出现滑移面，发生失稳破坏．
[1] ZHANG T, ZHENG H, SUN C. Global method for stability analysis of anchored slopes[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2019, 43(1): 124-137.
[2] ZHENG H, SUN G H, LI C G. Cauchy problem of three-dimensional critical slip surfaces of slopes[J]. International Journal for Numerical & Analytical Methods in Geomechanics, 2011, 35(4): 519-527.
[3] 李明霞, 董联杰. 层状反倾边坡变形特征及影响因素分析[J]. 计算力学学报, 2015, 32(6): 831-837.(LI Mingxia, DONG Lianjie. Analysis on influential factors and deformation characteristics of toppling slope[J]. Chinese Journal of Computational Mechanics, 2015, 32(6): 831-837.(in Chinese))
[4] 周健, 王家全, 曾远, 等. 土坡稳定分析的颗粒流模拟[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.(in Chinese))
[5] FAKHIMI A, CARVALHO F, ISHIDA T, et al. Simulation of failure around a circular opening in rock[J]. International Journal of Rock Mechanics & Mining Sciences, 2002, 39(4): 507-515.
[6] POTYONDY D O, CUNDALL P A. A bonded-particle model for rock[J]. International Journal of Rock Mechanics & Mining Sciences, 2004, 41(8): 1329-1364.
[7] BACKSTROM A, ANTIKAINEN J, BACKERS T, et al. Numerical modelling of uniaxial compressive failure of granite with and without saline porewater[J]. International Journal of Rock Mechanics & Mining Sciences, 2008, 45(7): 1126-1142.
[8] HSIEH Y M, LI H H, HUANG T H, et al. Interpretations on how the macroscopic mechanical behavior of sandstone affected by microscopic properties: revealed by bonded-particle model[J]. Engineering Geology, 2008, 99(1): 1-10.
[9] 杨冰. 边坡动力破坏过程及塌落范围的颗粒力学模拟[D]. 硕士学位论文. 北京: 清华大学, 2011.(YANG Bing. Particle dynamics simulation of slope dynamic failure process and collapse range[D]. Master Thesis. Beijing: Tsinghua University, 2011.(in Chinese))
[10] 张小雪. 基于颗粒流模拟的边坡稳定性分析[D]. 硕士学位论文. 哈尔滨: 哈尔滨工程大学, 2015.(ZHANG Xiaoxue. Slope stability analysis based on particle flow simulation[D]. Master Thesis. Harbin: Harbin Engineering University, 2015.(in Chinese))
[11] 杨玲. 黄土崩塌的离散元数值模拟[D]. 硕士学位论文. 西安: 西安科技大学, 2015.(YANG Ling. Discrete element numerical simulation of loess collapse[D]. Master Thesis. Xi’an: Xi’an University of Science and Technology, 2015.(in Chinese))
[12] 冯春, 李世海, 孙厚广, 等. 颗粒接触无网格法及其在边坡成灾范围模拟中的应用[J]. 岩土力学, 2016, 37(12): 3608-3617.(FENG Chun, LI Shihai, SUN Houguang, et al. Particle contact-based meshfree method and its application to slope disaster range simulation[J]. Rock and Soil Mechanics, 2016, 37(12): 3608-3617.(in Chinese))
[13] 邹育. 边坡稳定性的多尺度分析[D]. 硕士学位论文. 广州: 华南理工大学, 2017.(ZOU Yu. Multi-scale analysis of slope stability[D]. Master Thesis. Guangzhou: South China University of Technology, 2017.(in Chinese))
[14] WALKER D M, TORDESILLAS A, SMALL M, et al. A complex systems analysis of stick-slip dynamics of a laboratory fault[J]. Chaos: an Interdisciplinary Journal of Nonlinear Science, 2014, 24(1): 013132.
[15] TORDESILLAS A, PUCILOWSKI S, WALKER D M, et al. A complex network analysis of granular fabric evolution in three-dimensions[J]. Dynamics of Continuous, Discrete and Impulsive Systems(Series B): Applications and Algorithms, 2012, 19(4/5): 471-495.
[16] 宜晨虹, 苗天德, 慕青松. 颗粒介质力链的复杂网络法研究[C]//2014颗粒材料计算力学会议. 兰州, 2014.(YI Chenhong, MIAO Tiande, MU Qingsong. Research on complex network method of particle media force chain[C]//Conference of Computational Mechanics of Granular Materials(CMGM 2014). Lanzhou, 2014.(in Chinese))
[17] QUIST J, EVERTSSON C M. Cone crusher modelling and simulation using DEM[J]. Minerals Engineering, 2015, 85: 92-105.
[18] 汪儒鸿, 周海清, 彭国园. 堆积体边坡突变失稳特性的离散元模拟分析[J]. 兵器装备工程学报, 2018, 39(6): 192-196.(WANG Ruhong, ZHOU Haiqing, PENG Guoyuan. Dispersion element simulation analysis of mutant instability of piles slope[J]. Journal of Ordnance Equipment Engineering, 2018, 39(6): 192-196.(in Chinese))
[19] KATZ O, MORGAN J K, AHARONOV E, et al. Controls on the size and geometry of landslides: insights from discrete element numerical simulations[J]. Geomorphology, 2014, 220: 104-113.
[20] 廖静薇. 基于颗粒流强度折减法的粉质粘土边坡稳定性分析[D]. 硕士学位论文. 重庆: 重庆大学, 2014.(LIAO Jingwei. Slope stability analysis of silty clay particle flow based on strength reduction method[D]. Master Thesis. Chongqing: Chongqing University, 2014.(in Chinese))
[21] 宋占璞, 施斌, 汪义龙, 等. 削坡作用土质边坡变形分布式光纤监测试验研究[J]. 工程地质学报, 2016, 24(6): 1110-1117.(SONG Zhanpu, SHI Bin, WANG Yilong, et al. Analysis on the strain-field of soil cut slope based on DFOS technology monitoring[J]. Journal of Engineering Geology, 2016, 24(6): 1110-1117.(in Chinese))
[22] 史卜涛, 张云, 张巍. 边坡稳定性分析的物质点强度折减法[J]. 岩土工程学报, 2016, 38(9): 1678-1684.(SHI Butao, ZHANG Yun, ZHANG Wei. Strength reduction material point method for slope stability[J]. Journal of Geotechnical Engineering, 2016, 38(9): 1678-1684.(in Chinese))
[23] 陈永明, 滕光亮, 石玉成, 等. 地震作用下宝成铁路109隧道边坡失稳机理的离散元模拟[J]. 岩土工程学报, 2013, 35(S1): 23-32.(CHEN Yongming, TENG Guangliang, SHI Yucheng, et al. Dispersion element simulation of the instability mechanism of the 109 tunnel slope of Baocheng Railway under the influence of earthquake[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(S1): 23-32.(in Chinese))