基于两相双质点MPM的滑坡堵江灾害链生全过程分析

杜文杰，盛 谦，杨兴洪，魏鹏飞，李丽华，付晓东*

(1.中国科学院 武汉岩土力学研究所 岩土力学与工程国家重点实验室，湖北 武汉 430071；2.中国科学院大学，北京 100049；3.云南交投集团公路建设有限公司，云南 昆明 650100；4.湖北省地质局第八地质大队，湖北 襄阳 441000)

自然灾害事件的发生和发展并非彼此独立，相继发生的灾害在时间和空间上都存在着一定关联。由一种灾害引发或导致另外一种灾害发生，从而使得灾害的发生具有链式效应，这样的灾变现象称为灾害链。

滑坡—堵江—堰塞湖灾害链是高山峡谷地区典型的灾害链形式，例如：2018年10月11日和11月3日，西藏白格村金沙江右岸发生了两次大规模高位滑坡，先后堵塞金沙江形成堰塞湖，后期堰塞体溃决形成的下泄洪水给下游四川、云南造成了严重的洪涝灾害。相比于单一滑坡灾害而言，滑坡—堵江—堰塞湖灾害链的危险性通过时空上的延拓而大幅增加。然而，由于灾害之间链生关系的复杂性及当前研究方法的局限性，目前仍缺乏针对链生灾害的深入研究。

灾害成链的关键在于承灾体之间的能量传递，对于滑坡—堵江—堰塞湖这一灾害链形式，滑坡与水体间复杂的流-固耦合作用作为灾害链的衔接过程，往往决定了后续灾害的影响规模。其中，水以不同形式参与到灾害链的全过程，以自由水形式存在的江河、水库与海洋，以及入渗到滑坡体中的孔隙水，为模型的建立及灾害的精细化分析带来困难。

传统网格方法，如有限单元法，在模拟滑坡等大位移和大变形问题时存在网格畸变问题，水体的翻腾和破碎现象具有高度的非线性，基于网格的经典流体动力学方法很难进行建模与模拟。相比之下，无网格方法的优势明显，诸如：SPH为代表的无网格粒子类方法、DDA方法，DEM-SPH为代表的离散元-无网格耦合方法均得到了广泛的应用。其中，物质点法（MPM）作为一种兼具欧拉和拉格朗日描述的计算方法，综合了传统网格类方法和无网格方法的优势，二十多年来广泛应用于模拟滑坡失稳、流体流动及流-固耦合等复杂力学问题上。然而，传统的MPM只能独立地模拟滑坡或水体，不适用于模拟滑坡堵江过程中的复杂流-固耦合问题。近年来发展的两相双质点格式的MPM通过对固相与液相分开建模，可将自由水、孔隙水在统一的框架下进行显式表达。

本文针对滑坡堵江灾害链生过程，在Pradhana发展的两相双质点物质点法（TPDP-MPM）的基础上，将TPDP-MPM拓展应用至地质灾害（链）的分析与评估中。通过模拟再现模型试验尺度下的Lituya滑坡体-水体相互作用过程，验证了TPDP-MPM在处理复杂流-固耦合问题方面的能力。在此基础上，对2011年湖北十堰二荒村滑坡堵塞平渡河灾害链生全过程进行计算模拟与分析，进一步丰富了滑坡堵江链生灾害方面的认识。

1 两相双质点MPM

1.1 物质点法的基本原理

物质点法综合了网格类方法与无网格类方法的优势，将计算域视作一组物质点组成的连续体。每个质点（MPs）存储其代表区域的所有物理信息，如速度、质量、位移和其他量，通过形函数实现质点与背景网格节点（GNs）间信息的相互映射。如图1所示，允许代表实际计算域的 MPs（用红色圆点标记）自由移动，并在每一步累积映射自GNs的位移，而背景网格则在每一步结束后重置。

图1 物质点法计算流程Fig. 1 Computational cycles of MPM

在计算域 Ω中，质量守恒和线动量守恒是连续体的基本方程，可以表示为：

式（1）～（2）中， ρ 和 ρ˙分别为密度和密度的时间导数，

v

和

v

