基于PFC3D的滑坡堰塞坝堆积过程与形态模拟

吴 建 川,张 世 殊，吴 爽，冉 从 彦

(1.中国电建集团 成都勘测设计研究院有限公司，四川 成都 610072; 2.中国地质大学(武汉)工程学院，湖北 武汉 430074)

滑坡堰塞坝是指由于降雨、地震、火山喷发等原因引起山体滑坡而截堵山谷、河道形成的坝体[1]。依据成因的不同，堰塞坝可以划分为滑坡型堰塞坝、崩塌型堰塞坝以及滑坡、崩塌运动过程中形成的泥石流型堰塞坝，其中滑坡型堰塞坝约占全部堰塞坝的70%以上[2]。在西部高山峡谷区域的河道，滑坡失稳堵江极易形成较高的堰塞坝，堰塞坝的物质组成通常由松散堆积物组成，在河水漫顶、雨水冲刷的条件下极易发生溃坝并引发洪水，对水电工程的选址、建设、安全运行和下游居民的生命财产安全造成巨大威胁。如2018年10月10日、11月3日，西藏自治区境内的金沙江白格村在同一地点先后两次发生大规模的山体滑坡，堵塞金沙江形成堰塞湖，严重威胁着堰塞坝上下游人民群众的生命财产和下游已建及在建水电站的安全[3]。

滑坡型堰塞坝的形成过程和几何形态特征与河道纵剖面形状、滑坡体积以及滑坡体的运动形态有关[4]。为了探寻滑坡型堰塞坝的形成过程和几何形态特征，成本低、效果明确的数值模拟是一种有效的手段。连续力学数值模拟方法方面：原俊红[5]、罗刚[6]运用连续力学数值模拟方法，模拟边坡的破坏特征并预测堰塞坝高度；徐文杰[7]对肖家河滑坡稳定性做了分析，并用ABAQUS/Explicit动力有限元分析滑坡在地震工况下的堵江过程。非连续力学数值模拟方法方面：Chang和 Taboada[8]利用二维颗粒流离散元程序模拟了台湾九份二山滑坡力学过程及堆积行为；孟云伟和柴贺军[9]认为颗粒流离散元能充分考虑实际岩体的大位移、大变形，能够更有效、更真实地模拟滑坡运动的优点；张龙[10]认为PFC3D软件能较好地模拟高速远程滑坡的运动过程及堆积情况，对于灾害堆积形态与影响范围可以直接模拟得到；Poisel和Preh[11]利用PFC3D对滑坡运动导致的水库涌浪进行了模拟，并取得非常好的效果；梁承洋[12]用PFC2D模拟了川藏交通廊道冰磧物滑坡堵江堆积物的形态，同时研究了单因素变量对堵江的影响；李子隆[13-14]用PFC2D开展了覆盖层库岸边坡在静水位条件下的失稳堰塞模拟和滑坡失稳形成堰塞坝的高度预测模拟；Dal Sasso[15]采用地貌指数法对滑坡堰塞坝几何特征进行了研究；Zhao[16]等采用CFD-DEM耦合的方法模拟滑坡堵江形成堰塞坝的过程。针对滑坡堰塞坝几何特征的研究，Dong[17]等指出准确的数字地形数据在绝大多数情况下无法短时间内获得。

综上，连续力学数值模拟方法虽有模拟速度快等突出优势，但由于模拟的块体不能破裂，模拟的效果并不明显；而于非连续力学方法的颗粒流离散元方法能较好地模拟大变形、大位移，对模拟滑坡失稳堰塞过程及堰塞坝堆积形态具有较好的适用性。

针对堰塞坝的溃决，蔡耀军[18]研究了白格堰塞湖形成过程及成因，并分析了堰塞体结构及形态特征、溃决发展阶段、溃决特征值；陈祖煜[19]运用DB-IWHR[20]溃坝洪水分析程序反演分析了金沙江“10.10”白格堰塞湖漫顶自然泄流过程；王敏[21]采用MIKEl l模型和一维溃坝洪水模型(DBFM)，对金沙江“11·3”白格堰塞湖溃决洪水演进过程进行了对比计算。

