堆载诱发型土质滑坡运动特征物质点法模拟

吴方东，张 巍，史卜涛，施 斌，张 云，郑 帅

(南京大学地球科学与工程学院，江苏 南京 210046)

堆载诱发型土质滑坡运动特征物质点法模拟

吴方东，张 巍，史卜涛，施 斌，张 云，郑 帅

(南京大学地球科学与工程学院，江苏 南京 210046)

坡顶堆载是人类工程活动诱发滑坡的主因之一。物质点法(MPM)属于一种无网格数值计算方法，它能够有效模拟滑坡大变形全过程物质行为与运动特征。文章基于线性形函数离散方法、MUSL求解格式及Drucker-Prager屈服准则，建立了可用于滑坡全过程模拟的单套单相物质点模型；通过对比干燥铝棒堆积物模拟砂堆失稳过程的基准试验结果，对模型有效性进行了验证。对堆载诱发型土质滑坡典型工况进行了物质点法全过程模拟，获得了滑坡全过程中典型时刻坡体形态、塑性应变分布以及控制点滑速演化趋势。结果表明：算例堆载诱发型土质滑坡属推移型滑坡，具有渐进性破坏特征，可分为坡顶压缩、局部蠕滑、加速滑动与减速滑动等四个阶段。参数分析结果亦表明，堆载诱发型土质滑坡前缘物质运动特征量均与堆载量间存在强正相关性、而与土体黏聚力及内摩擦角存在强负相关性。统计29种典型工况，分别建立了峰值滑动加速度、最大滑速、最大滑距及坡体最大动能等运动特征量与堆载量、土体黏聚力及内摩擦角之间的线性回归方程，可用于堆载诱发型土质滑坡致灾行为预测。

土质滑坡；堆载；物质点法；运动参数

斜坡弃土堆载、坡顶兴建房屋或构筑物等人类工程活动是诱发滑坡的主因之一[1～3]。国内外学者已对堆载诱发型滑坡开展了一系列研究。试验方面，胡田飞等[4]基于室内土质边坡模型，对分级加载作用下滑坡演化过程进行了观察，发现堆载诱发型滑坡加速滑动阶段与剧滑启动的过渡时间极短，滑动面贯通后前缘坡体会迅速解体剪出。Li等[5]在室内对一个砂土边坡模型进行了加载破坏试验，并尝试使用离散元方法对斜坡破坏行为进行模拟。此外，Zhu等[6]还采用分布式光纤传感技术，表征土质边坡模型加载破坏试验过程中的应变场信息。理论方面，陈春利等[7]基于有限元模拟，分析了人工堆载黄土滑坡变形机制与破坏机理，结果表明，边坡堆载改变了坡体内的应力分布，致使坡体中、下部土体剪应力逐渐接近其抗剪强度，边坡稳定性随之降低直至失稳破坏。王洪兵等[8]采用有限差分法分析云南宝腾高速公路路堑边坡，发现坡体中部及坡脚堆载对边坡稳定有利，坡顶堆载对坡体稳定性最为不利。此外，董夫钱等[9]以浙江上三公路路堑1号滑坡为例，采用有限元方法分析发现堆载除直接改变边坡稳定性之外，还将间接影响到坡体内部地下水的渗流条件，并进一步降低坡体稳定性。

值得指出，既有研究多集中于堆载诱发型滑坡失稳机制的探讨，而研究斜坡从变形演化发展到失稳运动的全过程，能更全面地表征滑坡的地质力学行为及其动力演化机制[10～12]。数值模拟是研究滑坡的有力工具，但传统的极限平衡计算方法只能对斜坡体的稳定性进行分析[12]；有限元方法在处理大变形计算时则因网格畸变而出现数值困难[13～14]，这些方法都难以表征滑坡全过程中的物质运动特征。

物质点法属于一种无网格法，其在数值计算中不需要生成网格，而按任意分布的坐标点来构造插值函数实现离散的控制方程，可模拟各种复杂形态流场[14]，因而适用于滑坡大变形及超大变形模拟。该方法还可导入各种土体本构模型，考虑流固耦合[15～16]，而与SPH等其他无网格法相比较，物质点法具有更高的计算效率与数值稳定性[15]。目前国内外使用物质点法模拟土质滑坡运动全过程已取得一系列进展[17～21]，已被用于模拟降雨入渗诱发非饱和土滑坡等复杂物质运动问题[22]。本文采用物质点法模型模拟堆载诱发型土质滑坡的典型工况，重现了滑坡运动全过程，获取了典型时刻坡体的滑速、滑距、塑性应变等物质运动及变形特征。