˙分 别为速度和加速度矢量， σ为柯西应力张量，∇· 表 示散度算子，

b

表示体力。由于质量信息储存在MPs上，因此，物质点法质量守恒自动满足。动量守恒的弱形式可以通过加权残差法推导出来。将任意试函数

w

引入式（2），然后，在以 ∂Ω 为边界的整个计算域 Ω上进行积分，控制方程的弱形式可以改写为：

式中， τ表 示 d

S

面上的剪应力张量。

1.2 两相双质点MPM

通常，在流-固耦合力学问题中，水不仅作为孔隙中的耦合流体存在，而且作为自由流体存在，如河流、海洋和水库。两相双质点MPM通过对固相与液相分别建模，将自由水、孔隙水在统一的框架下进行显式地表达，如图2所示。

图2 两相双质点MPM示意图Fig. 2 Schematic diagram of two-phase double-point MPM approach

对TPDP-MPM的固相速度-液相速度格式进行简要介绍，详细内容参考文献[25-28]。

饱和多孔介质的动力学行为同样用动量守恒方程来表述：

式中： α=s,f 分别表示固相和液相； D /D

t

为物质导数算子； σ为α相的局部应力张量；

g

为重力加速度；

m

表示两相之间动量的交换，必须满足以确保动量守恒，在满足这一条件的前提下，有：

对式（5）求和，可得到混合物的动量守恒方程：

式中，表示混合物中的柯西应力，即各相应力之和。根据柯西应力的概念，各相的线性动量守恒意味着混合物的线性动量守恒。

Jassim等所提出的重力驱动流中液相与固相的动量交换项

m

、

m

为：

式中：为阻力系数，其中，

n

为 孔隙比，为滑坡体的固有渗透率，

g

为重力加速度标量；

p

为液相压强， ∇

φ

为水的体积分数。

c

(

v

−

v

)表示由固相质点与液相质点耦合产生的黏性力，由于其足以描述耦合系统的运动，忽略

p

∇

φ

，两相之间的动量交换被认为是纯粹的耗散过程。

针对滑坡堵江灾害链特点，为克服使用水体真实状态方程而导致的时间积分步长过小的问题，将水体视作微可压缩流体，采用Becker等提出的人工状态方程作为求解水体压强的控制方程：

式（8）～（11）中，

B

为水体压强相关项， γ为反映流体微可压缩特性的项，变形梯度

J

为液相质点当前体积与初始体积之比，

I

=[1,1,1,0,0,0]， ψ为势能函数。根据式（10），压强

p

由式（9）所给出的液相势能函数对变形梯度

J

进行求导运算得到。

此外，基于内接Mohr-Coulomb屈服面的Drucker-Prager准则建立岩土材料的弹塑性本构关系模拟滑坡体。

2 滑坡体-水体相互作用的方法验证

利用两相双质点MPM法对Fritz开展的模型试验进行建模与模拟，模拟结果进一步证明了TPDP-MPM在涉及耦合问题的灾害链评估中良好的应用前景。

1958年7月8日，Gilbert海湾一侧山体在里氏8.3级地震作用下失稳，约3×10m滑坡体滑入Gilbert海湾，激起的涌浪在对岸山体的爬升高度高达524 m，图3为滑坡区域航拍照片。Fritz等建立了模型试验尺度的Lituya滑坡概化模型，将实际滑坡按1∶675的比例进行了缩小。本文以模型试验尺度建立了数值计算模型，通过两套质点，将滑坡体离散为2 189个质点，将水体离散为5 877个质点，结构化背景网格尺寸为0.02 m，两相材料参数及其他输入参数见表1。

图3 Gilbert海湾Lituya滑坡航拍示意图Fig. 3 Aerial photo of Lituya landslide in Gilbert inlet

表1 计算参数
Tab. 1 Input parameters

滑坡体（砂土） 水体 滑床 时间步长/(10-5 s)密度/(kg·m-3)摩擦角/(°)黏聚力/kPa 泊松比 弹性模量/(107 Pa)(10-3 m·s-1) 初始孔隙比 体积模量/(107 Pa)渗透系数/密度/(kg·m-3) γ 摩擦系数2 400 39 0 0.26 1.6 5 0.2 2 1 000 7 0.2 2.5

