基于DL Breach模型的尾矿库溃坝影响数值分析

宁 伟 李 瑜 吴蒙蒙 杨 超

（1.洛阳栾川钼业集团股份有限公司；2.栾川县环境保护局；3.河海大学水利水电学院）

随着我国基础设施建设行业的大力发展，对各种金属和非金属矿产资源的需求也迅速增加。在矿产资源的开发和利用中，不可避免地产生大量以浆状形式排出的废弃矿渣，一般的矿山企业都会选择建设尾矿库堆放，来减少尾矿渣给环境带来的污染。尾矿库是一个重大的危险源，一旦发生溃坝，产生的后果相当严重［1］。因此，进行尾矿库溃坝数值模拟研究，分析尾砂演进与泥石流到达时间特性，可以为评价尾矿库建设方案的安全性和合理性提供依据。

尾矿库溃坝是一个复杂的水沙耦合过程，其研究成果丰硕。尹光志等［2］运用数值模拟和室内模拟的方法，对云南玉溪某尾矿库进行溃坝砂流演进规律模拟，探讨泥沙在下游的淹没范围。袁兵等［3］通过研究几个尾矿库的溃决案例，考虑尾矿渣的物理性质和流动过程中的变形来建立数学模型，对尾矿库溃坝灾害进行研究。林江等［4］等利用非恒定流输砂方程模拟溃坝沙流的演进过程，建立了一种新的模型来同时模拟下泄水流和尾矿砂的运动，并将其运用在云南某工程中。DL Breach模型在国内土石坝溃坝模拟分析中的应用愈来愈广泛，但鲜有人将其运用到尾矿库溃坝上。周清勇等［5］将Breach模型成功应用于七一水库风险分析；杨德玮等［6］将Breach模型与MIKE模型进行耦合分析，成功应用于郑州常庄水库；李相南等［7］同样对比分析了Breach模型与DBIWHR模型在唐家山堰塞湖溃决中的反演分析，认为2种模型均可较好地反演唐家山溃决洪水的过程。

本研究基于DL Breach模型，利用数值模拟的方法，探索沙流运动的规律，计算尾矿库下游的淹没范围、淹没深度，预测尾矿坝溃决可能引起的灾害，对矿区居民的安全保障和矿区环境的保护有十分重要的指导意义。

1 溃坝数学模型

1.1 模型分类与选择

参考土石坝的溃坝模型，尾矿库溃坝的数学模型也可以分为3类：第一类是基于参数的模型，该模型通过利用已有的尾矿库工程的溃坝数据，根据经验确定尾矿库溃坝的重要参数；第二类基于尾矿库溃坝机理的简化数学模型，该模型利用水动力学相关的推导，建立了一个简化的数学模型并用来模拟溃坝沙流的下泄流量过程线和溃坝沙流演进过程；第三类是基于尾矿库溃坝机理的精细数学模型，该模型利用水动力坝料冲蚀方程来建立溃坝的三维数学模型，可以更加精细地模拟尾矿库的溃坝过程［8］。DL Breach模型属于第二类溃坝模型，该模型基于数学机理，假定溃口形状为梯形，预测溃口发展趋势和溃坝导致的泄流水位过程线，通过输入大坝几何尺寸、大坝宽高值和材料内外属性值等重要参数，来求得含有坝体溃决的各类相关参数以及大坝溃口流量过程曲线。本研究以涩草湖尾矿库为例，基于DL Breach模型获取尾矿库溃坝的下泄沙流流量过程线，模拟下泄沙流的演进过程。

1.2 DL Breach模型

DL Breach模型是一个基于数学机理的预测溃口特征（尺寸、形成时间）和溃决土坝引起的泄流水位过程线的数学模型，该模型原理基于水力学、泥沙运动、土力学、大坝的几何属性和材料属性及水库属性［9］。DL Breach模型将漫顶溃决分为2个阶段：第一阶段为强烈决口或侵蚀阶段，决口水流为超临界，受上游控制。第二阶段为一般决口或进口演变阶段，其水流为亚临界流，受上游或上下游共同控制。在第一阶段，使用堰流方程计算决口流量；在第二阶段，使用Keulegan公式计算决口流量。该模型将漫顶溃坝模式下的溃口横截面近似为梯形，溃口纵截面近似为平顶。该模型已经通过使用50组大坝溃决数据和4组实验案例数据来对比，计算的峰值溃坝流量、溃口宽度等参数与实测数据基本一致。因此将DL Breach模型应用于尾矿库溃坝下泄沙流的计算中，同样具有很重要的参考意义。

1.2.1 溃口流量计算

溃口的流量计算可以分为2个阶段，第一个阶段可以采用宽顶堰公式来计算流量：

