抗滑桩加固边坡大变形破坏三维物质点模拟

Three-dimensional Simulation of Large Deformation Failure for Slope Reinforced with Anti-slide Piles Based on Material Point Method

摘要

目前，关于含加固物的边坡三维大变形的分析研究极为匮乏.本文利用物质点法构建了抗滑桩加固边坡三维数值模型，分析了抗滑桩加固条件下边坡大变形破坏模式，探讨了抗滑桩长度、布设位置以及布设间距对边坡大变形破坏的影响.结果表明，抗滑桩布设于坡脚时，抗滑桩上方土体易发生局部破坏并产生越顶滑移，桩长增加会缩减塑性区范围，但会加剧土体浅层滑动，增大滑移范围；抗滑桩布设于坡中时，边坡破坏后的位移量最小，表现出良好的加固效果，但当桩长与潜在滑动面深度大致相同时，桩体的存在可能会导致边坡破坏后的滑移距离更远；桩间距的增加会引起桩间土的塑性流动；当抗滑桩布设于坡顶时，抗滑桩下方土体易发生脱开滑移，加固效果较差.

Abstract

Current studies on the three-dimensional （3D） large deformation analysis of slopes with reinforcement structures remain limited. In this study， a 3D numerical model for a slope reinforced with anti-slide piles is constructed using the material point method. The failure modes of the slope with anti-slide piles are analyzed， and the effects of the pile length， location， and spacing on slope failure are also discussed. The results show that， when the anti-slide piles are located at the slope toe， the soil above the anti-slide piles is prone to local failure， inducing overtopping and sliding. Increasing the pile length reduces the area of the plastic zone， but exacerbates shallow sliding and increases runout. When the anti-slide piles are located in the middle of the slope， the displacement of the slope after failure is the smallest， and the reinforcement effect is the best； however， when the pile length is roughly the same as the depth of the potential slip surface， the presence of piles may lead to a longer runout distance after failure. The increase of the pile spacing causes the plastic flow of the soil between the piles. When the anti-slide piles are located at the slope crest， the soil beneath them is easily dislodged and slides and the reinforcement effect is poor.

Graphical abstract

关键词

抗滑桩 / 加固 / 大变形 / 物质点法 / 边坡稳定性

抗滑桩因具有加固效果突出、施工扰动小以及布置灵活等优点，在边坡防护中得到广泛应用. 目前，针对抗滑桩-边坡体系的物理评价方法主要可分为3类. 第一类是基于极限平衡法对抗滑桩-边坡稳定性进行理论分析，这类方法概念清晰、使用方便，但无法考虑边坡与抗滑桩的相互作用^［1-3］，且通常局限于二维平面应变分析，忽略了抗滑桩体积效应对边坡稳定性的影响^［4-5］；第二类是基于室内模型试验的评价方法^［6-8］，此类方法多聚焦于桩土系统破坏的力学规律，但受制于试验条件通常采用缩尺模型，因此无法评价真实尺寸下抗滑桩的作用效果，且费用昂贵；第三类是以有限元法为代表的模拟评价方法，因其使用成本较低且与真实工况相近，目前已被众多学者广泛应用.

彭文哲等^［9］基于抗滑桩边坡二维有限元分析，探讨了抗滑桩加固后最危险滑动面的变化规律；Han等^［10］针对张家窑滑坡降雨失稳机理对抗滑桩支护参数进行了优化设计，提高了抗滑桩的支护效果. 孙俊等^［11］分析了三维抗滑桩作用下土岩接触面滑坡下潜在的破坏模式与抗滑桩加护效应. 以上研究为抗滑桩加固边坡评价提供了良好借鉴，但受限于有限元网格畸变问题，这些研究均侧重于抗滑桩加固边坡稳定性的小变形分析，难以考虑抗滑坡在加固条件下边坡的大变形破坏情况，抗滑桩的存在可能改变边坡大变形破坏模式. 因此，亟须开展抗滑桩加固边坡破坏的三维大变形数值分析工作，全面探究边坡破坏工况下抗滑桩对土体大变形的影响规律.

近年来，无网格法（Mesh-free methods）的发展为岩土大变形问题分析提供了有力工具. 无网格法无须依赖网格对介质信息的近似，而是在任意分布的节点处构造插值函数离散控制方程进行求解^［12-13］. 常见的无网格算法有离散元法（DEM）、光滑粒子流体动力学（SPH）以及物质点法（MPM）等. 离散元法虽然可以较好模拟滑坡破坏模式演变等极端变形问题，但颗粒间微观参数与宏观参数对应机制不明，难以准确解释微观参数的物理意义^［14］；光滑粒子流体动力学基于土体本构模型较好地表征滑坡物质的运动全过程^［15］，但每一步需要对边界粒子进行搜索，计算量较大^［16］. 相较于前两种无网格方法，物质点法充分结合了拉格朗日和欧拉算法的优势，在保证连续介质力学计算效率的基础上，克服了网格畸变的缺点^［17-18］；此外，物质点法可以有效解决不同材料间的非滑动接触问题^［19-20］. 目前，物质点法已被广泛用于降雨滑坡^［21］、边坡管道相互作用^［22］、CPT贯入^［23］等问题中. 针对含加固物的边坡大变形问题，目前文献基于物质点法主要研究了非嵌入式挡墙与滑坡冲击间的相互作用^［24-25］，而对桩、锚加固物的研究十分匮乏. 此外，为保证数值稳定，现有含加固物的边坡大变形研究大多采用二维平面应变假设^［26-27］，无法考虑抗滑桩加固条件下的三维土拱效应与土体绕流现象.

