基于物质点法的白格滑坡堰塞坝形成过程数值模拟

钟启明 吴昊 单熠博 赵鲲鹏

摘要：为准确模拟堰塞坝的形成过程，在物质点方法中引入了滑动速度依赖型的摩擦系数模型，实现了滑面接触算法的求解，并与试验结果对比验证方法的可靠性。在此基础上，反演了白格滑坡堰塞坝的形成过程，探讨了滑床底摩擦对堰塞坝成坝模式的影响规律。研究结果表明：白格滑坡堰塞坝形成过程可分为滑坡体加速运动、河谷制动和堆积成坝三阶段。白格滑坡速度最大达到45 m/s且最大速度位于滑坡体后缘，而滑坡体前缘受到河床底部摩擦影响，速度不稳定，变化幅度大且比滑坡体后缘早20 s制动。不同滑面接触模式决定了滑坡最终堆积成坝的状态：若采用滑动速度依赖型的摩擦系数模型模拟，则能够真实反映白格滑坡堰塞坝堆积状态；采用滑面的峰值摩擦系数计算，则导致计算的白格滑坡运动距离减小了32%，不能形成完全堵塞金沙江的堰塞坝；而采用残余摩擦系数计算，则导致计算的白格滑坡运动距离增大了12%。研究成果对加深堰塞坝成坝模式的认识具有意义，对于滑坡堰塞坝链生灾害预测及防灾减灾应急抢险处置具有一定的参考价值。

关键词：物质点法； 数值模拟； 接触模型； 摩擦系数； 白格滑坡堰塞坝

中图法分类号： TV122.4

文献标志码： A

DOI：10.16232/j.cnki.1001-4179.2024.04.004

0引 言

在地震、降雨和冰川融雪等外部营力作用下河谷两岸斜坡发生失稳、大变形并运动堆积于谷底，形成自然的挡水体，即滑坡堰塞坝［1-3］。滑坡堰塞坝的形成可引发一连串的次生地质灾害［4-6］，如造成上游形成堰塞湖洪涝、下游溃坝洪水泥石流、河道失稳、堰塞湖区二次滑坡并产生涌浪及其级联灾害效应，这给公路、铁路工程、管线工程、水电工程等的建设和运营带来极大威胁，严重的，会对所在流域人民生命财产造成灾难性影响［7］。

滑坡堰塞坝的形成涉及复杂的多物理过程，特别是滑坡体与滑面的接触作用，直接影响滑坡体的运动距离，决定了其后所形成的堰塞坝状态，如是否完全堵塞河谷、坝体长度和高度等。数值模拟方法为研究这一复杂过程提供了高效的研究手段。目前数值方法主要分为基于网格的方法和无网格方法。网格方法，如有限单元法，在模拟滑坡等大位移、大变形、高速运动和碰撞等问题时存在网格畸变问题，故难以再现滑坡堰塞坝形成过程这类复杂的物理过程［8］。相比之下，无网格方法解决这类问题优势更加明显，如：以离散单元法（DEM）［9］为代表的无网格粒子类方法、光滑粒子流体动力学方法（SPH）［10］、非连续变形分析方法（DDA）［11-12］等均得到了广泛的应用。其中，物质点法（MPM）作为一种兼具歐拉和拉格朗日描述的计算方法，综合了传统网格类方法和无网格方法的优势，近年来广泛应用于模拟滑坡失稳大变形运动等复杂力学问题上［13-16］。然而，考虑滑动面摩擦系数随着滑坡体速度变化的物质点计算方法目前研究尚不足。

本文针对滑坡堰塞坝形成过程，在MPM方法中实现了滑面接触算法的求解，发展了滑动速度依赖型的摩擦系数模型，通过两个典型的物理模拟试验验证了方法的准确性，在此基础上，模拟了2018年第一次白格堰塞坝形成过程，并探讨了不同滑面接触摩擦系数条件下的堰塞坝成坝堆积状态，进一步丰富了对滑坡堰塞坝形成机制的认识。

1物质点法原理

MPM方法结合了固定的欧拉网格和移动的材料点（拉格朗日粒子）的优点，将材料视为连续体，并将其离散化为材料点［17］。粒子的状态，包括质量、位置、动量和变形梯度，通过欧拉网格进行更新，在此网格上计算粒子的运动。本文采用显式MPM算法，其步骤如下：在每个时间步的开始，将材料点的信息使用常见的形状函数传递到网格节点；然后，求解控制方程以获得节点加速度，通过这些加速度计算当前时间材料点的加速度、速度和位移；最后，在欧拉网格中更新材料点的位置，以准备计算下一个时间步。