1 物质点模型

1.1基本描述

物质点法采用拉格朗日质点系和欧拉网格双重描述的方法。如图1，离散的质点系携带介质的所有物质信息，计算网格仅用于动量方程的求解与空间导数求解。

图1 物质点法示意图Fig.1 Schematic diagram of the material point method

不考虑热交换，物质点法控制方程包括虚功方程、动量方程、几何方程、本构方程以及边界条件方程等，具体参见文[19]、[21]原理部分。

1.2物质点离散

采用具有线性形函数的背景网格进行离散，连续体的密度可近似表示为：

式中：np——质点总数；

mp——质点p的质量；

δ——Dirac Delta函数；

xi——空间坐标；

xip——质点p的坐标。

据文[13]可知，动量方程和给定面力边界条件的等效积分弱形式(即动量方程)为式(2)：

引入式(1)可将虚功方程式(2)简化为求和的形式：

式中，用带有下标p的量来表示物质点携带的相应物理量(下同)，同时引入假想边界层厚度h将虚功方程式左端最后一项转化为体积分。

用物质点法求解时，质点与背景网格不发生相对移动，故可以通过建立在背景网格结点上的有限元形函数NI(xi)实现质点和背景网格结点之间的信息映射，即：

其中，xip、uip、uip，j、δuip分别为质点p的坐标、位移、位移的导数以及虚位移，δuiI为结点的虚位移，用带有下标I的量来表示网格结点的变量。

式(5)～(7)结合式(3)，同时考虑结点虚位移δuiI在本质边界Γu上为零，且在其它结点处具有任意性，得到背景网格结点的运动方程：

∉

其中

为第I个结点在i方向上的动量；

为背景网格的质量矩阵；

为结点内力；

为结点外力。

当采用集中质量矩阵时，

式(9)简化为：

则背景网格结点的运动方程为：

∉

1.3求解格式

求解动量方程时，需要将某时刻质点的质量、动量等信息映射到背景网格，以计算背景网格结点的质量和动量等。动量方程采用显式时间积分求解。网格结点的速度变化量可以用于计算质点的应变增量，从而对质点进行应力更新。

根据应力更新时所采用的结点速度的不同，物质点法可分为USF、USL和MUSL等不同的求解格式。在时间步开始时进行应力更新的格式称为USF格式；在时间步结束时进行应力更新的格式称为USL格式；在时间步结束时的质点动量映射到背景网格后计算结点速度，并用此时的结点速度更新应力，称之为MUSL格式。本文模型采用MUSL格式，它具有较好的能量守恒特性[13]。

1.4本构模型

物质点法可分别采用单套、两套与三套物质点，并分别基于单相、两相与三相理论表征土体力学行为[22]，文[21]据此整理了单套单相、单套两相、单套三相、两套两相、两套三相与三套三相等6种土体物质点表征方法。为简化算法，提高计算效率，本文采用单套单相物质点表征模型，该模型适用于饱和土体，包括水饱和土与气饱和土(干燥土体)。

考虑到边坡初始应力场对坡体稳定性的影响，通常采用线弹性模型确定施加重力荷载后边坡的初始应力场。对于各向同性线弹性本构模型，焦曼应力率和变形率之间的关系为

边坡开始产生塑性变形直至失稳出现大变形阶段，则应采用弹塑性本构模型。Mohr-Coulomb强度准则虽广泛应用于岩土领域，但考虑其屈服面是由六个平面围成的锥形表面，容易导致数值求解困难。本文采用的Drucker-Prager屈服准则，其屈服面为一光滑圆锥面。屈服函数fs为：

式中：τ——等效剪应力；

σm——球应力；

qφ——摩擦系数；

kφ——纯剪切状态时的屈服应力。

qφ，kφ与材料的黏聚力c和摩擦角φ的关系由下式确定：

加号表示DP屈服面在π平面上内接MC屈服面，减号则表示外接。

2 模型验证

Bui等采用干燥铝棒堆积物进行试验，模拟了干燥砂堆的失稳过程[23]。铝棒被排列在一个长方体区域内，其长、高、宽分别为0.20 m、0.10 m、0.05 m。材料摩擦角约19.8°，泊松比0.3，平均体积模量0.7 MPa。