本文依托Anura3D平台（www.anura3d.com）建立抗滑桩加固边坡三维物质点模型. Anura3D采用了节点混合离散技术^［28］，有效增强了三维大变形数值分析稳定性. 基于该三维物质点模型，系统探究抗滑桩加固条件下边坡破坏大变形，总结桩土系统的大变形破坏模式，并对抗滑桩参数进行敏感性分析，研究结果可为基于风险的抗滑桩加固边坡设计提供技术参考.

1 物质点法原理

物质点法是由有限元法发展而来的无网格方法^［29-31］，1994年首次被Sulsky等学者引入到固体力学中^［30］. 该方法采用拉格朗日和欧拉两套描述^［32］；将连续空间介质离散为若干物质点，每个物质点携带所在区域的速度、应变和应力等信息，用于模拟变形（拉格朗日描述）；通过构建背景网格覆盖整个问题域，利用形函数将物质点信息映射到网格节点上求解控制方程（欧拉描述）. 基于连续介质力学质量及动量守恒原理，物质点法的控制方程可表示为下式：

\frac{\mathrm{d}}{\mathrm{d}t}\rho +\rho \nabla \mathit{v}=\mathrm{0}（1）

\rho \mathit{a}=\rho \mathit{b}+\nabla \mathit{\sigma }（2）

式中：ρ为介质密度； v 为速度； a 为加速度； σ 为对称应力张量； b 为体力，如重力. 为求解方便，将公式（2）通过虚功方程转化为弱形式，如下式所示：

\underset{\Omega }{\int }\rho \omega \mathit{a}d\Omega =

\underset{\partial {\Omega }_{\mathit{\sigma }}}{\int }\rho {\mathit{\tau }}^s\omega d\partial {\Omega }_{\mathit{\sigma }}-\underset{\Omega }{\int }\rho \mathit{\sigma }:\nabla \omega d\Omega +\underset{\Omega }{\int }\rho \omega \mathit{b}\mathrm{d}\Omega（3）

式中：Ω为研究对象的积分域；ω为测试函数，在积分域Ω边界处取值为0；∂Ω_σ 为面力 τ^s的积分域.

将连续介质离散为物质点，式（1）可变换为：

{\displaystyle \sum _{p=\mathrm{1}}^{N_p}}{m}_p={\displaystyle \sum _{p=\mathrm{1}}^{N_p}}\int \rho (\mathit{x},t)\delta (\mathit{x}-{\mathit{x}}_p)\mathrm{d}{V}_p=

\int \rho (\mathit{x},t)\mathrm{d}\Omega（4）

式中：m_p 为物质点p的质量；δ为Dirac Delta函数；N_p 为介质离散的物质点数； x 为物体介质中任意一点的空间位置； x_p 为物质点p当前的空间位置. 将式（4）代入式（3）中，可以得到动量守恒方程的离散形式，如下式所示：

{\displaystyle \sum _{p=\mathrm{1}}^{N_p}}{m}_p[\omega ({\mathit{x}}_{{}_p}^t){\mathit{\tau }}^s({\mathit{x}}_{{}_p}^t)h^{-\mathrm{1}}+\omega ({\mathit{x}}_{{}_p}^t)\mathit{b}({\mathit{x}}_{{}_p}^t)]=

{\displaystyle \sum _{p=\mathrm{1}}^{N_p}}{m}_p\mathit{\sigma }({\mathit{x}}_{{}_p}^t):\nabla \omega ({\mathit{x}}_{{}_p}^t)+{\displaystyle \sum _{p=\mathrm{1}}^{N_p}}{m}_p[\omega ({\mathit{x}}_{{}_p}^t)\mathit{a}({\mathit{x}}_{{}_p}^t)]（5）

式中：h为边界层厚度.

背景网格与物质点上的状态变量基于形函数进行相互映射：

{\mathit{x}}_p={\displaystyle \sum _{i=\mathrm{1}}^{N_g}}{N}_i({\mathit{x}}_p^t){\mathit{x}}_{{}_i}^t（6）