本文采用单点模式模拟干颗粒崩滑体的运动过程，不考虑水相对土体运动过程的影响，土体被视为单一介质，并通过一组材料点进行离散化，每个材料点代表固体和流体相［18］。

1.1平衡方程

由物质点求解颗粒的运动，整个系统满足质量和动量平衡方程。

1.2本构模型

对于土体材料，本文采用具有莫尔-库仑屈服准则的弹塑性本构模型。有效应力表达式为

1.3接触模型

本文MPM算法采用如下思路：首先，对每个物体基于动量方程分别求解其速度，称之为下一时间步的预测速度；然后，将接触对象视作整体，基于动量方程求解耦合体的整体速度和受力；最后，通过耦合体受力关系，根据接触关系判断，在此基础上确定是否对每个接触对象的预测速度进行校正。

若物体的接触关系是黏结接触，则不需要对速度预测值进行修正；若物体的接触关系是滑移接触，则需要对速度预测值进行校正。下一时间步的校正速度按式（9）

2模型验证

2.1颗粒柱崩塌试验验证

本节采用Lube等［20］开展的颗粒柱坍塌物理模拟试验，验证基于莫尔库伦强度准则的MPM模拟得到的颗粒柱坍塌速度和位移。颗粒柱初始堆积的长度和高度均为1 m，颗粒摩擦角为31°，黏聚力为零。试验开始前，滑坡体右侧被闸门限制移动，当打开闸门后，滑坡体开始坍塌、滑动。数值模拟的设置与物理模拟试验相同，MPM的网格尺寸为0.2 m。图1为MPM计算得到的颗粒柱坍塌演化过程。Fern和Soga［21］开展了同样的数值模拟研究，并在MPM模型的右上角设置了一个监测点记录颗粒柱崩塌的速度和位移。图2对比了本文计算得到的颗粒柱坍塌速度和位移结果与Fern和Soga［21］的数值模拟结果。可以看出，本文模拟得到的监测点速度和位移变化趋势与前人结果一致。0.5 s时刻，监测点滑动达到最大速度2.53 m/s，与前人结果相比误差在5%以内。但本文模拟得到的监测点运动时间略大于前人结果，监测点最大位移达到1.5 m，相较于前人结果高了7%。总体来说，颗粒柱崩塌试验结果定量误差在10%以内，验证了数值模拟结果的合理性。

2.2考虑底部摩擦作用的滑坡试验验证

本节采用Mangeney等［22］开展的滑坡侵蚀基底试验，验证考虑底部摩擦作用下滑坡演化模拟结果。试验槽长3 m，高1.2 m，整体呈22°倾角；滑坡体初始堆积高度为0.14 m，长度为0.2 m，滑坡底部有一层厚度为0.004 6 m的可侵蚀层铺满整个试验槽；试验槽、滑坡体和底部可侵蚀层宽度均为0.1 m。试验开始前，滑坡体右侧被闸门限制移动，当打开闸门后，滑坡体开始坍塌、滑动。

本节采用物质点法模拟再现了上述试验，分析不同时刻滑坡的滑动过程以及滑坡体在底部摩擦作用下的运动距离。滑坡体采用莫尔库伦强度准则。材料性质如下：弹性模量2.0×104 kPa，泊松比0.3，颗粒密度滑坡体的高度最初沿水槽壁方向为0.14 m；在重力作用下，滑坡体滑移，水槽壁方向的高度降低，沿水槽方向运动距离增加。在t=0.32 s时，本文数值模拟得到的滑坡体沿水槽壁方向的高度为0.10 m，低于试验结果16%，本文数值模拟得到的滑坡体沿着斜坡滑动距离为0.55 m，高于试验结果10%；在t=0.64 s时，数值模拟得到的滑坡体沿水槽壁方向高度降至约0.07 m，低于试验结果20%，数值模拟得到的滑坡体沿着斜坡滑动了0.9 m，基本与试验结果相同。总体来看，数值模拟计算得到的沿试验槽高度方向的位置低于试验结果20%以内，分析主要原因是试验槽侧壁与滑坡体之间存在一定的摩擦作用力，而数值模拟分析忽略了侧壁位置处的接触关系。但滑坡体整体的运动距离和断面与Mangeney等［22］试验结果相符，证明了该方法适用于模拟考虑底部摩擦作用的滑坡演化过程。

3白格滑坡堰塞坝案例与数值模型