本文以雅砻江上游水电规划、开发的重点关注对象唐古栋滑坡为研究对象，运用三维颗粒流离散元PFC3D建立三维滑坡模型，对唐古栋滑坡1967年失稳形成的堰塞坝堆积形态进行反演模拟，校核标定的细观参数。在此基础上，开展唐古栋滑坡变形强烈部位失稳形成堰塞坝的预测模拟，并采用DAMBRK溃坝洪水计算数学模型对堰塞坝溃坝造成的影响进行了探讨。

1 唐古栋滑坡

1.1 简 介

唐古栋滑坡位于雅砻江中游雅江县孜河乡雨日村，位于雅砻江两河口至卡拉中游规划河段力邱河口-蒙古山段的右岸，为一巨型高速岩质滑坡，主滑方向145°，与雅砻江流向近垂直[22]。滑坡体横河长约1.4 km，顺河宽约1.2 km，展布面积超过1.4 km2，前缘高程2 450 ～2 480 m，后缘高程3 500 m，高差约1 000 m。滑坡区地层岩性为三叠系上统侏倭组(T3zh)变质砂岩与印支-燕山期侵入的花岗伟晶岩脉(γρ)，岩层产状为N50°～75°W/NE∠40°～60°，倾向坡内[23]。

1.2 滑坡分区

唐古栋滑坡主要包括上游侧变形体和下游侧已滑坡区，根据岸坡变形破坏程度及残留物规模，由上游至下游可分为A、B、C、D四个区(见图1)，A区为未发生滑动的强变形区，B、C、D区为已发生滑动区，各区特征及规模见表1[24]。

1.3 滑坡体物理力学参数

唐古栋滑坡物质组成由外至内可分为滑坡堆积物、崩坡堆积体、强风化岩体、弱风化岩体、微-新岩体。基于滑坡岩土体物理力学试验成果，并通过抗剪强度参数反演，综合确定各层岩土体物理力学参数取值，见表2[25]。

表1 唐古栋滑坡分区特征及残余方量汇总

表2 唐古栋滑坡各层岩土体抗剪强度参数取值

1.4 滑坡失稳-堰塞-溃决

1967年6月8日，唐古栋滑坡发生大规模整体滑移失稳破坏，且在5 min之内堵断雅砻江并冲向对岸岸坡；滑坡体在高速滑动过程中充分解体，堆积形成顺河向呈梯形的堰塞坝，堰塞坝方量约6.41×107m3；堰塞坝底面沿江长1 500 m,顶面沿江长860 m，高225～270 m；对岸坝高335 m，本岸坝高175 m[26]。

堰塞坝阻断雅砻江达10 d，黄润秋[27]等调查发现，堰塞坝回水53 km，壅高水位高程达2 575 m，库容达6.8亿 m3；堰塞坝溃决后，洪峰流量达5 3000 m3/s，坝下游10 km处水位上涨达48 m，洪水沿雅砻江而下，冲毁房屋、田地、公路、桥梁等，给下游造成巨大损失。

2 PFC3D颗粒流参数标定

2.1 PFC3D软件

PFC3D是根据Cundall提出的细观离散元理论(粒子流理论)开发的一种数值计算软件。该软件模拟单元为三维球体颗粒(granular)，从微观结构的角度出发，研究颗粒集合体的破裂、破裂发展和大位移的颗粒流问题[28]。

模型的建立基于以下假设：① 颗粒单元被视为刚性体；② 颗粒间的接触发生在可忽略不计的点上；③ 颗粒间的接触视作柔性接触，颗粒允许在接触点上相互重叠；④ 重叠量与接触力大小有关，接触力通过力-位移定律计算，需要指出的是重叠量相比于颗粒尺寸非常小；⑤ 颗粒间的接触点可以存在连接(bonds)；⑥ 所有颗粒的形状都是球形，也可以用到clump生成任意的形状，在每个clump中颗粒相互重叠但没有相互作用力。

基本计算原理是运用显示中心差分法进行动态松弛，计算过程是在每个时步内交替利用力-位移定律与牛顿第二定律，时步迭代并遍历整个颗粒集合体(见图2)。

图2 PFC3D计算过程

材料的本构行为通过接触点上的本构模型来模拟，颗粒间的接触本构模型有3种：接触刚度模型、滑动模型和接触连接模型。