式中，H=z s-z b，m，z s是上游水头，m，z b是溃口的底部高程，m；b为溃口的底部宽度，m；m为下游溃口边坡的坡度；系数c1=1.7，c2=1.3；k sm为堰流尾水淹没修正系数；z t是尾水位的高程，m。

溃决的第二阶段采用Keulegan公式计算，该公式是简化的恒定非均匀流能量方程，可以考虑河道的收缩或扩张所引起的水头损失：

式中，L是下泄沙流方向的距离，m；Q是溃口处的流量，m3/s；A是尾矿库溃口的面积，m2；R是水力半径，m；B是溃口处的水面宽度，m；n是曼宁系数；g是重力加速度，m/s2；λe n和λex分别是下泄沙流入口和出口处的水头损失，m；ρ是水的密度，kg/m3；ρa是空气密度，kg/m3；C d是风阻系数；U w in是风速，m/s；θwin是溃口中心线与风向的夹角，（°）。

1.2.2 输砂模型

尾矿库溃坝的下泄沙流中存在着悬移质与河床质泥沙，总的输砂率qb可以表示为

式中，R为河槽水力半径，m；ωs为泥沙流速，m/s；U为水流速，m/s；q*为悬移质输砂率；q b*为河床质泥沙单位时间单宽输砂率，m2/s；d为泥沙粒径，mm；τc=为泥沙的饱和容重，N/m3，γ为泥沙的天然容重n为河床的曼宁糙率系数，τb为河床剪应力，Pa。

1.2.3 溃口扩展过程

尾矿库溃口边坡的稳定程度与尾矿砂的性质、堆积高度等因素密切相关，可以采用滑楔法来模拟溃口边坡的发展（图1）。

2 溃坝沙流运动计算原理

2.1 尾砂流变模型

尾矿库是一种特殊的水工建筑物，当库容较大、堆积坝较高时，库内尾砂具有较高的势能，一旦溃决水流挟带尾矿砂从溃口处冲出。尾矿库下泄的沙流是一种散粒体与水复合材料，组分复杂多样。从尾矿库排放的尾砂造成的下泄沙流性质类似泥石流，溃坝初期由于水量较多，其性质更接近稀性泥石流。由泥石流引起的土壤流动可以被认为是流体和散粒体之间特殊的运动形式。它可以用与流体流动类似的动态方程和连续方程来描述。当沙流流动时具有初始屈服应力，表现出非牛顿流体的特性，根据工程经验，尾砂下泄沙流符合Bingham流动模式，Bingham理想塑性体的切应力与变形速率成线性函数关系，但并未经过原点，存在屈服应力。不同流体的剪切应力与变形速率关系如图2所示。Bingham流体在运动的过程中，密度、黏度以及屈服应力等均不随时间和空间而改变，因此在计算过程中，把下泄沙流假设为Bingham流体，其运动控制方程仅需满足连续性定理、动量守恒和能量守恒，可以简化运动控制方程和数值模拟的过程。同时假设尾矿浆在流动过程中，尾矿砂与水充分混合，两者之间没有相对运动速度，这也符合下泄沙流的一般特征。

Bingham方程是根据Bingham塑性流体的流变曲线写出的，其形式为

式中，τ为剪切应力，Pa；τ0为屈服应力，Pa；ηp为塑性黏度为速度梯度，s-1。

符合Bingham方程的流体称为Bingham流体，塑性流体也称Bingham流体。应力τ0和塑性黏度ηp是公式中2个主要的流变参数。屈服应力是指Bingham流体为实现层流流动所需要克服的动态结构力，而塑性黏度则反映了层流状态下的Bingham流体在处于动平衡状态时内摩擦力的大小。

Bingham流体的表观黏度可表示为

2.2 控制方程

流动控制方程建立在基于Boussinesq和流体静压假定的不可压缩雷诺平均N-S方程的解决方案的基础上。考虑到上述针对尾矿坝溃坝沙流的基本假定，可近似成下述平面深度平均方程组来描述。

连续方程：

动量方程：

式中，t为时间，s；x,y为笛卡尔坐标系坐标，m；η为水面高度，m；d为静水深，m；h为总水深，m；f为科里奥利参数；g为重力加速度，m/s2；ρ为密度，kg/m3；s xx，s xy，s yx和s yy为辐射应力张量的分量，Pa；ρ0为水的密度，kg/m3；S为点源的流量大小，m/s；uˉ，vˉ为x，y方向上的速度分量，m/s。

3 工程实例

3.1 概况

涩草湖尾矿库位于三联村涩草湖下游的钻天道与大干沟内，初期坝为碾压堆石坝，设计总库容5 742.5万m3，尾矿库总坝高为199.0 m，平均堆积边坡1∶5.0。钻天道在初期坝下游约550 m处，原居住4处村宅，现房屋已拆除，人员已搬迁。初期坝下游1 680，1 960，2 560 m处有3座桥，分别定义为1#桥、2#桥、3#桥，这3座桥附近分布有17户村宅。