式中：N_g 为背景网格的结点总数；N_i （{\mathit{x}}_p^t）为节点i在物质点p上的形函数，表达为获取节点i信息时的加权系数； x^t_i 为t时刻节点i的坐标位置. 映射后，网格节点上的动量守恒方程为：

{\mathit{M}}_i{\mathit{a}}_i^t={\mathit{F}}_i^t=({\mathit{F}}_i^t{)}^{\mathrm{e}\mathrm{x}\mathrm{t}}-({\mathit{F}}_i^t{)}^{\mathrm{i}\mathrm{n}\mathrm{t}}（7）

式中： M_i 为节点i的质量矩阵， F^t_i 为节点i上的合力，（{\mathit{F}}_i^t）^ext和（{\mathit{F}}_i^t）^int分别为节点i上的外力以及内力.

物质点法数值求解采用显式积分算法，主要计算过程可分为4步：1）离散连续介质获得物质点，通过形函数将物质点信息映射到背景网格节点上；2）在网格节点上求解动量守恒方程，即公式（7）；3）基于动量守恒方程的求解结果更新物质点加速度、速度和位置等信息，并基于本构方程计算物质点应力、应变；4）重置背景网格，进入下一个时间步的计算.

2 抗滑桩加固边坡物质点模型

2.1 模型设置

为提高数值分析效率，采用Cai等^［33］提出的三维半无限对称模型分析方法，抗滑桩加固边坡模型设置为半桩及0.5倍桩间距内的边坡土体，该分析方法已在文献中得到论证^［34］. 参考Duan等^［35］的研究案例，基于GiD平台建立如图1所示的抗滑桩加固边坡三维几何模型，其中边坡坡度为1∶1.5，坡高为6 m. 抗滑桩采用方形桩，桩径W=1 m；桩间距为S，即桩体中心间的距离；桩长为L.

基于Anura3D建立抗滑桩加固边坡的三维物质点模型，如图2所示. 该物质点模型尺寸与图1一致，其中黑色网格代表背景网格，浅棕色和红色颗粒分别代表土和桩的物质点颗粒. 背景网格采用四面体网格，覆盖边坡破坏后的可能运动区域；非物质点区域的背景网格尺寸为3 m；桩土接触区域的背景网格尺寸为0.5 m，每个网格内包含4个物质点；物质点采用单套单相框架，仅模拟固相力学行为. 考虑单层匀质边坡，土体采用Mohr-Coulomb破坏准则的理想弹塑性本构，抗滑桩采用理想线弹性本构.

抗滑桩加固边坡通常处于稳定状态. 为模拟抗滑桩加固条件下边坡大变形破坏，采用强度折减的方式，将土体强度设置为残余强度. 参考超固结黏土的残余强度^［36-37］，本案例材料参数取值如表1所示. 可以看到，表1中抗滑桩杨氏模量取值小于混凝土材料杨氏模量但显著大于土体杨氏模量，这主要是为了在确保抗滑桩结构刚性行为的同时，提升物质点法的分析效率^{［19，22］}.

由于桩土系统破坏时桩土间相对位移通常较小，考虑桩土间接触模式为非滑动接触. 如引言所述，物质点法可以自动处理不同物体间的非滑动接触问题^［38-40］，模拟桩土材料间的黏结和脱开行为^{［30，41-42］}. 计算过程中将模型底部与背景网格顶部的X和Y方向位移、左右两侧的X方向位移、前后两侧的Z方向位移进行约束固定.

2.2 模型验证

目前关于抗滑桩加固边坡大变形分析的研究十分匮乏，难以直接获取相关模型试验或工程案例对本文建立的数值模型进行验证. 由于有限元法在工程领域被广泛认为具有较高的分析精度^{［19， 43］}，因此考虑通过对比有限元与物质点在小变形阶段的分析结果^{［14， 44-46］}，间接验证本文模型的有效性.

参照文献［47］中的案例，图3（a）（c）分别展示了桩间距与桩径的比值（S/W）分别为4和6时，基于物质点法与有限元法分析获得的抗滑桩加固边坡剪切应变. 可以看到，两种方法得到的边坡塑性区及滑动面基本一致，土体变形量也较为接近，表明该三维物质点模型具有较高的模拟精度.

3 抗滑桩加固边坡大变形破坏分析

过去有限元分析研究^［48］指出桩土系统在小变形条件下主要破坏模式包括桩间土塑性流动、边坡局部破坏和边坡整体破坏等. 本节利用上文建立的物质点模型系统研究不同抗滑桩设计参数下边坡的大变形破坏模式，探究抗滑桩对边坡大变形破坏的影响. 考虑3类抗滑桩设计参数^［34］，分别为桩中心距坡脚的水平距离L_x （桩位），桩长L以及桩间距S（图1）. 通过组合这3类抗滑桩设计参数，共设置27种工况，具体取值如表2所示. 后续分别对这27种工况下桩土系统的应变、位移依次进行分析，评价不同抗滑桩设计参数对边坡破坏后运动距离的影响. 本文采用台式机（处理器为AMD Ryzen 9 5950X 16-Core Processor 3.40 GHz）开展抗滑桩加固边坡物质点数值分析，预分析一组工况的平均物理时间为3.3 h，最大物理时间为7.5 h.