岩土体并非刚性体，系统的能量并不仅仅通过单元之间的摩擦方式进行耗散，颗粒与颗粒及颗粒与墙体之间的碰撞能耗也是系统能量损失的最重要途径之一。为了更真实地还原岩土体介质在运动过程中的碰撞问题，PFC3D提供了局部阻尼、黏性阻尼两类基本阻尼模型。局部阻尼的存在不适用于滑坡在重力作用下的动力学分析，模拟过程中将局部阻尼系数统一设为零；黏性阻尼通过模拟颗粒自由下落的“碰撞-反弹”，分析颗粒的速度恢复系数，确定出法向阻尼系数为0.4，切向阻尼系数为0.1。

2.2 参数标定

PFC3D模拟采用的是颗粒之间的微观参数，主要包括：球最小半径Rmin、球半径比Rmax/Rmin、球-球接触模量Ec、球刚度比kn/ks、球摩擦因数f、平行黏结半径乘子λ、平行黏结模量Ec、平行黏结刚度比kn/ks、平行黏结法向强度与切向强度σ和τ。试验获取的岩土体参数是宏观参数，岩土体的宏观参数与微观参数之间存在一定的联系，但目前尚没有成熟的理论和方法来界定这种联系。针对PFC3D微观参数的获取，笔者运用软件自带的虚拟三轴测试程序进行参数标定与反演，通过将反演参数与唐古栋滑坡综合物理力学参数进行对比，确定出微观物理力学参数。

参数标定按照以下过程进行：① 首先根据实际方量与计算机计算能力确定颗粒最小半径为3 m，最大颗粒与最小颗粒半径比为1.67；② 根据崩坡积物的变形参数,如弹性模量与泊松比反演出对应的法向刚度与切向刚度为1×106N/m；③ 采用平行连接模型模拟崩坡积物力学特性，根据强度参数凝聚力与内摩擦角反演出连接强度，分别为法向强度20 MPa，切向强度15 MPa；④ 颗粒摩擦系数会影响最后的堆积体形态，试算不同摩擦系数最终导致的堆积体形态，并与1967年唐古栋形成的堆积形态进行对比，摩擦系数的增加会导致最终堆积坝堆积坡度的增高，通过多次调整，最终得到相似堆积坝堆积坡度对应的摩擦系数为0.324，并将其用于本次数值计算中。模型参数见表3。

2.3 堰塞坝堆积反演分析

根据1967年唐古栋滑坡堆积堰塞坝形态反演，以验证标定细观参数的可靠性。

2.3.1构建模型

根据A区地形恢复唐古栋滑坡1967年以前的地形，1967年主要失稳范围B、C、D区的失稳方量约为9 260万m3。通过提取一定密度的高程点，运用三维地质建模，将每3个点组合成一个三角形墙面，最后将所有墙面有序连接。地形由18 409个三角形组成，滑体由半径3～5 m的171 343个球形单元组成，球与球之间设置平行连接，球与地表之间设置为接触刚度模型中的弹性模型(见图3)。

表3 唐古栋滑坡PFC3D数值模拟细观参数

图3 1967年失稳反演三维模型

2.3.2堰塞坝形态反演分析

在模拟进行到100 000时步时，滑坡运动达到收敛平衡，堆积形成堰塞坝(图4)，分别沿河流流向截取一条横剖面线，沿主滑方向截取一条纵剖面线，得到堰塞坝的横剖面形态图(图5)和纵剖面形态图(图6)。

图4 1967年失稳反演堰塞坝堆积形态

从堰塞坝横剖面形态图(图5)中可以看出，堰塞坝在河道内呈中部高两边低的梯形形态分布，坝底宽为1 483 m，坝顶宽为516 m，最大高度178 m，下游坡度为18°，上游坡度为22°。

图5 1967年失稳反演堰塞坝横剖面形态(单位：m)

从堰塞坝纵剖面形态图(图6)中可以看出，堰塞坝在滑坡对岸一侧高度明显较高，在靠近滑坡一侧高度较低，对岸坝高245m，本岸坝高178m。

图6 1967年失稳反演堰塞坝纵剖面形态(单位：m)

