DOI: 10.11817/j.issn.1672-7207.2017.07.022

基于离散单元法的溜砂坡堆积形态数值研究

来志强^{1, 2}，周伟^{1, 2}，杨利福^{1, 2}，马刚^{1, 2}，常晓林^{1, 2}，徐琨^{1, 2}

(1. 武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉，430072；

2. 武汉大学 水工岩石力学教育部重点实验室，湖北 武汉，430072)

摘要：采用离散单元法模拟了溜砂坡的堆积过程，研究粒间摩擦因数和颗粒形状对溜砂坡堆积形态的影响，建立粒间摩擦因数、颗粒形状与堆积形态的相关关系。研究结果表明：颗粒堆积体上、下部的粒径分布均具有分形特性，摩擦因数越大，堆积体上部分形维数和堆积角度越大。对比不同的颗粒形状，颗粒浑圆度越差，堆积体上部分形维数和堆积角度越大。随着粒间摩擦因数的增大，强力链数目减少，弱力链数目增多，力链结构的各向异性增强，颗粒间法向接触力的水平分量变小，接触法向的主方向向竖直方向偏转。圆形颗粒堆积体中力链分布呈树根状，考虑颗粒形状后的堆积体内强力链方向基本与斜坡平行。溜砂坡下部分形维数与其内部孔隙率呈较好的线性关系，下部分形维数越大，堆积体孔隙率越小。

关键词：离散元；溜砂坡；堆积形态；细观机制；颗粒形状；分形

中图分类号：P642        文献标志码：A         文章编号：1672-7207(2017)07-1839-10

Numerical study of accumulation state for sand-sliding slope based on distinct element method

LAI Zhiqiang^{1, 2}, ZHOU Wei^{1, 2}, YANG Lifu^{1, 2}, MA Gang^{1, 2}, CHANG Xiaolin^{1, 2}, XU Kun^{1, 2}

(1. State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China;

2. Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering, Ministry of Education, Wuhan University, Wuhan 430072, China)

Abstract: The discrete element method was employed to simulate the accumulation process of sand-sliding slope. The aim of this study was to explore the relationship among inter-particle friction coefficient, particle shape and accumulation characteristics of sand-sliding slope. The results indicate that the particle size distribution of both the upper and the lower sand-sliding slope has a fractal feature. The fractal dimension of the upper accumulation body and the repose angle of sand-sliding slope increase with inter-particle friction coefficient increasing. Comparing the simulation results with different particle shapes, it is found that the worse the roundness of the particles, the larger the fractal dimension of the upper accumulation body and the repose angle of sand-sliding slope. With the increase of friction coefficient, heterogeneous degree of contact force also increases, that is, the proportion of weak force chain increases while the proportion of strong force chain decreases. Furthermore, the horizontal component of normal contact force decreases and the main direction of the contact normal orientation between particles inclines to the vertical direction. The force distribution among circular particles is like tree roots while the direction of strong force with shaped particles is almost the same as the direction of the chute. There is a linear relationship between the fractal dimension of the lower accumulation body and the accumulation porosity, namely, the larger the former, the smaller the latter.

Key words: distinct element method; sand-sliding slope; accumulation state; meso-mechanism; particle shape; fractal