以往大变形分析中通常采用位移阈值作为边坡破坏判据^［49-51］，该阈值与边坡尺寸密切相关. 本案例中，采用1 m的位移阈值，即当边坡位移大于1 m时则认为边坡发生破坏. 同样，采用0.01 m的位移阈值作为平衡收敛准则，即当下一时间步的边坡位移变化量小于0.01 m，即判断该边坡破坏后再次处于平衡状态. 下文展示的边坡剪切应变或位移云图在未说明模拟时间情况下均是边坡失稳后再次建立平衡后的状态.

3.1 未加固边坡破坏分析

首先对未加固边坡的破坏情况进行分析. 图4（a）（b）分别展示了未加固边坡的剪切应变与位移分布情况.可以看到，边坡在自重作用下发生了整体剪切破坏，形成一条贯通的圆弧形滑动面，滑面深度约为4 m；坡脚处剪切带发育最为明显，坡面土体的最大位移约为4.0 m.

3.2 抗滑桩位于坡脚处边坡破坏分析

图5展示了抗滑桩布设在坡脚处（L_x =2 m），桩间距和桩长均为6 m时（S/W=6，L=6 m），边坡位移随模拟时间T的变化情况（图中红色框线代表抗滑桩所在位置）. 与图4进行对比可以看出，抗滑桩的存在有效抑制了坡脚处剪切带的发育，但对上部土体的加固效果较差，坡中和坡顶发生了局部浅层破坏，且滑体最终越过抗滑桩顶部发生滑移，形成了桩土破坏模式中典型的越顶破坏. 土体最大位移约为5.5 m，相较于未加固时更大，这主要是由于抗滑桩的存在使深层转动破坏转变为浅表滑移，此时更多势能转换成了动能，增大了滑塌距离. 类似现象在相关文献中也有过报道^{［49， 51-52］}.

当桩间距与桩径比值固定（S/W=6）时，图6展示了桩长变化对边坡越顶破坏程度的影响.一般来说，桩长的增加可以提升边坡的安全储备，但对比图6（b）（d）与图5（d）可以发现，抗滑桩加固条件下边坡发生破坏时，抗滑桩并未切穿滑动面致使边坡失效，此时长抗滑桩更易形成浅层越顶破坏，且滑体越过坡脚后的运动距离更远.

3.3 抗滑桩位于坡中处边坡破坏分析

图7展示了抗滑桩位于坡中时（L_x =5 m），小桩间距条件下（S/W=2）桩长变化对边坡剪切应变和位移的影响. 可以看到，当桩长较短时，边坡深层形成贯通的剪切滑移带；除了坡脚处，桩体附近沿桩身也形成了明显剪切带，坡体的滑动造成抗滑桩的整体倾覆，此时边坡破坏模式为整体深层破坏与桩前后局部浅层破坏的叠加.随着桩长增加，边坡深层剪切带被隔断，破坏模式逐渐完全转变为桩前后土体的局部浅层破坏；相较于4 m的桩长，6 m桩长下坡体浅层滑动的位移量更大.以上现象的产生主要是由于当抗滑桩的布设间距较小时，桩长增加虽然会减小边坡的整体位移使边坡安全性能得到提升，但同时会加剧边坡的局部浅层滑移，增大位移量.当抗滑桩布设间距较大时，桩长的增加对边坡局部浅层位移的影响将会减弱，这主要是由于桩间土发生塑性流动.以S/W=6、L=6 m为例，图8展示了较大桩间距条件下边坡破坏的位移演化过程. 可以看到，桩间土塑性绕流位移场由于桩体支护作用并不连续，近桩侧桩体拦挡效果较强，桩体前后的土体塑性变形较小.

图9和图10分别展示了S/W=6、L=4 m时边坡剪切应变和位移的变化. 可以发现，桩长与边坡整体剪切滑移带深度一致增大了边坡破坏后滑移. 以上现象主要是由于桩间土塑性流动导致边坡整体剪切破坏，造成桩体倾覆，抗滑桩的倾覆转动促使坡体前缘土体发生进一步滑移，坡面处桩间土体塑性区逐渐贯通并产生明显滑移，最终形成更大的滑体体积和影响范围. 因此，在工程实践中应通过增加桩长尽可能避免该类现象发生.

3.4 抗滑桩位于坡顶处边坡破坏分析

