考虑界面摩擦效应的土质滑坡物质点法模拟

Material Point Method Simulation of Soil Landslide Considering Interface Friction Effect

摘要

大型滑坡在高速远程滑动破坏过程中，滑带剧烈摩擦和生热，导致碎屑流和基岩之间存在界面摩擦效应，影响滑坡的最终堆积形态。该文以深圳杨宝地滑坡为例，采用物质点法建立二维滑坡模型，通过多重背景网格的方法引入接触算法，对滑坡的启动、滑动和堆积等过程进行了数值模拟。结果表明：滑坡启动后沿基岩呈流态状进行长距离运动，受基岩的地形地貌特征影响，在部分缓倾区域形成局部堆积，最终堆积形态与粒子有限单元法等方法的计算结果较为接近，初步验证了该文方法的有效性。改变碎屑流和基岩之间的摩擦系数进行计算，发现界面摩擦作用会对碎屑流的运动规律产生显著影响，最大滑动速度和最远滑动距离随着摩擦系数增大而减小，但碎屑流靠近基岩区域的土体变形受摩擦作用更为明显，等效塑性应变反而随着摩擦系数增大而增大，同时会形成大量不连续的塑性变形突触，影响碎屑流的最终堆积形态。研究成果可为分析界面摩擦效应作用下自然滑坡破坏机理提供有效的技术手段。

Abstract

During the high-speed and long-runout sliding failure process of large landslides, severe friction and heat generation in the sliding zone induce an interface friction effect between the debris flow and the bedrock, influencing the final accumulation morphology of the landslide. In this paper, a two-dimensional landslide model of the Yangbaodi landslide in Shenzhen was established using the material point method. A contact algorithm was introduced through a multi-background grid method to conduct numerical simulations on the initiation, sliding, and accumulation processes of the landslide. The results show that after initiation, the landslide moves over a long distance in a flow-like state along the bedrock. Influenced by the topographic and geomorphic characteristics of the bedrock, local accumulation is formed in some gently inclined areas. The final accumulation morphology is relatively close to the calculation results of methods such as the particle finite element method, which preliminarily verifies the effectiveness of the method in this paper. Calculations changing the friction coefficient between the debris flow and the bedrock reveal that the interface friction effect significantly influences the movement behavior of the debris flow. The maximum sliding velocity and farthest sliding distance decrease as the friction coefficient increases. However, the soil deformation of the debris flow near the bedrock area is more obviously affected by friction, and the equivalent plastic strain conversely increases as the friction coefficient increases. At the same time, a large number of discontinuous plastic deformation protrusions are formed, influencing the final accumulation morphology of the debris flow. The research results can provide an effective technical means for analyzing the failure mechanism of natural landslides under the interface friction effect.

Graphical abstract

关键词

滑坡 / 物质点法 / 界面摩擦效应 / 大变形破坏机理

Key words

landslide / material point method / interface friction effect / large deformation failure mechanism

0 引言

大型滑坡在高速远程破坏过程中，碎屑流和基岩之间存在界面摩擦作用，导致滑坡具有复杂的运动破坏机理^［1-7］。滑坡灾害通常会造成严重的生命财产损失，揭示界面摩擦作用对碎屑流运动规律和最终堆积形态的影响，有助于阐明滑坡破坏机理和提升防灾减灾水平。