边坡表面的岩体经物理风化而形成的碎石、砂砾等松散砂石，常常在重力作用下沿坡面溜动，最终以自然休止角堆积成锥状斜坡，即溜砂坡。此类堆积体结构松散，颗粒间的黏聚力几乎为零，完全依靠颗粒间的摩擦和咬合作用自稳。因此，溜砂坡的稳定性较差，极易因外界扰动而发生大面积的二次塌滑，给边坡防治增大了难度。目前，研究松散颗粒运动和堆积形态的方法主要有物理试验和数值模拟等方法。郝明辉等^[1]考虑了细颗粒含量、滑床糙率等因素对堆积物结构的影响，采用室内模型试验研究了颗粒在运动过程的分离问题；张元才等^[2]对天山公路的溜砂坡进行了动力学特性物理试验，认为细颗粒层在失稳过程中起控制性作用；ZHOU等^[3-4]采用室内模型试验与数值模拟相结合的方法研究了颗粒反序对颗粒流动特性的影响；何娜等^[5-6]在现场调查的基础上，通过大型物理试验研究了溜砂坡的运动及堆积特征、自组织特性等；JIANG等^[7]采用物理试验研究了溜砂坡的运动过程及其对挡墙的冲击力。上述研究采用模型试验方法可以直观地再现溜砂坡的形成过程和堆积形态，但是由于试验成本较高，且试验成果的代表性不强，具有一定的局限性^[8]。同时，受试验设备、测试手段的限制，物理模型试验很难从细观层面对溜砂坡的堆积特性进行深入研究。近年来，离散单元法的发展和计算机性能的飞速提高，使得通过细观数值方法研究砂石溜动过程和堆积形态成为可能。王忠福等^[8]利用颗粒流软件PFC^3D模拟了山体滑坡的整个过程，并分析了摩擦因数对砂石堆积特性的影响；王章琼等^[9]利用离散元软件UDEC研究了反倾岩坡发生崩塌灾害的控制性因素；ZHOU等^[10]采用有限差分法(FDM)和离散元方法(DEM)耦合方法分析了东河口滑坡的破坏过程；YANG等^[11]利用非连续变形分析(DDA)方法模拟了砂石体的运动过程。以上研究大多着重于对砂石运动过程的研究，涉及溜砂坡堆积形态的研究则较少。然而，溜砂坡的堆积形态、堆积范围以及内部结构的稳定性，直接影响其灾害的程度和防治工程的设计。因此，本文作者采用商业离散元软件PFC^2D对溜砂坡的堆积过程进行二维数值模拟，考虑了摩擦因数、颗粒形状对堆积体分形特性、内部细观组构分布规律的影响，分析了堆积体的密实度与其分形特性及内部力链结构的关联性，揭示了堆积体内部细观组构参数与其宏观堆积角度的关系，为进一步完善松散颗粒堆积理论提供参考依据。

1  溜砂坡堆积过程离散元模拟

1.1  数值计算模型与试样制备

数值计算模型参照文献[6]中的溜砂坡堆积试验平台的尺寸，如图1所示。图1中：H₁和H₂分别为挡板高度和滑动区高度；d为滑动区坡度；S为底板长度。根据图2所示的级配曲线在挡板区域内生成圆形颗粒1 426个。采用Clump方法生成6类代表性块体(类圆锯齿形、类梯形、类椭圆形、类三角形、类长方形)^[12]对圆形颗粒进行随机替换，并在自重作用下于挡板区域堆积至平衡状态。

图1  数值模拟模型

Fig. 1  Numerical simulation model

图2  颗粒级配曲线及粒径数分布图

Fig. 2  Grading curves and particle size distribution

1.2  数值阻尼与模拟参数

松散砂石在滑坡过程势能转化为动能，其中一部分动能通过颗粒间的摩擦耗散掉，另一部分动能在颗粒碰撞过程中耗散。摩擦因数、颗粒形状和阻尼系数对颗粒体堆积密度、堆积角等堆积特性有较大影响^[13-14]。在动力计算中，离散单元法通过设置阻尼来耗散系统动能，但计算结果对阻尼参数的取值较为敏感^{[3, 15]}。ZHOU等^[3]通过对岩土颗粒离散元参数模拟研究，发现当黏滞阻尼系数为0.2，局部阻尼系数为0~0.15时，数值模型所得颗粒恢复系数与球形岩块斜坡落体试验中所得的0.4~0.6的恢复系数最为接近。离散单元法中常用的阻尼有黏性阻尼和局部阻尼。图3所示为同时存在黏性阻尼和线性接触模型的颗粒—颗粒(墙体)接触处相互作用的示意图。黏性阻尼通过在每个接触上设置法向和切向阻尼来模拟颗粒碰撞耗散能量，黏性阻尼力D_i表达式为

 (i=n, s)             (1)

式中：下标i为n时，表示法向；s表示切向；C_i为阻尼系数；v_i为接触C_i的相对速度。

局部阻尼通过颗粒的运动方程对每个颗粒施加一个与速度相反的阻尼力，表达式为

(j=1, 2, 3)      (2)

式中：为阻尼系数；sign(x)为符号系数；为广义速度：当j=1, 2时表示平动速度，当j=3时表示转动速度。

图3  有黏性阻尼的线性接触受力模型

Fig. 3  Viscous damping activated at a contact with linear contact model

