改进DB-IWHR模型及其在尾矿库溃坝影响分析中的应用

2022-09-12沈鸿杰徐力群刘子茜

水利水电科技进展 2022年5期

沈鸿杰，徐力群，刘子茜

(河海大学水利水电学院，江苏南京 210098)

尾矿库作为一种蕴含着很高势能的人工泥石流危险源，当发生溃坝事故时，大量下泄尾砂在极短时间内冲向下游，威胁下游居民和生态环境[1]。据统计， 2001—2015年我国共发生尾矿库溃坝事故65起，其中洪水漫顶占事故总数的36.9%[2]。洪水漫顶是导致尾矿库发生溃坝的主要原因之一，较为准确地预测其溃坝过程及影响范围，是制定应急预案的基本条件[3]。

尾矿库的溃坝过程类似于散粒体土石坝及堰塞坝的溃坝，尾矿库溃坝计算理论的研究多参考水库溃坝实例[4]。近年来，不少学者基于尾矿库坝体材料的特殊性开展了相关的尾矿库溃坝模型研究。袁兵等[5]根据多个大坝溃决实测资料，建立了尾矿坝溃坝数学模型，提出了尾矿库溃坝沙流对下游影响的预测方法。刘磊等[6]通过建立尾矿坝溃决物理模型对溃坝过程进行观测，依据试验结果建立了尾矿库漫顶溃坝洪水预测数学模型。武立功等[7]采用不同粒径的尾砂作为筑坝材料，研究了尾砂粒径对尾矿库漫顶溃坝的影响。

基于唐家山堰塞湖实测溃决资料，Chen等[8]开发了DB-IWHR模型，提出了以双曲线形式代替指数形式的冲刷模型，引入了简化Bishop法模拟溃口扩展过程，并采用速度增量形式直接求解控制方程，避免了数值迭代。张强[9]运用DB-IWHR模型对“10·10”和“11·3”白格堰塞湖溃决过程进行了反演分析，结果表明模型可较为准确地计算洪水演进过程。张贵金等[10]采用DB-IWHR模型计算堆石坝溃坝洪水过程，利用MIKE21模拟洪水演进过程，评估了极端条件下堆石坝溃坝对下游的淹没风险。

虽然尾矿库溃坝过程与土石坝相似，但尾矿库建造过程和条件与土石坝不同，目前的溃坝溃口数值分析模型中很少考虑到尾矿库的自身特点，直接采用现有土石坝模型进行溃坝洪水演进计算，会产生较大误差。此外，据尾矿库溃坝调研实际情况，大量尾矿库发生溃坝后，总泄沙量约为实际库容的1/3，极少出现全部库容尾砂下泄的情况。为此，本文根据尾矿库结构和泄沙特点，提出了一种改进模型，以合理模拟尾矿库溃坝过程，准确预测尾矿库发生溃坝后对下游的影响，为尾矿库溃坝的防灾减灾提供理论依据。

1 改进DB-IWHR模型

1.1 尾矿库结构特点

尾矿库筑坝方法可分为上游式筑坝法、下游式筑坝法、中线筑坝法、高浓度堆积法和水库式堆积法5种。我国目前修筑的尾矿库多采用上游式筑坝法，约占金属矿山尾矿库的95%。上游式筑坝法修筑的尾矿库一般由初期坝和多级子坝构成，初期坝一般为透水堆石坝，而后期堆积坝主要利用库内尾矿料填筑而成，随着生产的进行，尾矿坝不断加高形成多级子坝。在结构特点及溃坝分析方面，尾矿坝与土石坝相比具有明显差别，主要表现在以下几个方面：

a.尾矿坝的安全性普遍低于土石坝。现阶段我国的尾矿坝多采用上游式筑坝法，其筑坝方法简单、生产管理方便，但存在抗地震液化能力差、排渗系统易淤堵导致浸润线抬高等问题，影响尾矿库安全性。

b.尾矿库内沉积尾矿料物理力学性质的时空变化规律更为复杂。土石坝坝体通常采用分区设计，同一分区内材料的级配、压实度和物理力学特性基本一致；而尾矿料在水力沉积作用下尾矿颗粒的空间差异性和沉积尾矿在自重荷载和上覆荷载作用下的固结变形和强度变化产生的时间变异性，使其物理力学特性的时空变化规律与土石坝相比更为复杂[11]。