许强等^［8］总结分析了大型崩滑灾害动力特征，认为滑坡通常以碎屑流、碎裂‒滑移崩塌、碎裂‒俯冲滑移等多种形式发生高速运动。滑坡属于典型的大变形破坏问题，采用传统有限单元法（Finite Element Method， FEM）模拟该过程存在网格畸变或计算不收敛等局限，因此有学者建议采用新型数值计算方法分析滑坡致灾机理。例如，张卫杰等^［9］采用滑粒子流体动力学（Smoothed Particle Hydrodynamics，SPH），模拟了四川省盐源县玻璃村“7·19”滑坡，分析了不同强度折减系数对滑坡滑动距离和滑动速率的影响；沈佳轶等^［10］采用耦合欧拉‒拉格朗日（Coupled Eulerian-Lagrangian，CEL）框架下的欧拉方法揭示了加拿大Saint-Jude滑坡的渐进性破坏机制，计算结果表明，坡脚侵蚀导致的局部剪切带的形成和拓展是灵敏性黏土斜坡发生渐进性破坏的主要诱因；Li等^［11］和Zhang等^［12］分别采用离散单元法（Discrete Element Method，DEM）和粒子有限单元法（Particle Finite Element Method，PFEM）针对深圳杨宝地滑坡的破坏模式展开细致研究，分析了不同材料参数对其运动规律的影响。近些年发展起来的物质点法（Material Point Method，MPM）同时具有拉格朗日和欧拉方法的优势，在模拟固体材料的大变形破坏和界面演化等问题时具有良好的计算效率和精度^［13-15］。该方法同样被广泛用于探究滑坡的致灾机制。Li等^［16］和Fernández等^［17］采用物质点法模拟了汶川地震诱发的大光包巨型滑坡的破坏过程，分别从二维和三维的角度揭示了地质条件和地貌特点等因素对滑坡运动规律的影响。但是，上述研究并未充分考虑碎屑流与下部基岩之间的界面摩擦作用。碎屑流在高速运动过程中与基岩表面发生剧烈摩擦会产生大量热量，由此产生的气垫层将降低摩擦系数。此外，碎屑流会改变基岩表面轮廓特征，造成接触面出现摩擦弱化现象。Zhang等^［18］的研究成果表明，受碎屑流与基岩之间的摩擦生热影响，大光包巨型滑坡的界面摩擦系数显著下降，致灾范围显著扩大；朱晨光等^［19］采用自主开发的MatDEM离散元软件分析了滑坡滑带摩擦生热过程，结果显示滑带附近会形成高热量区，引发的气垫作用将降低土体抗剪强度。因此，在采用物质点法分析滑坡破坏机理时，需要考虑界面摩擦效应，以得到更为准确的计算结果。

本文在传统的物质点法中引入接触算法，用于探究界面摩擦效应对滑坡运动规律的影响。针对深圳杨宝地滑坡，首先，将物质点法计算得到的最终堆积形态和特征点处的速度变化规律与现有的研究成果进行对比，初步验证了该方法的有效性。随后，探究了界面摩擦效应对该滑坡流滑破坏规律的影响，建立了界面摩擦系数与最远滑动距离和最大滑动速度之间的联系。

1 物质点法理论基础和接触算法

1.1 滑坡的界面摩擦效应

以砂泥岩互层的红层地区斜坡为例，简要论述滑坡的破坏机制^［20］。根据图1，泥岩层在风化剥蚀后将形成砂岩顺层斜坡，表面同时覆盖有一定厚度的残坡积层。人工开挖或强降雨作用可能导致残坡积层发生失稳破坏，并沿砂岩顺层斜坡呈流态状进行长距离运动。碎屑流在高速运动过程中会破坏砂岩顺层斜坡的表面形态并导致局部温度上升，致使其摩擦系数发生改变，形成界面摩擦效应。

滑坡的界面摩擦效应会对其运动规律和最终堆积形态产生重要影响。因此，在采用物质点法模拟滑坡时需要引入新的接触算法，将有助于揭示滑坡的破坏机理，为滑坡灾害防治工作提供科学的指导。

1.2 物质点法计算原理

物质点法的计算原理如图2所示，采用拉格朗日和欧拉双重描述，将材料区域离散为携带所有物理信息的质点集，并采用欧拉背景网格求解空间导数和运动方程。物质点可以独立于背景网格发生运动，能够有效避免传统拉格朗日计算方法存在的网格畸变问题，因此，在处理固体材料的大变形破坏问题时具有显著优势。在实际计算过程中，首先，采用插值函数将物质点的质量和动量等物理信息映射至背景网格结点上；然后，施加边界条件并求解运动方程；最后，将计算得到的网格结点的速度和加速度映射回物质点，用于更新其空间位置和受力状态。

在计算域Ω中，物质点法的动量方程的弱形式可以表示为：

{\int }_{\Omega }\rho {\ddot{u}}_i\delta u_{i}d\Omega +{\int }_{\Omega }{\mathit{\sigma }}_{ij}\delta u_{i,j}d\Omega -{\int }_{\Omega }\rho b_i\delta u_{i}d\Omega -{\int }_{{\Gamma }_t}{\overline{t}}_i\delta u_{i}d\Gamma =\mathrm{0}（1）