数值计算过程中采用线性刚度接触模型来模拟颗粒间及颗粒与墙体的接触行为，采用黏滞阻尼和局部阻尼来吸收系统能量。颗粒间的摩擦因数一般在0.5~0.8之间，本文根据文献[5]取值为0.7，阻尼系数通过试凑法^[3]获得，即采用Clump方法生成的6种颗粒形状，将1 426个圆形颗粒进行随机等质量替换，在此基础上，模拟了多组在不同阻尼系数下颗粒混合体的堆积过程。将坡体沿高度方向等分为上、中、下3部分。量取坡体每个部位中能够形成长坡面线的堆积角度为该部位的休止角度。当局部阻尼取0.15，法向和切向黏性阻尼均取0.20时，颗粒堆积体的上、中、下3个部位的堆积角度分别为39.0°，40.0°和38.0°，与文献[6]中相应的堆积角度相同，取该组阻尼系数为文章数值模拟的计算参数。本文最终数值模拟参数选取如表1所示。

表1  数值模型参数

Table 1  Parameters of numerical model

1.3  堆积过程模拟

图4所示为不同时刻松散颗粒沿坡面向下运动过程图。在去除挡板后的瞬间，挡板底部(红色)颗粒首先加速向下运动，随后临空面上部(黄色和绿色)颗粒开始向下倾倒运动，而顶部靠近墙体(灰色)颗粒则几乎保持静止的状态。临空面颗粒首先到达坡底，是堆积体中、下层的主要来源。顶部靠近墙体(灰色)颗粒对之前形成的堆积体表层进行冲击与补给，使得堆积体各部位坡度发生改变，这与文献[6]描述的堆积过程基本一致。当堆积体的坡度增加至休止角上限时，之前形成的坡体发生崩塌，系统调整至平衡状态，最终形成稳定的堆积形态。数值模拟的颗粒堆积体中部的倾角最大，上部次之，下部最小。坡体上部主要以细颗粒层面为主，下部则粗颗粒居多，少量粗颗粒滚离坡脚，这与文献[5-6]中溜砂坡的堆积特征相一致。

图4  离散元模拟的颗粒堆积过程

Fig. 4  Numerical simulated process of granular accumulation

2  摩擦因数对溜砂坡堆积特性的影响

为研究摩擦因数对溜砂坡堆积特性的影响，分别取摩擦因数为0.2，0.4，0.6，0.8和1.0，相应的试验编号为A1，A2，A3，A4和A5，对颗粒堆积过程进行模拟，其他计算参数均按表1选取。

2.1  摩擦因数对堆积体分形特性的影响

颗粒堆积体的粒度组成可以采用质量-粒径分布分维模型进行描述^[5-6]，表示为

式中：R为颗粒粒径；M(r＜R)为颗粒粒径r小于R的颗粒质量；M₀为颗粒集合体的总质量；R₀为最大颗粒尺寸；D为粒度的分布维数。

若统计得到的与具有线性关系，根据斜率即可求出分维数D，反之则说明堆积体不具有明显的分形结构。

将堆积体沿高度方向平均分为上、下2部分，分别统计出不同摩擦因数下各粒径区间的颗粒质量占总质量的比例(见表2)，并对分形维数进行拟合，如图5所示。由拟合结果可知：不同摩擦因数下所形成的堆积体上、下2部分均具有良好的分形特性。随着摩擦因数的增大，细颗粒在堆积体上部的比重加大，坡体上部分形维数增大，下部分形维数减小。在摩擦因数一定的情况下，堆积体上部细颗粒的相对密度总是大于下部细颗粒的相对密度，上部相应的分形维数总是大于下部相应的分形维数，这与室内试验结果一致。以上分析结果表明，摩擦因数的增大使得堆积体粗细颗粒的分选现象更加明显。

表2  不同摩擦因数下堆积体颗粒分析成果

Table 2  Particle analytical results of accumulation body under different friction coefficient

图5  不同摩擦因数下堆积体-线性拟合

Fig. 5  Linear fitting relationships between and of accumulation body under different friction coefficients

2.2  摩擦因数对宏观堆积特性的影响

图6所示为不同摩擦因数下颗粒堆积体的坡面线。由图6可见：摩擦因数越大，堆积体的坡面线越陡，坡面线从光滑逐渐变得凹凸不平(如图6箭头所示)。这是因为摩擦因数越大，颗粒间摩擦力也越大，坡体局部颗粒堆积使得坡面线凹凸明显。