在计算过程中实时记录滑坡诱发涌浪的运动，如图4所示。通过与Fritz模型试验结果对比发现，TPDP-MPM可以准确模拟滑坡体与水体相互作用及耦合过程，从而验证了TPDP-MPM在模拟复杂流-固耦合问题方面的能力。进一步地，可以将TPDP-MPM应用于高山峡谷地区广泛分布的滑坡堵江灾害链的力学分析与风险评估。

图4 滑坡诱发涌浪过程数值模拟与试验结果对比Fig. 4 Comparing the numerical and experimental results of landslide induced surge process

3 典型滑坡堵江灾害案例与仿真模型

二荒村滑坡位于湖北省西北部房县南部山区，海拔高程大多在1 100～1 500 m之间，地形切割深度在400～1 000 m之间，属构造侵蚀剥蚀中低山地形。滑体所在平渡河河段自南东向北西径流，构成斜向河谷。河谷横断面呈“V”字形，河谷宽20～40 m，两岸地势陡峻，坡角35°～65°，局部呈陡崖。滑坡纵长约300 m，后缘宽50 m，前缘宽 150 m，滑体厚10～60 m，平均厚40 m，总滑动量15×10m。

2011年6月14日16时40分，湖北省十堰市房县上龛乡二荒村二组平渡河右岸山体发生滑坡，滑坡堆积体堵塞平渡河，形成堰塞湖。图5（a）为二荒村滑坡最终堆积并堵塞平渡河的现场照片。滑坡造成百户沟电站（装机容量7 200 kW）前池和部分压力管道毁坏，致使百户沟电站停止发电。二荒村滑坡地质剖面图如图5（b）所示，滑坡堆积体由巨石、碎块石、碎石土及粉质黏土等组成，杂乱无序堆积。巨石块度达3.0 m×2.0 m×1.5 m，块石直径一般为0.4～1.0 m，黏土所占比例较少，土石比约为1∶9。堆积体结构松散，空隙大，可见堰塞湖内水流从空隙中下泄。

图5 二荒村滑坡概况Fig. 5 Overall situation of Erhuang Village landslide

通过两套质点分别对滑坡体与水体进行建模，如图6所示。地层自下而上依次为灰岩、泥质粉砂岩和粉质黏土夹杂碎块石。基于泊松采样技术，将滑坡体离散为2 491个质点，平渡河水体离散为1 172个质点，结构化背景网格尺寸为1.0 m。采用显示更新格式，参考材料经验与文献取值，确定两相材料参数，见表2。

表2 计算参数
Tab. 2 Input parameters

滑坡体（粉质黏土夹杂碎块石） 水体 滑床 时间步长/(10-5 s)密度/(kg·m-3)摩擦角/(°)黏聚力/kPa 泊松比 弹性模量/(107 Pa)(10-3 m·s-1) 初始孔隙比 体积模量/(107 Pa)渗透系数/密度/(kg·m-3) γ 摩擦系数2 300 19 19 0.26 6 75 0.2 2 1 000 7 0.6 2.5

图6 二荒村滑坡仿真模型Fig. 6 Numerical model of Erhuang Village landslide

通过解除初始地应力约束实现二荒村滑坡的启动：根据地质勘察资料，二荒村滑坡的滑动面是明确的，以滑坡体的初始形状作为计算边界对滑坡体进行约束，输入偏大的强度参数计算滑坡体初始地应力；随后解除初始约束，给滑坡体赋以实际强度参数，使滑坡体进入不稳定状态后开始下滑。

4 滑坡堵江灾害链生全过程分析

图7为TPDP-MPM模拟得到的二荒村滑坡失稳堵塞平渡河的全过程持续约32.3 s间从滑坡启动

t

=0到最终堆积稳定

t