式中：\rho为物质密度；{\ddot{u}}_i为加速度；δ为Dirac符号；u_i为位移；{\mathit{\sigma }}_{ij}为柯西应力张量；b_i为物体体力；{\mathrm{\Gamma }}_{\mathrm{t}}为自然边界；{\overline{t}}_i为边界力。下标i和j分别表示物理量的空间分量，满足爱因斯坦求和约定，下标中的逗号表示求导。

为提升计算精度，采用Bardenhagen等^［21-22］提出的广义插值物质点法（Generalized Interpolation Material Point Method，GIMP），将计算域离散为由特征函数表示的质点，根据Petro-Calerkin方法建立新的求解格式。特征函数{\chi }_p(x)定义了质点占据的空间区域，在初始构型中满足单位分解条件\sum _p{\chi }_p(x)=\mathrm{1}，因此，可以将质点的体积表示为：

V_p={\int }_{{\Omega }_p\cap \Omega }{\chi }_p(x)d\Omega （2）

式中：\mathrm{\Omega }为物体现时构形所占据的区域；{\mathrm{\Omega }}_p为质点p的特征函数在现时构形中的支撑域，可以视为质点p在现时构形中占据的区域。

在广义插值物质点法中，某个物理量f可以用其在质点处的值f_p近似表示为：

f(x)=\sum _p{f}_p{\chi }_p(x) （3）

采用上述方法对质点的密度、应力和加速度等物理量进行光滑处理，代入式（1）给出的动量方程的弱形式，可以得到：

{\dot{p}}_{iI}=f_{iI}^{\mathrm{i}\mathrm{n}\mathrm{t}}+f_{iI}^{\mathrm{e}\mathrm{x}\mathrm{t}},I=\mathrm{1,2},\dots ,n_g （4）

其中：

{\dot{p}}_{iI}={\displaystyle \sum _{p=\mathrm{1}}^{n_p}}{S}_{Ip}{\dot{p}}_{ip} （5）

f_{iI}^{\mathrm{i}\mathrm{n}\mathrm{t}}=-{\displaystyle \sum _{p=\mathrm{1}}^{n_p}}{\sigma }_{ijp}{S}_{Ip,j}{V}_p （6）

f_{iI}^{\mathrm{e}\mathrm{x}\mathrm{t}}={\displaystyle \sum _{p=\mathrm{1}}^{n_p}}{m}_p{S}_{Ip}{f}_{ip}+{\int }_{{\Gamma }_t}{N}_I(x){\overline{t}}_i\mathrm{d}\mathrm{\Gamma } （7）

式中：p_{iI}为第I个网格结点在i方向的动量，上标符号“\bullet”表示对时间进行求导；p_{ip}和{\sigma }_{ijp}分别为质点动量和应力张量；f_{iI}^{\mathrm{i}\mathrm{n}\mathrm{t}}和f_{iI}^{\mathrm{e}\mathrm{x}\mathrm{t}}分别为第I个网格结点在i方向的内力和外力；I代表网格节点；n_g为背景计算网格的结点总数；n_p为质点数目；f_{ip}为质点的体力密度；m_p为质点质量；V_p为质点体积；形函数S_{Ip}及其导数S_{Ip,j}为关于背景网格形函数N_I的函数。

1.3 接触算法

为了考虑滑坡的界面摩擦效应，可采用多重背景网格的方法，在物质点法中引入接触算法。Bardenhagen等^［21］建议在接触区域的背景网格结点处引入多重速度场，即两个相互接触的物体s和r，在相同的网格结点I处具有不同的结点速度v_{iI}^r和v_{iI}^s。根据图3，若两个物体对该网格结点的动量均有贡献，则认为两者接触，满足条件：

\left(v_{iI}^{r,k-\mathrm{1}/\mathrm{2}}-v_{iI}^{s,k-\mathrm{1}/\mathrm{2}}\right)n_{iI}^{r,k-\mathrm{1}/\mathrm{2}}>\mathrm{0}（8）

式中：n_{iI}^{r,k}为物体r的边界在网格结点I处的外法向线单位矢量。

两个物体接触后，接触结点处产生的接触力f_{iI}^{\mathrm{g},\mathrm{c},k} 由法向接触力f_{iI}^{\mathrm{g},\mathrm{n}\mathrm{o}\mathrm{r},k}和切向接触力f_{iI}^{\mathrm{g},\mathrm{t}\mathrm{a}\mathrm{n},k}两部分组成，计算方法分别为：