建立以上铝棒堆积物的单套单相物质点模型，示意如图2。

图2 物质点模型(80 000个物质点)Fig.2 MPM model (80 000 particles)

X=0 m,X=0.2 m取对称边界条件。Y=0 m,Y=0.05 m取对称边界条件。Z=0 m取固定边界条件；Z=0.1 m取自由边界条件。计算网格间距为1 mm，物质点间距为0.5 mm，物质点总数为80 000个。典型计算时刻铝棒堆积物失稳形态与塑性应变云图见图3。

图3 失稳形态与塑性应变云图Fig.3 Instability morphology and plastic strain contour map

图4为试验过程中铝棒堆积物的最终失稳形态。图5为分别采用SPH模型[23]和本文MPM模型所得最终失稳形态界面与试验失稳形态界面的比较结果。对比三者结果发现，堆积物滑动迹线与稳定后的坡面迹线吻合程度较高，验证了前文所建立的单套单相物质点模型的合理性与有效性。

图4 铝棒堆积物失稳最终形态(Bui，2008)Fig.4 Final form of the aluminum rods debris

图5 模拟与实验结果比较Fig.5 Comparison of the results of simulation and experiment

3 土质滑坡物质点法模拟

建立一坡高25 m、坡角45°、地基厚10 m的理想均匀土质边坡模型。采用规则的八结点六面体单元，网格间距1 m，质点间距0.5 m，每个计算网格内在水平和竖直方向各有2个物质点，垂直纸面方向有1个物质点。

根据背景网格三维坐标值确定每个物质点的几何位置。物质点的质量为其所在网格区域所代表土体单元的质量。如图6所示，边坡离散体设置时，由于其坡面线切割六面体背景网格，故而在坡面附近的物质点将偏离背景网格的体心，需按照三角形区域的形心布置，坡面物质点的质量也是三角形区域代表的土体质量。

根据以上方法所建立的物质点离散模型如图7所示。模型底部采用固定边界，两侧施加对称边界条件。此外，通过约束三维模型的y方向位移来模拟平面应变状态，即模型的前后面亦施加对称边界条件。在坡肩、坡中与坡趾等位置坡面上分别设置A、B、C三个计算监测点。

土体的主要物理力学参数取值：密度ρ=2 000 kg/m3，弹性模量E=70 MPa，泊松比u=0.3，黏聚力c=20 kPa，内摩擦角φ=25°。

图6 物质点及网格布置示意图Fig.6 Schematic diagram of the material point and grids layout

图7 边坡物质点模型(10 975个物质点)Fig.7 Slope MPM model (10 975 particles)

4 滑坡物质运动特征

4.1荷载设置

在坡顶设置高度为5 m的堆载体,分布范围自坡肩A点向内20 m。为了消除参数分析过程中由于堆载体体积变化而对坡体运动可能产生的影响，文章通过改变堆载体的重力加速度的方式来模拟不同的坡顶荷载。另外，为了减小计算过程中的数值震荡，施加重力荷载时均采用线性加载的方式。此外，考虑到边坡初始应力场对坡体稳定性的影响，0～3 s时刻间采用线弹性本构模型，3 s后采用Drucker-Prager弹塑性本构[19]，如图8所示。图9为重力施加完成后边坡内部初始应力场分布。

图8 加载模式Fig.8 Loading mode

图9 初始应力场Fig.9 Initial stress field

4.2滑坡全过程分析

基于图9的初始应力场进行加载，并对堆载体赋3倍的重力加速度值,即相当于294 kPa的坡顶荷载，计算时间为15 s。图10绘出了滑坡各典型时刻坡体位移形态与塑性应变区演化过程。

图11绘出了坡体失稳滑动前后A、B、C三个测点的水平滑速曲线。由图可见，从第3秒坡体失稳后3个监测点均出现加速段，随后在第6～7秒内依次转变为减速段，计算时刻10 s左右滑动停止，其中，坡脚C点在起滑后迅速达到0.45 m/s滑速，随后基本以匀速滑动至4 s时刻，才继续加速滑动。3个测点中，坡中B点最大滑速最大，达到5.8 m/s；坡肩A点最大滑速达到5 m/s；坡角C点最大滑速最小，仅1.5 m/s。分别采用式(20)、(21)指数衰减函数对A点滑速曲线的上升段和下降段进行分段拟合；而B、C两点滑速曲线则可由式(22)、(23)表示的Gauss函数进行全段拟合，拟合函数的表达式具体如下：