图11展示了抗滑桩位于坡顶（L_x=8 m）且小桩间距（S/W=2）条件下，桩长变化对边坡剪切应变和位移的影响. 由图11可知，由于桩体完全位于边坡滑面以上，虽然桩后土体塑性变形较小，但桩底土体仍形成了完整的塑性剪切带，桩前土体脱开桩体形成明显滑移，此时边坡破坏模式为桩前土体局部浅层破坏和坡身整体深层破坏的叠加. 随着桩长增加，桩后缘靠近坡肩处的土体变形量不断减小，边坡破坏模式逐渐完全转变为桩前土体的局部浅层破坏. 图12展示了抗滑桩位于坡顶且桩间距较大（S/W=6）时边坡破坏的位移演化过程. 可以看到，与抗滑桩位于坡中规律类似，桩间距较大诱发桩间土的塑性流动导致桩长对边坡局部浅层破坏的影响变弱.

3.5 抗滑桩-边坡潜在大变形模式总结

表3总结了不同工况下抗滑桩-边坡体系破坏后的最大位移和潜在大变形模式. 抗滑桩位于坡脚时，边坡易发生浅层越顶滑移，边坡破坏后产生的位移量最大；桩间距较小时，土体侧向绕流变形困难，更倾向于沿滑动方向变形，此时桩长的增加会显著加剧桩后非稳定土体的浅层破坏，导致更大的土体位移量；当S/W=2、L=6 m时，最大位移量达到8.5 m. 抗滑桩位于坡中时，边坡位移量最小；小桩间距条件下，坡体深层滑动面随桩长的增加被隔断，形成更加明显的桩前后坡面的局部浅层破坏；当桩间距增大时，桩间土发生塑性流动破坏. 抗滑桩位于坡顶时，主要发生桩前土脱开桩体的局部浅层破坏.

4 抗滑桩参数敏感性分析

4.1 桩位对边坡最大位移的影响

为直观量化不同桩位抗滑桩的加固效果，图13展示了S/W=6时桩位对边坡最大位移的影响. 基于最大位移指标对比可以发现，当抗滑桩布设于坡中位置时能对较好地限制边坡变形，其次将抗滑桩布设于坡顶时对边坡变形的控制效果略优于坡脚处布设，该结论与基于有限元的抗滑桩加固边坡稳定性分析结论一致^{［6， 53］}. 当桩长L=6 m时，边坡的最大位移对于桩位的变化最敏感.

4.2 桩长对边坡最大位移的影响

以坡中为例，图14分析了不同桩间距条件下桩长对边坡最大位移的影响. 由图14可知，当桩间距较小时，随着桩长增加，边坡变形逐渐得到控制，加固效果明显；但当桩长过长时，边坡整体位移会得到控制，但会加剧土体局部的浅层滑动. 当桩间距较大时，边坡破坏主要为桩间土塑性流动，此时桩长的增加对边坡加固效果提升有限. 图14中S/W=6且L= 4 m时的边坡位移过大是由于桩长与滑动面深度相当，抗滑桩的倾覆加剧了坡体大变形破坏，如前文3.3小节所述.

4.3 桩间距对边坡最大位移的影响

同样以坡中为例，分析不同桩长条件下桩间距对边坡最大位移的影响，结果如图15所示. 随着桩间距的增大，抗滑桩的加固效果逐渐减弱，其原因在于随着桩间距扩大，桩间土发生塑性流动，难以形成土拱对土体进行有效拦挡^［54］. 当桩较短，桩的加固效果更依赖桩间土拱的拦挡作用，因此当桩长L=2 m时，边坡的最大位移对于桩间距的变化最敏感. 当桩长L=4 m时，桩间距增大使得边坡位移骤增的现象也是由桩长与滑动面深度一致引起的，如前文3.3小节所述.

5 结 论

Anura3D为建立抗滑桩加固边坡三维物质点模型提供了有力工具，有效改善了三维数值分析不稳定问题. 本文基于建立的三维物质点模型探明了抗滑桩加固边坡破坏大变形模式，揭示了桩位、桩长和桩间距对抗滑桩加固边坡破坏大变形的影响，得到如下主要结论：

1）抗滑桩位于坡脚或坡顶附近时难以对边坡滑动面进行有效加固，会使边坡破坏后产生更远的浅表滑移，但桩靠近坡顶处能对高势能易滑动土体进行拦挡，其加固效果略优于坡脚.

2）桩长变化影响边坡破坏后大变形特性的演变. 随着桩长增加，边坡破坏模式由深层转动转变为浅表滑移；边坡破坏后，桩身过长会增加坡体的表层位移. 值得注意的是，当桩长与滑动面深度相当时，桩体转动会促进塑性区向坡前缘发育，加剧土体破坏.

3）桩间距增加致使近桩侧与远桩侧土体的破坏模式存在差异. 随着桩间距增加，桩间土的塑性流动趋势逐渐增强，抗滑桩对远桩侧土体的支护效果逐渐减弱，此时边坡破坏模式呈现出多种破坏模式的叠加.