f_{iI}^{\mathrm{g},\mathrm{n}\mathrm{o}\mathrm{r},k}=[f_{I,\mathrm{1}}^{\mathrm{n}\mathrm{o}\mathrm{r},k}+\mathrm{m}\mathrm{i}\mathrm{n}(f_{I,\mathrm{2}}^{\mathrm{n}\mathrm{o}\mathrm{r},k},\mathrm{0})]n_{iI}^{\mathrm{g},k}（9）

f_{iI}^{\mathrm{g},\mathrm{t}\mathrm{a}\mathrm{n},k}=\mathrm{m}\mathrm{i}\mathrm{n}(\mu \left|f_I^{\mathrm{n}\mathrm{o}\mathrm{r},k}\right|,f_I^{\mathrm{t}\mathrm{a}\mathrm{n},k})t_{iI}^{\mathrm{g},k}（10）

式中：f_{I,\mathrm{1}}^{\mathrm{n}\mathrm{o}\mathrm{r},k}为将背景网格接触结点的初始速度场调整为满足非嵌入条件的速度场所需施加的力；f_{I,\mathrm{2}}^{\mathrm{n}\mathrm{o}\mathrm{r},k}为满足加速度连续条件所需施加的力；f_I^{\mathrm{n}\mathrm{o}\mathrm{r},k}和f_I^{\mathrm{t}\mathrm{a}\mathrm{n},k}分别为法向和切向接触力；t_{iI}^{\mathrm{g},k}为切向单位矢量；\mu为摩擦系数。

该接触算法的详细论述可以参见文献［21］、［23］。

2 深圳杨宝地滑坡物质点法模拟

2.1 滑坡的地质条件

以深圳某填土滑坡为案例，采用物质点法进行模拟，并基于接触算法探究界面接触效应对该滑坡大变形破坏规律的影响。据相关资料，该填土滑坡发生在深圳市梅林关口杨宝地山，地质条件如图4所示，在强降雨作用下约2.5×10⁴ m³的松散填土发生滑塌，造成重大人员伤亡和经济损失^［11］。

根据图4，该填土滑坡主要由燕山期花岗岩、花岗岩残积土和人工填土组成，其中填土层的土质较为松散，厚度3~10 m。滑坡在平面上整体呈勺状，主轴走向大致为NE45°，总长超过315 m。发生破坏后，滑体沿NE57°向前运动约140 m，形成舌状堆积，属于典型的流滑破坏。

赵春宏等^［24］指出降雨入渗引起地下水位上升是导致该土质滑坡发生破坏的主要原因，根据现场取样的松散击实填土的等压固结不排水剪切试验和常应力排水剪切试验，测得该土体的内摩擦角\phi =\mathrm{28}°，黏聚力c=0 kPa。Li等^［11］和Zhang等^［12］分别采用离散单元法和粒子有限单元法模拟了该滑坡的大变形破坏过程，深入分析了不同材料参数对最远滑动距离和破坏后的堆积形态的影响，得到了较为相似的计算结果。

上述研究成果能够为本文提供重要参考。后文将在相同材料参数条件下采用物质点法模拟深圳杨宝地滑坡，并且将计算结果与Li等^［11］和Zhang等^［12］的研究成果进行对比，用于验证该方法的有效性；随后，探究界面摩擦效应对滑坡流滑破坏规律的影响。

2.2 物质点法计算模型

该填土滑坡的物质点法计算模型如图5所示。由基岩和填土层两部分组成，分别采用线弹性模型和Mohr-Coulomb模型来描述其力学行为，对应的材料参数如表1所示。在计算过程中，组成基岩和填土层的物质点数目分别为470 472和17 676，背景网格间距为0.5 m，总的计算时间为20 s。根据朱晨光等^［19］的研究成果，土质边坡的摩擦系数为0.1~0.4，为考虑碎屑流和基岩之间的相互摩擦作用，假定界面摩擦系数分别为0.20、0.25、0.30、0.35、0.40，模拟该土质滑坡的高速远程破坏过程。

2.3 计算结果分析

图6给出了不同计算时刻滑坡的界面演变规律及其最终堆积形态，并与粒子有限单元法的计算结果进行对比^［9］。

