某堰塞坝冲刷溃决数值模拟分析

吕 霞，杨兴国，2，范 刚，2，林子钰，刘大瑞

(1.四川大学水利水电学院，四川 成都 610065；2.四川大学水力学与山区河流开发保护国家重点实验室，四川 成都 610065)

0 引言

堰塞坝(体)是在一定的地形地貌条件下，由滑坡、泥石流、熔岩流等堵塞山谷和河道形成的天然坝体[1-3]。堰塞坝是自然作用的产物，具有坝体形态不规则、物质组成随机性强、内部结构及空间分布复杂等特点，受上游来流、渗流侵蚀等外载荷的影响，容易在短期内发生失稳破坏，对下游社会和自然环境造成重大危害[4-7]。1974年秘鲁Mayunmarca村发生滑坡-碎屑流灾害，导致450人死亡，随后碎屑物质堵塞河流形成堰塞坝，溃决后造成下游巨大灾难[8]。2008年汶川8.0级大地震诱发多处堰塞坝，其中唐家山是规模最大、危险性最高的堰塞湖，严重威胁绵阳市城区及下游130万人的生命安全[9]。2018年8月6日，汉源县富泉镇西沟发生150万m3的大规模滑坡，并在主沟内形成体积约40万m3的堰塞坝[10]。2018年10月11日及11月3日，金沙江白格滑坡两次堵塞金沙江，形成的堰塞湖淹没了上游村镇，堰塞坝溃决后洪水冲击下游造成了巨大损失[11]。2020年6月17日，四川省丹巴县梅龙沟发生特大泥石流灾害，泥石流局部堵塞小金川河河道，形成堰塞坝。坝体溃决后，导致电站大坝和大量房屋被淹，直接经济损失约2000万元[12]。

美国地质调查局研究发现，漫顶导致的溃坝案例高达90%，而渗流破坏导致的溃坝数量仅占10%[13]。因此有必要深入研究堰塞坝漫顶溃决过程以及溃决机理，为溃坝洪水风险分析和应急预案的编制提供理论和技术支持，以期削弱溃决洪水对下游的破坏程度。目前，已有众多学者对堰塞坝溃决过程和机理分析做了大量研究。赵天龙等通过不同来流条件下粗粒级配土料的冲刷性能实验，揭示坝体宏观溃口的形成机制[14]。邓明枫等通过坝体漫顶溃决模型实验，研究了堰塞坝漫顶溃决的动力机制[15]。傅旭东等采用Osman和Thorne模型发展了堰塞坝溃决过程模拟模型，并将其应用于唐家山堰塞湖溃决过程的模拟，再现了溃决洪水流量过程线、溃口展宽和下切过程[16]。Zhang等在天然河道中开展了大型物理模型，提出了堰塞坝纵向以及横向溃决演化模型，并利用白格滑坡的两次溃坝事件，验证了该演化模型的准确性[17]。刘若星等将NWS BREACH模型应用于唐家山堰塞湖泄流过程的模拟，验证了改模型的合理性[18]。杨兴国等构建了考虑水流侵蚀与溃口间歇性崩塌的堰塞坝溃决演化模型，初步揭示了冰碛土滑坡-泥石流运移与堵江机制[19]。目前，针对堰塞坝溃决机理和溃决过程的研究大多未使用堰塞坝的实地模型，而对坝体做了相应的概化，其研究成果可为堰塞坝溃坝机制研究提供重要借鉴，基于概化模型的研究成果需要实际堰塞坝模拟结果的检验与验证。

2020年8月31日四川省甘洛县黑西洛沟发生滑坡-泥石流灾害，沟内冲出大量固体物源堵塞尼日河并形成堰塞坝。本论文通过等高线地形图，建立了黑西洛堰塞坝实地模型，使用Flow-3D软件对堰塞坝的溃决过程进行数值模拟，得到坝体泄流时的流速特征和溃口变化过程，以进一步了解溃口的发展规律，为堰塞坝风险分析和应急预案的编制提供技术支撑。

1 黑西洛滑坡-堰塞湖灾情概况