图6  不同摩擦因数下堆积体坡面线变化图

Fig. 6  Slope lines of accumulation body under different friction coefficients

不同摩擦因数下试样上、中、下部位的堆积角度如图7所示，其中平均坡度为3个部位坡度的加权平均值。由图7可以看出：每组堆积体的中部倾角最大，上部次之，下部最小，这与文献[5-6]中物理试验得出的结论一致。堆积体3个部位的坡体角度均随着摩擦因数的增大而增加，当摩擦因数增大到一定程度后，坡面线上升的幅度变小，这与图5中堆积体坡面线变化的幅度一致。当摩擦因数增大时，堆积体内部孔隙率增大，密实度降低。当摩擦因数增大到一定程度后，堆积角度和内部孔隙率增幅降低，密实度变化不大。

2.3  摩擦因数对堆积组构的影响

溜砂坡的堆积特性与其细观组构的演化密切相关。ROTHENBURG等^[16]采用傅里叶函数来拟合颗粒间接触法向的角域分布，其表达式为

图7  不同摩擦因数下坡体角度和孔隙率变化曲线

Fig. 7  Repose angles and porosity curves of accumulation body under different friction coefficients

式中：为颗粒接触法向与试样水平面的夹角; 分别为接触法向各向异性的主方向；为反映接触法向各向异性程度的傅里叶系数。

图8所示为不同摩擦因数下堆积体粒间接触法向各向异性分布图和相应的傅里叶函数拟合结果。图8中每9°为1个区间，统计颗粒间接触法向落入角度区间内的个数占总接触数的比例(见图8中左边数据)。随着摩擦因数增大，颗粒间接触法向各向异性程度逐渐增强，主方向逐渐向竖直方向偏转，颗粒体堆积角度增加(坡面线向竖直方向偏转)。

堆积体在重力作用下，颗粒间接触形成力链。定义接触力大于平均接触力的为强力链，反之则为弱力链^[17]。强力链构成了颗粒集合体承载结构的骨架，对系统稳定起决定性作用。图9所示为不同摩擦因数下堆积体法向力和切向力的概率分布图，图9横坐标为力链强度除以接触力的平均值并取对数。由图9可以看出：颗粒间的法向接触力具有“双峰”式的分布规律，随着摩擦因数的增大，堆积体中力链结构发生调整，法向和切向接触力中强力链数目减少，弱力链数目增多，“双峰”式分布规律减弱。

3  颗粒形状对溜砂坡堆积特性的影响

为研究颗粒形状对溜砂坡堆积特性的影响，分别采用圆形、类圆锯齿形、类梯形和类椭圆形、类长方形对初始圆形颗粒进行替换，试验编号为B1，B2，B3，B4和B5。生成颗粒集合体时通过提高颗粒簇的密度来保证初始纯圆颗粒转换为颗粒簇块体时各组试样的质量相同，其他计算参数均按表1取得。

图8  不同摩擦因数下堆积体接触法向各向异性分布图

Fig. 8  Contact orientation anisotropy diagrams of accumulation body under different friction coefficients

图9  法向力和切向力分布图对比

Fig. 9  Distribution comparison between normal force and tangential force

3.1  颗粒形状对堆积体分形特性的影响

将堆积体沿高度方向平均分为上、下2部分。不同颗粒形状下堆积体颗粒分析成果如表3所示，不同颗粒形状下堆积体-线性拟合如图10所示。由表3及图10可知：不同颗粒形状组成的堆积体上、下2部分均具有良好的分形特性，细颗粒主要堆积于坡体上部，大粗颗粒主要堆积于坡体下部。圆形颗粒堆积体上部细颗粒比例最低，其相应的分形维数最小；对于异形颗粒，类圆锯齿形堆积体上部的分形维数最大，类长方形堆积体上部的分形维数最小。

表3  不同颗粒形状下堆积体颗粒分析成果

Table 3  Particle analytical results of accumulation body under different particle shapes

图10  不同颗粒形状下堆积体-线性拟合

Fig. 10  Linear fitting relationships between and of accumulation body under different particle shapes