对比分析反演堰塞坝与1967年堆积堰塞坝的形态，结果见表4。沿河流流向的横剖面，堰塞坝底宽偏差1.13%，顶宽偏差40%，最大坝高偏差20.89%～34.07%；沿滑坡主滑方向的纵剖面，对岸坝高偏差26.87%，本岸坝高偏差1.71%。

通过上述对比，验证出标定的细观参数是合理的。1967年堆积堰塞坝的形态是通过调查残余痕迹获得，且坝顶宽残余形态保留较差，模拟堰塞坝与堆积堰塞坝上述3个数据的对比偏差达到20%～40%。

表4 堰塞坝堆积形态对比

3 堰塞坝堆积过程与形态分析

根据对唐古栋滑坡的地质调查及稳定性分析，滑坡A区(见图7)蠕滑体在地震工况下可能发生整体失稳，失稳方量达1 700万m3，规模巨大，故需对A区失稳后堆积形成的堰塞坝开展预测分析，为溃决分析、应急治理等提供依据。

3.1 构建模型

构建三维模型，地形由19 572个三角形组成，滑体由半径为3～5 m的35 109个球形单元组成，球与球之间设置为平行连接，球与地表之间设置为接触刚度模型中的弹性模型(见图8)。参数选取反演标定的细观参数。

图7 滑坡A区平面示意

图8 A区三维模型

3.2 失稳堆积过程分析

根据失稳过程中不同时刻滑坡形态的改变及运动的特点，分别选取6 000,10 000,20 000时步截取位移云图和纵剖面图进行堆积过程分析。

3.2.16000时步堆积分析

6 000时步滑坡整体下滑约500～600 m，部分后缘滑体绕过基岩沿着滑坡边界冲沟滑动，前缘部分岩土体有部分已滑入河谷处但是堆积方量不大。由于地形原因前缘滑体滑动速率较后缘滑体慢，大部分后缘及中部滑体集中于整个滑坡中部。滑动前期由于动能较大，此时岩土体移动速率较大(见图9,10)。

图9 6000时步位移云图(单位：m)

图10 6000时步纵剖面图(单位：m)

3.2.210000时步堆积分析

10 000时步滑坡整体下滑约1 000 m，此时滑体整体下滑速率依旧很快。中后部滑体在前缘地形较陡处发生飞跃的现象。下部河谷处堆积体方量仍不大，河谷处堆积体明显沿着河谷方向向两边扩展，堆积体主要由前缘及中部岩土体堆积而成，堆积高度约20～30 m左右(见图11,12)。

图11 10000时步位移云图(单位：m)

图12 10000时步纵剖面示意(单位：m)

3.2.320000时步堆积分析

20 000时步滑体运动速度开始明显减小，逐步趋于稳定状态；滑坡整体下滑约1 500～1 600 m左右。滑体主要集中在前缘及河谷处，在河谷处形成较为稳定的堆积形态，堆积体高度约为60 m。堆积体主要由滑坡前缘及中部的滑体形成，有少量后缘岩土体。由于受到深冲沟的限制，堆积体主要集中于A区、B区正下方，沿着河谷的扩展宽度并不是很大(见图13,14)。

图14 20000时步纵剖面示意(单位：m)

3.3 堆积形态分析

滑坡失稳经过80 000时步后达到稳定状态，在河谷处形成稳定的堰塞坝。岩土体整体滑移1 600 m左右。分别沿河流流向截取一条横剖面线，沿主滑方向截取一条纵剖面线(图15)，得到堰塞坝的横剖面形态图(图16)和纵剖面形态图(图17)。

图15 堰塞坝堆积形态

从堰塞坝横剖面形态图(图16)中看出:堰塞坝形态呈中部高两边低的梯形形态分布，坝顶宽437.91 m，坝底宽994.39 m，坝高92.65 m，下游坡度为18.6°，上游坡度为18.7°。

图16 堰塞坝横剖面形态(单位：m)

从堰塞坝纵剖面形态图(见图17)中看出:堰塞坝在滑坡对岸一侧堆积高度较高，在靠近滑坡一侧堆积高度较低，对岸坝高100.80 m，本岸坝高92.65 m。

图17 堰塞坝纵剖面形态(单位:m)

4 溃坝影响分析