=32.3 s中的6个时刻的滑坡形态。由图7可以看出：初始时刻，在重力作用下，二荒村滑坡体启动下滑，滑坡体中后缘拉裂缝开始发育并向深部发展；12.5 s左右，滑坡体前缘到达平渡河水面后开始侵入河体，15.75 s时的模拟结果如图7（c）所示；15.8 s左右，滑坡体前缘撞击河谷后开始沿岸坡堆积；32.3 s左右，滑坡运动停止，如图7（f）所示，堆积体在靠近对岸岸坡的最低高程为697.2 m，平渡河水面高程为700 m，即二荒村滑坡几乎完全堵塞平渡河，并形成堰塞体。图8为TPDP-MPM计算得到的二荒村滑坡最终堆积形态与现场实际情况（图8中虚线）的对比，两者在轮廓上吻合良好。

图7 二荒村滑坡失稳后运移—堵江全过程Fig. 7 Whole process of migration and river blocking after failure of landslide

图8 滑坡最终堆积形态结果对比Fig. 8 Final accumulation profile of landslide

4.1 速度分布特征

不同时刻滑坡体速度场如图9所示。

根据滑坡失稳过程呈现的模式，将滑坡体沿滑动方向等分为前缘部分、中间部分与后缘部分（标记如图7（b）所示）。失稳滑动前期（

t

=0～9.8 s），在重力作用下，滑坡体中部及前缘开始滑动，并逐渐产生比后缘更高的速度；

t

=9.8 s时，前缘部分速度就达到了21.0 m/s。在这一阶段，滑坡后缘的运动始终滞后于滑坡前缘。随着时间的推移，滑坡后缘与中前缘的速度差距进一步拉大，

t

=9.8 s后，后缘与中前缘连接处拉裂缝开始发展，如图9（b）所示。随着拉裂缝逐步向深部扩展，后缘“坡脚处”产生自由表面，由于约束减少，后缘开始加速下滑，此时滑坡后缘运移速度高于上一阶段的速度，滑坡趋于整体的高速运动。

t

=12.3 s时，滑坡体前缘开始由滑床抛出，如图9（c）所示，滑坡体质点速度进一步提高，在16.0 s左右，滑坡体入水速度达到了39.6 m/s。此后，随着滑坡体前缘开始侵入水体，对比平渡河水面上下的滑坡体质点速度可以发现：入水部分的滑坡体在水体制动作用下速度开始大幅降低，如图9（e）所示。此后，滑坡整体向河谷内减速、堆积，直至

t

=32.3 s末，滑坡整体速度趋于0，此时滑坡体与平渡河水体均接近于稳定状态，堆积过程结束，如图9（h）所示。

4.2 能量演化规律

图10为滑坡—堵江灾害全过程中的滑坡体质点动能、滑坡体总动能随时间的变化情况。由图10可知：滑坡由启动到前缘抵达河面前（0～11.9 s），滑坡动能快速提高，图10（a）中的两条曲线可以说明滑坡体前缘和后缘的动能演化规律，前期滑坡后缘动能始终低于前缘，后缘动能在6.2 s后开始与滑坡前缘以一致的增速增长；9.68 s左右，滑坡前缘到达临河缓坡区域后，前缘速度有所降低，如图9（b）所示，滑坡体总动能增长放缓，如图10（b）的范围Ⅲ所示；11.9 s滑坡前缘入江后，如图10（a）范围Ⅰ所示，质点速度的离散性增强，在水体制动作用下，滑坡体前缘入水部分的质点速度明显降低。随着入水滑坡体积的增大，水体制动作用越发显著，滑坡体总动能保持增长但趋势进一步放缓，如图10（b）中范围Ⅰ所示；15.8 s左右滑坡体前缘到达河谷前，滑坡体总动能达到峰值，滑坡体前缘撞击河谷后速度趋于0，滑坡体开始在河谷底部向上堆积，由于重力势能的减小，滑坡体总动能开始锐减，如图10（b）范围Ⅱ所示，直至32.3 s时完全堆积、静止。

4.3 滑坡—堵江全过程分析

从TPDP-MPM模拟结果来看，二荒村滑坡在重力作用下启动，经历下滑过程中的加速、破裂解体后，滑坡体前缘侵入平渡河水体，并在撞击平渡河河谷后开始在谷底堆积，最终完全堵塞平渡河，形成方量巨大的堰塞体。滑坡—堵江全过程可分为4个阶段。

