Stability and Dynamic hazard Analysis of Dumpsites Based on MPM-DEM Coupling Algorithm
-
摘要:
随着山区高速公路建设的迅速推进，弃渣场地的稳定性及潜在失稳灾害评估日益受到重视。带有巨大能量的滑坡体可能会对沿程的结构物等造成冲击破坏，危害生命财产安全。充分利用物质点法（MPM）对连续介质大变形过程的模拟及离散元法（DEM）准确接触判断的优势，MPM-DEM 耦合算法可解决滑坡体与复杂地形和沿程结构物之间的相互作用问题。本文基于GPU并行高性能计算软件CoSim中的MPM-DEM耦合算法，实现了弃渣场边坡的稳定性及潜在失稳灾害动力学分析。首先以散粒体冲击结构物的算例验证了该算法的合理性和准确性，在此基础上，以云南某高速公路弃渣场为案例，计算了其稳定性系数，预测了其潜在失稳灾害的影响范围和危害程度。结果表明，该弃渣场边坡处于稳定状态，若发生失稳，会对下游高速公路桥桩产生巨大的冲击作用。数值模拟结果证明了该耦合算法在弃渣场边坡稳定性与失稳灾害动力学分析中具有强大的优势，可实现边坡“稳定性→大变形→流动→堆积”的全过程分析。
-
关键词:
- 物质点法(MPM) /
- 离散元法(DEM) /
- MPM-DEM耦合算法 /
- 稳定性分析 /
- 失稳评估
Abstract:The rapid expansion of highway construction in mountainous areas has heightened the focus on the stability of dumpsites and the assessment of potential failure disasters. High-energy landslides can severely impact structures along the way, posing significant risks to life and property. Leveraging the material point method (MPM) for simulating large deformations in continuous media and the discrete element method (DEM) for precise contact detection, the MPM-DEM coupling algorithm can effectively models interactions between landslides and complex terrains or structures. This study employs the MPM-DEM coupling algorithm within the GPU-accelerated high-performance computing software CoSim to analyze the stability and dynamic hazard potential of dumpsite slopes. The validity and accuracy of the proposed algorithm were initially verified through a case study involving granular material impacting a structure. Subsequently, a highway dumpsite in Yunnan Province was selected as an example to evaluate the stability factor and the predict the potential failure’s impact range and severity. Results indicate that the dumpsite slope is currently stable; however, should instability occur, it would exert significant impact on downstream highway bridge piles. The numerical simulation results demonstrate the strong capabilities of the coupling algorithm in the stability and dynamic failure analysis of dumpsite slopes, which can simulate the entire process of "stability-large deformation-flow-accumulation" for slopes.
-
0. 引言
山区高速公路建设不断推进，大量余泥渣土堆填形成弃渣场。由于渣土较为松散，其稳定性问题需要给予高度重视，一旦失稳会对周围的结构物产生破坏，对人民生命财产安全造成严重威胁，因此开展高速公路弃渣场稳定性及潜在失稳灾害分析具有重要意义。滑体失稳后的运动过程中，会受到复杂地形以及周围结构物的影响，这涉及到岩土体与相对刚性结构物之间的多体动力相互作用。为准确模拟滑体失稳后运动过程及评估灾害风险，需要采用一种有效的方法模拟刚柔多体间的相互作用问题。
随着计算机技术的发展，数值方法成为了边坡稳定性及灾害分析的重要手段[1]。常用的有限元法（FEM）将连续体离散化为若干个有限大小的单元体集合，以求解连续介质力学问题[2]。然而，FEM通常用于分析小变形问题，对于岩土体的大变形问题通常由于网格变形而造成收敛困难[3 − 4]。基于无网格的物质点法（MPM）[5]在模拟大变形问题时具有显著优势。谢艳芳等[6]基于MPM对新磨村滑坡进行了模拟，并通过分析滑体碰撞破碎等特征揭示了其动力演化过程；Li 等[7]利用MPM模拟了柔性散粒体材料冲击刚性挡板的物理试验过程，并开展了泥石流冲击挡土墙的数值计算；Cuomo 等[8]基于MPM接触算法进行了滑坡冲击建筑物结构的模拟，并分析了不同滑体材料和冲击速度造成结构破坏的响应程度。虽然MPM已经在滑坡灾害模拟中成功应用，但是基于背景网格结点的接触判断无法保证接触精度，存在虚假接触现象，而且针对地形的物质点离散会增加计算规模，占用大量内存并降低计算效率[9]。
与MPM相比，离散元法（DEM）[10]作为一种典型的非连续介质方法，其优势在于可以描述复杂模型边界并计算多体系统接触力[11 − 13]。刘广煜等[14]提出了一种基于块体DEM的黏结断裂模型，并系统研究了西藏易贡滑坡中崩塌区崩解启动到堰塞体形成的全过程；Shen等[15]基于DEM 研究了泥石流对刚性屏障的影响，并分析了撞击过程中颗粒流能量和动量的变化过程；Basinskas等[16]采用DEM模拟了带有复杂边界的搅拌机中颗粒的流动过程，并研究了颗粒数量和转速对搅拌程度的影响。然而，由于DEM接触检测计算成本较高，难以模拟具有大量不规则形状颗粒组成的大规模真实颗粒材料[17]。
结合MPM处理连续材料大变形行为的能力和DEM接触搜索在处理多体相互作用中的优势，MPM-DEM耦合算法可准确模拟可变形体与刚性复杂边界的相互作用机理。而目前的耦合算法大多基于CPU架构，三维真实尺度下的滑坡精细化大规模模拟将受到计算规模和效率的限制。本文将基于GPU加速的耦合模拟软件CoSim中的MPM-DEM耦合算法[18]实现了高速公路弃渣场的稳定性及潜在失稳风险分析。首先对MPM和DEM基本理论及基于GPU加速的MPM-DEM耦合算法进行介绍，然后通过散粒体冲击结构物的典型算例，证明了采用MPM-DEM耦合算法进行刚柔接触分析的可行性。最后以云南某高速公路弃渣场地边坡为例，研究了其稳定性及潜在的失稳灾害动力过程，从而论证了本文方法在弃渣场地边坡稳定性及灾害分析中的准确性，同时为类似工程的选址设计及风险分析防治提供依据。
1. MPM-DEM耦合算法
1.1 MPM
作为拉格朗日-欧拉混合计算方法，物质点法MPM将连续体离散为一系列拉格朗日物质点，这些点在欧拉背景网格中移动。每个物质点代表一块子域并存储所有状态信息，包括质量、速度、应力、应变等，背景网格用于求解动量方程，通过形函数实现与物质点的信息交互[19 − 20]。
图1显示了MPM计算流程。在单步MPM计算循环中，首先将物质点信息映射到网格结点上（图1a）：
![]()
图 1 物质点法计算步骤：（a）物质点映射网格P2G；（b）网格结点更新；（c）网格映射物质点G2P；（d）物质点更新Figure 1. Calculation procedure of the MPM: (a) particle to grid (P2G); (b) grid node update; (c) grid to particle (G2P); (d) material particle update$$ {m_n} = \sum\limits_p {{m_p}{N_{p - n}}} $$ (1) $$ {v_n} = \frac{1}{{{m_n}}}\sum\limits_p {{m_p}{v_p}{N_{p - n}}} $$ (2) 式中：下标p、n——物质点和网格结点；
${m_p}$、${m_n}$——质量；
${v_p}$、${v_n}$——速度；
${N_{p - n}}$——形函数。
之后由结点速度计算物质点的应变增量$d{\varepsilon _p}$：
$$ d{\varepsilon _p} = \frac{1}{2}dt\sum\limits_n {[{v_n}\nabla {N_{n - p}} + {{\left( {{v_n}\nabla {N_{n - p}}} \right)}^T}]} $$ (3) 式中：$\nabla {N_{n - p}}$——形函数梯度；
$dt$——时间步长。
根据材料本构模型更新物质点应力${\sigma _p}$，并将其插值到网格上得到结点内力$f_n^{int}$，同时根据物质点的体力$ f_p^{body} $插值得到结点外力$f_n^{ext}$：
$$ f_n^{int} = \sum\limits_p {{\sigma _p}{V_p}\nabla {N_{p - n}}} $$ (4) $$ f_n^{ext} = \sum\limits_p {f_p^{body}{N_{p - n}}} $$ (5) 式中：${V_p}$——物质点体积。
计算总结点力$f_n^{tot}$，并积分计算结点动量${P_n}$（图1b）：
$$ f_n^{tot} = f_n^{int} + f_n^{ext} $$ (6) $$ P_n^{t + dt} = P_n^t + f_n^{tot}dt $$ (7) 然后将更新后的结点信息映射回物质点（图1c），并更新物质点的速度${v_p}$和位置${x_p}$（图1d）：
$$ v_p^{t + dt/2} = v_p^{t - dt/2} + dt\mathop \sum \limits_n \frac{{f_n^{tot}}}{{{m_n}}}{N_{p - n}} $$ (8) $$ x_p^{t + dt} = x_p^t + dt\mathop \sum \limits_n \frac{{P_n^{t + 1/2}}}{{{m_n}}}{N_{p - n}} $$ (9) 单步计算完毕，将背景网格重置为初始状态：
$$ {P_n} = 0,f_n^{tot} = 0 $$ (10) 1.2 MPM-DEM耦合算法
MPM-DEM耦合算法可充分利用MPM模拟连续介质大变形问题的能力及DEM在复杂系统中准确计算接触力的优势，因而在模拟颗粒和刚体相互作用方面具有很高的应用价值。与MPM传统接触算法[21]相比，MPM-DEM没有额外的内存开销，并且采用了更精确的接触模型来处理多体交互系统。
对于物质点与复杂边界的相互作用问题，采用物质点-三角面局部接触计算方法，即将DEM中球颗粒-三角面之间的接触检测和接触力计算引入到 MPM中。图2显示了基于MPM-DEM耦合算法的物质点与复杂边界相互作用示意图。
![]()
图 2 基于MPM-DEM耦合算法的物质点与复杂边界相互作用示意图Figure 2. Schematic diagram of interaction between material points and complex boundaries based on the MPM-DEM coupling algorithm该耦合算法的基本思想可以概括如下：
（1）定义局部边界背景网格
根据刚性边界模型的三角面单元，与MPM全局背景网格进行相交检测判断，将被三角形切割的网格单元标记为局部边界网格单元。如果物质点运动至该单元内，则将其转化为DEM球颗粒，并将物质点的位置和速度传递给球颗粒。
（2）基于DEM的球-三角面接触检测
计算球颗粒球心与边界三角面的距离d，球半径为r，定义接触深度$ {u^n} = r - d $，以此检查球颗粒与边界三角面之间的接触情况。若$ {u^n} > 0 $，则证明已经发生接触，需要计算球颗粒与刚性边界之间的接触力。
（3）基于DEM的球-三角面接触力计算
如图3所示，球颗粒受到的接触力合力$ \overrightarrow {{F^C}} $可分解为法向接触力$ \overrightarrow {{F^{{C_{\text{n}}}}}} $和切向接触力 $ \overrightarrow {{F^{{C_{\text{s}}}}}} $：
![]()
图 3 球与三角面接触示意图Figure 3. Schematic Diagram of contact between a sphere and a triangular surface$$ \overrightarrow {{F^C}} = \overrightarrow {{F^{{C_{\text{n}}}}}} + \overrightarrow {{F^{{C_{\text{s}}}}}} $$ (11) $$ \overrightarrow {{F^{{C_{\text{n}}}}}} = {K_n}{u^n}\overrightarrow n $$ (12) $$ \overrightarrow {{F^{{C_{\text{s}}}}}} = {K_s}\overrightarrow {\Delta {u^s}} + {\{ \overrightarrow {{F^{{C_{\text{s}}}}}} \} _{update}} $$ (13) 式中：$ {u^n} $——法向接触深度；
$ \overrightarrow {\Delta {u^s}} $——剪切位移增量；
$ \overrightarrow n $——接触法线方向；
$ {\{ \overrightarrow {{F^{{C_{\text{s}}}}}} \} _{update}} $——上一步的切向力；
$ {K_n} $——法向刚度；
$ {K_s} $——切向刚度，由材料和几何特性计算：
$$ {K_n} = 1/\left(\frac{1}{{{E_A}{r_A}}} + \frac{1}{{{E_B}{r_B}}}\right) $$ (14) $$ {K_s} = 1/\left(\frac{1}{{{E_A}{\nu _A}{r_A}}} + \frac{1}{{{E_B}{\nu _B}{r_B}}}\right) $$ (15) 式中：$ {E_{A/B}} $、$ {\nu _{A/B}} $、$ {r_{A/B}} $——两相交单元A和B的杨氏模 量、接触剪切与法向刚度之 比、颗粒半径。
对于无半径的三角面，由与其相交的球半径替代。
球与三角面发生接触时，二者会存在相交平面截面圆，不同的接触关系产生的截面圆形态不同。基于CoSim-DEM已有的接触模型[22]，共考虑了11种不同的接触情况。因此，引入面积系数$ {S _{ratio}} $：
$$ {S _{ratio}} = {S _{in}}/{S _c} $$ (16) 式中：$ {S _{in}} $——截面圆与当前三角形重叠部分的面积；
$ {S _c} $—截面圆面积，图3所示情况为完全接触， $ {S _{ratio}} $=1。
接触力$ {F^C} $由每个分力乘以面积系数$ {S _{ratio}} $求和计算：
$$ \overrightarrow {{F^C}} = \sum {{S _{ratio}}} (\overrightarrow {{F^{{C_{\text{n}}}}}} + \overrightarrow {{F^{{C_{\text{s}}}}}} ) $$ (17) （4）MPM求解计算
球颗粒接触力更新后，将球颗粒转回物质点并将接触力作为体力施加到物质点上，以此完成MPM-DEM的耦合信息传输，之后回到MPM计算循环中进行MPM求解计算。
$$ f_{MPM}^{body} = F_{DEM}^C $$ (18) 2. MPM-DEM算例验证
为验证MPM-DEM耦合算法在颗粒材料与复杂边界或刚体系统相互作用中的可行性和准确性，模拟了散粒体材料对结构物的冲击影响，如图4所示。在L=0.8 m，W=0.1 m，H=0.5 m的容器左侧放置砂样，砂样模型的尺寸为0.2 m×0.1 m×0.3 m，被离散为392,818个物质点，背景网格设定为0.005 m。在距离砂样右侧d=0.3 m的位置放置一个圆柱形结构物，尺寸为h=0.2 m，r=0.025 m，并利用三角面进行剖分。材料参数设置见表1，基于砂土料参数设置MPM材料参数，MPM采用Drucker-Prager（DP）弹塑性本构模型[23]，基于参数反演设置DEM边界参数，接触类型设置为粘弹性接触模型[22]。物质点系统的阻尼系数设为0.02。
![]()
图 4 散粒体冲击结构物模型示意图Figure 4. schematic diagram of the granular material impact model on structure表 1 散粒体冲击结构物模拟材料参数Table 1. Material parameters for granular material impact simulation on structureDEM边界参数 MPM材料参数 弹性模量/MPa 100 弹性模量/MPa 10 泊松比 0.2 泊松比 0.2 静摩擦角/（°） 30 内摩擦角/（°） 38 动摩擦角/（°） 20 粘聚力/kPa 0 回弹系数 0.1 密度/（kg·m−3） 2000 模拟开始后，砂样首先在重力作用下形成初始应力场，然后迅速释放右侧挡板，砂样开始坍塌，流动下滑并冲击结构物。图5显示了颗粒流动过程中与结构物的相互作用情况及其速度分布结果。模型右上侧颗粒因重力而向下坍塌，并不断向前流动；在t=0.3 s时，流动颗粒的最前缘与结构物发生接触，而后颗粒流的前缘中间部分被结构筒壁阻挡，颗粒流绕开结构物沿两侧通道继续滑动行进，t=0.4 s时流速已明显减慢；由于结构物的阻挡作用，结构物左侧颗粒堆积高度不断增加，部分颗粒继续以较小速度绕结构物两侧向前运动，并逐渐在结构物右侧汇集，整个流动过程在t=2.0 s时停止。
![]()
图 5 不同时刻颗粒流形态和速度分布Figure 5. Particle flow morphology and velocity distribution at different time steps不同时刻颗粒流速度可以直观地反映结构物的缓冲和阻挡作用，同时颗粒流对结构的冲击力分析对于评估防冲击结构的稳定性和安全性也很重要。图6显示了结构物所受冲击力的演化曲线。整个冲击过程可分为3个阶段：（1）加速阶段：由于颗粒流前缘速度较高，结构物受冲击效果较强，冲击力迅速增大。点A时刻，冲击效应集中分布在结构底部，随着堆积高度增加，冲击力分布不断向上扩展，自下而上缓慢减弱，点B时刻的冲击力达到峰值21.6 N。（2）冲击阶段：在结构物的阻碍作用下，颗粒流动速度减慢，冲击作用减弱。结构物冲击力的分布范围和数值开始减小。但曲线上有一个小的反弹（点C），这是由于初始模型右上部位的颗粒到达并再次撞击结构体（图5 t=0.5 s），产生了微弱的二次冲击作用。（3）减速阶段：t=0.8 s后，颗粒流速度已经很小，冲击力曲线趋于平缓，结构体受力来自左侧堆积物的土压力，受力分布呈现随着高度增加而逐渐减小的规律（点D）。
![]()
图 6 结构物冲击力演化曲线及其受力分布Figure 6. Evolution curve and stress distribution of impact force on the structure3. 弃渣场工程地质概况
该弃渣场位于云南省玉溪市，位于山坳沟谷地带，属于沟谷型弃渣场，于2017年9月开始堆渣，2019年堆渣完成。堆渣区域沟谷平面上呈“C”形，沟心高程
1802 m～1810 m。沟谷左侧山顶高程2102 m，山坡坡度20°～40°；右侧山体低矮，顶部高程1820 m，山坡坡度30°。弃渣场堆土量约30×104 m3，堆渣高度23.10 m，顶部平台高程1824 m。渣体北东侧沿山坡堆砌，其余三个侧面按1∶2放坡至原地面，堆渣边坡脚设置高度为8 m的混凝土挡墙。挡墙下游7 m为澄川高速桥梁，桥梁下部结构为柱式墩、桩基础；弃土场顶部南东侧约20 m为养白牛社区居民房，周边环境相对复杂（图7）。弃渣场上覆地层主要为第四系全新统人工堆积层（${\mathrm{Q}}_4^{{\mathrm{ml}}} $）渣土，下伏基岩为侏罗系中统张家河组（J2z）泥岩（图8）。渣土主要为碎石土，以全、强风化泥岩为主，含少量黏性土，平均厚度15.38 m。根据弃渣场实地勘察及室内试验，两类岩土的物理力学参数如表2所示。
![]()
图 8 弃渣场地工程地质特征：（a）工程地质平面图；（b）1-1’工程地质剖面图Figure 8. Engineering geological characteristics of the dumpsite: (a) Engineering geological map; (b) Engineering geological cross-section of 1-1’表 2 弃渣场各地层岩土物理力学参数Table 2. Physical and mechanical parameters of geomaterials in different layers岩土体类型 杨氏模量
/MPa泊松比 容重
/（kN·m−3）粘聚力
/kPa内摩擦角
/（°）渣土 200 0.2 19.5 6 26 强风化泥岩 2000 0.2 22.0 50 40 4. 弃渣场稳定性及灾害动力学分析
4.1 数值计算模型
为了研究弃渣场稳定性及其失稳后危害程度，为相关工程选址及设计提供借鉴，基于MPM-DEM耦合算法模拟分析了上述弃渣场的稳定性及潜在失稳灾害。根据弃渣场工程地质特征，建立了该边坡的数值计算模型（图9）。弃渣土体采用MPM，共被划分为870,707个物质点，背景网格设置为2.5 m×2.5 m×2.5 m；地形面及与下覆的基岩接触面采用三角网表征，共划分为
92019 个三角面；并根据场地周边的高速公路实际情况，建立了32根桥桩和一个水池构筑物，桩体和水池构筑物均采用三角网来表征。计算中弃渣土体与地形、基岩接触面、桥桩基及水池构筑物间的接触采用MPM-DEM耦合方法；MPM采用Drucker-Prager（DP）弹塑性本构模型；弃渣土体与下覆基岩面间的接触采用DEM的粘结摩擦模型。基于场地地质资料设置模拟所用参数如表3所示。表 3 弃渣场数值计算材料参数Table 3. Material parameters for numerical simulation of the dumpsiteDEM边界参数 MPM材料参数 弹性模量/MPa 2000 弹性模量/MPa 200 泊松比 0.2 泊松比 0.2 静摩擦角/（°） 30 内摩擦角/（°） 26 动摩擦角/（°） 20 粘聚力/kPa 6 法向粘聚力/kPa 50 密度/（kg·m−3） 1950 切向粘聚力/kPa 50 — — 4.2 稳定性分析
为了对弃渣土体的稳定性进行定量分析，在计算中引入了强度折减法。强度折减法计算边坡稳定性的基本原理是通过不断降低岩土体的抗剪强度参数（粘聚力和内摩擦角），使岩土体达到破坏状态，此时的折减系数为边坡的稳定性系数[23]。图10显示了边坡表面三个典型监测点P1、P2、P3（图9）的位移随折减系数Rf增加的变化曲线。由图得知，随着Rf增加，位移变化可划分为两个阶段：初始缓慢变形阶段，位移随Rf缓慢变化，基本呈线性增长趋势；当Rf大于1.25后，位移进入快速发展阶段，监测点位移出现明显拐点，呈现急剧增加趋势。由上述分析得知，当Rf=1.25时，弃渣场地边坡位移发生突变，达到破坏状态，因而其稳定性系数Fos=1.25。
![]()
图 10 监测点 P1, P2, P3的位移随Rf的变化曲线Figure 10. Displacement curves of monitoring points P1, P2, and P3 with respect to changes in Rf由数值计算得到了不同Rf下的边坡变形破坏状态，图11展示了不同Rf下弃渣土体A-A’剖面的等效塑形应变云图和三维弃渣场的位移云图。当Rf<Fos=1.25时，边坡位移值较小，随着Rf的增加，边坡位移缓慢增加，边坡塑性区率先在顶部平台边缘处发展；当Rf =Fos=1.25时，边坡整体位移增大，位移增速变快，坡顶塑性区范围不断扩大，同时坡脚处也出现较大塑性应变，塑性区即将贯通整个坡体；当Rf >Fos=1.25时，坡体岩土体位移急剧增加，坡顶和坡脚处均产生较大位移，部分岩土体滑塌，同时塑性剪切破坏带贯通，滑动面近似呈圆弧形，滑体前缘位于坡脚处，后缘延伸至顶部平台，此时坡体已进入失稳破坏状态。由上述分析得知，该堆填场地边坡在天然状态下处于稳定状态；随着强度不断折减，边坡抗滑力减小，位移和塑性应变均不断增大，在重力作用下顶部平台边缘处和前缘坡脚处均产生较大变形，表现为后缘推移和前缘牵引相结合的复合型滑坡破坏模式。
![]()
图 11 不同Rf对应的边坡等效塑性应变与位移结果Figure 11. Equivalent plastic strain and displacement results of the slope under different Rf values4.3 灾害动力学分析
作为一种人工松散堆积体，弃渣场受到降雨或其他因素时强度会降低，易产生滑坡、泥石流等灾害，而对下游人民群众的生命财产安全构成威胁。为了预测弃渣场滑坡的影响范围和危害程度，基于上述算法模拟分析了弃渣场潜在的失稳滑动过程。在Rf >Fos=1.25时边坡失稳，破坏方式为渐进破坏，逐渐失稳形成了大规模滑坡，并对下游的高速公路桥桩产生冲击。图12显示了滑坡体的失稳运动过程及其影响范围。
![]()
图 12 滑坡体失稳滑动过程：（a）5 s；（b）15 s；（c）30 s；（d）60 sFigure 12. Landslide Body Instability and sliding process: (a) 5 s; (b) 15 s; (c) 30 s; (d) 60 s在5 s后（图12 a），边坡进入大变形发展阶段，边坡失稳，沿滑动面整体向下滑动；在t=15 s时（图12 b），处于临空面的滑体前缘已经剪出，顺着地形下滑，并对靠近渣场一侧的公路桥桩产生冲击，此时前缘已形成碎屑流，表现出一定的流动性，同时滑体后缘失去支撑后变形加剧形成牵引拉伸破坏；在t=30 s时（图12 c），碎屑流沿沟道流动，影响范围进一步扩大，并冲击到远离渣场一侧的桥桩和沟道下游的桥桩，同时与附近的构筑物产生接触撞击，此时最大滑距已达67.8 m；在t=60 s时（图12 d），碎屑流绕过桥桩向下游流动，构筑物已基本被掩埋，受影响桥桩增多，此后滑体滑动速度逐渐减缓并堆积到沟道中，最大滑距达104.8 m。
滑坡运动过程中，高速公路的多根桥桩受到了较大影响。在图9中选取了三根桥桩对其受到的冲击力进行了实时监测，图13显示了这三根桥桩各自受到的总撞击力随时间的变化曲线。由图可知，桥桩均受到了较大的冲击力。以1号桩为例，在10 s左右，滑坡体已经下滑到该桩附近，桩开始受到岩土体的冲击碰撞作用；此后岩土体不断下滑，大量颗粒连续冲击该桩，导致所受冲击力有所波动，同时从整体趋势来看，冲击力值不断增加，在70 s左右达到峰值；此后碎屑流流速变慢，冲击力有下降趋势，冲击作用减弱，到80 s后趋于稳定，整个流动过程已基本完成，总冲击力值趋于静止土压力值。此外，比较三根桥桩的冲击力曲线，可发现2号桩的冲击力相比1号桩有一定的滞后性，而且3号桩更加滞后，这是因为边坡失稳下滑后，岩土体率先撞击到相对靠近的1号桩，其次是2号桩和3号桩。同时，1号桩和2号桩的冲击力峰值均在1.1×106 N附近，3号桩冲击力的峰值约为8×105 N，低于1号桩和2号桩，这是因为3号桩在远离边坡的一侧，靠近边坡一侧的桩对下滑岩土体有一定的阻挡作用，致使岩土体对远离边坡桥桩的冲击作用相对弱化。
![]()
图 13 桥桩冲击力演化曲线：（a）1号桩；（b）2号桩；（c）3号桩Figure 13. Evolution curves of total impact force on bridge piles: (a) pile 1; (b) pile 2; (c) pile 3综合上述分析，弃渣场若失稳发生滑坡灾害，影响范围涉及高速公路桥梁的多根桥桩及周边的构筑物。滑坡体会对桥桩产生巨大的冲击作用，足以使桥桩断裂破坏，从而使高速公路桥梁坍塌，同时碎屑流在流动过程中，会掩埋影响范围内的构筑物及河道。因此，弃渣场的潜在失稳灾害会对人民群众的生命财产安全造成巨大的威胁。
5. 结 论
随着山区基础建设的不断发展，弃渣场作为一种特殊的边坡，其结构松散、稳定性差，易形成滑坡灾害，因而对周边结构物、居民等存在威胁。实现弃渣场地边坡的稳定性分析，及极端工况下的灾害动力学预测与影响范围分析，对于弃渣场地的选择与设计具有重要的意义。物质点法（MPM）可高效模拟连续介质的大变形问题，离散元法（DEM）在处理多体接触作用中优势显著，因而MPM-DEM 耦合算法可较好地解决连续介质与复杂边界之间的多体接触问题。本文基于GPU加速的耦合模拟软件CoSim中的MPM-DEM耦合算法开展相关研究，实现了弃渣场的稳定性及失稳动力学分析，得到如下结论：
（1）MPM-DEM算法的关键点在于将边界背景网格单元中的物质点激活为DEM球颗粒实现耦合信息交互，从而克服了传统MPM接触算法的精细化模型边界描述困难及接触检测不准确的缺点。通过典型的散粒体冲击结构物的算例验证了该耦合算法在进行颗粒物质大变形及其与复杂边界之间的碰撞接触机理研究和动力学特性分析方面的可靠性。
（2）针对云南某高速公路弃渣场工程，本文构建了其三维精细化数值计算模型，使用MPM-DEM耦合算法模拟分析了其稳定性和潜在失稳滑动过程。结果表明，该边坡处于稳定状态，稳定性系数为1.25，随着强度折减，表现为后缘推移和前缘牵引相结合的复合型滑坡破坏模式。在Rf >Fos=1.25时边坡失稳，破坏方式为渐进破坏，对下游的高速公路桥桩产生巨大的冲击作用，冲击力最高达百万牛顿；与此同时，对极端工况下堆填场地的颗粒流动和淹没范围进行了分析和评价。
（3）MPM-DEM耦合算法有效地分析了边坡稳定性，预测了弃渣场滑坡的影响范围和危害程度，因此该算法在边坡稳定性及失稳动力学研究中具有较好优势，可实现边坡的“稳定性→大变形→流动→堆积”的全过程分析，从而为工程选址设计及灾害风险防治提供理论和技术支持。
-
![]()
图 1 物质点法计算步骤：（a）物质点映射网格P2G；（b）网格结点更新；（c）网格映射物质点G2P；（d）物质点更新
Figure 1. Calculation procedure of the MPM: (a) particle to grid (P2G); (b) grid node update; (c) grid to particle (G2P); (d) material particle update
![]()
图 2 基于MPM-DEM耦合算法的物质点与复杂边界相互作用示意图
Figure 2. Schematic diagram of interaction between material points and complex boundaries based on the MPM-DEM coupling algorithm
![]()
图 3 球与三角面接触示意图
Figure 3. Schematic Diagram of contact between a sphere and a triangular surface
![]()
图 4 散粒体冲击结构物模型示意图
Figure 4. schematic diagram of the granular material impact model on structure
![]()
图 5 不同时刻颗粒流形态和速度分布
Figure 5. Particle flow morphology and velocity distribution at different time steps
![]()
图 6 结构物冲击力演化曲线及其受力分布
Figure 6. Evolution curve and stress distribution of impact force on the structure
![]()
图 8 弃渣场地工程地质特征：（a）工程地质平面图；（b）1-1’工程地质剖面图
Figure 8. Engineering geological characteristics of the dumpsite: (a) Engineering geological map; (b) Engineering geological cross-section of 1-1’
![]()
图 10 监测点 P1, P2, P3的位移随Rf的变化曲线
Figure 10. Displacement curves of monitoring points P1, P2, and P3 with respect to changes in Rf
![]()
图 11 不同Rf对应的边坡等效塑性应变与位移结果
Figure 11. Equivalent plastic strain and displacement results of the slope under different Rf values
![]()
图 12 滑坡体失稳滑动过程：（a）5 s；（b）15 s；（c）30 s；（d）60 s
Figure 12. Landslide Body Instability and sliding process: (a) 5 s; (b) 15 s; (c) 30 s; (d) 60 s
![]()
图 13 桥桩冲击力演化曲线：（a）1号桩；（b）2号桩；（c）3号桩
Figure 13. Evolution curves of total impact force on bridge piles: (a) pile 1; (b) pile 2; (c) pile 3
表 1 散粒体冲击结构物模拟材料参数
Table 1 Material parameters for granular material impact simulation on structure
DEM边界参数 MPM材料参数 弹性模量/MPa 100 弹性模量/MPa 10 泊松比 0.2 泊松比 0.2 静摩擦角/（°） 30 内摩擦角/（°） 38 动摩擦角/（°） 20 粘聚力/kPa 0 回弹系数 0.1 密度/（kg·m−3） 2000 
下载: 导出CSV
表 2 弃渣场各地层岩土物理力学参数
Table 2 Physical and mechanical parameters of geomaterials in different layers
岩土体类型 杨氏模量
/MPa泊松比 容重
/（kN·m−3）粘聚力
/kPa内摩擦角
/（°）渣土 200 0.2 19.5 6 26 强风化泥岩 2000 0.2 22.0 50 40
下载: 导出CSV
表 3 弃渣场数值计算材料参数
Table 3 Material parameters for numerical simulation of the dumpsite
DEM边界参数 MPM材料参数 弹性模量/MPa 2000 弹性模量/MPa 200 泊松比 0.2 泊松比 0.2 静摩擦角/（°） 30 内摩擦角/（°） 26 动摩擦角/（°） 20 粘聚力/kPa 6 法向粘聚力/kPa 50 密度/（kg·m−3） 1950 切向粘聚力/kPa 50 — —
下载: 导出CSV
-
[1] 王崇敬,张龙,刘国伟. 模糊数学评判和数值模拟相结合的土质边坡稳定性综合评价[J]. 中国地质灾害与防治学报,2023,34(6):69 − 76. [WANG Chongjing,ZHANG Long,LIU Guowei. Integrated assessment of soil cutting slope stability using fuzzy mathematics and numerical simulation[J]. The Chinese Journal of Geological Hazard and Control,2023,34(6):69 − 76. (in Chinese with English abstract)] WANG Chongjing, ZHANG Long, LIU Guowei. Integrated assessment of soil cutting slope stability using fuzzy mathematics and numerical simulation［J］. The Chinese Journal of Geological Hazard and Control, 2023, 34（6）: 69 − 76. （in Chinese with English abstract）
[2] FARIAS M M,NAYLOR D J. Safety analysis using finite elements[J]. Computers and Geotechnics,1998,22(2):165 − 181. DOI: 10.1016/S0266-352X(98)00005-6
[3] BHANDARI T,HAMAD F,MOORMANN C,et al. Numerical modelling of seismic slope failure using MPM[J]. Computers and Geotechnics,2016,75:126 − 134. DOI: 10.1016/j.compgeo.2016.01.017
[4] 孙玉进,宋二祥. “12•20” 深圳滑坡动态模拟[J]. 岩土工程学报,2018,40(3):441 − 448. [SUN Yujin,SONG Erxiang. Dynamic simulation of “12•20” Shenzhen landslide[J]. Chinese Journal of Geotechnical Engineering,2018,40(3):441 − 448. (in Chinese with English abstract)] DOI: 10.11779/CJGE201803007 SUN Yujin, SONG Erxiang. Dynamic simulation of “12•20” Shenzhen landslide［J］. Chinese Journal of Geotechnical Engineering, 2018, 40（3）: 441 − 448. （in Chinese with English abstract） DOI: 10.11779/CJGE201803007
[5] 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.
[6] 谢艳芳,李新坡,赵曙熙,等. 基于物质点法的新磨村滑坡动力特性分析[J]. 山地学报,2018,36(4):589 − 597. [XIE Yanfang,LI Xinpo,ZHAO Shuxi,et al. MPM-based numerical analysis of the kinematic characteristics of xinmo landslide in Maoxian county,Sichuan,China[J]. Mountain Research,2018,36(4):589 − 597. (in Chinese with English abstract)] XIE Yanfang, LI Xinpo, ZHAO Shuxi, et al. MPM-based numerical analysis of the kinematic characteristics of xinmo landslide in Maoxian county, Sichuan, China［J］. Mountain Research, 2018, 36（4）: 589 − 597. （in Chinese with English abstract）
[7] LI Xinpo,XIE Yanfang,GUTIERREZ M. A soft–rigid contact model of MPM for granular flow impact on retaining structures[J]. Computational Particle Mechanics,2018,5(4):529 − 537. DOI: 10.1007/s40571-018-0188-5
[8] CUOMO S,DI PERNA A,MARTINELLI M. MPM modelling of buildings impacted by landslides[M]//ICL Contribution to Landslide Disaster Risk Reduction. Cham:Springer International Publishing,2020:245-266.
[9] 廉艳平,张帆,刘岩,等. 物质点法的理论和应用[J]. 力学进展,2013,43(2):237 − 264. [LIAN Yanping,ZHANG Fan,LIU Yan,et al. Material point method and its applications[J]. Advances in Mechanics,2013,43(2):237 − 264. (in Chinese with English abstract)] DOI: 10.6052/1000-0992-12-122 LIAN Yanping, ZHANG Fan, LIU Yan, et al. Material point method and its applications［J］. Advances in Mechanics, 2013, 43（2）: 237 − 264. （in Chinese with English abstract） DOI: 10.6052/1000-0992-12-122
[10] CUNDALL P A,STRACK O D L. A discrete numerical model for granular assemblies[J]. Geotechnique,1979,29(1):47 − 65. DOI: 10.1680/geot.1979.29.1.47
[11] 陈宇,沈位刚,宋忠友,等. 土垫层缓冲落石冲击力特性离散元数值模拟分析[J]. 中国地质灾害与防治学报,2024,35(2):90 − 97. [CHEN Yu,SHEN Weigang,SONG Zhongyou,et al. Analysis of soil cushion buffering characteristic for rockfall impact force through discrete element numerical simulation[J]. The Chinese Journal of Geological Hazard and Control,2024,35(2):90 − 97. (in Chinese with English abstract)] CHEN Yu, SHEN Weigang, SONG Zhongyou, et al. Analysis of soil cushion buffering characteristic for rockfall impact force through discrete element numerical simulation［J］. The Chinese Journal of Geological Hazard and Control, 2024, 35（2）: 90 − 97. （in Chinese with English abstract）
[12] 臧佳园,常文斌,邢爱国,等. 含构造节理的崩塌体动力破碎特征[J]. 中国地质灾害与防治学报,2024,35(3):1 − 11. [ZANG Jiayuan,CHANG Wenbin,XING Aiguo,et al. Dynamic fragmentation characteristics of rock avalanche with tectonic joints[J]. The Chinese Journal of Geological Hazard and Control,2024,35(3):1 − 11. (in Chinese with English abstract)] ZANG Jiayuan, CHANG Wenbin, XING Aiguo, et al. Dynamic fragmentation characteristics of rock avalanche with tectonic joints［J］. The Chinese Journal of Geological Hazard and Control, 2024, 35（3）: 1 − 11. （in Chinese with English abstract）
[13] SINGER V,SAUTTER K B,LARESE A,et al. A partitioned material point method and discrete element method coupling scheme[J]. Advanced Modeling and Simulation in Engineering Sciences,2022,9(1):16. DOI: 10.1186/s40323-022-00229-5
[14] 刘广煜,徐文杰,佟彬,等. 基于块体离散元的高速远程滑坡灾害动力学研究[J]. 岩石力学与工程学报,2019,38(8):1557 − 1566. [LIU Guangyu,XU Wenjie,TONG Bin,et al. Study on dynamics of high-speed and long Run-out landslide hazards based on block discrete element method[J]. Chinese Journal of Rock Mechanics and Engineering,2019,38(8):1557 − 1566. (in Chinese with English abstract)] LIU Guangyu, XU Wenjie, TONG Bin, et al. Study on dynamics of high-speed and long Run-out landslide hazards based on block discrete element method［J］. Chinese Journal of Rock Mechanics and Engineering, 2019, 38（8）: 1557 − 1566. （in Chinese with English abstract）
[15] SHEN Weigang,ZHAO Tao,ZHAO Jidong,et al. Quantifying the impact of dry debris flow against a rigid barrier by DEM analyses[J]. Engineering Geology,2018,241:86 − 96. DOI: 10.1016/j.enggeo.2018.05.011
[16] BASINSKAS G,SAKAI M. Numerical study of the mixing efficiency of a ribbon mixer using the discrete element method[J]. Powder Technology,2016,287:380 − 394. DOI: 10.1016/j.powtec.2015.10.017
[17] JIANG Yupeng,ZHAO Yidong,CHOI C E,et al. Hybrid continuum–discrete simulation of granular impact dynamics[J]. Acta Geotechnica,2022,17(12):5597 − 5612. DOI: 10.1007/s11440-022-01598-2
[18] FENG Zekang,XU Wenjie,ULLAH JAN KHAN K. A GPU based Hybrid Material point and Discrete element method (MPDEM) algorithm and validation[J]. Computers and Geotechnics,2023,159:105462. DOI: 10.1016/j.compgeo.2023.105462
[19] SULSKY D,ZHOU Shijian,SCHREYER H L. Application of a particle-in-cell method to solid mechanics[J]. Computer Physics Communications,1995,87(1/2):236 − 252.
[20] SULSKY D,SCHREYER H L. Axisymmetric form of the material point method with applications to upsetting and Taylor impact problems[J]. Computer Methods in Applied Mechanics and Engineering,1996,139(1/2/3/4):409 − 429.
[21] 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.
[22] ZHOU Qian,XU Wenjie,LIU Guangyu. A contact detection algorithm for triangle boundary in GPU-based DEM and its application in a large-scale landslide[J]. Computers and Geotechnics,2021,138:104371. DOI: 10.1016/j.compgeo.2021.104371
[23] FENG Zekang,XU Wenjie. GPU material point method (MPM) and its application on slope stability analysis[J]. Bulletin of Engineering Geology and the Environment,2021,80(7):5437 − 5449. DOI: 10.1007/s10064-021-02265-8
计量
- 文章访问数: 32
- HTML全文浏览量: 21
- PDF下载量: 7
邮件订阅
RSS