应当指出，本文抗滑桩加固边坡体系破坏研究仅通过数值比对初步证明了本文模型建立的合理性，但本文方法的正确性和有效性需通过大变形实际工程分析案例进行进一步证明. 由于相关研究较为匮乏，未来仍需开展抗滑桩加固边坡大变形破坏的室内模型试验或现场试验进行验证.

参考文献

[1] VIGGIANI C ．Ultimate lateral load on piles used to stabilize landslides［C］//10th International Conference on Soil Mechanics and Foundation Engineering （Stockholm），1981，3：555-560．

[2] ITO T， MATSUI T ．Methods to estimate lateral force acting on stabilizing piles［J］．Soils and Foundations，1975，15（4）：43-59．

[3] 戴自航，沈蒲生 ．抗滑桩内力计算悬臂桩法的改进［J］．湖南大学学报（自然科学版），2003，30（3）：81-85．

[4] DAI Z H， SHEN P S ．Improvements on calculation of internal forces of cantilever anti-sliding piles［J］．Journal of Hunan University （Natural Sciences），2003，30（3）：81-85．（in Chinese）

[5] 陈昌富， 杜翠翠，张根宝. 基于双参数法刚性抗滑桩嵌固深度可靠性设计［J］. 湖南大学学报（自然科学版），2014，41（2）：40-46．

[6] CHEN C F， DU C C， ZHANG G B ．Reliability design of the embedded depth of the rigid anti-slide piles based on bi-parameter method［J］．Journal of Hunan University （Natural Sciences），2014，41（2）： 40-46．（in Chinese）

[7] 侯超群，丁莹，孙志彬，等 ．抗滑桩加固边坡三维稳定因素的敏感性分析［J］．合肥工业大学学报（自然科学版），2022，45（8）：1092-1099．

[8] HOU C Q， DING Y， SUN Z B，et al ．Sensitivity analysis of 3D stability factors of slopes reinforced with piles［J］．Journal of Hefei University of Technology （Natural Science），2022，45（8）：1092-1099．（in Chinese）

[9] 张泽蔚，刘素嘉，张嘎 ．降水条件下抗滑桩-锚杆加固边坡的离心模型试验研究［J］．岩土工程学报，2023，45（）：206-209．

[10] ZHANG Z W， LIU S J， ZHANG G ．Centrifuge modeling of anti-slide piles with anchored bolt-reinforced slopes under drawdown conditions［J］. Chinese Journal of Geotechnical Engineering，2023，45（Sup.1）： 206-209．（in Chinese）

[11] 宋晓东 ．怀邵衡铁路DK324边坡坍滑分析及处理方案优化［J］．湖南大学学报（自然科学版），2018，45（）： 30-34．

[12] SONG X D. Analysis of DK324 slope slump in Huaishaoheng railway and treatment optimization［J］. Journal of Hunan University （Natural Sciences）， 2018， 45（Sup.1）： 30-34．（in Chinese）

[13] ZHANG S F， LI C， QI H，et al. Soil arch evolution characteristics and parametric analysis of slope anchored anti-slide pile［J］. KSCE Journal of Civil Engineering， 2021， 25（11）：4121-4132．

[14] 彭文哲，赵明华，肖尧，等 ．抗滑桩加固边坡的稳定性分析及最优桩位的确定［J］．湖南大学学报（自然科学版），2020，47（5）：23-30．

[15] PENG W Z， ZHAO M H， XIAO Y，et al ．Stability analysis of anti-slide pile reinforced slope and determination of optimal pile position［J］．Journal of Hunan University （Natural Sciences），2020，47（5）：23-30．（in Chinese）

[16] HAN C X， HOU T S， CHEN Y ．Three-dimensional slope stability and anti-slide pile treatment of Zhangjiayao landslide under rainfall［J］．Engineering Computations，2023，40（9/10）：2667-2692．

[17] 孙俊，刘亦威，邢爱国 ．土-岩接触面滑坡孕灾机理与加护效应研究［J/OL］．武汉大学学报（工学版），2023：1-14．［2024-04-30］．

[18] SUN J， LIU Y W， XING A G ．Research on the disaster pregnant mechanism and reinforcement effect of a soil-rock contact surface landslide［J/OL］．Engineering Journal of Wuhan University，2023：1-14．［2024-04-30］．（in Chinese）

[19] 孙玉进， 宋二祥 ．大位移滑坡形态的物质点法模拟［J］．岩土工程学报，2015，37（7）：1218-1225．

[20] 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-1225．（in Chinese）

[21] ZHANG M Q， ZAINAL ABIDIN A R， TAN C S ．State-of-the-art review on meshless methods in the application of crack problems［J］. Theoretical and Applied Fracture Mechanics，2024，131：104348．