c.尾矿库溃坝下泄流体性质更为复杂。土石坝溃决后下泄至河道的流体主要为库水，尾矿库溃决后主要下泄物为水砂混合物，其流变性质与泥石流具有一定相似性，目前对尾砂浆体流变性的研究很少，且尾砂浆体的流体参数选取和数值模型没有统一的标准[12]。

1.2 改进DB-IWHR模型的溃口流量计算

1.2.1溃口流量

溃口流量按宽顶堰公式[8]进行计算，即：

式中：Q为溃口流量，m3/s；C为综合流量系数，m1/2/s；B为溃口宽度，m；z为溃口底高程，m；H为水库水位，m；mb为宽顶堰流量系数；mq为侧收缩系数；g为重力加速度，m/s2。

1.2.2溃口冲刷

溃口冲刷数值分析采用双曲线模型[8]，公式如下：

(2)

其中τ′=k(τ-τc)

式中：r为侵蚀率，mm/s；τ为剪应力，Pa；a、b为冲刷侵蚀参数；k为在剪应力范围内允许r接近其极值的单位变换因子；τc为临界剪应力，Pa。

临界剪应力是判断溃坝溃口开始与终止发展的重要参数，DB-IWHR模型中一个重要假定就是坝体组成物质较均匀。对于上游式筑坝法修筑的尾矿库，初期透水坝一般为碾压堆石坝，在现代施工条件下，孔隙率较小，稳定性较好，不易被水流冲刷破坏，而尾矿堆积坝主要利用库内尾矿料填筑而成，粒径较小，临界剪应力也较小，遇洪水漫顶时易被冲刷发生溃坝。因此式(2)中τc的取值不符合尾矿坝实际情况。

为使模型更适用于尾矿库漫顶溃坝的实际情况，对DB-IWHR模型进行改进，即在不完全溃坝模式下，根据尾矿坝剖面进行材料分区，并采用修正Shields曲线[13]计算临界剪应力τc。当τ>τc时，溃口开始冲蚀；当τ<τc时，溃口完全发展，冲蚀结束。临界剪应力τc的计算公式为

(3)

式中：D*为无量纲泥沙尺寸参数；Re*为颗粒雷诺数；θc为希尔兹数；γs为尾矿砂容重，N/m3；γ为尾矿砂和水混合物容重，N/m3；D为尾矿颗粒粒径，m；υ为水的运动黏度，m2/s。

根据曼宁公式，可以得到临界剪应力与临界流速的关系：

(4)

式中：R为水力半径，m；J为水力坡度；h为水深，m；vc临界流速, m/s；n为曼宁系数，s/m1/3。

修正的Shields曲线是由钱宁等[14]在原始Shields曲线基础上，在Re*<2的层流区域及Re*很大的紊流区域内对θc的取值进行修正所得。修正的Shields曲线根据流体和颗粒的特性引入了无量纲泥沙尺寸参数D*，避免了数值迭代，简化了求解过程。

表1 改进DB-IWHR溃口模型参数

1.2.3溃口扩展

溃口的边坡稳定性分析采用简化的毕肖普圆弧滑动面分析方法[8]。在实际应用中将计算所得图形假设为一系列梯形，采用双曲线模型逐级计算溃口宽度和倾角[4]，溃口宽度计算公式为

(5)

式中：B0为溃口初始宽度，m；Δz为溃口底高程增量，m；β为溃口边坡倾角，(°)。

1.2.4初始条件

溃口流量可以通过单位时间内水库库容的损失来确定，根据质量守恒定律，可得水量平衡方程：

(6)

式中：W为库容，m3；t为时间，s；q为入库流量，m3/s。

在单溃情况下，假设溃坝起始条件为

(7)

可得初始溃口宽度及初始溃口底高程，计算公式为

(8)

(9)

式中：z0为初始溃口底高程，m；H0为初始库水位，m；q0为天然入流量，m3/s；m为跌落系数。

1.3 尾矿库溃坝下泄沙流演进模拟方法