圆形颗粒堆积体上部的细颗粒比例低，主要是因为圆形颗粒间的咬合作用较弱，粗颗粒之间容易形成孔隙，细颗粒穿过粗颗粒之间的缝隙滚入坡体下部。相对于纯圆形颗粒，异形颗粒由于颗粒间的咬合作用使得细颗粒容易停留在坡体上部，坡体上部的分形维数较大。

3.2  颗粒形状对堆积组构的影响

将坡体沿高度方向三等分为上、中、下3部分，不同颗粒形状组成的堆积体上、中、下3个部位的堆积角度变化曲线如图11所示。由图11可知：各组堆积体均表现为中部倾角最大，上部次之，下部最小。其中，圆形颗粒堆积体平均坡度最小，仅为23.0°，异形颗粒堆积体平均坡度均较大，分别达35.7°，36.7°，37.5°和40.0°(类圆锯齿形、类梯形、类椭圆形和类长方形)。

圆形颗粒在挡板区域内初始孔隙率为0.151，异形颗粒在挡板区域内初始孔隙率分别为0.247，0.245，0.246和0.245(类圆锯齿形、类梯形、类椭圆形和类长方形)。圆形颗粒堆积后孔隙率增大至0.176。类圆锯齿形、类梯形、类椭圆形堆积后孔隙率均增大，分别达0.300，0.281和0.272。类长方形堆积后孔隙率则变化不大，为0.249。圆形堆积体密实度最大，这是因为圆形颗粒之间的咬合作用较弱，在堆积过程当中细颗粒容易填充于粗颗粒之间而使得其孔隙率最低。异形堆积体的孔隙率较圆形堆积体高，这是因为其内部拱架结构形成的空隙所致。相比于类圆锯齿形等其他异形堆积体，类长方形堆积体中较多的细颗粒填充了拱架结构所形成的空隙，使得其整体结构较为密实。

图12所示为不同颗粒形状下堆积形态图和局部力链结构图。圆形颗粒堆积体的坡面最为平滑，异形颗粒堆积体局部凹凸不平的现象明显。圆形堆积体内部力链分布呈网状，几乎没有拱效应，而异形堆积体下部粗颗粒相互咬合，形成明显的拱架结构。

图11  不同颗粒形状下坡体角度和孔隙率变化曲线

Fig. 11  Repose angles and porosity curves of accumulation body under different particle shapes

3.3  颗粒形状对堆积组构的影响

图13所示为不同颗粒形状组成的堆积体内部力链结构分布。圆形堆积体中力链分布呈明显的树根状；考虑颗粒形状后，堆积体中力链分布稀疏且清晰可见，强力链方向基本与斜坡平行，颗粒间形成的拱架结构增多。颗粒的浑圆度越差，颗粒体内部沿斜坡方向的强力链长度增加，在宏观层面上表现为颗粒体的堆积高度和角度增加。

图12  不同颗粒形状的堆积形态图

Fig. 12  Accumulation state of accumulation body under different particle shapes

图13  不同颗粒形状下堆积体力链结构图

Fig. 13  Force chain structure of accumulation body under different particle shapes

4  溜砂坡宏细观堆积特性的联系

颗粒的堆积形态实质是介观尺度颗粒间的咬合作用以及宏观尺度粗细颗粒位置分布的综合作用结果。孔亮等^[18]认为异形颗粒间由于咬合作用而形成的拱架结构是导致堆积体孔隙率增大的主要原因。拱架结构主要由坡体下部粗颗粒形成，细颗粒受该结构的阻挡作用容易堆积于坡体上部。当坡体下部拱架结构较多时，坡体下部的细颗粒比例较低，较多的空隙形成于坡体下部，使得堆积体整体的密实度降低。利用坡体下部分形维数表征堆积体下部细颗粒比例，建立其与堆积体孔隙率线性拟合关系，分别如图14所示。由图14可以发现：坡体下部分形维数越大，细颗粒比例越大，空隙的填充率越高，相应的孔隙率越小。与不同摩擦因数的堆积体相比，不同颗粒形状下的堆积体孔隙率受坡体下部分形维数的影响较为敏感，这是因为颗粒形状的改变对堆积体内部结拱率影响较大。

图14  孔隙率与坡体下部分形维数关系

Fig. 14  Relationship between porosity and fractal dimension under slope body

