工程地质学报  2018, Vol. 26 Issue (5): 1139-1154   (#KB#)    
Article Options
  • PDF (#KB#)
  • Full Text HTML
  • Abstract
  • Figures
  • References
  • History
  • 收稿日期:2018-04-10
  • 接受日期:2018-07-19
  • 扩展功能
    把本文推荐给朋友
    加入引用管理器
    Email Alert
    文章反馈
    浏览反馈信息
    本文作者相关文章
    王颖
    庄建琦
    李威
    赵勇
    贾艳军

    引用本文  

    王颖, 庄建琦, 李威, 等. 2018. 地震作用下黄土斜坡失稳及运动过程的离散元模拟[J]. 工程地质学报, 26(5): 1139-1154. doi: 10.13544/j.cnki.jeg.2018026.
    WANG Ying, ZHUANG Jianqi, LI Wei, et al. 2018. Discrete element simulation of instability and movement process of loess slope under seismic loads[J]. Journal of Engineering Geology, 26(5): 1139-1154. doi: 10.13544/j.cnki.jeg.2018026.

    地震作用下黄土斜坡失稳及运动过程的离散元模拟
    王颖, 庄建琦①②, 李威, 赵勇, 贾艳军    
    ① 长安大学地质工程与测绘学院 西安 710054;
    ② 长安大学西部矿产资源与地质工程教育部重点实验室 西安 710054
    摘要:在山区,地震诱发滑坡造成的人员伤亡往往占到地震总伤亡人数的一半以上,尤其是在黄土地区,由于地震造成的黄土滑坡具有运动距离远等特征,其灾难性更严重。针对地震诱发黄土斜坡的失稳峰值加速度、临界位移和运动距离等问题,本文利用离散元(PFC)方法,通过对室内三轴实验的应力-应变曲线进行标定和高精度航拍三维地形数据进行转化,以宁夏海口村黄土斜坡为研究对象,开展三维地震作用下斜坡失稳破坏和运动过程的数值模拟研究。通过监测不同位置在地震作用下的应力分量,计算监测点的pq值并与室内三轴实验所得到p-q破坏线进行比较,获得斜坡破坏过程中的应力路径。结合颗粒的监测位移,获得斜坡失稳破坏的临界位移,得出该黄土斜坡失稳峰值加速度为0.135 g,临界位移为50 cm。同时根据不同地面环境条件,预测了在不同摩擦系数下斜坡失稳后的危害范围,为黄土地区边坡工程的抗震设计及防震减灾工作提供一种新的可视化方法。
    关键词地震滑坡    离散元模拟    室内三轴实验    峰值加速度    临界位移    可视化    
    DISCRETE ELEMENT SIMULATION OF INSTABILITY AND MOVEMENT PROCESS OF LOESS SLOPE UNDER SEISMIC LOADS
    WANG Ying, ZHUANG Jianqi①②, LI Wei, ZHAO Yong, JIA Yanjun    
    ① School of Geological Engineering and Geomatics, Chang'an University, Xi'an 710054;
    ② Key Laboratory of Western China's Mineral Resources and Geological Engineering, Xi'an 710054
    Abstract: In mountainous areas, casualties caused by earthquake-induced landslides often account for more than half of the total earthquakes. Especially in the loess area, due to the far-distance loess landslide caused by the earthquake, the disastrousness is even more serious. This paper attempts to solve the problems of earthquake-induced loess slope instability such as peak acceleration, critical displacement and moving distance. It uses discrete element(PFC)method, calibrates the stress-strain curve of the indoor three-axis experiment, and transforms high-precision aerial three-dimensional topographic data. It takes the loess slope in Haikou Village, Ningxia Province as the research object to carry out numerical simulation of slope unstable failure and movement process under the three-dimensional earthquake. It monitors the stress components at different locations under the earthquake, compares the p and q values of the monitoring points with the p-q failure lines obtained from the indoor three-axis experiments. It obtains the stress paths and peak accelerations in the process of slope failure. Combining with the maximum displacement of the monitored particles, the critical maximum displacement of slope failure is obtained. The peak acceleration of the slope failure is 0.135g and the critical displacement is 50 cm. At the same time, according to different ground environmental conditions, the damage range of the slope after sliding failure under different friction coefficients is predicted. The paper provides a new visualization method for the seismic design and earthquake prevention and mitigation of slope engineering in the loess area.
    Key words: Earthquake landslides    Discrete element simulation    Indoor triaxial experiment    Peak acceleration    Critical displacement    Visualization    

    0 引言

    我国是一个地质构造活动剧烈,地震频发的国家,一次强震往往诱发大量的滑坡灾害并造成大量的人员死亡和财产损失。根据历史地震诱发滑坡的研究,当地震震级(MS≥4)时便会触发地震滑坡灾害,在强震(MS≥6)时,地震触发滑坡地质灾害更为明显(Keefer,1984张铎等,2013)。一次大型地震可触发数以万计的滑坡且地震滑坡造成的生命财产损失达到地震所造成的损失的50%,甚至更多(张铎等,2013)。1920年宁夏海原MS8.5级特大地震,诱发滑坡657处,地震滑坡造成死伤达10万人之多,约占地震死亡人数的一半之余(孔纪明等,2009)。1933年四川叠溪MS7.5级地震,造成叠溪古城被毁,500人直接丧生,引起数处大型滑坡堵塞岷江后溃堤,造成下游2500余人死亡(柴贺军等,1995)。2005年Kashmir MW7.6级地震中,25 500人死于地震滑坡(Dunning et al., 2007)。2008年四川汶川MS8.0级地震,滑坡崩塌等地质灾害估计造成2万人死亡,约占地震造成死亡人数的1/3(黄润秋,2009)。如何评价斜坡在地震作用下的稳定性及其失稳后造成的危害范围是预测地震诱发滑坡的难点之一,也是研究地震诱发滑坡的一个重要方向。

    现阶段对地震滑坡的机理还有许多未知之处,因此对地震滑坡的预测有一定的困难。目前针对地震诱发滑坡的预测主要集中在对滑坡的稳定性评价、诱发滑坡峰值加速度、滑坡地震动响应、滑坡的滑距和滑速(张铎等,2013)。对于地震滑坡稳定性评价,Terzaghi et al.(1948)将拟静力法应用到斜坡稳定性评价中,将地震作用简化为水平和竖直加速度,根据极限平衡原理,将加速度作用下的潜在滑体分解,得出斜坡的稳定系数。Seed et al. (1968)利用岩体动三轴实验并考虑地震过程中的惯性力,利用力多边形计算土条块的动力稳定系数。Xu(1991)认为黏土斜坡在某些情况下采用拟静力法,对于给定圆心的圆形滑动面没有最小安全系数,得不到临界滑动面,因此无法利用拟静力法计算斜坡稳定性。拟静力法在计算斜坡稳定系数时,没有考虑斜坡的变形对斜坡稳定性的影响,因此Newmark(1965)提出Newmark位移分析方法,该方法是基于位移评价地震斜坡稳定性最常用的工具之一,利用永久累计位移是否大于一定值(临界位移)时来判断滑坡是否失稳,但是在区域评价上,斜坡的临界加速度与岩土的类型、孔隙水压力、坡体几何形态有关,因此其在均质区域土地区应用较广,但是在地质结构复杂地区很难推广(王涛等,2013)。Wieczorek(1985)基于Newmark方法和以往地震滑坡位移记录,以不同地下水条件、不同临界加速度和滑坡位移为5 cm作为判别界限对圣马特奥地区的斜坡进行敏感性划分。Jibson(1993)通过Newmark方法计算厄尔士山区的黏性斜坡在地震作用下永久位移大于10 cm的斜坡位置,计算结果与真实地震诱发的滑坡位置相吻合。对于峰值加速度研究,Wilson et al. (1985)根据美国历史地震记录数据和诱发滑坡距离发现触发滑坡的峰值加速度下限为0.05 g,并给出了震级与滑坡分布范围公式。黄台丰(1999)对1998年台湾瑞里地震诱发滑坡研究发现地震诱发滑坡最大峰值加速度为2.5 g以上。在地震诱发滑坡滑动距离的研究上,刘忠玉等(2000)假定滑体运动形式是连续可变的,建立了预测高速远程滑坡的块体运动模型,该模型可预测滑坡的最大滑速和最大滑距。郝明辉等(2015)以文家沟滑坡为依据,进行沟槽实验,表明随着下垫面的粗糙度增加,滑体的滑动距离明显减少。Zhou et al.(2015)结合室内模型实验对易贡滑坡的滑动距离进行分析,认为下垫面的干湿对滑距有着很大的影响。综上所述,针对地震诱发滑坡的诱发条件预测研究一直处于探索阶段,目前针对单体滑坡的失稳峰值加速度和临界位移等研究较少。

    随着计算机的发展,数值模拟成为研究地震诱发滑坡的主要方法。夏敏等(2010)利用FLAC3D软件对黑龙江上某滑坡的复活机制进行研究分析,得出滑坡对地震波的敏感度高于周围基岩。崔芳鹏等(2010)利用UDEC软件对唐家山滑坡在地震纵横波时差耦合作用下产生的崩滑破坏的动力全过程模拟认为地震纵波产生的竖向加速度起到优势破坏作用。崔铁军等(2016)基于PFC研究不同地震加速度情况下尾矿库内部结构的变形和稳定性,得出了尾矿库发生位移的颗粒都集中在初期坝承受土压力的部分的结论。孙新坡等(2017)采用PFC软件对汶川地震诱发的牛圈沟滑坡动力过程进行数值模拟,并通过与实际滑坡形态和运动过程特点相对比反演了离散元参数。以上数值模拟大多数基于二维模拟情况下进行研究分析,对于反映真实的启动和运动过程存在困难。三维离散元对滑坡的运动过程有着天然的优势,张龙等(2012)利用三维颗粒流软件PFC3D对鸡尾山高速远程滑坡堆积过程进行模拟。Lu et al.(2014)基于三维离散元对台湾芦山温泉区不稳定斜坡进行危险性预测,并获得了很好的效果,为预测地震滑坡启动条件和堆积范围提供了良好途径。

    本文以宁夏固原市海口村北侧不稳定斜坡为对象,结合无人机获取高精度三维地形,采用PFC3D模拟在地震作用下斜坡失稳破坏和运动过程,并预测斜坡发生的峰值加速度、临界位移以及斜坡失稳后危害范围,为黄土地区边坡工程设计和防灾减灾工作服务。

    1 研究区地质环境

    研究区位于宁夏固原市彭阳县海口村北侧,处于黄土高原的西北区域,具备黄土地区特殊的自然地理和地质条件,也表现出黄土地区的典型地质环境问题。该区域跨鄂尔多斯西缘坳陷带和鄂尔多斯台坳两个二级构造单元,位于六盘山断裂和云雾山断裂之间(图 1),根据中国区域烈度表,该区域烈度为Ⅷ度,未来50 a内该区域具有发生6级以上地震的可能,超越概率10%对应的峰值加速度0.2~0.3 g。

    图 1 研究区位置图 Fig. 1 Location map of the study area

    在地貌上,该区域属于黄土丘陵沟壑区和残塬沟壑区的过渡带。地形复杂,地貌破碎,沟壑纵横,梁峁起伏,地貌形态主要表现为构造剥蚀的中低山、构造剥蚀的残山丘陵、剥蚀堆积的黄土地貌和堆积地貌,海拔高度1900~2000 m,山势起伏较缓,坡度为20°~30°。研究区域出露地层主要为第四系的堆积物和保山群(K1zd)山麓河流沉积物。上覆土层第四纪堆积物主要为风成黄土,颜色为褐黄、灰黄色,土质均匀,含植物根系,厚度为3~50 m。下部蓝灰、灰绿色泥岩、砂质泥岩,易风化(图 2)。该区属于半湿润、半干旱气候,具有雨热同期、干旱少雨等特点,年平均气温6.3 ℃。地区降雨量为400~500 mm。一年内降水分布极不均匀,4~9月份降水占全年降水量的84%,而7~9月占55%,分为明显的雨季和旱季。地区降雨常以阴雨和暴雨形式出现,也是该区境内滑坡、崩塌和泥石流形成的主要诱发因素之一。

    图 2 斜坡不同视角的照片 Fig. 2 Images of the slope from different perspectives a.研究区航拍影像;b.斜坡航拍影像;c.斜坡主视照片;d.斜坡右视照片

    通过对该区域的调查发现,该区域一处斜坡处于欠稳定状态,在强降雨或地震作用下有失稳的可能,该斜坡位于山体中部,高程为1812~1846 m。斜坡上覆黄土平均厚度5~6 m,长约85 m,宽约42 m,体积方量约为1.6×104 m3,下覆为白垩纪六盘山群(K1lp)的泥质粉砂岩,斜坡为岩质顺向坡,岩层产状SW170°∠30°。该区域曾经出现过剧烈的构造活动,受构造活动的影响斜坡后缘出现了大量的裂缝,航拍影像(图 2a图 2b)和斜坡正面和侧面照片(图 2c图 2d)均有裂缝存在(图 2b~图 2d),因此该斜坡未来极有可能在地震作用下发生失稳。

    2 离散元参数标定
    2.1 室内三轴实验

    在PFC3D模型中,没有与土体宏观性质直接对应的模拟参数,只能给出颗粒和黏结的微观参数。Potyondy et al.(2004), Cundall et al.(1979)提出了一种校准过程,通过采用单轴压缩试验结果,获得了PFC3D的合适的微观参数,因此本文采用室内三轴实验来获取与宏观土体所对应的颗粒和黏结的微观参数。室内三轴实验土样取自固原市彭阳县海口村模拟斜坡位置,土样基本物理参数(表 1)。室内三轴实验采用TFB-1型非饱和土应力-应变控制式三轴仪。由于斜坡土层平均厚度5~6 m,故围压设为0、50 kPa、100 kPa、150 kPa。为模拟在地震作用下黄土快速剪切破坏,故实验采用剪切速率为1 mm · min-1进行等压固结不排水剪(ICU),得到不同围压下应力-应变曲线(图 3)。通过应力-应变曲线可以看出黄土在应变为20%前呈硬化状态,但应变在20%左右应力有趋于平稳趋势(图 3)。通过广义剪应力q和平均主应力p绘制p-q破坏线(图 4),其直线方程为q=1.168 8p+83.534。

    表 1 重塑黄土的基本物理性质参数 Table 1 The basic physical property parameters of the reshaped loess

    图 3 剪切速率为1 mm · min-1的应力-应变曲线 Fig. 3 Stress-strain curve with shear rate of 1 mm · min-1

    图 4 剪切速率为1 mm · min-1p-q强度破坏线 Fig. 4 p-q strength failure curve with shear rate of 1 mm · min-1

    2.2 室内三轴实验标定
    2.2.1 颗粒流原理介绍

    PFC软件是以1979年由Cundall et al. (1979)提出离散单元法为理论的颗粒流模拟软件,可以模拟刚体颗粒的运动及互相作用,允许颗粒产生转动和有限的位移。PFC软件在计算中表现为颗粒与颗粒之间不平衡状态向平衡状态转化的动态过程,该动态过程求解采用显式差分算法。刚体颗粒具有一定质量,颗粒的运动规律遵循牛顿第二定律,颗粒在接触力和体力作用下产生新的位移,通过力-位移定律更新接触部分的接触力,在颗粒模拟过程中交替应用力-位移定律和牛顿第二定律直到颗粒达到新的平衡。离散单元法的核心是颗粒接触特征及接触本构,颗粒流的参数决定岩土体的宏观表现,因此选择正确的岩土参数,才能确保较为准确的模拟结果。

    为选择正确微观参数,不同的学者开展了大量的实验研究,Yang et al.(2006)对黏性颗粒材料采用平行黏结模型进行了二维单轴压缩实验,研究黏结模型中微观参数与宏观样本的弹性参数及单轴抗压强度的关系。Huang(1999)采用contact-bond接触模型,对黏性颗粒材料微观参数和宏观参数之间的相似关系进行了研究。徐小敏等(2010)通过室内三轴试验和PFC3D模拟的结果进行回归分析,基于线性接触模型建立了颗粒材料初始杨氏模量、初始泊松比宏观弹性常数与颗粒法向刚度、颗粒刚度比等细观弹性常数间的经验公式。李识博等(2013)采用PFC3D软件对陇西地区黄土三轴固结不排水剪切试验进行了数值模拟,并与室内试验结果进行对比,得到了宏观的应力-应变曲线的变化规律。以上是对离散元微观参数规律研究和微观参数与室内试验进行对比的研究成果。本次模拟基于前人经验基础上,利用PFC3D软件,对室内黄土三轴实验应力-应变曲线结果进行微观参数标定和选择。

    2.2.2 模型建立及模拟三轴原理

    PFC3D三轴实验模拟土样尺寸大小采用与室内实验相同,高度80 mm,半径40 mm,由6750个颗粒构成黄土试样模型(图 5a图 5b),该数量可以反映颗粒材料基本属性(周健等,2008)。散体颗粒半径主要为0.0011~0.0015 m,最大最小半径之比为1.36,颗粒密度2600 kg · m-3,颗粒接触采用线性平行黏结。颗粒外部采用墙体约束,侧面墙体约束法向刚度取颗粒法向刚度的1/10。围压通过伺服控制墙体速度来模拟三轴实验下围压,通过上下圆盘墙体施加速度来实现竖直方向应变控制加载方式(图 5b)。

    图 5 室内三轴实验与模拟三轴实验 Fig. 5 Indoor triaxial experiment and simulation triaxial experiment

    对于圆柱形试样计算竖向应变(εH)和径向应变(εR),计算公式如下:

    $ {\varepsilon _H} = \frac{{H - {H_0}}}{{\frac{1}{2}\left( {H + {H_0}} \right)}} $ (1)

    $ {\varepsilon _R} = \frac{{R - {R_0}}}{{\frac{1}{2}\left( {R + {R_0}} \right)}} $ (2)

    其中,H为试样当前竖直方向的高度;H0为试样竖直方向的原始高度;R为试样当前径向的尺寸;R0为试样径向的原始尺寸,一个时步长内约束墙体位移所引起的约束力增量为:

    $ \Delta {F^{\left( w \right)}} = k_n^{\left( w \right)}{N_c}{{\dot u}^{\left( w \right)}}\Delta t $ (3)

    其中,ΔF(w)为墙体位移运动引起的约束力增量;kn(w)为墙体与颗粒接触面的平均刚度;Nc墙体约束上接触面个数;$\mathit{\dot u} $(w)约束墙体的移动速度;Δt为一个时间步长。因此该区域的平均约束应力为:

    $ \Delta {\sigma ^{\left( w \right)}} = \frac{{k_n^{\left( w \right)}{N_c}{{\dot u}^{\left( w \right)}}\Delta t}}{A} $ (4)

    其中,A为约束区域的面积,为该区域约束应力。为了保持稳定约束应力的绝对值必须小于实测应力与目标应力差值的绝对值,在实际中,使用松弛因子(α)满足稳定要求变为下列表达式:

    $ \left| {\Delta {\sigma ^{\left( w \right)}}} \right| < \alpha \left| {\Delta \sigma } \right| = \alpha \left| {{\sigma ^{measured}} - {\alpha ^{required}}} \right| $ (5)

    墙体运动速度$\left({{{\mathit{\dot u}}^{\left(\mathit{w} \right)}}} \right) $用实测应力与目标应力差值来表示:

    $ \Delta {\sigma ^{\left( w \right)}} = G\left( {{\sigma ^{measured}} - {\alpha ^{required}}} \right) $ (6)

    其中,G为实现目标应力的增量参数,将式(4)带入式(5)中,其中得:kn(w)Nc,ΔtA都大于0,因此得:

    $ \frac{{k_n^{\left( w \right)}{N_c}\Delta t}}{A}\left| G {\Delta \sigma } \right| < \alpha \left| {\Delta \sigma } \right| $ (7)

    Gα大于0,消去|Δσ|得:

    $ \frac{{A\alpha }}{{k_n^{\left( w \right)}{N_c}\Delta t}} < G $ (8)

    模拟中α取常数,确定G来实现用速度实现应力控制,则G取:

    $ \frac{{A\alpha }}{{k_n^{\left( w \right)}{N_c}\Delta t}} = G $ (9)

    2.2.3 参数标定分析与选择

    根据室内三轴实验条件和斜坡特征,采用围压为0、50 kPa、100 kPa、150 kPa模拟三轴实验,对室内实验所得到应力-应变曲线进行标定。模拟所需要的微观参数主要包括颗粒的微观参数(Eckn/ksμ)和黏结的微观参数($\overline {{E_c}}, \overline {{\mathit{k}_\mathit{n}}} /\overline {{k_s}}, \overline {{\mathit{\sigma }_\mathit{b}}}, \overline {{\mathit{\tau }_\mathit{b}}}, \mathit{\overline \lambda } $)。Ec$\overline {{E_c}} $分别是颗粒和黏结的杨氏模量,kn/ks$\overline {{k_n}} /\overline {{k_s}} $分别为颗粒和黏结的法向刚度和切线刚度之比,$\overline {{\mathit{\sigma }_\mathit{b}}} $为平行法向黏结应力,$\overline {{\mathit{\tau }_\mathit{b}}} $为平行切向黏结应力,$\mathit{\overline \lambda } $为颗粒黏结半径系数。在对滑坡演化进行计算时,滑体不仅向下滑动消耗能量还有土体之间断裂和碰撞问题,黏滞阻尼是在PFC3D模型反映黏结断裂和碰撞的一个重要参数,用来表明当一个颗粒与另一个颗粒接触时的能量耗散。Giani(1992)在研究黏滞阻尼时,给出恢复系数与阻尼系数比的关系,对于岩土恢复系数为0.5对应的阻尼比为0.22。

    如何让宏观参数与微观参数对应,许多学者对此作了研究,Rodríguez-Sastre et al.(2008)通过单轴抗压强度计算弹性模量:E=429.56(UCS)0.9122(UCS为单轴抗压强度)。Potyondy et al.(2004)认为:kn=4REc$\overline {{E_\mathit{c}}} = \overline {{k_\mathit{n}}} \left({{R_1} + R2} \right)$(R为颗粒平均半径,R1R2分别为两个接触颗粒的半径)。Pande et al.(1990)认为由于材料的泊松比在0~5之间,因此kn/ks取值范围在2~3之间。一般弹性阶段曲线的斜率,kn/ks其值越大,直线斜率越大,颗粒法向刚度(kn)变大时,峰值应变变小。影响峰值强度的主要参数(kn),峰值强度随着颗粒法向刚度(kn)增大而增大。综上所述调节微观参数的总体思路:先计算弹性模量再算颗粒刚度,改变kn/ks$\overline {{\mathit{k}_\mathit{n}}} /\overline {{k_s}} $值调整曲线斜率,改变颗粒摩擦调整峰值应变,最后改变颗粒黏结和颗粒摩擦校正峰值应力。

    模拟过程曲线分为3个阶段,弹性阶段,宏观弹性模量主要由颗粒接触杨氏模量和颗粒黏结杨氏模量决定,且呈线性关系。第2阶段为塑性破坏,曲线的斜率受μ和黏结半径系数($\mathit{\overline \lambda } $)控制,曲线的斜率随μλ增大而增大。颗粒黏结力、μ的大小总体控制总强度的变化。第3阶段为破坏阶段,内摩擦角(φ)决定曲线的硬化程度,φ值越大曲线硬化越明显。从图 6可以看出,实验和模拟的应力-应变曲线都表现出应变硬化现象且曲线吻合度较好。取不同围压下的微观参数的平均值为最终的微观参数。

    图 6 不同围压下标定应力-应变曲线 Fig. 6 Calibration of stress-strain curves under different confining pressures

    3 地震作用三维斜坡分析
    3.1 三维模型建立
    3.1.1 三维模型生成

    PFC3D对研究崩塌、滑坡、泥石流等自然灾害后破坏程度及灾害影响范围有着明显的可视效果。本文通过无人机航拍获取该研究区域高精度的航拍影像,通过Pix4 d影像处理软件获得该区域的数字高程模型,再利用Arcgis将dem生成带有数据的三维等高线,通过Sketchup软件和Rhino软件处理得到PFC3D模型,对于上部黄土用颗粒代替,下部基岩统一简化为墙体(图 7)。

    图 7 三维模型建立流程 Fig. 7 Establishment process of three-dimensional model

    3.1.2 模型参数选择及地震波输入

    对于该斜坡取样黄土的数值模拟实验与室内三轴实验已经非常接近了,然而,在实际的滑坡模型中使用的颗粒的大小范围与模拟三轴实验粒径不同,这是由于数值三轴试验与大滑坡模型间存在尺度问题且受限于计算机的计算效率与存储能力,在保证宏观参数和微观参数的一致性前提下,粒径的更改使其他微观参数也在改变(Chang et al., 2009Lin et al., 2015),随着粒径尺度的改变,保证计算精度的要求,颗粒的刚度和黏结刚度也随其变化(Li et al., 2012邓益兵等,2017)(表 2表 3)。为了确定临界滑面墙体摩擦因数,通过为墙面赋予一个值使斜坡在重力作用下稳定系数大于1.0,通过fish语言将墙体摩擦系数以0.01不断减小迭代,以最大不平衡力与颗粒接触力的比值小于求解值为判定依据,确定墙体临界摩擦因数为0.38,因此采用滑面墙体摩擦因数为0.4。

    表 2 三轴模拟结果参数 Table 2 Parameters of triaxial simulation result

    表 3 模拟的参数 Table 3 Parameters of simulation

    根据中国地震烈度区划图和地震反应谱,未来50 a该区可能遭遇超越概率为10%的烈度值对应的峰值加速度为0.2~0.3 g,因此地震波采用最大峰值加速度0.28 g的阪神地震波(图 8b),将地震加速度积分为速度形式加载在墙体上,选取模型的底部边界作为地震施加源。Sato et al.(2007)认为地震诱发的滑坡大多数滑动方向与震时地壳主运动方向一致,因此本文输入地震波水平方向为炭山—嵩店大断裂的倾向北西30°左右和垂直方向为重力加速度方向(图 8c)。对斜坡9个位置的颗粒进行速度、位移、不平衡力监测,对斜坡4个位置应力分量进行监测记录(图 8a)。

    图 8 模型条件示意图 Fig. 8 Model conditions diagram

    3.2 三维地震作用下斜坡启动分析
    3.2.1 斜坡破坏应力过程分析

    分别以峰值加速度为0.01 g、0.05 g、0.1 g、0.12 g、0.15 g、0.17 g、0.2 g进行地震波加载,监测应力分量随时间变化情况。以监测x方向的应力(xx-stress)为例进行分析,从图 9中可以看到,当峰值加速度为0.01 g时,斜坡应力基本没有波动且处于平稳状态。当峰值加速度为0.05~0.12 g,斜坡应力在2.5~7.5 s内出现的波动,其波动幅度随着峰值加速度增加而增加,但斜坡依然处于稳定阶段,土体应力过程总结为:稳定-波动-稳定。当峰值加速度为0.12~0.15 g时,随着峰值加速度增加,斜坡应力分量在9 s后出现了微小波动,斜坡可能发生失稳现象,其具体情况还应结合位移来判断,土体应力过程总结为:稳定-波动-稳定-失稳。当峰值加速度大于0.15 g时,随着峰值加速度的增大,斜坡应力在失稳前波动幅度变化不大但频率增大,土体发生失稳的时间随着峰值加速度增加而提前,土体应力过程总结为:稳定-波动-失稳。

    图 9 不同峰值加速度下应力变化曲线 Fig. 9 Stress curves under different peak accelerations a.峰值加速度0.01 g;b.峰值加速度0.05 g;c.峰值加速度0.1 g;d.峰值加速度0.12 g;e.峰值加速度0.15 g;f.峰值加速度0.17g;g.峰值加速度0.2 g;h.峰值加速度0.25 g

    将监测点的应力分量式(10)通过式(11)、式(12)(徐秉业,1988)可将应力分量转化为平均正应力(p)、广义剪应力(q),得到监测点在地震作用下的应力路径图(图 10)。

    图 10 不同峰值加速度下应力路径 Fig. 10 Stress path under different peak accelerations

    $ {\sigma _{ij}} = \left[ {\begin{array}{*{20}{c}} {{\sigma _{xx}}}&{{\sigma _{xy}}}&{{\sigma _{xz}}}\\ {{\sigma _{yx}}}&{{\sigma _{yy}}}&{{\sigma _{yz}}}\\ {{\sigma _{zx}}}&{{\sigma _{zy}}}&{{\sigma _{zx}}} \end{array}} \right] $ (10)

    $ p = \frac{1}{3}\left( {{\sigma _{xx}} + {\sigma _{yy}} + {\sigma _{zz}}} \right) $ (11)

    $ q = \sqrt {\frac{1}{2}\left[ \begin{array}{l} {\left( {{\sigma _{xx}} - {\sigma _{yy}}} \right)^2} + {\left( {{\sigma _{xx}} - {\sigma _{yy}}} \right)^2}\\ + {\left( {{\sigma _{xx}} - {\sigma _{yy}}} \right)^2} + 6\left( {{\sigma _{xx}}^2 + {\sigma _{xz}}^2 + {\sigma _{yz}}^2} \right) \end{array} \right]} $ (12)

    图 10可以看出,当峰值加速度小于0.1 g时,随着峰值加速度的增加,坡体内的平均正应力(p)基本波动不大,广义剪应力(q)随着峰值加速度增大而增大。当峰值加速度大于0.15 g后,斜坡应力坡脚(1、2号应力点)破坏主要是由广义剪应力增大而引起的剪切破坏,斜坡坡顶(3号应力点)的破坏主要是由平均主应力较小而造成的拉张破坏。因此可以通过应力路径来判断斜坡在地震作用下坡体不同位置的破坏形式。

    3.2.2 斜坡启动过程分析

    为了获得该斜坡的临界位移,以0.03 g的增量加速度,将峰值加速度从0.03 g增大到0.45 g分别进行地震波加载模拟,绘制不同加速度峰值条件下加载时间与颗粒位移曲线图(图 11)。将不同峰值加速度的颗粒的位移曲线分为2类,斜坡失稳前的状态(图 11a)和斜坡失稳后的状态(图 11b)。从图 11a可以看出,位移放大幅度从坡底(2号颗粒)到坡腰(颗粒3,5号)再到坡顶(8号颗粒)先增大再减小,在坡腰处最大。随着峰值加速度的增大,斜坡位移放大幅度也在增大,但坡腰部(颗粒3,5号)与坡底部(2号颗粒)的位移比值不变(3倍)。从图 11a可以看出,当最大峰值加速度小于0.125 g时,在25 s地震波结束后,斜坡处于稳定状态。当峰值加速度大于0.135 g时,地震波结束后,斜坡底部开始发生局部的失稳,随后斜坡在28 s时开始失稳,因此可以确定在地震作用下该斜坡的失稳峰值加速度约为0.135 g。

    图 11 不同峰值加速度下颗粒位移曲线 Fig. 11 Particles displacement curves under different peak accelerations a.失稳前;b.失稳后

    为了准确判断地震作用下该斜坡失稳的临界位移,将监测的所有颗粒位移取平均值来表示斜坡在不同峰值加速度下的位移历时曲线(图 12)。斜坡进入加速变形阶段是滑坡发生的基础和前提,因此,斜坡的加速变形阶段对于滑坡的预测预报具有十分重要的意义(许强等,2008)。从图 12a可以看出,当地震峰值加速度小于0.135 g和最大永久位移小于等于9.5 cm时,斜坡在地震作用后都处于稳定状态。当峰值加速度等于0.135 g时,斜坡在25 s地震波结束后进入不稳定阶段,直到斜坡失稳,因此可以确定斜坡进入蠕变加速阶段的永久位移约为9.5 cm。

    图 12 不同峰值加速度下斜坡位移曲线 Fig. 12 The slope displacement curve under different peak accelerations a.失稳前;b.失稳后

    当峰值加速度大于0.125 g时,斜坡位移曲线可以分为两种类型(图 12b):(1)当峰值加速度在0.135~0.24 g之间时,斜坡颗粒运动具有明显的匀加速到加加速的过程,具有和土样蠕变曲线相似的规律。(2)当峰值加速度大于0.24 g时,斜坡位移曲线表现为直接进入加速破坏阶段,斜坡发生失稳。钱海涛等(2012)认为临界位移值实质上对应着系统势能的突变点,与滑带的几何特征、力学特性和滑体自重相关,与具体的地震波输入等无关,因此该斜坡在不同峰值加速度下,斜坡的临界位移是一个定值。许强等(2009)认为某个滑坡的等速变形阶段的速率为定值,通过位移除以v的方法,将位移-时间曲线改为纵横坐标具有相同量纲的曲线,所确定的改进切线角具有唯一性。本文按其方法确定改进切线角约89°时所对应的累计位移为该边坡的临界位移,从图 13中可以看到不同峰值加速度下切线角约为89°对应的临界位移为50 cm。综上所述,确定该斜坡加速变形阶段的起点位移约为9.5 cm,临界位移约为50 cm。

    图 13 不同峰值加速度下位移-切线角曲线 Fig. 13 The slope displacement-tangent angle curves under different peak accelerations

    3.3 地震作用下斜坡运动分析
    3.3.1 斜坡运动过程分析

    斜坡破坏后的位移和影响范围是对斜坡危险性评价的重要依据,根据中国地震烈度区划图和地震反应谱,该区域未来50 a具有发生6.0级地震的可能,当地震峰值加速大于0.125 g时,该斜坡失稳可能性很大,因此本文以地震峰值加速度为0.3 g、地震作用时长为25 s、墙体摩擦系数为0.1、模拟时间为45 s为条件来模拟该斜坡运动破坏后的运动过程。从图 14上看在5 s左右颗粒已经开始运动,在10 s左右颗粒运动最大速度达到12 m · s-1,颗粒整体平均速度8 m · s-1。在地震作用的12 s后颗粒速度开始降低,颗粒速度在40 s时逐渐趋于0,颗粒的运动过程受地形和地震作用共同影响。

    图 14 滑坡运动过程 Fig. 14 Landslide movement process

    3.3.2 不同摩擦力下滑体的堆积范围

    Reid et al.(2011)在水槽泥石流实验中发现不同底面环境对滑动距离影响较大,尤其是在底面湿润条件下坡体滑动产生的滑动距离比干燥时更大。通过模拟计算,当墙体摩擦因数为0.38时,斜坡处于临界稳定状态。Lin et al.(2015)利用PFC3D软件对台湾南部暴雨引发的大型滑坡进行堆积范围标定,标定在暴雨情况下对应的墙体滑面的摩擦因数为0.09。因此本文采用0.1、0.2、0.3、0.4摩擦因数来模拟不同下底面特征,预测在不同地面环境下滑体堆积范围和村庄危害范围。从图 15中看出随着墙体摩擦力增大,滑体堆积范围缩小。当摩擦因数为0.1、0.2时,颗粒堆积主要为两个位置,一部分主要堆积在坡脚到山口之间,另一部分主要堆积在山口前60 m内,堆积体在山口前的滑体危害到村民房屋处,南北长70 m,东西最大范围为80 m,危害范围内村庄房屋面积为1805~2502 m2,覆盖体积3602~5502 m3(表 4)。当摩擦系数为0.3、0.4时,靠近山的第一排村民房屋受少量滑体威胁,向山口东侧堆积危害范围40 m左右,南北20 m左右,危害范围内村庄房屋面积为1024~1357 m2,覆盖体积3602~5502 m3(表 4),不同摩擦因数下滑坡的危害几何范围和滑距见表 4。若地震作用后发生强降雨,在持续降雨条件下,坡脚到山口之间滑体易发生泥石流,将会造成更大的危害。

    表 4 不同摩擦因数下滑体对村庄危害信息 Table 4 Slipping mass to village hazard information under different friction factors

    图 15 不同摩擦因数下滑坡堆积范围和厚度图 Fig. 15 Range of accumulation and thickness of landslides with different friction factors

    4 结论

    本文利用离散元(PFC)方法,通过室内三轴实验得到力学参数及黄土p-q破坏线,并对室内三轴实验力学参数进行标定,再利用高精度航拍三维地形和标定好的三轴微观参数,以宁夏海口村黄土斜坡为研究对象,开展三维地震作用斜坡失稳破坏和运动过程数值模拟研究,得到如下结论:

    (1) 以不同峰值加速度对该斜坡进模拟得出其应力过程分为3类:当加速度为0.5~0.12 g,土体应力过程可总结为:稳定-波动-稳定。当加速度为0.12~0.15 g,土体应力过程总结为:稳定-波动-稳定-失稳。当峰值加速大于0.15 g土体应力过程总结为:稳定-波动-失稳。

    (2) 本文通过模拟分析不同峰值加速度下该斜坡稳定情况,确定该斜坡的失稳峰值加速度为0.135 g,斜坡位移大于9.5 cm时,斜坡将进入不稳定阶段,斜坡临界位移为50 cm。随着临界加速度增大,斜坡启动时间提前,随着时间提前所需的峰值加速增量也不断增大。

    (3) 斜坡滑体在地震作用下失稳后滑体平均速度为8 m · s-1,最大速度达到12 m · s-1。通过模拟不同摩擦系数,预测了在环境条件下斜坡滑动破坏后的危害范围,当摩擦因数为0.1、0.2时,颗粒堆积主要为两个位置,一部分主要堆积在坡脚到山口之间,另一部分主要堆积在山口前60 m内,堆积在山口前的滑体危害到村民居住区,南北长70 m,东西最大范围为80 m,危害范围中村庄房屋面积为1805~2502 m2,覆盖体积3602~5502 m3。当摩擦系数为0.3、0.4时,靠近山的第一排村民房屋受少量滑体威胁,向山口东侧的堆积体危害范围40 m左右,南北20 m左右,危害范围中村庄房屋面积为1024~1357 m2,覆盖体积3602~5502 m3

    (4) 本文通过室内三轴实验提供模拟所需的力学参数和无人机提供精准三维地形图为基础利用离散元在模拟大变形和运动过程的优势确定该斜坡的失稳临界峰值加速度和失稳临界位移并预测不同地面环境下斜坡失稳的堆积危害范围,为黄土地区斜坡工程的抗震设计及防震减灾工作提供一种新的可视化方法。

    参考文献
    Cai H J, Liu H C, Zhang Z Y. 1995. The catalog of Chinese landslide dam events[J]. Journal of Geological Hazards and Environmental Preservation, 6(4): 1~9.
    Chang K J, Alfredo T. 2009. Discrete element simulation of the Jiufengershan rock-and -soil avalanche triggered by the 1999 Chi-Chi earthquake, Taiwan[J]. Journal of Geophysical Research-Earth Surface, 114(F4).
    Cui F P, Hu R L, Yin Y P, et al. 2010. Discrete element analysis of collapsing and sliding response of slope triggered by time difference coupling effects of P and S seismic waves:Taking Tangjiashan landslide in Beichuan County for Example[J]. Chinese Journal of Rock Mechanics and Engineering, 29(2): 319~327.
    Cui T J, Ma Y D, Wang L G. 2016. Simulation and analysis of the earthquake stability of the tailing reservoir based on PFC3D[J]. Journal of Safety And Environment, 16(1): 95~99.
    Cundall P A, Strack O D L. 1979. A discrete numerical model for granular assemblies[J]. Géotechnique, 29(1): 47~65. DOI:10.1680/geot.1979.29.1.47
    Deng Y B, Yang Y P, Shi D D, et al. 2017. Refinement and application of variable particle-size methods in 3D discrete element modelling for large-scale problems[J]. Chinese Journal of Geotechnical Engineering, 32(7): 62~70.
    Dunning S A, Mitchell W A, Rosser N J, et al. 2007. The Hattian Bala rock avalanche and associated landslides triggered by the Kashmir Earthquake of 8 October 2005[J]. Engineering Geology, 93(3-4): 130~144. DOI:10.1016/j.enggeo.2007.07.003
    Giani G P. 1992. Rock slope stability analysis[M]. Rotterdam: A. A. Balkema..
    Hao M H, Xu Q, Yang X G, et al. 2015. Physical model tests on inverse grading of particles in high speed landslide debris[J]. Chinese Journal of Rock Mechanics and Engineering, 34(3): 472~479.
    Huang H Y. 1999. Discrete element modeling of tool-rock interaction[D]. Minneapolis: University of Minnesota.
    Huang R Q. 2009. Mechanism and geomechanical modes of landslide hazards triggered by Wenchuan 8.0 earthquake[J]. Chinese Journal of Rock Mechanics and Engineering, 28(6): 1239~1249.
    Huang T F. 1999. Landslides triggered by Ruli earthquake[D]. Taibei: Central University.
    Jibson R W, Keefer D K. 1993. Analysis of the seismic origin of landslide:Examples from the New Madrid seismic zone[J]. Geological Society of American Bulletin, 105(4): 521~536. DOI:10.1130/0016-7606(1993)105<0521:AOTSOO>2.3.CO;2
    Keefer D K. 1984. Landslides caused by earthquakes[J]. Geological Society of America Bulletin, 95(4): 406~421. DOI:10.1130/0016-7606(1984)95<406:LCBE>2.0.CO;2
    Kong J M, A F Y, Wu W P. 2009. Typical examples analysis the types of Wenchuan Earthquake landslide[J]. Journal of Soil and Water Conservation, 23(6): 66~70.
    Li S B, Wang C M, Ma J Q, et al. 2013. Microscopic changes of Longxi loess during triaxial shear process[J]. Rock and Soil Mechanics, 34(11): 3299~3305.
    Li W C, Li H J, Dai F C, et al. 2012. Discrete element modeling of a rainfall-induced flow slide[J]. Engineering Geology, 149-150(4): 22~34.
    Lin C H, Lin M L. 2015. Evolution of the large landslide induced by Typhoon Morakot:A case study in the Butangbunasi River, southern Taiwan using the discrete element method[J]. Engineering Geology, 197(30): 172~187.
    Liu Z Y, Ma C W, Miao T D, et al. 2000. Kinematic block model of long run-out prediction for high-speed landslides[J]. Chinese Journal of Rock Mechanics and Engineering, 19(6): 742~746.
    Lu C Y, Tang C L, Chan Y C, et al. 2014. Forecasting landslide hazard by the 3D discrete element method:A case study of the unstable slope in the Lushan hot spring district, central Taiwan[J]. Engineering Geology, 183: 14~30. DOI:10.1016/j.enggeo.2014.09.007
    Newmark N M. 1965. Effects of earthquakes on dams and embankment[J]. Geotechinque, 15(2): 139~160. DOI:10.1680/geot.1965.15.2.139
    Pande G N, Sorriso-Valvo M. 1990. Numerical methods in rock mechanics[M]. New York: John Wiley and Sons.
    Potyondy D O, Cundall P A. 2004. A bonded-particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences, 41(8): 1329~1364. DOI:10.1016/j.ijrmms.2004.09.011
    Qian H T, Zhang L F, Jing J Y, et al. 2012. Analyses of critical sliding displacement of landslide induced by strong earthquake in mountainous area[J]. Chinese Journal of Rock Mechanics and Engineering, 31(S1): 2619~2628.
    Reid M E, Iverson R M, Logan M, et al. 2011. Entrainment of bed sediment by debris flows:results from large-scale experiments[J]. Journal of Engineering Geology and Environment, B(42): 367~373.
    Rodríguez-Sastre M A, Gutiérrez-Claverol M, Torres-Alonso M. 2008. Relationship between cleavage orientation, uniaxial compressive strength and Young's modulus for slates in NW Spain[J]. Bulletin of Engineering Geology & the Environment, 67(2): 181~186.
    Sato H P, Hasegawa H, Fujiwara S, et al. 2007. Interpretation of landslide distribution triggered by the 2005 Northern Pakistan earthquake using SPOT 5 imagery[J]. Landslides, 4(2): 113~122. DOI:10.1007/s10346-006-0069-5
    Seed H B. 1968. Landslides during earthquakes due to soil liquefaction[J]. Journal of the Soil Mechanics and Foundations. Division, 94(SM5): 1053~1122.
    Sun X P, He S M, Gao C F, et al. 2017. Discrete element numerical analysis of Niujuangou landslide[J]. Journal of Lanzhou University(Natural Sciences, 53(1): 48~53.
    Terzaghi k, Peek R B. 1948. Soil mechanics in engineering practice[M]. New York: John Wile.
    Wang T, Wu S R, Shi J S, et al. 2013. Case study on rapid assessment of regional seismic landslide hazard based on simplified newmark displacement model:Wenchuan MS 8.0 Earthquake[J]. Journal of Engineering Geology, 24(1): 16~24.
    Wieczorek G F, Wilson R C, Harp E L. 1985. Map showing slope stability during earthquakes in San Mateo County, California[J]. Annals of the New York Academy of Sciences, 681(6431): 155~64.
    Wilson R C, Keefer D K. 1985. Predicting a real limits of earthquake induced landslides[C]//Ziony J I. Evaluating earthquake hazards in the Los Angeles Region-an Earth-science Perspective.[S. l]: [s. n.]: 317-345.
    Xia M, Ren G M, Guo Y S, et al. 2010. Flac3D numerical simulation of resuscitation mechanism of landslide earthquake loading[J]. Journal of Engineering Geology, 18(3): 305~311.
    Xu B Y. 1998. Plasticity mechanics[M]. Beijing: Higher Education Press.
    Xu Q, Tang G M, Xu K X, et al. 2008. Research on space-time evolution laws and early warning-prediction of landslides[J]. Chinese Journal of Rock Mechanics and Engineering, 27(6): 1104~1112.
    Xu Q, Zeng Y P, Qian J P, et al. 2009. Study on a improved tangential angle and the corresponding landslide prewarning criteria[J]. Geological Bulletin of China, 28(4): 501~505.
    Xu X M, Ling D S, Chen Y M, et al. 2010. Correlation of microscopic and macroscopic elastic constants of granular materials based on linear contact model[J]. Chinese Journal of Geotechnical Engineering, 32(7): 991~999.
    Xu X Y, Law K T, Selvadurai A P S. 1991. Examination of the Examination of the pseudo-static limit equilibrium method for dynamic stability analysis of slopes[C]//Proceedings of 44th Canadian Geotechnical Conference. Calgary: [s.n.].
    Yang B D, Jiao Y, Lei S T. 2006. A study on the effects of microparameters on macro properties for specimens created by particles[J]. International Journal for Computer-Aided Engineering and Software, 23(6): 607~631. DOI:10.1108/02644400610680333
    Zhang D, Wu Z H, LI J C, et al. 2013. An overview on earthquake-induced landslide research[J]. Journal of Geomechanics, 19(3): 225~226.
    Zhang L, Tang H M, Xiong C R, et al. 2011. Movement process simulation of high-speed long-distance Jiweishan Landslide with PFC3D[J]. Chinese Journal of Rock Mechanics and Engineering, 30(S1): 2601~2611.
    Zhao J W, Cui P, Hao M H. 2015. Comprehensive analyses of the initiation and entrainment processes of the 2000 Yigong catastrophic landslide in Tibet, China[J]. Landslide, 13(1): 39~54.
    Zhou J, Jia M C. 2008. Microscopicmodel test and numerical simulation for soil engineering[M]. Beijing: Science Press.
    柴贺军, 刘汉超, 张倬元. 1995. 中国滑坡堵江事件目录[J]. 地质灾害与环境保护, 6(4): 1~9.
    崔芳鹏, 胡瑞林, 殷跃平, 等. 2010. 纵横波时差耦合作用的斜坡崩滑效应离散元分析——以北川唐家山滑坡为例[J]. 岩石力学与工程学报, 29(2): 319~327.
    崔铁军, 马云东, 王来贵. 2016. 基于PFC3D的尾矿库地震稳定性模拟于分析[J]. 安全与环境学报, 16(1): 95~99.
    邓益兵, 杨彦聘, 史旦达, 等. 2017. 三维离散元大尺度模拟中变粒径方法的优化及其应用[J]. 岩土工程学报, 39(1): 62~70.
    郝明辉, 许强, 杨兴国, 等. 2015. 高速滑坡-碎屑流颗粒反序试验及其成因机制探讨[J]. 岩石力学与工程学报, 34(3): 472~479.
    黄润秋. 2009. 汶川8.0级地震触发崩滑灾害机制及其地质力学模式[J]. 岩石力学与工程学报, 28(6): 1239~1249. DOI:10.3321/j.issn:1000-6915.2009.06.021
    黄台丰. 1999.瑞里地震诱发之山崩[D].台北: 中央大学.
    孔纪名, 阿发友, 吴文平. 2009. 汶川地震滑坡类型及典型实例分析[J]. 水土保持学报, 23(6): 66~70. DOI:10.3321/j.issn:1009-2242.2009.06.016
    李识博, 王常明, 马建全, 等. 2013. 陇西黄土三轴剪切过程微观变化研究[J]. 岩土力学, 34(11): 3299~3305.
    刘忠玉, 马崇武, 苗天德, 等. 2000. 高速滑坡远程预测的块体运动模型[J]. 岩石力学与工程学报, 19(6): 742~746. DOI:10.3321/j.issn:1000-6915.2000.06.012
    钱海涛, 张力方, 兰景岩, 等. 2012. 强震作用下山区滑坡稳定临界位移分析[J]. 岩石力学与工程学报, 31(增1): 2619~2628.
    孙新坡, 何思明, 高成凤, 等. 2017. 牛圈沟滑坡离散元数值分析[J]. 兰州大学学报(自然科学版), 53(1): 48~53.
    王涛, 吴树仁, 石菊松, 等. 2013. 基于简化Newmark位移模型的区域地震滑坡危险性快速评估——以汶川MS 8.0级地震为例[J]. 工程地质学报, 21(1): 16~24. DOI:10.3969/j.issn.1004-9665.2013.01.003
    夏敏, 任光明, 郭亚莎, 等. 2010. 地震诱发滑坡复活机制的FLAC3D数值模拟分析[J]. 工程地质学报, 18(3): 305~311. DOI:10.3969/j.issn.1004-9665.2010.03.003
    徐秉业. 1988. 塑性力学[M]. 北京: 高等教育出版社.
    徐小敏, 凌道盛, 陈云敏, 等. 2010. 基于线性接触模型的颗粒材料细-宏观弹性常数相关关系研究[J]. 岩土工程学报, 32(7): 991~999.
    许强, 曾裕平, 钱江澎, 等. 2009. 一种改进的切线角及对应的滑坡预警判据[J]. 地质通报, 28(4): 501~505. DOI:10.3969/j.issn.1671-2552.2009.04.011
    许强, 唐明高, 徐开祥, 等. 2008. 滑坡时空演化规律及预警预报研究[J]. 岩石力学与工程学报, 27(6): 1104~1112. DOI:10.3321/j.issn:1000-6915.2008.06.003
    张铎, 吴中海, 李家存, 等. 2013. 国内外地震滑坡研究综述[J]. 地质力学学报, 19(3): 225~226. DOI:10.3969/j.issn.1006-6616.2013.03.001
    张龙, 唐辉明, 熊承仁, 等. 2012. 鸡尾山高速远程滑坡运动过程PFC3D模拟[J]. 岩石力学与工程学报, 31(增1): 2601~2611.
    周健, 贾敏才. 2008. 土工细观模型试验与数值模拟[M]. 北京: 科学出版社.
    Cai H J, Liu H C, Zhang Z Y. 1995. The catalog of Chinese landslide dam events[J]. Journal of Geological Hazards and Environmental Preservation, 6(4): 1~9.
    Chang K J, Alfredo T. 2009. Discrete element simulation of the Jiufengershan rock-and -soil avalanche triggered by the 1999 Chi-Chi earthquake, Taiwan[J]. Journal of Geophysical Research-Earth Surface, 114(F4).
    Cui F P, Hu R L, Yin Y P, et al. 2010. Discrete element analysis of collapsing and sliding response of slope triggered by time difference coupling effects of P and S seismic waves:Taking Tangjiashan landslide in Beichuan County for Example[J]. Chinese Journal of Rock Mechanics and Engineering, 29(2): 319~327.
    Cui T J, Ma Y D, Wang L G. 2016. Simulation and analysis of the earthquake stability of the tailing reservoir based on PFC3D[J]. Journal of Safety And Environment, 16(1): 95~99.
    Cundall P A, Strack O D L. 1979. A discrete numerical model for granular assemblies[J]. Géotechnique, 29(1): 47~65. DOI:10.1680/geot.1979.29.1.47
    Deng Y B, Yang Y P, Shi D D, et al. 2017. Refinement and application of variable particle-size methods in 3D discrete element modelling for large-scale problems[J]. Chinese Journal of Geotechnical Engineering, 32(7): 62~70.
    Dunning S A, Mitchell W A, Rosser N J, et al. 2007. The Hattian Bala rock avalanche and associated landslides triggered by the Kashmir Earthquake of 8 October 2005[J]. Engineering Geology, 93(3-4): 130~144. DOI:10.1016/j.enggeo.2007.07.003
    Giani G P. 1992. Rock slope stability analysis[M]. Rotterdam: A. A. Balkema..
    Hao M H, Xu Q, Yang X G, et al. 2015. Physical model tests on inverse grading of particles in high speed landslide debris[J]. Chinese Journal of Rock Mechanics and Engineering, 34(3): 472~479.
    Huang H Y. 1999. Discrete element modeling of tool-rock interaction[D]. Minneapolis: University of Minnesota.
    Huang R Q. 2009. Mechanism and geomechanical modes of landslide hazards triggered by Wenchuan 8.0 earthquake[J]. Chinese Journal of Rock Mechanics and Engineering, 28(6): 1239~1249.
    Huang T F. 1999. Landslides triggered by Ruli earthquake[D]. Taibei: Central University.
    Jibson R W, Keefer D K. 1993. Analysis of the seismic origin of landslide:Examples from the New Madrid seismic zone[J]. Geological Society of American Bulletin, 105(4): 521~536. DOI:10.1130/0016-7606(1993)105<0521:AOTSOO>2.3.CO;2
    Keefer D K. 1984. Landslides caused by earthquakes[J]. Geological Society of America Bulletin, 95(4): 406~421. DOI:10.1130/0016-7606(1984)95<406:LCBE>2.0.CO;2
    Kong J M, A F Y, Wu W P. 2009. Typical examples analysis the types of Wenchuan Earthquake landslide[J]. Journal of Soil and Water Conservation, 23(6): 66~70.
    Li S B, Wang C M, Ma J Q, et al. 2013. Microscopic changes of Longxi loess during triaxial shear process[J]. Rock and Soil Mechanics, 34(11): 3299~3305.
    Li W C, Li H J, Dai F C, et al. 2012. Discrete element modeling of a rainfall-induced flow slide[J]. Engineering Geology, 149-150(4): 22~34.
    Lin C H, Lin M L. 2015. Evolution of the large landslide induced by Typhoon Morakot:A case study in the Butangbunasi River, southern Taiwan using the discrete element method[J]. Engineering Geology, 197(30): 172~187.
    Liu Z Y, Ma C W, Miao T D, et al. 2000. Kinematic block model of long run-out prediction for high-speed landslides[J]. Chinese Journal of Rock Mechanics and Engineering, 19(6): 742~746.
    Lu C Y, Tang C L, Chan Y C, et al. 2014. Forecasting landslide hazard by the 3D discrete element method:A case study of the unstable slope in the Lushan hot spring district, central Taiwan[J]. Engineering Geology, 183: 14~30. DOI:10.1016/j.enggeo.2014.09.007
    Newmark N M. 1965. Effects of earthquakes on dams and embankment[J]. Geotechinque, 15(2): 139~160. DOI:10.1680/geot.1965.15.2.139
    Pande G N, Sorriso-Valvo M. 1990. Numerical methods in rock mechanics[M]. New York: John Wiley and Sons.
    Potyondy D O, Cundall P A. 2004. A bonded-particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences, 41(8): 1329~1364. DOI:10.1016/j.ijrmms.2004.09.011
    Qian H T, Zhang L F, Jing J Y, et al. 2012. Analyses of critical sliding displacement of landslide induced by strong earthquake in mountainous area[J]. Chinese Journal of Rock Mechanics and Engineering, 31(S1): 2619~2628.
    Reid M E, Iverson R M, Logan M, et al. 2011. Entrainment of bed sediment by debris flows:results from large-scale experiments[J]. Journal of Engineering Geology and Environment, B(42): 367~373.
    Rodríguez-Sastre M A, Gutiérrez-Claverol M, Torres-Alonso M. 2008. Relationship between cleavage orientation, uniaxial compressive strength and Young's modulus for slates in NW Spain[J]. Bulletin of Engineering Geology & the Environment, 67(2): 181~186.
    Sato H P, Hasegawa H, Fujiwara S, et al. 2007. Interpretation of landslide distribution triggered by the 2005 Northern Pakistan earthquake using SPOT 5 imagery[J]. Landslides, 4(2): 113~122. DOI:10.1007/s10346-006-0069-5
    Seed H B. 1968. Landslides during earthquakes due to soil liquefaction[J]. Journal of the Soil Mechanics and Foundations. Division, 94(SM5): 1053~1122.
    Sun X P, He S M, Gao C F, et al. 2017. Discrete element numerical analysis of Niujuangou landslide[J]. Journal of Lanzhou University(Natural Sciences, 53(1): 48~53.
    Terzaghi k, Peek R B. 1948. Soil mechanics in engineering practice[M]. New York: John Wile.
    Wang T, Wu S R, Shi J S, et al. 2013. Case study on rapid assessment of regional seismic landslide hazard based on simplified newmark displacement model:Wenchuan MS 8.0 Earthquake[J]. Journal of Engineering Geology, 24(1): 16~24.
    Wieczorek G F, Wilson R C, Harp E L. 1985. Map showing slope stability during earthquakes in San Mateo County, California[J]. Annals of the New York Academy of Sciences, 681(6431): 155~64.
    Wilson R C, Keefer D K. 1985. Predicting a real limits of earthquake induced landslides[C]//Ziony J I. Evaluating earthquake hazards in the Los Angeles Region-an Earth-science Perspective.[S. l]: [s. n.]: 317-345.
    Xia M, Ren G M, Guo Y S, et al. 2010. Flac3D numerical simulation of resuscitation mechanism of landslide earthquake loading[J]. Journal of Engineering Geology, 18(3): 305~311.
    Xu B Y. 1998. Plasticity mechanics[M]. Beijing: Higher Education Press.
    Xu Q, Tang G M, Xu K X, et al. 2008. Research on space-time evolution laws and early warning-prediction of landslides[J]. Chinese Journal of Rock Mechanics and Engineering, 27(6): 1104~1112.
    Xu Q, Zeng Y P, Qian J P, et al. 2009. Study on a improved tangential angle and the corresponding landslide prewarning criteria[J]. Geological Bulletin of China, 28(4): 501~505.
    Xu X M, Ling D S, Chen Y M, et al. 2010. Correlation of microscopic and macroscopic elastic constants of granular materials based on linear contact model[J]. Chinese Journal of Geotechnical Engineering, 32(7): 991~999.
    Xu X Y, Law K T, Selvadurai A P S. 1991. Examination of the Examination of the pseudo-static limit equilibrium method for dynamic stability analysis of slopes[C]//Proceedings of 44th Canadian Geotechnical Conference. Calgary: [s.n.].
    Yang B D, Jiao Y, Lei S T. 2006. A study on the effects of microparameters on macro properties for specimens created by particles[J]. International Journal for Computer-Aided Engineering and Software, 23(6): 607~631. DOI:10.1108/02644400610680333
    Zhang D, Wu Z H, LI J C, et al. 2013. An overview on earthquake-induced landslide research[J]. Journal of Geomechanics, 19(3): 225~226.
    Zhang L, Tang H M, Xiong C R, et al. 2011. Movement process simulation of high-speed long-distance Jiweishan Landslide with PFC3D[J]. Chinese Journal of Rock Mechanics and Engineering, 30(S1): 2601~2611.
    Zhao J W, Cui P, Hao M H. 2015. Comprehensive analyses of the initiation and entrainment processes of the 2000 Yigong catastrophic landslide in Tibet, China[J]. Landslide, 13(1): 39~54.
    Zhou J, Jia M C. 2008. Microscopicmodel test and numerical simulation for soil engineering[M]. Beijing: Science Press.
    柴贺军, 刘汉超, 张倬元. 1995. 中国滑坡堵江事件目录[J]. 地质灾害与环境保护, 6(4): 1~9.
    崔芳鹏, 胡瑞林, 殷跃平, 等. 2010. 纵横波时差耦合作用的斜坡崩滑效应离散元分析——以北川唐家山滑坡为例[J]. 岩石力学与工程学报, 29(2): 319~327.
    崔铁军, 马云东, 王来贵. 2016. 基于PFC3D的尾矿库地震稳定性模拟于分析[J]. 安全与环境学报, 16(1): 95~99.
    邓益兵, 杨彦聘, 史旦达, 等. 2017. 三维离散元大尺度模拟中变粒径方法的优化及其应用[J]. 岩土工程学报, 39(1): 62~70.
    郝明辉, 许强, 杨兴国, 等. 2015. 高速滑坡-碎屑流颗粒反序试验及其成因机制探讨[J]. 岩石力学与工程学报, 34(3): 472~479.
    黄润秋. 2009. 汶川8.0级地震触发崩滑灾害机制及其地质力学模式[J]. 岩石力学与工程学报, 28(6): 1239~1249. DOI:10.3321/j.issn:1000-6915.2009.06.021
    黄台丰. 1999.瑞里地震诱发之山崩[D].台北: 中央大学.
    孔纪名, 阿发友, 吴文平. 2009. 汶川地震滑坡类型及典型实例分析[J]. 水土保持学报, 23(6): 66~70. DOI:10.3321/j.issn:1009-2242.2009.06.016
    李识博, 王常明, 马建全, 等. 2013. 陇西黄土三轴剪切过程微观变化研究[J]. 岩土力学, 34(11): 3299~3305.
    刘忠玉, 马崇武, 苗天德, 等. 2000. 高速滑坡远程预测的块体运动模型[J]. 岩石力学与工程学报, 19(6): 742~746. DOI:10.3321/j.issn:1000-6915.2000.06.012
    钱海涛, 张力方, 兰景岩, 等. 2012. 强震作用下山区滑坡稳定临界位移分析[J]. 岩石力学与工程学报, 31(增1): 2619~2628.
    孙新坡, 何思明, 高成凤, 等. 2017. 牛圈沟滑坡离散元数值分析[J]. 兰州大学学报(自然科学版), 53(1): 48~53.
    王涛, 吴树仁, 石菊松, 等. 2013. 基于简化Newmark位移模型的区域地震滑坡危险性快速评估——以汶川MS 8.0级地震为例[J]. 工程地质学报, 21(1): 16~24. DOI:10.3969/j.issn.1004-9665.2013.01.003
    夏敏, 任光明, 郭亚莎, 等. 2010. 地震诱发滑坡复活机制的FLAC3D数值模拟分析[J]. 工程地质学报, 18(3): 305~311. DOI:10.3969/j.issn.1004-9665.2010.03.003
    徐秉业. 1988. 塑性力学[M]. 北京: 高等教育出版社.
    徐小敏, 凌道盛, 陈云敏, 等. 2010. 基于线性接触模型的颗粒材料细-宏观弹性常数相关关系研究[J]. 岩土工程学报, 32(7): 991~999.
    许强, 曾裕平, 钱江澎, 等. 2009. 一种改进的切线角及对应的滑坡预警判据[J]. 地质通报, 28(4): 501~505. DOI:10.3969/j.issn.1671-2552.2009.04.011
    许强, 唐明高, 徐开祥, 等. 2008. 滑坡时空演化规律及预警预报研究[J]. 岩石力学与工程学报, 27(6): 1104~1112. DOI:10.3321/j.issn:1000-6915.2008.06.003
    张铎, 吴中海, 李家存, 等. 2013. 国内外地震滑坡研究综述[J]. 地质力学学报, 19(3): 225~226. DOI:10.3969/j.issn.1006-6616.2013.03.001
    张龙, 唐辉明, 熊承仁, 等. 2012. 鸡尾山高速远程滑坡运动过程PFC3D模拟[J]. 岩石力学与工程学报, 31(增1): 2601~2611.
    周健, 贾敏才. 2008. 土工细观模型试验与数值模拟[M]. 北京: 科学出版社.