2020年8月31日，四川省甘洛县阿兹觉乡发生山体滑坡-泥石流-堰塞湖灾害链，尼日河被大量泥石流物质堵塞，形成堰塞坝，严重威胁当地人民生命财产安全。黑西洛沟道两侧松散物源分布广泛，植被不发育。灾害发生时，沟内的松散物质被山洪裹挟带走，并在搬运过程中不断铲刮底层物源，从而演变为滑坡和崩塌灾害。在沟内通道持续下切的同时，两侧岸坡持续垮塌，坍塌体进入沟道后，滑坡规模急速扩大，最终演变为泥石流灾害。大量的泥石流碎屑以近乎垂直的方式冲入尼日河，阻塞尼日河，形成堰塞坝。经现场实测，堰塞坝坝体纵向方向长度约200m，顺河向长度约为400m，高度约30m，堰塞坝体积约100万m3。堰塞坝形成约20min之后，尼日河上游水位迅速上涨，达到了772.0m，上游河水漫过堰塞坝，并在坝体中偏左侧低洼地带形成天然的泄流槽，由于坝体材料松软，堰塞坝在10min左右即溃决完毕。堰塞坝溃决后，该段河道截弯取直，近乎顺流，如图1所示。坝体溃决后冲毁下游场镇、学校、房屋和公路，造成阿兹觉乡受灾人口1730人，3人失联。沟口的成昆铁路大桥被冲毁，坝体下游约1.2km的道路以及多处桥梁被冲毁，经济损失严重[20]。

图1 2020年黑西洛沟滑坡-泥石流-堰塞湖灾害情况

通过现场调查和资料收集，本次滑坡-泥石流-堰塞湖灾害链是典型的“小水大灾”灾害，灾害发生时大量松散堆积物被冲出沟道，堵塞尼日河，形成具有松散结构的堰塞坝，并且堰塞坝在短时间内溃决，坝体上游库区内的洪水快速下泄，造成下游严重的生命财产损失。本文采用现场实测地形，使用数值模拟方法重现和分析黑西洛堰塞坝冲刷溃决过程，以期为今后类似堰塞坝的应急处置和编制避险预案提供技术支持。

2 数值计算模型

本文以黑西洛沟堰塞湖为研究对象，利用Flow-3D数值方法，开展堰塞湖溃决过程数值模拟研究。Flow-3D是基于CFD解算技术的仿真模拟软件，其计算原理是有限差分法，采用质量连续性方程和Navier-Stokes方程作为控制方程。其FAVOR技术既能进行网格内的流场计算，又能进行固态流场的数值模拟，特别是对自由表面流的模拟具有显著优势。另外，该技术还可以通过简单的矩形网格来构建任何复杂的形状，并通过VOF对自由面进行追踪，从而更接近于真实的流体运动，并能够对流体的多种特性进行模拟。堰塞坝溃决后，洪水迅速冲出，对坝体物质产生强烈的冲蚀作用，考虑到下游弯曲河道处存在紊流现象和剧烈变形等特点，采用RNG湍流模型对溃坝洪水进行模拟，该模型能很好地模拟溃决过程中洪水的复杂流态。此外，将坝体设置成泥沙冲刷模式，以模拟其在水流冲击作用下逐渐溃决的过程。Flow-3D软件中的泥沙冲刷模型，通过对泥沙起动、堆积、推移质输移等计算来描述泥沙运动，与坝体冲刷破坏机理一致。

本文根据获取的等高线地形图，建立了黑西洛堰塞坝的三维实体模型，并将其导入Flow-3D软件中。模型的计算范围长900m，宽约460m，高约125m，如图2所示。本文使用正交的矩形网格对模型的计算区域进行划分，网格大小为2.5m×2.5m×2.5m，计算网格总数为331万个。进一步对计算区域进行边界条件的设定，将模型入口边界设为压力边界，出口边界设为自由出流边界，底面和两侧边界设为Wall(墙)边界，自由液面上方设一个标准大气压。并将坝体上游初始水位条件设为772.0m，为上游库区水位达到的最大值。堰塞体物质来源于沟内松散堆积体，颗粒级配曲线如图3所示[20]，将6个监测点颗粒中值粒径的平均值作为本次模拟计算的典型坝体粒径，即d50=11.75mm。本次设定堰塞坝堆积区域泥沙区域，泥沙密度为2000kg/m3，临界Shields系数为0.32，挟带系数设为0.013。根据现场调查和模型试算，坝体在600s后基本稳定，故设定求解时间为600s，最小步长为1×10-7s。