基于模拟分析得到堰塞坝坝高为92.65 m，堰塞坝坝顶最低高程约为2 483 m，当上游持续来流导致水位漫顶时，堰塞坝前最大水深约为93 m，回水长度约为33 km。

采用DAMBRK溃坝洪水计算数学模型[29]，进行楞古水电站施工期和运行期的唐古栋滑坡A区失稳堆积堰塞坝溃坝影响预测分析。溃决起始水位取为堆积体坝顶高程2 483 m，溃口尺寸参照1967年唐古栋垮山滑坡溃坝，考虑为局部溃决，溃口深度按相同比例推算为55.8 m，溃口顶宽155 m，溃口底宽50 m，溃口边坡坡度取1967年溃坝自然发展的左岸边坡坡度，坡比为1∶0.94。计算结果见表5。

(1)施工期溃坝洪水计算结果表明：在各种流量方案下，强变形A区失稳堆积堰塞坝溃坝洪峰流量均远大于楞古水电站拟选坝址校核洪水流量；上游不同来流方案对溃口最大洪峰的影响较小，各方案差值均在3 700 m3/s以内；当遭遇5a一遇洪水流量时，溃口洪峰流量为25 200 m3/s，约为楞古水电站拟选坝址校核洪水流量8 690 m3/s的2.9倍。

表5 A区失稳堆积堰塞坝溃坝计算方案及参数

(2)运行期溃坝洪水计算结果表明：在各种流量方案下，强变形A区失稳堆积堰塞坝溃坝洪峰流量均远大于楞古水电站拟选坝址校核洪水流量；上游不同来流方案对溃口最大洪峰的影响较小，方案间差值均在3 700 m3/s以内；当遭遇5a一遇洪水流量时溃口洪峰流量为13 000 m3/s，约为楞古水电站拟选坝址校核洪水流量8 690 m3/s的1.5倍。

5 结 论

本文运用三维离散元软件PFC3D，以唐古栋滑坡为研究对象，通过反演模拟1967年滑坡失稳堰塞坝堆积形态，验证标定了离散元细观参数。并基于此细观参数进行强变形A区失稳堰塞坝的预测模拟。在此基础上，开展堰塞坝溃坝影响预测分析，主要研究结论如下。

(1)基于滑坡室内物理力学试验和抗剪强度反演得到的唐古栋滑坡体综合物理力学参数，运用PFC3D软件进行虚拟三轴测试，标定了一组细观力学参数。

(2)通过对比分析1967年失稳堆积堰塞坝与PFC3D反演堰塞坝的形态，校核标定的细观参数是合理的。堰塞坝底宽和本岸坝高偏差较小，顶宽、最大坝高、对岸坝高的偏差20%～40%，是由于1967年堆积堰塞坝的形态是通过调查残余痕迹获得，且坝顶宽残余形态保留较差。

(3)强变形A区失稳堆积过程如下：A区整体下滑约1 600 m；滑动前期，岩土体移动速率较大，前缘滑体滑动速率较后缘滑体慢；滑动后期，堰塞坝主要由滑坡前缘及中部的滑体形成，有少量后缘岩土体；受冲沟限制，堰塞坝沿河谷的扩展宽度并不大。

(4)强变形A区失稳堆积堰塞坝的形态，顺河向呈中部高两边低的梯形形态分布，坝顶宽437.91 m，坝底宽994.39 m，坝高92.65 m，下游坡度为18.6°，上游坡度为18.7°；横河向对岸一侧堆积高度较高，对岸坝高100.8 m，本岸坝高92.65 m。

(5)采用DAMBRK溃坝洪水计算数学模型，选取年平均流量、汛期平均流量和5 a一遇洪水流量，进行楞古水电站施工期和运行期的溃坝影响分析。在各种流量方案下，堰塞坝溃坝洪峰流量均远大于楞古水电站拟选坝址校核洪水流量，上游不同来流方案对溃口最大洪峰的影响较小；强变形A区失稳导致的堰塞是影响楞古水电站坝址选择的重要因素。

(6)通过模拟发现，PFC3D软件用于模拟滑坡堰塞坝堆积过程及堆积形态有较好的适用性，尤其可以获取堰塞坝堆积的三维形态，且堰塞坝高度较一般经验公式更为准确，可有效运用于溃坝影响分析。