通过上述模型可以模拟尾矿库溃坝溃口发展情况，预测下泄流量过程，但无法定量分析溃坝对库区下游的影响。因此需要根据下游库区地形，选取合适的计算模型，模拟尾矿库溃坝下泄沙流的演进过程，为尾矿库风险预警提供数值参考依据。

本文基于分段模型法模拟尾矿库漫顶溃坝下泄沙流演进过程，具体步骤为： ①根据尾矿库相关参数，利用改进DB-IWHR模型计算尾矿库漫顶溃坝溃口演化情况和下泄流量过程线；②基于尾矿库下游地形资料，确定模型范围，建立非结构化网格，并生成尾矿库下游影响区地形模型；③设置边界位置、初始条件、糙率等相关参数，将计算得到的下泄流量过程线作为二维水动力模型的入口边界条件；④通过二维水动力模型计算求解，模拟尾矿库溃坝下泄沙流向下游的演进情况，并对模拟结果进行分析，得到尾矿库漫顶溃坝对下游的影响，评估尾矿库的溃坝风险。

1.4 改进模型溃口演化模拟精度验证

以唐家山堰塞坝的泄流过程模拟为例验证改进DB-IWHR溃口模型的合理性，模型的参数选择参考相关资料[15]，具体参数如表1所示，原始模型和改进模型模拟结果对比如图1所示。

图1 唐家山堰塞坝溃坝原始模型与改进模型模拟结果对比

通过下泄流量过程(图1(a))可以得出，改进模型与原始模型对洪峰流量和洪峰流量到达时间的模拟结果差别不大，改进模型模拟的洪峰流量为7 423.33 m3/s，洪峰到达时间为5.7 h。在到达洪峰流量之后，原始模型模拟的流量过程线随着时间变化与实测值误差逐渐变大，从溃坝开始至14 h时，实测总下泄量为12 434万m3，原始模型模拟的总下泄量为20 565万m3，比实测值大65%，而改进模型模拟的总下泄量为13 891万m3，仅比实测值大12%，因此，模拟的下泄流量过程更为合理。通过库水位过程(图1(a))可以得出，对于实测数据，库水位从742.5 m下降至719.48 m，原始模型库水位模拟值下降至698.76 m，远小于实测值，改进模型库水位模拟值下降至719.10 m，误差为0.5%，模拟精度明显提高。

通过溃口底高程和溃口宽度变化过程(图1(b))可以得出，原始模型模拟的溃口底高程从740 m下降至700 m，属于完全溃坝，与实际明显不符。改进模型通过控制最终溃口底高程，在底高程下降至719 m时停止纵向冲刷。由于最终的溃口宽度是根据边坡稳定计算得出，因此原始模型和改进模型模拟结果较为一致，与实测数据基本符合。当溃口发展到一定深度时，其纵向下切和横向扩展均已停止，溃口不再发展，因此导致原始模型模拟值与实测值误差较大。而改进模型通过修正临界剪应力控制最终溃口形态，模拟结果更符合实际情况，采用改进模型进行溃坝沙流演进模拟能更准确地预测尾矿库漫顶溃坝对下游的影响。

2 工程应用

2.1 工程概况

某拟建尾矿库位于湖北省宜都市丘陵地带，库区东、西、南三面环山，在北部地形狭窄处修建尾矿库。初期坝为碾压堆石坝，坝高58 m，坝顶高程340 m；堆积坝采用上游式筑坝法，最终堆积标高为360 m，总坝高为78 m，坝顶轴线长约422 m，设计最终总库容为1 222.4万m3，尾矿库等级为三等。堆积坝分为尾粉砂和尾粉土两个区，尾矿库材料分区见图2。

图2 尾矿库剖面及材料分区

2.2 模型建立

2.2.1二维水动力模型

根据尾矿库下游周边环境情况确定模型模拟的范围，依据地势下降最快方向，选取自初期坝坝址向下游延伸2 km作为模型范围。在初期坝下游165 m、660 m、1 000 m、1 275 m、1 555 m处共有5处居民区，分别定义为居民区1、居民区2、居民区3、居民区4和居民区5，居民区建筑物最低高程分别为239.4 m、207.2 m、194.1 m、189.5 m、和190.8 m。用模拟范围内的等高线地形图对模拟区域进行网格划分和地形插值，本次建模依据为该尾矿库下游沟谷真实地形，按1∶1比例建立地形模型，如图3所示。