由图6可知：当计算时间t=5 s时，位于测点P₁和P₂处的土体首先发生破坏，具有向下运动的趋势。在t=10 s时刻，该滑坡将发生流滑破坏，整体呈流体状沿着基岩面向下作长距离运动，水平方向最大运动速度约为15 m/s，与粒子有限单元法的计算结果较为接近。当t=15 s时，碎屑流将在水平承台上形成堆积，运动速度也逐渐减小，并最终在t=20 s时刻处于相对静止状态。

物质点法和粒子有限单元法计算得到的滑坡最远滑动距离均为300 m左右，接近现场勘查结果。不过两者在界面演化规律方面仍存在差异，具体表现为：在粒子有限单元法计算过程中，滑坡在流滑破坏后会全部堆积在坡脚平坦区域，而物质点法计算结果则表明，受基岩的地形地貌特点影响，填土会在基岩较为缓倾的区域形成局部堆积。

为进一步验证本文计算结果的有效性，分别将位于图5中不同测点处土体的水平方向速度v_x 和竖直方向速度v_y 与Li等^［11］给出的离散单元法计算结果进行对比分析，结果如图7~图9所示。

由图7~图9可以看出：在测点P₂和P₃处物质点法计算得到的v_x 和v_y 与Li等^［11］的计算结果较为接近，在滑坡启动后碎屑流的水平方向速度随计算时间快速增大，分别在12 s和6 s左右达到最大值15.2 m/s和12.6 m/s，由于P₃点位于坡角，会先于P₂点破坏。竖直方向速度具有不同的变化特点，在流滑破坏初始阶段并不会显著增加，直至分别在12 s和8 s时迅速达到峰值状态，随后又在极短时间内快速减小。不过，在测点P₁处物质点法计算得到的v_x 和v_y 的大小与离散单元法计算结果差异较为明显，分析原因主要在于受基岩表面几何轮廓影响，在碎屑流向下滑动过程中位于测点P₁处的土体尽管在开始时具有较大的运动速度，但在10 s的计算时刻，其周围土体即达到局部稳定状态，运动速度会显著减小。对比图6中给出的物质点法和离散单元法计算结果发现，即使在滑坡破坏结束以后，仍有部分物质点堆积在基岩表面，而后者则主要位于坡角。

2.4 界面摩擦效应对碎屑流运动规律的影响

在前文研究基础上，通过改变基岩的摩擦系数来探究界面摩擦效应对该滑坡破坏规律的影响，尤其是对最远滑动距离、最大运动速度和堆积形态的影响。

图10分别给出摩擦系数μ取0.20、0.25、0.30、0.35、0.40时碎屑流的最终堆积形态。

由图10可以看出：增大摩擦系数μ会导致碎屑流的最远滑动距离L_s减小，两者之间具有较好的线性相关性（图11）。该计算结果初步表明：界面摩擦效应会改变滑坡的动能损耗规律，并影响其最终堆积形态。需要补充的是，在本文计算过程中摩擦系数μ始终为定值，而最新的研究结果^［18］表明在碎屑流高速运动过程中会与基岩发生摩擦生热，导致摩擦系数μ逐渐减小，这也会对碎屑流的最终堆积形态产生影响。

图12给出了滑坡在休止状态下的等效塑性应变{\overline{\epsilon }}_p的分布规律。

可以看出：尽管增大摩擦系数\mu会导致碎屑流的最远滑动距离减小，但位于测点P₂和P₃处的土体的最大等效塑性应变{\overline{\epsilon }}_p却明显增大，表明界面摩擦效应会加剧滑坡的表面演变过程，与基岩接触的土体在变形之后将形成触手状的塑性剪切带，表现出更为复杂的内在作用机理。在碎屑流高速运动过程中，在界面摩擦效应作用下，具有不规则几何边界条件的基岩对其表面接触的土体产生更大的扰动，导致该区域土体的运动规律变得更为复杂，其塑性变形突触会充分发育，能量耗散速率也会随之增大，导致碎屑流的最远滑动距离减小。

图13~图15给出了不同摩擦系数条件下位于测点P₁、P₂和P₃处土体的水平方向速度v_x 和竖直方向速度v_y 与时间t的变化规律。