图10 滑坡位移形态与塑性应变区演化Fig.10 Landslide displacement morphology and plastic strain zone evolution

图11 滑速演化趋势Fig.11 Slide speed evolution trends

结合图10、11分析可知，本算例堆载诱发型土质滑坡全过程可划分为以下4个阶段:

(1)Ⅰ坡顶压缩

如图10(b)所示，坡顶土体在堆载的作用下发生压缩变形，堆载体两侧边缘以及坡脚产生应力集中，开始出现塑性应变。

(2)Ⅱ局部蠕滑

如图10(c)、(d)所示，荷载的持续作用使得弧形塑性区(潜在滑带)首先从坡脚开始形成并逐渐向坡体内部发育。滑体前缘开始鼓胀隆起，坡脚处土体的剪应力超过了抗剪强度并开始向前蠕滑。

(3)Ⅲ加速滑动

如图10(e)所示，持续的蠕滑使得弧形塑性区完全贯通，滑体前缘进一步鼓胀形成滑舌，最终沿滑带从坡脚处剪出产生整体滑动，滑速迅速增大进入巨滑阶段，如图11所示。

(4)Ⅳ减速滑动

如图11所示，从6～7 s阶段，坡体先后由加速滑动阶段转入减速滑动阶段，滑坡体继续前移，滑舌下方出现高塑性应变区，滑坡前缘最大滑距达到17 m，如图10(f)所示。

可见，本算例斜坡在荷载作用初期以坡顶压缩变形为主，随坡体内塑性应变的累积，逐渐向临空面发展，沿弧形滑带发生整体滑动，表现出推移式滑坡的变形破坏特点，即以压剪破坏为主[24]，这与文[4]所观察到的试验现象相吻合。边坡破坏一般表现为滑带形成、发育直至整体贯通后滑动的渐进破坏过程[24-25]，算例堆载诱发型滑坡除具有这种阶段性演化特征之外，还表现出剧滑前的蠕动变形时间短，滑带形成后迅速解体破坏的特点[4]。

5 参数分析

高速运动的滑体携带巨大的动能，具有极高的冲击力。因此，运动过程中动能大小可以用来评估滑坡的运动灾害性[26]。通过改变不同坡顶荷载和土体物理力学参数，分析了滑坡动能的变化情况，结果如图12所示。图12(a)为考虑不同坡顶荷载影响，且保持黏聚力为20 kPa，内摩擦角为25°的各工况；图12(b)为考虑不同土体内摩擦角影响，且保持坡顶堆载为294 kPa，土体黏聚力20 kPa的各工况；图12(c)为考虑不同土体黏聚力影响，且保持坡顶堆载为294 kPa，土体内摩擦角为25°的各工况。

从图12可见，滑坡动能峰值随坡顶堆载增大而接近线性增加，随内摩擦角以及黏聚力的增大而接近线性的减小，且动能变化曲线基本都呈现出上升段与下降段基本对称的特征。图12(a)表明，随坡顶荷载的增加，系统动能峰值随之增加，而达到峰值所需要的时间略有减小，坡体达到变形稳定的时间随之缩短。图12(b)表明，内摩擦角的增大使得系统动能峰值明显减小。在高摩擦角情况下，坡体滑动达到动能峰值的时间也有较明显的推迟，滑坡体需要更长的时间达到变形稳定阶段。图12(c)表明，黏聚力的变化只对系统动能峰值有影响，即使在高黏聚力状态下，滑坡达到变形稳定所需要的时间也基本不变。需指出， Drucker-Prager屈服准则未考虑材料的软化特性，因此在时间尺度上，计算时间与实际滑坡时间尺度的对应关系需另文研究。

总体而言，坡顶荷载的增加与土体强度的降低，将使得滑坡动能峰值显著增加。同时，高荷载、低摩擦角还将缩短坡体动能达到峰值的时间，使得滑坡的突发性变得更加显著。

滑速、滑距以及动能等滑坡运动特征参数对于评估滑坡致灾范围和致灾程度具有重要意义。为了进一步研究上覆荷载、黏聚力和内摩擦角对滑速、滑距以及坡体动能等参数的影响程度，建立了29种不同堆载量、土体黏聚力与内摩擦角组合工况下边坡的物质点模型，对滑坡前缘C点的峰值加速度、滑速、最大滑距以及坡体最大动能进行了统计(表1)。