图2 堰塞坝的数值模拟模型

图3 堰塞体颗粒级配测点位置及级配曲线

3 模拟结果分析

3.1 流速

图4是堰塞坝溃决过程中水流流速的分布情况。堰塞坝形成泄流通道后流速较高，流道中后段流速最大为17.5m/s。因流道内的水流流速较大，水流的侵蚀能力较强，溃口随时间不断扩展。坝体溃决持续76s后，堰塞体形成贯通的流道，流道内流速高达17.5m/s；溃决持续203s后，流速在流道内分布不均匀，流速相对较高，水流的侵蚀作用强烈，流道不断下切变宽；溃决持续时间417s后，坝体泄流过程中，上游水位在逐渐降低，导致泄流槽内流速开始缓慢下降，此时流速最大达到15.3m/s。如图4所示，泄流过程中流速分布有一定相似性，水流进入泄流槽之前，流速基本为0，上游库区处于静压平衡状态；水流流入泄流槽后，进口处水流流速明显上升，整个溃决过程中，泄流槽中后段水流流速一直较快。泄流槽末端，槽内斜坡道上的水流流速急速上升，水流侵蚀作用增强，不断侵蚀坝体，导致泄流槽不断扩展，且末端斜坡道上的跌坎不断向上游推进，泄流后期，跌坎消失。

图4 堰塞坝冲刷溃决过程水面流速动态演进模拟结果

3.2 溃口扩展过程

3.2.1冲淤形态

水流运动时会对泥沙颗粒产生上举力和拖曳力，当水流强度达到泥沙起动条件时，能使泥沙颗粒跃起，并被水流带走，坝体产生冲刷现象，导致泄流槽底部逐渐下切。图5是坝体溃决过程中泥沙的冲淤演进变化过程，其中正值表示淤积，负值表示冲刷。如图5所示，溃决持续发展过程中，坝体的淤积高度和冲刷深度不断增加。持续冲刷76s后，泄流槽整个连通并不断扩展，水流流速高达17.5m/s，水流的挟沙能力强，冲刷速度快，此时堰塞体的最大冲刷深度达22.7m；坝体溃决持续203s后，水流不断带动泥沙向下运动，堰塞坝的冲刷深度持续增大，泄流槽底部下切和两侧侵蚀明显，此时堰塞体的冲刷深度达25.2m；溃决持续417s后，上游水位降低，泄流槽内流量变小，流道中水面宽度和水深逐渐降低，泄流槽底部下切和侧壁扩展速度变慢，坝体的冲刷深度增长变缓，此时坝体最大冲刷深度达25.7m。在泄流槽内，各位置的冲蚀状况也有差别，泄流槽进口处以及坝体中上游区域的冲刷深度变化缓慢，是由于该位置的水流流速较小，而泄流槽下游段流速变高，水流携带大量泥沙运动，导致下游坝体的冲刷深度和溃口展宽都较上游面更大，泄流槽明显下切，并呈向下的“喇叭口”形状。

图5 堰塞坝冲刷溃决过程泥沙冲淤动态演进模拟结果

堰塞坝下游，受重力以及床面摩擦等因素的影响，被水流带走的泥沙在此处沉积，出现淤积现象。泄流初始阶段，水流携带着流道内的泥沙向下游运输，由于水流流速不大而堆积在近坝体河段。如图5所示，堰塞坝持续泄流76s后，河道内出现淤积现象，最大淤积高度达到4.5m；随着水流持续下泄以及溃口不断发展，靠近坝体的河段也被水流冲刷，初期沉积在此处的泥沙被冲向下游。泄流后期，流道内流速基本达到稳定，流道内冲刷深度变化不明显，泥沙不断淤积在坝体下游，堆积高度最大值为8.5m。泥沙在坝体下游左岸堆积，是因为黑西洛堰塞坝处在弯曲河段，溃决洪水产生横向环流对岸坡形成顶冲之势，使得凹岸冲刷而凸岸淤积。水流流经弯曲河道，由于离心作用，导致下游左岸水深变浅，流速降低，水流的搬运能力下降，泥沙在河道左侧逐渐堆积。