3.2 溃坝计算参数确定

涩草湖尾矿库初期坝上、下游边坡为1∶2.0，坝高79 m，坝顶宽4 m，坝轴线长约165 m，后期堆积坝平均堆积边坡1∶5.0，尾矿堆高120 m，总坝高199 m，设计总库容5 742.5万m3，有效库容4 708.8万m3。

根据提供的工程地质勘察资料和设计资料，结合现状尾矿库渗流场反演分析成果，在尾矿库剖面上，可以将其划分为5个区域，具体分区和参数如图3和表1所示。

?

3.3 边界条件及初始条件

边界条件包括上游边界和下游边界。上游的边界条件为下泄沙流的流量过程线，由DL Breach模型算出；下游边界条件为开边界，位置在距离初期坝坝顶3 400 m处，其余模型边界作为陆地边界，取法线方向速度分量为零。

尾矿库下游河道中水流的流量很小，水深很浅，因此忽略其影响，认为初始水深为0 m，尾矿库下游的初始条件取天然的地形高程。

3.4 计算工况与流量过程

根据相关国家规范，并参考类似工程经验，涩草湖尾矿库可能发生的溃坝情况为洪水漫顶溃坝工况。漫顶溃坝是大坝溃坝的最主要形式，溃口流量过程线由DL Breach模型计算得出，由图4可知，DL Breach模型得出的溃口最大流量为246.33 m3/s。

3.5 溃坝影响预测

3.5.1 淹没范围预测

本模型范围内沿程沟谷狭长且较平缓，断面基本呈V字形，沟谷地表多树木，两侧山体高。由DL Breach模型计算出的溃口流量过程线可知，整个溃坝流程为3 600 s，选取几个典型时刻：T=600，1 200，1 800，2 800 s，下泄沙流不同时刻的淹没范围见图5。

3.5.2 下游淹没深度

选取下游区域有代表性的几个敏感点，计算下泄沙流在这些地方的淹没深度，并绘制了如图6的溃坝沙流淹没深度随时间变化的曲线。

从计算结果来看，尾砂会积聚在下游钻天道内，下游空地、溪沟、灌木等会被淹没；钻天道附近居民已经搬迁，一直到初期坝下游1.5 km以内都没有重要建筑物，该范围内沙流淤积面积大，扩散范围广。截渗坝和拦沙坝上游（钻天道沟谷内）会有大量的泥石流淤积，而其下游淹没深度小，树木、植被也对泥石流的演进具有明显的阻碍作用，危害较小。

在整个沙流运动的过程中，尾砂在下游的堆积体的深度也随着距离的增大而减小，距离初期坝越远，淹没深度越小。溃坝后水深总体上先增大后减小，大致可分为3个阶段：增长阶段、稳定阶段、衰减阶段，这也反映了沙流运动的过程。尾砂行进过程中，堆积深度由0 m开始逐渐增加，增大至最大深度后，继续向下游流动，此后堆积深度逐渐减小，并趋于平缓，这符合沙流演进的一般规律。钻天道原居民点处的淹没深度最大可达到2.49 m，根据计算结果可以发现，距离初期坝越近，尾矿砂的堆积深度越深，而不会随着下泄沙流的运动往下游推进从而堆积得更高，这也与尾砂这一特殊材料的性质有关。对于漫顶溃坝，洪水漫过坝顶冲刷坝体表面，这种情况下尾矿坝溃口缓慢发展，按照DL Breach模型计算出的流量过程曲线，下泄沙流的流量先增大再减小，数值模拟计算结果表明，该模型不仅适用于普通的土石坝，也适用于尾矿坝这一特殊的水工建筑物。利用DL Breach模型来预测尾矿库溃坝后下泄沙流的演进过程是合理可靠的。

4 结语

（1）结合具体的工程案例，利用数值模拟的方法，计算某尾矿库溃坝后下泄沙流的演进过程，得到了下游淹没范围、各敏感点淹没深度以及流量、溃坝沙流抵达时间等关键数据，符合尾矿沙流运动的一般规律。

（2）实例计算结果：溃坝下泄沙流的堆积深度总体上先增大后减小，增大到最大深度后，沙流继续向下游流动，此后堆积深度慢慢减小；下泄沙流在钻天道沟谷内淹没深度较大，淤积较多，最大堆积深度达2.48 m，但此处没有民宅，因此不涉及下游群众的安问题；随着与初期坝的距离逐渐变大，下游的3#桥处的最大淹没深度不超过1 m，沙流对下游主要建筑物影响较小，整体安全风险可控，不会威胁到村民的生命安全。

（3）根据利用DL Breach模型进行计算的尾矿库溃坝数值模拟，可以预测尾矿坝溃决可能引起的灾害，对尾矿库的安全运行具有指导意义。