图12 不同条件下系统动能变化曲线Fig.12 System kinetic energy variation curves under different conditions

黏聚力/kPa内摩擦角/(°)总荷载/kPa滑速峰值/(m·s-1)加速度峰值/(m·s-2)坡体滑距/m动能峰值/kJ202519613505212153820252451570681525592025294186082173596202534321110721463415252942551102349291625294236110224654172529421510421436718252942011002041041925294192094193841212529418208717335322252941760831631122325294166073152878242529415406414265425252941520581324202625294149058122180272529414405211196928252941310451017622925294120038915633025294099032813712020294430167286814202129439612926602620222943031322353052023294261109214677202429422110419411220262941600721631052027294150062142607202829412005012215520292940980381117902030294072020101382

对以上计算结果进行回归统计分析，并依次用函数式(24)～(27)建立C点速度峰值v、加速度峰值a、滑距s以及最大动能k与黏聚力c、内摩擦角φ和坡顶堆载p三者的线性表达式，四式的相关性系数分别达到0.92、0.98、0.99和0.99，表明三个影响因素与滑速、加速度和滑距之间存在强线性相关性。

观察上式可以发现，黏聚力和内摩擦角与以上4个滑坡运动特征量呈负线性相关性，即二者减小都会使得峰值加速度、最大滑速、滑距及动能显著增加，且内摩擦角的影响程度远甚于黏聚力。而堆载量则与以上4个滑坡运动特征量呈正相关性，其增大将导致峰值加速度、最大滑速、最大滑距及最大动能增大。

6 结论

基于物质点法基本理论，采用数值算例，模拟了堆载诱发型土质滑坡运动全过程，分析了物质运动特征，得出以下结论：

(1)物质点法是一种能够有效模拟滑坡大变形的无网格数值方法。该方法可对滑坡运动全过程进行重现，可提取出任意计算时刻坡体的滑速、滑距、塑性应变等物质运动及变形特征。

(2)算例堆载诱发型土质滑坡属推移式滑坡，其力学机制为压剪破坏，以沿弧形滑带的整体滑动为主，并具有渐进式破坏特征，可划分为坡顶压缩、局部蠕滑、加速滑动与减速滑动等4个阶段。

(3)堆载诱发型土质滑坡运动特征与上覆堆载量和土体物理力学性质密切相关。上覆堆载量越大，土体内摩擦角越小，滑坡越趋向于在短时间内达到较高滑速，突发性更为明显，体现出“突然滑动，迅速稳定”的运动特征，而土体黏聚力对此影响则相对较小。

(4)统计了29种典型工况下不同堆载量与土体物理力学指标数据，建立了滑坡峰值加速度、最大滑速、最大滑距与最大动能的线性回归方程，可用于此类滑坡致灾行为预测。

[1] 张茂省, 李同录. 黄土滑坡诱发因素及其形成机理研究[J]. 工程地质学报, 2011, 19(4): 530-540. [ZHANG M S, LI T L. Triggering factors and forming mechanism of loess landslides[J]. Journal of Engineering Geology, 2011, 19(4): 530-540. (in Chinese)]

[3] 刘传正. 南昆铁路八渡滑坡成因机理新认识[J]. 水文地质工程地质, 2007, 34(5): 1-5. [LIU C Z. A new discussion about genesis and failure mechanism of Badu landslides in Nanning-Kunming railway[J]. Hydrogeology amp; Engineering Geology, 2007,34(5): 1-5. (in Chinese)]

[4] 胡田飞. 堆载诱发型滑坡变形演化机理的模型试验研究[J]. 石家庄铁道大学学报(自然科学版), 2016,29(1):38-42. [HU T F. Model test study of deformation mechanism of landslide induced by loading[J].Journal of Shijiazhuang Railway Institute,2016, 29(1):38-42.(in Chinese)]

[5] LI N, CHENG Y M. Laboratory and 3-D distinct element analysis of the failure mechanism of a slope under external surcharge[J]. Natural Hazards and Earth System Sciences, 2015, 15:35-43.

[6] ZHU H H, SHI B, ZHANG J,etal. Distributed fiber optic monitoring and stability[J]. Journal of Mountain Science, 2014, 11(4): 979-989.