3.2.2溃口断面历时变化

对黑西洛堰塞坝溃决的横、纵剖面进行截取，进而揭示堰塞坝漫顶溃决的破坏规律。截面具体位置如图6所示，其中1-1′截面是纵剖面，坐标为x=200m，2-2′截面、3-3′截面、4-4′截面为泄流槽不同位置的横断面，坐标分别为y=250m、y=430m、y=600m。

图6 截面示意图

堰塞坝纵剖面(1-1′截面)不同时间的冲淤变化如图7所示，色条表示堰塞坝的冲淤情况，其中正值代表淤积，负值代表冲刷。由图7可知，在堰塞坝泄流槽斜坡道上，冲刷现象最显著，泄流槽末端水流流速急速上升，冲蚀能力明显增强。泄流槽末端不断被水流侵蚀破坏，最大冲刷深度达到25.0m。漫顶水流侵蚀斜坡段，使其产生跌坎，水流流经跌坎，流速陡然增加，挟沙能力变强，致使跌坎处的溃决速度明显更快。在跌坎的上游，水流较稳定，流速小，其冲刷强度较低；跌坎的下游，水流湍急，水流流速较大，水流对坝体产生强烈的冲蚀作用，使跌坎附近的泥沙不断被冲走，泄流过程中跌坎不断向上游推进，产生明显的溯源冲刷。溯源侵蚀会在短期内对堰塞坝产生强烈的冲蚀作用，由图7可知，溃决持续381s后，溯源侵蚀已经发展到坝体上部，跌坎消失。

图7 1-1′截面冲淤演变过程

图8—10分别展示了堰塞坝不同横断面(2-2′截面、3-3′截面、4-4′截面)在不同时间的冲淤演进过程。如图8所示，该断面侵蚀程度较弱，由于溃口上游段流速较小，水流的侵蚀能力不强，因此流道发展比较缓慢，坝体的冲刷深度最大为4.7m；图9为3-3′横断面溃口的演变过程，该断面流速增大，水流冲刷能力变强，水流携带更多的泥沙向下运动，持续冲刷417s之后，此截面最大冲刷深度达到了21.3m；如图10所示，水流在泄流槽末端侵蚀作用最强烈，此断面位置的流速明显增大，冲刷能力强，最大冲刷深度达25.7m，溃口扩展最为迅速，且泄流槽宽度不断拓宽。

图8 2-2′截面冲淤演变过程

图9 3-3′截面冲淤演变过程

图10 4-4′截面冲淤演变过程

整个泄流过程中溃口的发展过程如下，泄流槽下游端冲刷最强烈，溃口向底部下切和两侧侵蚀，之后逐渐向上游发生侵蚀，当侵蚀扩展到坝体上游后，又开始向坝体两侧缓慢侵蚀。水流挟带泥沙向下游运动不断侵蚀坡脚，侧壁发生破坏，导致泄流槽的宽度逐渐加宽。总的说来初始溃口发展较为缓慢，然后逐渐加快，最后流道内发展达到平稳状态。

4 结论

本文采用Flow-3D软件对黑西洛沟堰塞湖自然泄流溃决过程进行模拟和分析，得出以下结论：

(1)堰塞坝漫顶冲刷是一个复杂过程，坝体溃决过程中流速变化大且分布不均匀，在泄流槽的下游段流速明显增大，水流的冲刷能力增强。

(2)泄流过程中，泄流槽末端斜坡道上流速最大，坝体的侵蚀破坏最明显，流道下切深度达25.7m，并在堰塞坝下游出现淤积现象，最大淤积高度达8.5m。

(3)堰塞坝泄流过程中出现了溯源侵蚀现象，水流不断使溃口加宽、加深，泄流跌坎不断向坝体上游移动。

(4)黑西洛堰塞坝自然泄流过程中，初期溃口发展缓慢，水流对坝体的冲刷破坏主要发生在泄流槽末端，随后向其两侧和上游侵蚀。