[22] 史卜涛，张云，张巍 ．边坡稳定性分析的物质点强度折减法［J］．岩土工程学报，2016，38（9）：1678-1684．

[23] 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）

[24] PASTOR M， BLANC T， HADDAD B，et al ．Application of a SPH depth-integrated model to landslide run-out analysis［J］. Landslides，2014，11（5）：793-812．

[25] SUN Z， LI H Q， GAN Y，et al. Material point method and smoothed particle hydrodynamics simulations of fluid flow problems： a comparative study［J］. Progress in Computational Fluid Dynamics，an International Journal，2018，18（1）：1-18．

[26] 李光耀，卡里鲁 ．弹塑性大变形畸变问题的无网格分析［J］．湖南大学学报（自然科学版），2003，30（1）：47-49．

[27] LI G Y， KA L L ．Element free Galerkin method for elasto-plastic problems with large deformation distortion［J］．Journal of Hunan University （Natural Science），2003，30（1）：47-49．（in Chinese）

[28] YERRO A， ALONSO E E， PINYOL N M ．The material point method for unsaturated soils［J］．Géotechnique，2015，65（3）：201-217．

[29] SOGA K， ALONSO E， YERRO A，et al ．Trends in large-deformation analysis of landslide mass movements with particular emphasis on the material point method［J］．Géotechnique，2016，66（3）：248-273．

[30] YERRO A ．MPM modelling of landslides in brittle and unsaturated soils［D］．Universitat Politècnica de Catalunya，Departament d'Enginyeria Civil i Ambiental，2015．

[31] LIU X， WANG Y ．Probabilistic simulation of entire process of rainfall-induced landslides using random finite element and material point methods with hydro-mechanical coupling［J］．Computers and Geotechnics，2021，132：103989．

[32] ZHU H H， XIE T C， ZHANG W，et al ．Numerical simulations of a strip footing on the soil slope with a buried pipe using the material point method［J］．International Journal of Geomechanics，2023，23（11）：04023190．

[33] GAO L， GUO N， YANG Z X，et al ．MPM modeling of pile installation in sand：Contact improvement and quantitative analysis［J］．Computers and Geotechnics，2022，151：104943．

[34] WANG J， CHAN D ．Frictional contact algorithms in SPH for the simulation of soil-structure interaction［J］．International Journal for Numerical and Analytical Methods in Geomechanics， 2014， 38（7）： 747-770．

[35] CHANDRA B， LARESE A， IACONETA I，et al ． Soil-structure interaction simulation of landslides impacting a structure using an implicit material point method［C］//2nd International Conference on the Material Point Method for Modelling Soil-water-structure Interaction，2019： 72-78．

[36] BUI H， SAKO K， FUKAGAWA R，et al ．SPH-based numerical simulations for large deformation of geomaterial considering soil-structure interaction［C］//The 12th International Conference of International Association for Computer Methods and Advances in Geomechanics （IACMAG）， Goa， India， 2008， 1： 570-578．

[37] GONZÁLEZ ACOSTA J L， VARDON P J， HICKS M A ．Study of landslides and soil-structure interaction problems using the implicit material point method［J］．Engineering Geology，2021，285：106043．

[38] LIAN Y P， ZHANG X， LIU Y ．An adaptive finite element material point method and its application in extreme deformation problems［J］．Computer Methods in Applied Mechanics and Engineering，2012，241：275-285．

[39] HARLOW F H ．PIC and its progeny［J］．Computer Physics Communications，1988，48（1）：1-10．

[40] SULSKY D， CHEN Z， SCHREYER H L ．A particle method for history-dependent materials［J］．Computer Methods in Applied Mechanics and Engineering，1994，118（1/2）：179-196．

[41] YORK 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（6）：901-924．

[42] WANG B， VARDON P J， HICKS M A，et al ．Development of an implicit material point method for geotechnical applications［J］．Computers and Geotechnics，2016，71：159-167．

[43] CAI F， UGAI K ．Numerical analysis of the stability of a slope reinforced with piles［J］．Soils and Foundations，2000，40（1）：73-84．

[44] 年廷凯，徐海洋，刘红帅 ．抗滑桩加固边坡三维数值分析中的几个问题［J］．岩土力学，2012，33（8）：2521-2526．

[45] NIAN T K， XU H Y， LIU H S ．Several issues in three-dimensional numerical analysis of slopes reinforced with anti-slide piles［J］．Rock and Soil Mechanics，2012，33（8）：2521-2526．（in Chinese）

[46] DUAN X R， ZHANG J ．Numerical method for slope stabilization design with piles［J］．IOP Conference Series：Earth and Environmental Science，2021，861（3）：032080．

[47] SKEMPTON A W ．Long-term stability of clay slopes［J］．Géotechnique，1964，14（2）：77-102．

[48] STARK T D， EID H T. Drained residual strength of cohesive soils［J］. Journal of Geotechnical Engineering，1994，120（5）：856-871．