图3 尾矿库下游影响区地形

2.2.2溃口下泄流量过程

采用湖北省某尾矿库设计资料作为溃口流量过程模拟的输入参数，应用改进DB-IWHR溃口模型对溃口演化过程进行分析，将得到的溃口流量变化过程作为后续模拟的边界条件，根据GB 50863―2013《尾矿设施设计规范》及设计报告确定冲刷侵蚀参数。由于构筑尾矿坝的尾矿砂粒径较小，计算得到的临界剪应力远小于唐家山堰塞湖，因此尾矿库的临界流速小于唐家山堰塞湖模拟中的取值。模型具体参数见表1。

计算得出的溃坝水位和下泄流量过程如图4所示。模拟初期溃口侵蚀较慢，溃口流量增长较为缓慢，水位变化小；约15 min时，溃口发生崩塌，溃口流量迅速增大，直至29.5 min时，溃口流量达到峰值1 420 m3/s；随后下泄流量逐渐减小， 60 min时，水位稳定在340 m，溃口发展达到稳定状态，溃坝过程结束。整个溃坝过程共60 min，峰值发生在29.5 min。溃坝过程主要分为3个阶段，即初始缓慢泄流、流量快速增长至峰值、下泄沙流逐渐落峰。

图4 溃坝下泄流量过程及库水位过程模拟结果

2.2.3模型参数设定及边界条件

根据溃坝下泄流量过程确定模型模拟时间为120 min，重点分析溃坝下泄沙流对下游的影响。尾矿库下游主要为农田及灌木，根据天然河道糙率表，取平均糙率0.04。由于下游无河道，因此将初始水深设置为0。

图5 溃坝下泄沙流演进过程

边界条件包括入口边界和出口边界，选取初期坝轴线位置作为模型的入口边界，采用模拟得到的下泄流量过程线作为边界条件；选取尾矿坝下游2 km处作为模型出口边界并设置为自由出流。

2.3 模拟结果与敏感性分析

2.3.1溃坝下泄沙流演进过程

不同时刻下泄沙流的淹没范围和淹没深度如图5所示。溃坝开始初期溃口流量较小，淹没深度及淹没范围较小，在29.5 min时溃口流量达到峰值，溃坝沙流开始快速向下游演进并迅速堆积，在41.5 min时下泄尾砂演进至下游边界处，由于溃口流量的下降及下泄沙流从下游出口不断流出，溃坝淹没深度逐渐下降。

2.3.2溃坝下泄沙流对下游影响分析

根据溃坝下泄尾砂演进路径，绘制各敏感点淹没深度随溃坝演进的时间变化曲线，如图6所示。从下游各敏感点淹没深度随演进时间的变化曲线来看，随着演进时间的增加，淹没深度逐渐增大，其间某一时刻流量达到最大，淹没深度达到峰值之后逐渐减小，最后趋于稳定。

图6 各敏感点淹没深度及流速随演进时间变化曲线

从溃坝沙流演进淹没模拟结果来看，溃坝沙流在坝址处流量大，流速快，淹没深度大，破坏性强；距离坝址越远，溃坝沙流到达时间越晚，最大淹没深度越小，最大流速也越小。对于距初期坝165 m的居民区1，其淹没深度最大达到9.21 m，而对于距初期坝较远的居民区5，其最大淹没深度降低至2.77 m，未达到建筑物最低高程。下泄尾砂大多堆积在坝址附近，对1 km外的主要建筑物影响较小，符合沙流演进的一般规律。

2.3.3改进模型敏感性分析

考虑溃口初始宽度的影响，分别取B0=3.29 m、6.58 m、13.16 m；考虑冲刷侵蚀参数的影响，分别取a=1.0、b=0.000 1，a=1.1、b=0.000 3，a=1.1、b=0.000 7，根据下泄流量过程曲线判断模型对其的敏感性。参数敏感性分析曲线如图7所示。