2018年10月10日22：06，西藏自治区昌都市江达县波罗乡白格村和四川省甘孜藏族自治州白玉县交界处（31°04′56.41″ N，98°42′17.98″ E）的金沙江河道右岸发生大规模山体滑坡（后文简称为“10.10”白格滑坡），堵塞金沙江干流河道，形成白格堰塞坝［23-24］。

“10.10”白格滑坡发生时，金沙江江水面高程约为2 880 m，滑坡后缘高程约3 680 m，前缘高程约2 900m，高差达780 m。白格滑坡滑动面倾角在海拔3 400 m以上为25°～38°，平均约为31°，在海拔3 400 m以下较陡，坡度范围为34°～50°，平均约为39°。滑坡物质主要为岩石和砾石，岩石主要有片麻岩和蛇纹岩两种类型，两者具有不同程度的风化［25］。根据现场调查［26］，金沙江白格滑坡及堰塞坝典型纵断面如图5所示，白格滑坡体的主要物理力学参数列于表1［27］，其中参与计算的基岩参数按照工程类比法选定［15］。

白格滑坡堰塞坝形成过程的物质点法数值计算模型如图6所示。采用三角形网格，基岩网格尺寸为60 m，滑坡体网格为20 m，共生成了8 186个网格和4 202个节点。白格滑坡体与基岩之间的摩擦作用是本文研究的一个重要方面。Hu等［28］试验研究发现，当滑坡体剪切速度超过1 m/s时，滑面的摩擦系数非常低，甚至能达到0.05。因此在本研究中，参考前人对大光包滑坡的研究成果［15］，如图6所示定义了两个滑床接触表面。接触表面A位于滑坡体下方，而接触表面B是原始的山体表面，由于滑坡体滑移，该表面摩擦由峰值摩擦系数迅速降低为残余的摩擦系数。本研究中，为使滑坡体触发移动，参考前人方法［29］，取滑坡体内摩擦角的70%作为接触表面A的残余摩擦系数（即0.4）；而接触表面B的摩擦系数根据基岩摩擦系数设定为初始值0.6。Li等［15］研究表明，滑坡体启动速度增加至2 m/s后，基岩接触面摩擦系数降低，本文参考前人研究取25%的初始值作为接触表面B的残余摩擦系数（即0.15）。如图6所示，分别在白格滑坡体的前缘、中部和后缘设置了一个速度监测点，以分析在滑坡演化过程中的速度变化规律。数值计算时间采用显示积分，时间步长为0.001 s，共计算70 s。

4计算结果分析

白格滑坡堰塞坝形成过程的速度场演化如图7所示，图中虚线表示白格堰塞坝实际断面。白格堰塞坝形成过程大约经历70 s，可分为3个基本阶段：t=0～20 s是滑坡体加速运动阶段；t=20～30 s，滑坡体运动至河谷中，由于河谷具有较高的粗糙度以及对岸岸坡的空间限制，滑坡体进入河谷制动阶段；t=30～70 s，此阶段，滑坡体前缘运动至最远距离，后缘滑坡体逐渐堆积于前缘滑坡体上，速度逐渐减为零，滑坡体在此阶段堆积成坝。t=70 s时，数值模拟得到的白格堰塞坝最终断面与实际白格堰塞坝断面对比，可以发现数值模拟得到的滑坡体运动至河谷最远端的距离，即堰塞坝长度，与实际情况较为一致，但数值模拟得到的白格堰塞坝高度高于实际值。这主要是由于本文采用的是二维的物质点法模拟，而滑坡体运动实际是三维的，在堆积成坝阶段，滑坡体由于碰撞对岸斜坡，会向着河谷上下游展宽，目前二维的数值模拟无法考虑这一情况。此外，白格滑坡堰塞坝的形成实际上存在滑坡体冲向对岸刮铲对岸崩坡积体后崩回过程，但这一过程涉及新的物质加入运动，需要在数值分析的质量守恒方程和动量方程中加以考虑，而本文模型尚未考虑这一过程，这也是计算结果存在差异的原因之一。