[7] 陈春利, 李同录, 贺凯, 等. 人工堆载诱发黄土滑坡失稳机制分析[J]. 中国地质灾害与防治学报, 2014, 25(1): 1-5. [CHEN C L, LI T L, HE K,etal. Mechanism of loess landslide induced by loading[J]. The Chinese Journal of Geological Hazard and Control, 2014,25(1): 1-5. (in Chinese)]

[8] 王洪兵, 姚磊华, 文海. 基于FLAC3D进行堆载作用下边坡稳定分析[J]. 勘察科学技术, 2012 (4): 5-10. [WANG H B, YAO L H, WEN H. Stability analysis of slope under preloading based on FLAC3D[J].Site Investigation Science and Technology, 2012 (4): 5-10. (in Chinese)]

[9] 董夫钱, 缪志顺, 吕庆, 等. 公路堆载诱发型滑坡稳定性分析[J]. 岩石力学与工程学报, 2004 (增刊): 4517-4520. [DONG F Q, MIAO Z S, LV Q,etal. Mechanism analysis of landslide induced by accumulation loading of highway[J]. Chinese Journal of Rock Mechanics and Engineering, 2004 (Sup1): 4517-4520. (in Chinese)]

[10] 张伟锋, 黄润秋, 裴向军. 大光包滑坡运动特征及其过程分析[J]. 工程地质学报, 2015, 23(5): 866-885. [ZHANG W F, HUANG R Q, PEI X J. Analysis on kinematics characteristics and movement process of Daguangbao landslide[J].Journal of Engineering Geology, 2015, 23(5): 866-885. (in Chinese)]

[11] 罗锋, 柴波, 方恒, 等. 顺层岩质滑坡运动过程数值模拟研究[J]. 水文地质工程地质, 2016, 43(3): 124-130. [LUO F,CHAI B,FANG H,etal. Numerical simulation of the bedding rock landslide motion process[J]. Hydrogeology amp; Engineering Geology, 2016, 43(3): 124-130. (in Chinese)]

[12] 殷坤龙, 姜清辉, 汪洋. 滑坡运动过程仿真分析[J]. 地球科学(中国地质大学学报), 2002,27(5): 632-636. [YIN K L, JIANG Q H, WANG Y. Simulation of landslide movement process by discontinuous deformation analysis[J]. Earth Science(Journal of China University of Geosciences), 2002, 27(5):632-636. (in Chinese)]

[13] 张雄, 廉艳平, 刘岩,等. 物质点法[M]. 北京：清华大学出版社, 2013. [ZHANG X, LIAN Y P, LIU Y,etal. Material point method[M]. Beijing: Tsinghua University Press, 2013. (in Chinese)]

[14] 张雄, 刘岩, 马上. 无网格法的理论及应用[J]. 力学进展, 2009,39(1): 1-36. [ZHANG X, LIU Y, MA S. Meshfree methods and their applications[J]. Advances in Mechanics, 2009,39(1) : 1-36. (in Chinese)]

[15] York II A R, Sulsky D, Schreyer H L. The material point method for simulation of thin membranes[J]. International Journal for Numerical Methods in Engineering,1999,44: 1429-1456.

[16] York II A R, Sulsky D, Schreyer H L. Fluid membrane interaction based on the material point method[J]. International Journal for Numerical Methods in Engineering, 2000,48: 901-924.

[17] Beuth L, Więcknoski Z, Vermeer P. Solution of quasi-static large-strain problem by the material point method[J]. International Journal for Numerical and Analytical method in Geomechanics, 2011,35(13): 1451-1465.

[18] Andersen S, Andersen L. Modelling of landslides with the material-point method[J]. Computational Geosciences, 2010, 14(1): 137-147.

[19] 史卜涛, 张云, 张巍. 边坡稳定性分析的物质点强度折减法[J]. 岩土工程学报, 2016, 38(9): 1678-1684. [SHI B T, ZHANG Y, ZHANG W. Strength reduction material point method for slope stability[J]. Chinese Journal of Geotechnical Engineering, 2016,38(9): 1678-1684. (in Chinese)]