[49] KAFAJI I K J. Formulation of a dynamic material point method （MPM） for geomechanical problems［M］．Stuttgart：Institut für Geotechnik der Universität Stuttgart，2013．

[50] BARDENHAGEN S G， BRACKBILL J U， SULSKY D ．Numerical study of stress distribution in sheared granular material in two dimensions［J］．Physical Review E，2000，62（3）： 3882-3890．

[51] ANDERSEN S. Material-point analysis of large-strain problems： modelling of landslides［D］. Aalborg University，Department of Civil Engineering，2009．

[52] YORK A R， SULSKY D， SCHREYER H L. The material point method for simulation of thin membranes［J］. International Journal for Numerical Nethods in Engineering，1999，44（10）：1429-1456．

[53] HU W， CHEN Z ．A multi-mesh MPM for simulating the meshing process of spur gears［J］．Computers & Structures，2003，81（20）：1991-2002．

[54] AUGARDE C E， LEE S J， LOUKIDIS D. Numerical modelling of large deformation problems in geotechnical engineering：a state-of-the-art review［J］. Soils and Foundations，2021，61（6）：1718-1735．

[55] LU M， CECCATO F， ZHOU M L，et al ．Evaluating the exceedance probability of the runout distance of rainfall-induced landslides using a two-stage FEM-MPM approach［J］. Acta Geotechnica， 2024， 19（6）： 3691-3706．

[56] ALSARDI A， YERRO A. Coseismic site response and slope instability using periodic boundary conditions in the material point method［J］. Journal of Rock Mechanics and Geotechnical Engineering， 2023， 15（3）： 641-658．

[57] ALSARDI A， COPANA J， YERRO A ．Modelling earthquake-triggered landslide runout with the material point method［J］．Proceedings of the Institution of Civil Engineers-Geotechnical Engineering， 2021， 174（5）： 563-576．

[58] DUAN X R， ZHANG J， LIU L L，et al ．Hybrid response surface method for system reliability analysis of pile-reinforced slopes［J］．Journal of Rock Mechanics and Geotechnical Engineering，2024，16（9）：3395-3406．

[59] ZHANG J， WU C G， TAN X H，et al ．Hierarchical response surface method for reliability analysis of a pile-slope system［J］．Canadian Geotechnical Journal，2023，60（4）： 397-409．

[60] LIU X， WANG Y， LI D Q ．Investigation of slope failure mode evolution during large deformation in spatially variable soils by random limit equilibrium and material point methods［J］．Computers and Geotechnics，2019，111：301-312．

[61] 梅伟，顾世祥，刘鑫，等 ．基于滑坡大变形过程的滑坡定量风险评估方法［J］．武汉大学学报（工学版），2022，55（5）：443-453．

[62] MEI W， GU S X， LIU X，et al ．Quantitative risk assessment of landslides based on large deformation process of landslide［J］．Engineering Journal of Wuhan University，2022，55（5）：443-453．（in Chinese）

[63] 刘鑫，王宇，李典庆 ．考虑土体参数空间变异性的边坡大变形破坏模式研究［J］．工程地质学报，2019，27（5）：1078-1084．

[64] LIU X， WANG Y， LI D Q. Slope failure modes at large deformation in spatially variable soils［J］. Journal of Engineering Geology，2019，27（5）：1078-1084．（in Chinese）

[65] 戴自航，徐祥 ．边坡抗滑桩设计计算的三维有限元法［J］．岩石力学与工程学报， 2012， 31（12）： 2572-2578．

[66] DAI Z H， XU X. 3D finite element method for design computations of anti-slide piles［J］. Chinese Journal of Rock Mechanics and Engineering， 2012， 31（12）： 2572-2578．（in Chinese）

[67] 冯文娟， 琚晓冬. 基于FLAC^3D的抗滑桩设计方法研究［J］．岩土工程学报， 2011， 33（）： 256-259．

[68] FENG W J， JU X D ．Design method for anti-slide piles based on FLAC^3D ［J］. Chinese Journal of Geotechnical Engineering， 2011，33（Sup.2）： 256-259．（in Chinese）

[69] 崔晓艳， 庄妍， 张希栋， 等. 循环荷载下桩承式路堤中土拱效应动力折减系数离散元研究［J］. 湖南大学学报（自然科学版），2022， 49（9）： 164-172．

[70] CUI X Y， ZHUANG Y， ZHANG X D，et al ．Discrete element analysis on dynamic reduction coefficient of soil arching effect in pile supported embankment under cyclic load［J］．Journal of Hunan University （Natural Sciences），2022，49（9）：164-172．（in Chinese）

基金资助

国家自然科学基金资助项目（42072302）; 国家重点研发计划项目（2021YFB2600504）; 国家资助博士后研究人员计划（GZB20240533）