白格滑坡体前缘、中间及后缘位置处速度特征如图8所示，最大速度约45 m/s，出现在滑坡体后缘，这一结果基本与Li等［30］基于离散元模拟得到的白格滑坡最大滑速50 m/s吻合。滑坡体前缘受到底部摩擦影响，速度变化幅度较大；而滑坡体后缘及中间速度基本呈现三个阶段，即快速增大后维持一段较快的运动速度，最后逐渐减小。受滑坡体前缘牵引与后缘挤压作用，滑坡体中间位置在滑坡体加速运动阶段的速度增长幅度最大。滑坡体前缘与中部速度在t=50 s前降为零，而滑坡体后缘速度在此后20 s逐渐降低。Hu等［31］基于计算流体动力学研究也揭示出白格滑坡体运动41 s后抵达对岸山体，与本文模擬得到的滑坡体前缘50 s后停止运动的认识基本一致。

进一步对比分析不同滑面接触模式下白格滑坡最终堆积成坝状态。分别采用了峰值摩擦系数0.4和残余摩擦系数0.15模拟白格滑坡成坝过程，计算结果如图9所示。图9（a）是考虑速度状态相关的滑面接触模式下白格滑坡最终堆积成坝状态，图9（b）是滑面接触为峰值摩擦系数的白格滑坡最终堆积成坝状态，图9（c）是滑面接触为残余摩擦系数的白格滑坡最终堆积成坝状态。

由图可知，若采用考虑速度状态相关的滑面接触模式计算，白格滑坡运动至水平距离2 129.4 m后停止，最大运动距离为841.5 m，这与实际情况基本一致；若采用滑面接触为峰值摩擦系数计算，白格滑坡运动至水平距离1 858.2 m后停止，最大运动距离为570.3 m，与图9（a）计算结果相比，减小了32%；若采用滑面接触为残余摩擦系数计算，白格滑坡运动至水平距离2 231.5 m后停止，最大运动距离为943.6 m，与图9（a）计算结果对比，增大了12%。故精确考虑速度状态相关的滑面接触模式，能更加精确地计算滑坡体的运动状态以及模拟其堵江成坝的模式。

5结 论

本文在MPM方法中实现了滑面接触算法的求解，提出了与滑动速度相关的摩擦系数模型，在此基础上反演了白格滑坡堰塞坝形成过程，得到如下结论：

（1） 通过颗粒柱崩塌试验和考虑底部摩擦作用的滑坡试验，验证了所提出的数值模型的正确性，MPM方法对滑坡大变形产生的位移和速度的计算均与试验结果一致。

（2） 白格滑坡堰塞坝形成过程可分为三个阶段，即滑坡体加速运动阶段、河谷制动阶段和堆积成坝阶段。滑坡体速度最大达到45 m/s且位于滑坡体后缘，而滑坡体前缘受到底部摩擦影响，速度不稳定，变化幅度大比滑坡体后缘早20 s制动。

（3） 不同滑面接触模式决定了滑坡最终堆积成坝状态。若采用考虑速度状态相关的滑面接触模式模拟，则能够真实反映白格滑坡堰塞坝堆积状态；采用滑面接触为峰值摩擦系数计算，则计算的白格滑坡运动距离减小了32%，不能形成完全堵塞金沙江的堰塞坝；采用滑面接触为残余摩擦系数计算，则导致计算的白格滑坡运动距离增大了12%。

本文研究还存在一些局限性，未来值得进一步深入研究。如本研究采用的是二维物质点法，而沟梁相间真实地形条件下考虑滑坡体基底侵蚀与碰撞刮产增容机制的三维物质点法值得进一步研究。

参考文献：

［1］石振明，张公鼎，彭铭，等.非均质结构堰塞坝溃决机理模型试验研究［J］.工程科学与技术，2023，55（1）：129-140.

［2］蔡耀军，杨兴国，张利民，等.堰塞湖风险评估快速检测与应急抢险技术和装备研发研究构想与成果展望［J］.工程科学与技术，2020，52（2）：10-18.

［3］吴昊，年廷凯，单治钢.滑坡堵江成坝的形成演进机制及危险性预测方法研究进展［J］.岩石力学与工程学报，2023，42（增1）：3192-3205.

［4］ZHONG Q，WANG L，CHEN S，et al.Breaches of embankment and landslide dams-State of the art review［J］.Earth-Science Reviews，2021，216：103597.

［5］崔鵬，郭剑.沟谷灾害链演化模式与风险防控对策［J］.工程科学与技术，2021，53（3）：5-18.

［6］钟启明，钱亚俊，单熠博.崩滑堰塞湖的形成—孕灾—致灾机理与模拟方法［J］.人民长江，2021，52（2）：90-99.

［7］吴昊.滑坡堵江成坝过程模拟及危险性预测方法研究［D］.大连：大连理工大学，2021.