颗粒体的堆积特性受宏观力学特性的影响显著，而宏观力学特性与细观组构关联性密切^[19]。因此，有必要对颗粒体的宏观堆积特性与其细观组构指标之间的关系进行进一步研究。图15所示为堆积体粒间接触法向傅里叶系数a与堆积体平均坡度的关系。由图15可以看出：堆积体的平均坡度越高，a越大，接触法向各向异性程度越大，两者的关系可以用直线进行拟合。拟合结果表明溜砂坡的宏观堆积特性与接触法向的各向异性规律保持了良好的一致性。

图15  接触法向傅里叶系数与堆积体平均坡度关系

Fig. 15  Relationship between contact orientation anisotropy coefficient and average dip angle

5  结论

1) 溜砂坡上、下2部分均具有良好的分形结构，可以采用分形维数表征细颗粒比例。随着摩擦因数的增大，坡体上部细颗粒比例增大，堆积体比较松散。考虑颗粒形状后，堆积体上部细颗粒比例不同程度地增大，堆积体比较松散。

2) 在堆积形态方面，堆积体中部倾角最大，下部最小，平均坡度随着摩擦因数的增大而增加，坡面线凹凸度加剧。由于粒间咬合作用异形颗粒堆积体平均堆积坡度比圆形颗粒平均堆积坡度大，坡面线凹凸显著。

3) 随着摩擦因数的增大，颗粒间法向接触力的水平分量减小，颗粒间接触法向的各向异性程度增强，堆积体内强力链数目减少，弱力链数目增加。圆形颗粒堆积体内部力链分布呈树根状，接触法向的各向异性程度较低。考虑颗粒形状后，堆积体内部存在较多的拱架结构，强力链方向基本与斜坡平行。

4) 堆积体下部分形维数与其内部孔隙率呈较好的线性关系，分形维数越大，孔隙率越小，堆积体越密实。颗粒体宏观堆积特性与细观组构关联密切，粒间接触法向傅里叶系数a与堆积体的平均坡度呈线性关系，随平均坡度的增加而增大。

参考文献：

[1] 郝明辉, 许强, 杨兴国, 等. 高速滑坡–碎屑流颗粒反序试验及其成因机制探讨[J]. 岩石力学与工程学报, 2015, 34(3): 472-479.

HAO Minghui, XU Qiang, YANG Xingguo, et al. Physical modeling tests on inverse grading of particles in high speed landslide debris[J]. Chinese Journal of Rock Mechanics and Engineering, 2015, 34(3): 472-479.

[2] 张元才, 黄润秋, 傅荣华, 等. 溜砂坡大规模失稳动力学机制试验研究[J]. 岩石力学与工程学报, 2010, 29(1): 65-72.

ZHANG Yuancai, HUANG Runqiu, FU Ronghua, et al. Experimental research on dynamic failure mechanism of large-scale talus slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(1): 65-72.

[3] ZHOU G G D, SUN Q C. Three-dimensional numerical study on flow regimes of dry granular flows by DEM[J]. Powder Technology, 2013, 239: 115-127.

[4] ZHOU G G D, NG C W W. Numerical investigation of reverse segregation in debris flows by DEM[J]. Granular Matter, 2010, 12(5): 507-516.

[5] 何娜. 撒落型散粒体斜坡变形破坏机制试验研究[D]. 成都: 成都理工大学环境与土木工程学院, 2013: 22-38.

HE Na. Deformation failure mechanism test studying on the scattered sand-sliding slope[D]. Chengdu: Chengdu University of Technology, College of Environment and Civil Engineering, 2013: 22-38.

[6] 何娜, 傅荣华, 卜祥航, 等. 散粒体斜坡运动堆积特征试验研究[J]. 地质灾害与环境保护, 2012, 23(1): 50-53.

HE Na, FU Ronghua, BU Xianghang, et al. An experimental study on movement and accumulation characteristics of the sand-sliding slope[J]. Journal of Geological Hazards and Environment Preservation, 2012, 23(1): 50-53.

[7] JIANG Y J, TOWHATA I. Experimental study of dry granular flow and impact behavior against a rigid retaining wall[J]. Rock Mechanics and Rock Engineering, 2013, 46(4): 713-729.

[8] 王忠福, 何思明, 李秀珍. 西藏樟木后山危岩崩塌颗粒离散元数值分析[J]. 岩土力学, 2014, 35(Supp.1): 399-406.