1）失稳启动

t

=0～9.8 s时，在重力作用下，滑坡体的动能开始累积，滑坡体中前缘很快加速到10 m/s以上。在这一阶段，受制于中前缘的约束，滑坡后缘滑动速度滞后于中前缘，如图9（b）所示。

2）高速滑动

t

=9.8～12.5 s时，随着中前缘与后缘速度差增加，在滑坡体后缘与中前缘连接处形成了拉裂缝，并向深部发展，滑坡体中前缘开始脱离后缘，约束减弱，后缘开始加速下滑，滑坡体趋于整体滑动。由于临河岸坡地形逐渐平缓，滑坡体前缘由滑床向河谷中央抛射，进一步提高了滑坡体前缘的入水速度。在这一阶段，滑坡表现为整体的高速滑动，如图9（c）所示。

3）入江制动

t

=12.5～15.8 s时，滑坡前缘开始侵入水体，对比图9（d）平渡河水面上下的质点速度，以及图10（a）入水前后质点动能变化，可以发现：在水体撞击制动及流-固耦合耗散作用下，入水部分滑坡体动能显著降低；这一阶段，在中后缘持续的推动作用下，滑坡体依然整体向河谷高速滑移；但随着入水滑坡体积的增大，水体的制动作用越发明显，滑坡体总动能增长速度放缓。

图10 滑坡体能量演化Fig. 10 Energy evolution of landslide mass

4）堆积成坝

t

=15.8～32.3 s时，随着滑坡体前缘开始碰撞河谷，质点速度急剧衰减并趋近于0，如图9（f）～（h）所示。滑坡体开始在河谷谷底沿岸坡向上堆积，随着重力势能减小，滑坡整体运移速度开始衰减，并在32.3 s左右滑坡整体堆积，形成堰塞坝堵江。

图9 滑坡体内部速度场演化Fig. 9 Velocity field of the landslide mass

5 结论与展望

为再现典型滑坡堵江链生灾害演化全过程，将两相双质点MPM扩展应用于地质灾害链的分析与评估中。以2011发生的湖北十堰二荒村滑坡为研究对象，利用两套质点分别对滑坡体与水体分别进行建模，模拟了二荒村滑坡堵江全过程，对滑坡启动、运移、入水制动及堆积堵江过程进行了运动学分析，得出以下结论：

1）通过两相双质点MPM模拟再现了二荒村滑坡堵塞平渡河的全过程，整个过程持续了32.3 s，计算得到的二荒村滑坡最终堆积形态与现场勘查形态大致相同，这也证明了TPDP-MPM在分析地质灾害链问题上的可靠性。

2）二荒村滑坡在重力作用下启动，经历下滑过程中的加速、破裂解体后，滑坡体前缘侵入平渡河水体，并在撞击平渡河河谷后在谷底堆积，最终完全堵塞平渡河形成堰塞坝。从TPDP-MPM模拟结果来看，二荒村滑坡失稳—堵江全过程可分为失稳启动、高速滑动、入江制动和堆积成坝4个阶段。

3）基于速度场分布揭示了滑坡体不同阶段的运移模式：在滑坡初始启动阶段，滑坡后缘的运动始终滞后于滑坡前缘；随着中前缘与后缘速度差的拉大，连接处拉裂缝开始发展；此后滑坡体后缘开始加速，滑坡趋于整体的高速运动。

4）质点动能及总动能演化规律揭示了水体制动效应对滑坡体堆积过程的影响：滑坡体入水后质点速度的离散性增强，滑坡体总动能保持增长但趋势放缓。随着入水滑坡体积的增大，水体的制动作用越发明显；随后滑坡体前缘撞击河谷后开始沿岸坡堆积，由于重力势能的减小，滑坡体总动能开始锐减直至最终堆积成坝。

受限于计算能力不足，本文仅进行了2维滑坡堵江灾害链生过程模拟，后续将引入针对性的岩土失效准则，针对滑坡启动机制、3维全过程模拟、高效并行计算等方面开展进一步研究。