可以看出：改变摩擦系数μ并不会改变其整体演变规律，但增大摩擦系数μ会使土体的最大水平方向速度v_xmax和最大竖直方向速度v_ymax减小，滑坡整体运动速度会逐渐减慢，其最远滑动距离L_s也随之减小。以P₂测点为例，当摩擦系数μ由0.2增长至0.4时，最大水平方向速度v_xmax 分别为18.7 m/s和12.6 m/s，减小幅度为32.6%，在模拟土质滑坡运动过程时若不考虑界面摩擦效应可能导致计算结果与现场监测结果存在明显差异，不利于其破坏机理分析。

图16给出不同测点处土体的最大水平方向速度v_xmax和最大竖直方向速度v_ymax随摩擦系数μ的变化规律，具体数值见表2。整体上，v_xmax和v_ymax的绝对值均随摩擦系数μ增大而减小，表明碎屑流和基岩之间的界面摩擦效应随着摩擦系数增大而增强，滑坡在运动过程中能量耗散速率逐渐增大，导致其最大水平方向和竖直方向的速度逐渐减小。

滑坡流滑破坏过程中，作用于基岩的平均法向摩擦力F_x 和切向摩擦力F_y 随时间的变化规律如图17所示。由图17可知：两者均随着计算时间增大呈周期性变化，并且水平方向的摩擦力会随着时间增大逐渐减小，而竖直方向的摩擦力则几乎不发生变化。

3 结论

（1） 碎屑流在启动后会发生加速破坏，速度在极短时间内达到最大值，之后随着动能耗散逐渐减小，其最终堆积形态受基岩表面几何特征影响。

（2） 碎屑流的运动规律受界面摩擦效应影响，最大运动速度和最远滑动距离会随着摩擦系数增大而减小，但由于碎屑流和基岩之间的相互摩擦作用加剧，周围土体会产生更大变形，其等效塑性应变反而会增大。

（3） 界面摩擦效应会影响碎屑流中土体的剪切破坏规律，靠近基岩的区域会形成局部触手状的等效塑性破坏区，表明在碎屑流呈流状作长距离运动时会不断发生剧烈界面演化。本研究暂未考虑摩擦生热作用，研究工作后续将会展开。

参考文献

[1] 中国自然资源报.去年我国共成功预报地质灾害 534起[EB/OL] 北京：中国自然资源，2021-01-18.

[2] China Natural Resources News. China successfully predicted 534 geological disasters last year [EB/OL], 2021-01-18.

[3] HE X, LIANG D, BOLTON M D. Run-out of cut-slope landslides: mesh-free simulations[J]. Géotechnique, 2018, 68(1): 50-63.

[4] LIN C, PASTOR M, YAGUE A, et al. A depth-integrated SPH model for debris floods: Application to Lo Wai(Hong Kong)debris flood of August 2005[J]. Géotechnique, 2019, 69(12): 1035-1055.

[5] 吴兵, 梁瑶, 李杰, 等. 注浆微型钢管桩组合结构在公路滑坡处治中的应用[J]. 中外公路, 2025, 45(4): 17-24.

[6] WU Bing, LIANG Yao, LI Jie, et al. Application of grouted steel pipe micropile composite structure in highway landslide treatment[J]. Journal of China & Foreign Highway, 2025, 45(4): 17-24.

[7] 芮国荣, 许建聪, 孙钧. 软岩公路隧道下穿古滑坡稳定性与风险分析[J]. 公路, 2024, 69(3): 399-404.

[8] RUI Guorong, XU Jiancong, SUN Jun. Stability and risk analysis of soft rock highway tunnel passing through ancient landslide[J]. Highway, 2024, 69(3): 399-404.

[9] ZHANG W J, ZHENG H, JIANG F Y, et al. Stability analysis of soil slope based on a water-soil-coupled and parallelized smoothed particle hydrodynamics model[J]. Computers and Geotechnics, 2019, 108: 212-225.

[10] 黄洪超. 某高速公路巨型滑坡稳定性分析及处置方案研究[J]. 中外公路, 2023, 43(4): 24-29.

[11] HUANG Hongchao. Study on stability analysis and treatment of a expressway landslide[J]. Journal of China & Foreign Highway, 2023, 43(4): 24-29.

[12] 许强, 黄润秋. 5⋅12汶川大地震诱发大型崩滑灾害动力特征初探[J]. 工程地质学报, 2008, 16(6): 721-729.

[13] XU Qiang, HUANG Runqiu. Kinetics charateristics of large landlides triggered by may 12th Wenchuan earthquake[J]. Journal of Engineering Geology, 2008, 16(6): 721-729.