［8］杜文杰，盛谦，杨兴洪，等.基于两相双质点 MPM 的滑坡堵江灾害链生全过程分析［J］.工程科学与技术，2022，54（3）：36-45.

［9］LI D，NIAN T，WU H，et al.A predictive model for the geometry of landslide dams in V-shaped valleys［J］.Bulletin of Engineering Geology and the Environment，2020，79：4595-4608.

［10］WU H，NIAN T，SHAN Z，et al.Rapid prediction models for 3D geometry of landslide dam considering the damming process［J］.Journal of Mountain Science，2023，20：928-942.

［11］ZHANG Y，WANG J，XU Q，et al.DDA validation of the mobility of earthquake-induced landslides［J］.Engineering Geology，2015，194：38-51.

［12］NIAN T，ZHANG Y，WU H，et al.Runout simulation of seismic landslides using discontinuous deformation analysis （DDA） with state-dependent shear strength model［J］.Canadian Geotechnical Journal，2020，57：1183-1196.

［13］ZHAO K L，QIU L C，LIU Y.Two-layer two-phase material point method simulation of granular landslides and generated tsunami waves［J］.Physics of Fluids，2022，34：123312.

［14］ZHAO K L，QIU L C，LIU Y.Numerical study of water wave generation by granular-liquid mixture collapse using two-phase material point method［J］.Applied Ocean Research，2023，137：103608.

［15］LI X，TANG X，ZHAO S，et al.MPM evaluation of the dynamic runout process of the giant Daguangbao landslide［J］.Landslides，2020，18：1509-1518.

［16］ZHU Y，ISHIKAWA T，ZHANG Y，et al.A FEM-MPM hybrid coupled framework based on local shear strength method for simulating rainfall/runoff-induced landslide runout［J］.Landslides，2022，19：2021-2032.

［17］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：248-273.

［18］TRONCONE A，PUGLIESE L，CONTE E.Analysis of an excavation-induced landslide in stiff clay using the material point method［J］.Engineering Geology，2022，296：106479.

［19］ABE K，SOGAOGA K，BANDARA S.Material point method for coupled hydromechanical problems［J］.Journal of Geotechnical Geoenvironmental Engineering，2014，140：04013033.

［20］LUBE G，HUPPERT H E，SPARKS R S J，et al.Collapses of two-dimensional granular columns［J］.Physical Review E，2005，72：041301.

［21］FERN E J，SOGA K.The role of constitutive models in MPM simulations of granular column collapses［J］.Acta Geotechnica，2016，11：659-678.

［22］MANGENEY A，ROCHE O，HUNGR O，et al.Erosion and mobility in granular collapse over sloping beds［J］.Journal of Geophysical Research Earth Surface，2010，115（F3）.

［23］FAN X，YANG F，SIVA S S，et al.Prediction of a multi-hazard chain by an integrated numerical simulation approach：the Baige landslide，Jinsha River，China［J］.Landslides，2019，17：147-164.

［24］ZHONG Q，CHEN S，WANG L，et al.Back analysis of breaching process ofBaige landslide dam［J］.Landslides，2020，17：1681-1692.

［25］蔡耀軍，栾约生，杨启贵，等.金沙江白格堰塞体结构形态与溃决特征研究［J］.人民长江，2019，50（3）：15-22.

［26］陈祖煜，张强，侯精明，等.金沙江 “10·10” 白格堰塞湖溃坝洪水反演分析［J］.人民长江，2019，50（5）：1-4.

［27］陈祖煜，陈生水，王琳，等.金沙江上游 “11.03” 白格堰塞湖溃决洪水反演分析［J］.中国科学：技术科学，2020，50：763-774.

［28］HU W，HUANG R，MCSAVENEY M，et al.Superheated steam，hot CO2 and dynamic recrystallization from frictional heat jointly lubricated a giant landslide：Field and experimental evidence［J］.Earth and Planetary Science Letters，2019，510：85-93.

［29］WANG J，WANG S，SU A，et al.Simulating landslide-induced tsunamis in the Yangtze River at the Three Gorges in China［J］.ActaGeotechnica，2021，16：2487-2503.

［30］LI D，NIAN T，TIONG RLK，et al.River blockage and impulse wave evolution of the Baige landslide in October 2018：Insights from coupled DEM-CFD analyses［J］.Engineering Geology，2023，321：107169.

［31］HU Y，YU Z，ZHOU J.Numerical simulation of landslide-generated waves during the 11 October 2018Baige landslide at the Jinsha River［J］.Landslides，2020，17：2317-2328.

（编辑：郑 毅）