[20] 孙玉进, 宋二祥. 大位移滑坡形态的物质点法模拟[J]. 岩土工程学报, 2015, 37(7): 1218-1225. [SUN Y J, SONG E X. Simulation of large-displacement landslide by material point method[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(7): 1218-1325. (in Chinese)]

[21] 张巍,史卜涛,施斌,等. 土质滑坡运动全过程的物质点法模拟及其应用[J]. 工程地质学报, 2017, 25(3)：815-823. [ZHANG W, SHI B T, SHI B,etal. Run-out process simulation of soil landslide using material point mehtod and its application[J]. Journal of Engineering Geology, 2017，25(3)：815-823. (in Chinese)]

[22] Yerro A, Alonso E, Pinyol N. The material point method for unsaturated soils[J]. Géotechnique,2015, 61(3): 201-217.

[23] BUI H H, FUKAGAWA R, SAKO K,etal. Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geomaterial using elastic-plastic soil constitutive model[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2008, 32(12):1537-1570.

[24] 王宝亮, 彭盛恩, 陈洪凯. 推移式滑坡形成机制的力学演绎[J]. 地质灾害与环境保护, 2010 (2): 74-77. [WANG B L,PENG S E,CHEN H K. Mechanical deduction of formation mechanism for the thrust load caused landslide[J].Journal of Geological Hazards and Environment Preservation, 2010 (2): 74-77.(in Chinese)]

[25] 卢应发, 黄学斌, 刘德富. 推移式滑坡渐进破坏机制及稳定性分析[J]. 岩石力学与工程学报, 2016, 35(2): 333-345. [LU Y F,HUANG X B,LIU D F. Mechanism and stability analyses of progressive failure of thrust-type landslides[J].Chinese Journal of Rock Mechanics and Engineering, 2016, 35(2): 333-345. (in Chinese)]

[26] 王效宁, 杜永廉, 巫锡勇. 岩质顺层滑坡的运动灾害性模拟试验研究[J]. 水文地质工程地质, 1988, 15(6): 1-5. [WANG X N, DU Y L, WU X Y. Simulation and experimental study on movement severe of bedding rock landslide[J].Hydrogeology amp; Engineering Geology, 1988, 15(6):1-5. (in Chinese)]

责任编辑

：汪美华

Run-outcharacteristicsimulationofasurcharge-inducedsoillandslideusingthematerialpointmethod

WU Fangdong, ZHANG Wei, SHI Butao, SHI Bin, ZHANG Yun, ZHENG Shuai

(SchoolofEarthSciencesandEngineering,NanjingUniversity,Nanjing,Jiangsu210046,China)

Spoil ground surcharge is one of the main factors inducing landslides in hilly grounds triggered by human engineering activities. The material point method (MPM) belongs to one of the meshless numerical analysis method, which can effectively simulate the material behavior and the run-out characteristics of landslides. Based on the discrete method using the linear shape function, the MUSL solving format and the Drucker-Prager yield criterion, we develop a single-layer and single-phase MPM model to simulate the run-out process of landslides. In comparison with the benchmark experiment results of the process of a sand pile losing its stability and simulated with accumulated dry aluminum bars, the proposed MPM model is verified. The run-out process of the typical scenario of the surcharged-induced soil landslide is simulated using MPM. The related results during the representative moments of the run-out process, including the morphology of slope mass, the plastic strain distribution and sliding velocity evolution trends of reference points, are obtained. The results of the numerical example show that the surcharge-induced soil landslide belongs to the thrust-type landslide, and is of the progressive failure. The whole run-out process can be divided into 4 stages, namely, the slope crest compression, local creep sliding, accelerated sliding and decelerated sliding. Concerning the sliding front, the results of parametric analysis also show that a strong positive correlation property exists between both of the kinetic representative parameters and the surcharge amount, and a strong negative correlation property exits between the kinetic representative parameters and the cohesion or internal friction angle of the soil. It is noted that the kinetic representative parameters include the maximal sliding acceleration, velocity, distance and kinetic energy of the slope mass. Based on 29 typical scenarios, the linear regression equations of all of the above 4 kinetic representative parameters, denoted by the surcharge amount, cohesion and internal friction angel of the soil, are established to predict the disastrous behavior of the surcharge-induced soil landslides.

soil landslide; surcharge; material point method; motion parameters

10.16030/j.cnki.issn.1000-3665.2017.06.19

P642.22

A

1000-3665(2017)06-0126-09

2016-12-27;

2017-04-05

国家自然科学基金重点项目(41230636)；江苏省自然科学基金项目(BK20160366)；苏州市科技计划项目(SYG201613)

吴方东(1993-)，男，硕士研究生，主要从事地质工程大变形数值模拟研究。E-mail: njuwfd@163.com

张巍(1974-)，男，副教授，从事工程地质数值方法研究。E-mail: wzhang@nju.edu.cn