WANG Zhongfu, HE Siming, LI Xiuzhen. Numerical analysis of Tibet Zhangmu dangerous rock collapse by discrete element method[J]. Rock and Soil Mechanics, 2014, 35(Supp.1): 399-406.

[9] 王章琼, 晏鄂川, 尹晓萌, 等. 层状反倾岩质边坡崩塌机理研究: 以湖北鹤峰红莲池铁矿边坡为例[J]. 中南大学学报 (自然科学版), 2014, 45(7): 2295-2302.

WANG Zhangqiong, YAN Echuan, YIN Xiaomeng, et al. Study on collapse mechanism of anti inclined rock slope:A case study of Honglianchi Iron Mine slope in Hefeng, Hubei province[J]. Journal of Central South University (Science and Technology), 2014, 45(7): 2295-2302.

[10] ZHOU Jiawen, CUI Peng, YANG Xingguo. Dynamic process analysis for the initiation and movement of the Donghekou landslide-debris flow triggered by the Wenchuan earthquake[J]. Journal of Asian Earth Sciences, 2013, 76: 70-84.

[11] YANG Qingqing, CAI Fei, SU Zhiman, et al. Numerical simulation of granular flows in a large flume using discontinuous deformation analysis[J]. Rock Mechanics and Rock Engineering, 2014, 47(6): 2299-2306.

[12] 刘洋, 李晓柱, 吴顺川. 多块体形状堆石体碾压颗粒破碎数值模拟[J]. 岩土力学, 2014, 35(11): 3269-3280.

LIU Yang, LI Xiaozhu, WU Shunchuan. Numerical simulation of particle crushing for rockfill of different particles shape under rolling compaction[J]. Rock and Soil Mechanics, 2014, 35(11): 3269-3280.

[13] NG T T, ZHOU Wei, MA Gang, et al. Damping and particle mass in DEM simulations under gravity[J]. Journal of Engineering Mechanics, 2015, 141(6): 04014167-1-04014167-9.

[14] NAKASHIMA H, SHIOJI Y, KOBAYASHI T, et al. Determining the angle of repose of sand under low-gravity conditions using discrete element method[J]. Journal of Terramechanics, 2011, 48(1): 17-26.

[15] 周国庆, 周杰, 陆勇,等. 颗粒流程序(PFC^2D)中阻尼参数的适用性研究[J]. 中国矿业大学学报, 2011, 40(5): 667-672.

ZHOU Guoqing, ZHOU Jie, LU Yong, et al. Selection of damping parameters used in a particle flow code (PFC^2D)[J]. Journal of China University of Mining & Technology, 2011, 40(5): 667-672.

[16] ROTHENBURG L, BATHURST R J. Analytical study of induced anisotropy in idealized granular materials[J]. Geotechnique, 1989, 39(4): 601-614.

[17] 孙其诚, 辛海丽, 刘建国, 等. 颗粒体系中的骨架及力链网络[J]. 岩土力学, 2009, 30(Suppl 1): 83-87.

SUN Qicheng, XIN Haili, LIU Jianguo, et al. Skeleton and force chain network in static granular material[J]. Rock and Soil Mechanics, 2009, 30(Suppl 1): 83-87.

[18] 孔亮, 彭仁. 颗粒形状对类砂土力学性质影响的颗粒流模拟[J]. 岩石力学与工程学报, 2011, 30(10): 2112-2119.

KONG Liang, PENG Ren. Particle flow simulation of influence of particle shape on mechanical properties of quasi-sands[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(10): 2112-2119.

[19] 常晓林, 马刚, 周伟, 等. 颗粒形状及粒间摩擦角对堆石体宏观力学行为的影响[J]. 岩土工程学报, 2012, 34(4): 646-653.

CHANG Xiaolin, MA Gang, ZHOU Wei, et al. Influences of particle shape and inter-particle friction angle on macroscopic response of rockfill[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(4): 646-653.

(编辑  杨幼平)

收稿日期：2016-08-07；修回日期：2016-10-15

基金项目(Foundation item)：国家优秀青年科学基金资助项目(51322905)；中国博士后科学基金面上资助项目(2015M572195) (Project(51322905) supported by the National Science Foundation for Excellent Young Scholars of China; Project(2015M572195) supported by China Postdoctoral Science Foundation)

通信作者：周伟，博士，教授，从事高坝结构数值仿真方面的研究；E-mail: zw_mxx@163.com