[14] 张卫杰, 余瑞华, 陈宇, 等. 强度指标影响下滑坡运动特征及参数反分析[J]. 岩土工程学报, 2022, 44(12): 2304-2311, F0004.

[15] ZHANG Weijie, YU Ruihua, CHEN Yu, et al. Post-failure movement characteristics and parameter back-analysis for landslides considering effect of strength parameters[J]. Chinese Journal of Geotechnical Engineering, 2022, 44(12): 2304-2311, F0004.

[16] 沈佳轶, 陈前, 库猛, 等. 基于CEL的灵敏性黏土斜坡渐进破坏数值模拟研究[J]. 岩土工程学报, 2022, 44(12): 2297-2303.

[17] SHEN Jiayi, CHEN Qian, KU Meng, et al. Numerical simulation of progressive failure of sensitive clay slopes using CEL method[J]. Chinese Journal of Geotechnical Engineering, 2022, 44(12): 2297-2303.

[18] LI W C, LI H J, DAI F C, et al. Discrete element modeling of a rainfall-induced flowslide[J]. Engineering Geology, 2012, 149: 22-34.

[19] ZHANG X, KRABBENHOFT K, SHENG D C, et al. Numerical simulation of a flow-like landslide using the particle finite element method[J]. Computational Mechanics, 2015, 55(1): 167-177.

[20] SULSKY D, ZHOU S J, SCHREYER H L. Application of a particle-in-cell method to solid mechanics[J]. Computer Physics Communications, 1995, 87(1/2): 236-252.

[21] SULSKY D, KAUL A. Implicit dynamics in the material-point method[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193(12/13/14): 1137-1170.

[22] SULSKY D, GONG M. Improving the material-point method[M/OL]//WEINBERG K, WALLIN M. Innovative numerical approaches for multi-field and multi-scale problems. Cham: Springer International Publishing, 2016: 217-240.

[23] LI X P, TANG X, ZHAO S X, et al. MPM evaluation of the dynamic runout process of the giant Daguangbao landslide[J]. Landslides, 2021, 18(4): 1509-1518.

[24] FERNÁNDEZ F, VARGAS E, MULLER A L, et al. Material point method modeling in 3D of the failure and run-out processes of the Daguangbao landslide[J]. Acta Geotechnica, 2024, 19(7): 4277-4296.

[25] ZHANG Y F, ZHANG W G, WANG L Q, et al. Mechanism of the high-speed and long-run-out landslide considering the evolution of the frictional heat in the sliding zone[J]. Natural Hazards, 2024, 120(4): 3299-3317.

[26] 朱晨光, 刘春, 许强, 等. 滑坡滑带摩擦热离散元数值模拟研究[J]. 工程地质学报, 2019, 27(3): 651-658.

[27] ZHU Chenguang, LIU Chun, XU Qiang, et al. Discrete element numerical simulation research on friction heat in sliding zone of the landslide[J]. Journal of Engineering Geology, 2019, 27(3): 651-658.

[28] 周翠英, 刘镇, 薛翊国, 等. 关于红层灾变基础研究的若干思考[J]. 工程地质学报, 2023, 31(3): 689-705.

[29] ZHOU Cuiying, LIU Zhen, XUE Yiguo, et al. Some thoughts on basic research of red beds disaster[J]. Journal of Engineering Geology, 2023, 31(3): 689-705.

[30] BARDENHAGEN S G, KOBER E M. The generalized interpolation material point method [J]. Computer Modeling in Engineering and Sciences, 2004, 5(6): 477-496.

[31] BARDENHAGEN S G, BRACKBILL J U, SULSKY D. The material-point method for granular materials[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 187(3/4): 529-541.

[32] MA J, WANG D, RANDOLPH M F. A new contact algorithm in the material point method for geotechnical simulations[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2014, 38(11): 1197-1210.

[33] 赵春宏, 戴福初. 深圳某填土滑坡破坏机理研究[J]. 中国地质灾害与防治学报, 2007, 18(2): 1-8.

[34] ZHAO Chunhong, DAI Fuchu. Study on failure mechanism of a fill slope in Shenzhen[J]. The Chinese Journal of Geological Hazard and Control, 2007, 18(2): 1-8.

基金资助

天津市教委科研计划项目(XJ2022009601)