漫顶溃坝尾砂流演进模拟及避险转移方案分析

彭浩，朱光源，张邵祥，蒋水华，蔡木良，黄劲松

(1.南昌大学工程建设学院，江西 南昌 330031;2.国网江西省电力有限公司电力科学研究院，江西 南昌 310058)

尾矿库是矿山企业安全生产过程中的重要基础设施，但同时也是矿山企业最大的危险源。一旦尾矿库发生溃坝事故，高势能下泄尾砂流将对下游人民生命财产安全和生态环境造成巨大影响[1]。洪水漫顶溃坝是尾矿库较为常见的溃坝模式。据统计，2001—2013年我国发生尾矿库溃坝事故共计65起，其中由洪水漫顶引起的溃坝事故占溃坝事故总数的36.9%[2]。因此，研究尾矿库漫顶溃坝的发展过程，分析下泄尾砂流的演进规律，科学预测尾矿库失事后果和制定科学合理的避险转移方案，具有非常重要的工程应用价值。

随着计算机技术的日益成熟，数值模拟方法因成本低、耗时少等优势越来越多地被应用于解决尾矿库溃坝过程仿真难题。郑欣[3]利用Fluent软件模拟了金山尾矿库溃坝尾砂演进过程，得到了下泄砂流的影响范围。阮德修等[4]结合Flo2D软件和3DMine数字矿山软件实现了湖南某尾矿库溃坝过程数值模拟。Mahdi等[5]发展了考虑下泄尾砂流非牛顿性质和黏塑性关系的尾矿库溃坝数值模拟方法，实现了对加拿大Alberta尾矿库溃坝尾砂演进过程的合理预测。

尽管这些研究工作为解决尾矿库溃坝模拟问题提供了有效途径，但是目前的方法通常将溃坝尾砂流视作物理性质恒定不变的非牛顿流体，并没有考虑尾矿库在发生漫顶溃坝过程中水对尾砂的冲刷作用以及下泄尾砂流含沙量的时变特征。为此，本文依托江西省永平铜矿燕仓尾矿库5号副坝，利用Civil-3D软件建立尾矿库及下游区域三维几何模型，采用RNGk-ε湍流模型和泥沙冲刷模型构建尾矿库漫顶溃坝水动力学数值模型，对尾矿库漫顶溃坝过程进行精细化模拟，预测溃坝淹没范围及下游村庄的淹没水深、泥沙淤积厚度、流速等关键信息。在此基础上，制定避险转移方案，为尾矿库溃坝风险管理与决策提供理论依据。

1 尾矿库漫顶溃坝数值模拟

本文在Flow-3D软件平台上，采用RNGk-ε湍流模型和泥沙冲刷模型构建尾矿库漫顶溃坝水动力学数值模型[6-8]。

1.1 水流控制方程

在笛卡尔坐标系中，溃坝水流可视做不可压缩流体，采用连续方程和Navier-Stokes动量方程描述流体运动[9]，其中连续方程为

(1)

Navier-Stokes动量方程[9]为

(2)

式中：u，v，w分别为流体在x，y，z3个方向的流速分量；Ax，Ay，Az分别为在计算单元中各截面流体所占面积与截面总面积的比值；VF为计算单元中的流体体积分数；p为压强；ρ为流体密度；gz，gy，gz分别为3个方向的重力加速度；fx，fy，fz分别为3个方向的黏滞力加速度。

1.2 RNG k-ε湍流模型

RNGk-ε湍流模型由Yakhot等[10]于1986年提出。该模型基于严格的统计技术，在传统k-ε湍流模型的基础上改进而来，能够有效地解决低雷诺数和具有强剪切区的湍流流动问题，计算精度高，适用范围广。其控制方程包括紊动能k方程和耗散率ε方程如下：

(3)

(4)

式中：k为紊动能；ε为紊动能耗散率；μ为流体动力黏滞系数；μt为流体湍动黏度；Gk为紊动能k的产生项；αk，αε，C1ε，C2ε为无量纲系数，根据文献[6]，这些参数取值为αk=αε=1.39，C1ε=1.42，C2ε=1.68。

1.3 泥沙冲刷模型

泥沙冲刷模型[9,11]通过以下4个方面的计算来实现对泥沙冲刷、淤积和运移的精细化模拟：(1)泥沙起动计算；(2)泥沙沉降计算；(3)悬浮泥沙随流体运移计算；(4)推移质泥沙运移计算。

泥沙起动指的是水流的拖拽使得沉积在河床上的泥沙被水流带走，成为悬浮泥沙；沉降指的是水流中的悬浮泥沙在重力作用下沉积。泥沙的挟带和沉降需依据水流速度而定，当水流速度大于泥沙起动速度时，泥沙将悬浮在水流中运动；当水流速度小于泥沙沉降速度时，泥沙将逐渐沉积或沿河床推移运动。

1.4 边界条件

数值模型所涉及的边界条件主要包括入流边界、壁面边界、出流边界和对称边界四大类。

(1) 入流边界。主要包括压力边界、流速边界、流量边界。由于水深和压力具有相关性，压力边界主要约束了边界处的压力和水深；流速边界和流量边界主要通过约束边界处的流速和流量来控制入流条件。

(2) 壁面边界。Flow-3D中常用固壁边界来模拟壁面约束，在边界处应用无滑移条件以及垂直于边界的零速度条件来模拟在流体运动问题中的固体边壁约束。

(3) 出流边界。出流边界包括连续边界和自由出流边界。连续边界流速梯度和压力梯度为0，即流体流经边界时各个物理参数恒定不变，该边界对计算结果干扰较大。自由出流边界假设边界符合Sommerfeld Radiation条件，流体在流经该边界时的扰动和波动是平顺的，较为符合实际情况。

(4) 对称边界。对称边界条件在边界处应用零梯度条件以及垂直于边界的零速度条件。对于具有对称性的计算模型，可运用对称边界减小计算量。另外，若在建模时划分了多个网格区域，则网格区域的边界重合处需要设置为对称边界条件以进行网格区域间的数据交换，实现网格区域的连接。

2 模型验证

为了验证所建立的尾矿库漫顶溃坝水动力学数值模型的有效性，首先选取文献[12]中漫顶溃坝水槽试验为典型算例进行模型验证。

2.1 计算模型

该试验在一长17 m，宽1.2 m的水槽中进行，如图1所示[12]。槽顶距地面2.4 m，槽深1.5 m，坡度为0。水槽上端为一个5.0 m×5.0 m蓄水池，池深2.4 m。试验中坝高设置为1.2 m，顶宽0.2 m，上游坝坡为1:1.25，下游坝坡为1:1.2。为保证下游坡面冲刷均匀，试验开始前在坝顶布置一块高0.1 m的挡板，之后供水水箱持续向蓄水池内供水。当坝前水位与挡板上缘平齐时，迅速提起挡板，溃坝开始。溃坝试验过程中供水水箱的供水流量恒定为20 L·s-1。

图1 水槽模型试验平面布置图[12]

根据漫顶溃坝水槽试验建立三维数值模型并对其进行网格划分，数值模型尺寸与室内试验装置尺寸相同，如图2所示。整个坝体由无黏性砂组成，中值粒径为3.77 mm，密度为2 650 kg·m-3，临界希尔兹系数采用Soulsby-Whitehouse公式计算[11]。为了模拟室内溃坝试验初始条件，上游水深设为1.3 m(高出坝顶0.1 m)。

图2 水槽模型网格剖分

采用正六面体结构化网格对模型进行网格剖分。根据该模型几何特征，设置了两个网格区域，这两个网格区域的网格尺寸均为0.1 m×0.1 m×0.1m。网格区域2的上游侧设置为流量边界(Q)，用于模拟试验过程中的上游来水；网格区域1和区域2重合边界处定义为对称边界(S)，实现网格区域连接；模型上方设置为压力边界(P)，模拟大气压条件，压力设置为101.325 kPa；网格区域1的下游侧设置为自由出流边界(O)；其余部分定义为无滑移固壁边界(W)，如表1所示。为与文献[12]中物理试验时间保持一致，数值模拟时间也设为250 s。

表1 水槽模型边界条件

2.2 溃口流量分析

图3给出了水槽模型试验和数值模拟获得的溃口流量qV随时间t的变化关系曲线。可知，溃坝发生以后，溃口流量短时间内迅速增大，到达峰值后开始下降，并逐渐趋近于上游来流流量，这显然符合漫顶溃坝溃口流量发展的一般规律[13]。本文通过数值模拟得到的各个时刻溃口流量变化规律与实测值基本吻合，表明本文所建立的尾矿库漫顶溃坝水动力学数值模型有效、可行。

t/s

3 工程应用

3.1 工程概况

本文以江西省永平铜矿燕仓尾矿库5号副坝[7]为例，进行漫顶溃坝尾矿演进过程模拟。该副坝初期坝为均质黏土坝，轴线底部高程为129.5 m，初期坝顶高程为135.0 m，设计外坡比为1:2，设计内坡比为1:2.5，坝顶宽为3 m；堆积坝采用上游法进行分散管放矿堆筑，最终堆积坝顶高程为150 m。正常蓄水位为146.5 m，设计外坡比为1:5，设计干滩坡度为1:167。尾矿坝下游约2 km范围内有两个村庄以及多条公路和大面积农田，如图4所示。

图4 尾矿库下游卫星图

3.2 溃坝模型

根据尾矿库实际尺寸以及下游地形图，基于地形等高线和高程点信息，采用Civil-3D软件建立尾矿库及下游2 km范围内地形和房屋的实体几何模型。为提高计算效率，本数值模拟洪水重现期取规范上限1 000年，根据《江西省暴雨洪水查算手册》推荐方法，计算出千年一遇洪水总量为6.15×105m3，与库区内能容纳水量(6.19×105m3)相差不大。故可在模拟开始时预先在库区内充满水体，保持库水位与坝顶高程相等。尾矿库遭遇洪水漫顶时，坝顶最薄弱部位一般先形成主溃口，但是数值模拟无法模拟坝顶最薄弱部位。故本文按照文献[8]的做法，在尾矿坝中部位置预先设置一个宽30 m和深3 m的矩形断面初始溃口，将漫顶过程逐渐导向最危险工况发展。

为了分析溃坝尾砂流演进过程及其对下游村庄的影响，在尾矿库下游设置了若干监测点以获取下游关键区域的淹没信息。如图5所示，在新岩前村内设置监测点1，其高程为107.12 m；因塘棣源村面积较大，房屋分布较广，设置了2个监测点，分别是位于村口的监测点2和位于村内的监测点3，其高程分别为90，92.69 m。

图5 尾矿库及下游地形三维模型

为了保证数值模拟精度和效率，将研究区域划分为两个网格区域，如图6所示。其计算网格尺寸均为3 m×3 m×3 m，总有效网格数约为450万。模型上方设置为压力边界(P)，模拟大气压条件，压力值为101.325 kPa；网格区域1和区域2重合边界处定义为对称边界(S)，实现网格区域连接；模型下游出口处设置为自由出流边界(O)；其余部分有山体阻挡或者为地面，定义为无滑移固壁边界(W)。具体边界条件设置见表2。此外，泥沙模型参数设置如下：平均粒径为0.112 mm，密度为2 940 kg·m-3，临界希尔兹系数为0.05；挟带系数为0.01，推移质系数为8，休止角为34.7°。经多次试算确定漫顶溃坝数值模拟时间为3 200 s。

表2 尾矿库三维模型边界条件设定

图6 三维数值模型网格剖分

3.3 溃坝尾砂流演进过程及流速

为分析溃坝尾砂流演进过程及流速变化规律，图7给出了不同时刻的溃坝尾砂流流速变化云图。

漫顶发生后100 s，见图7(a)，溃口水流通过初始溃口向下游演进并逐渐冲刷坝体，形成的挟砂水流演进至尾矿库下游约200 m处，溃口流速约为7 m·s-1。溃坝发生后500 s，见图7(b)，溃口继续扩宽并向上游发展，溃坝尾砂流已经演进至尾矿库下游约1 000 m处，其流速约为5 m·s-1，即将对塘棣源村造成影响。相比之下，新岩前村地势较高，只有少量溃坝尾砂流越过山地，进入新岩前村。溃坝发生后1 000 s，见图7(c)，在水流的冲刷作用下溃口持续向上游发展，此时溃坝尾砂流一部分已经进入新岩前村内，村内尾砂流流速在0～3 m·s-1范围内变化，但是由于新岩前村内建筑物地势较高，溃坝尾砂流并没有对其产生影响。另一方面，溃坝尾砂流已经演进至塘棣源村，靠近房屋处尾砂流流速约为5 m·s-1。溃坝发生后1 500 s，见图7(d)，由于库水较少，库内尾砂的冲刷和侵蚀作用较小，因而下泄尾砂流也较少，但是沟谷中和新岩前村内的尾砂流仍然持续向下游演进，导致下游地势较低的塘棣源村的尾砂流流速并未明显减弱，维持在4 m·s-1左右。溃坝发生后2 000 s，见图7(e)，库水已经下泄完毕，下游区域仅剩塘棣源村前还有溃坝尾砂流以较低的流速演进。溃坝发生后3 200 s，见图7(f)，尾矿库下游区域已无尾砂流演进，大量泥沙已经沿程淤积。

(a) 100 s

3.4 淹没深度

图8给出了3个监测点的淹没水深h随时间t的变化关系曲线。可知，溃坝尾砂流在漫顶后约750 s到达新岩前村口监测点1处，约1 190 s达到峰值，峰值水深为3.42 m；溃坝尾砂流在漫顶后约570 s到达位于塘棣源村口的监测点2处，约1 330 s达到峰值，峰值水深为4.56 m；溃坝尾砂流在漫顶后约860 s到达位于塘棣源村内的监测点3处，约970 s达到峰值，峰值水深为1.77 m。由于溃坝尾砂流在沿下游沟谷演进过程中遭遇了山体地形阻拦被分为了两股(见图7)，其中一股在分岔口跨过地形直接演进至塘棣源村，导致布置于该村的监测点2和监测点3处的水位上升，随着该股尾砂流持续向下游演进，监测点2和监测点3的水位逐渐达到局部峰值后开始下降；另一股尾砂流越过地形演进至新岩前村，随后继续演进至塘棣源村，到达塘棣源村的时间较晚，这便造成了塘棣源村监测点2的水位第2次上升至局部峰值而后再次下降。此外，监测点3由于地势较高，当第1股尾砂流消散后，第2股尾砂流还未到达，故该监测点水深在溃坝后1 000 s骤减至0。

t/s

结合布置于村庄内3个监测点的高程和峰值水深信息可知，新岩前村内底部高程低于110.5 m的房屋建筑物将会受到溃坝尾砂流影响；而塘棣源村内底部高程低于94.5 m的房屋建筑物将会受到溃坝尾砂流影响。

3.5 最大淹没范围及尾砂淤积情况

由上述溃坝尾砂流演进分析可知，新岩前村和塘棣源村出现最大淹没的时刻分别是在漫顶发生后1 160，1 280 s，见图9。溃坝尾砂流进入新岩前村后主要聚集于村内低洼处，但村内房屋总体地势较高，房屋最低高程为112.12 m，高于110.5 m，故该村房屋均不受溃坝尾砂流的影响。塘棣源村内部分房屋低于94.5 m，存在被尾砂流冲击和淹没的可能性，故该村受到溃坝的影响较大。

(a) 新岩前村

下泄尾砂总量约为4.82×105m3，大量淤积于尾矿库下游沟谷和两个村庄中，掩埋或损坏村内的房屋、道路等设施。图10为溃坝结束后下游尾砂淤积云图。可以清楚看出，溃坝后的泥沙淤积分布情况。尾砂在新岩前村淤积较轻，并且由于该村房屋建筑物地势较高，淤积尾砂并未对其造成影响。大量尾砂淤积于塘棣源村以及塘棣源村前农田处，地势低洼处尾砂淤积厚度达到10 m左右，部分地势较低的房屋建筑物将受到较为严重的淤积尾砂的影响，甚至存在被尾砂完全掩埋的可能性。

图10 尾砂淤积云图

4 避险转移方案

由上述尾砂流演进分析结果可知，永平铜矿燕仓尾矿库发生漫顶溃坝后，新岩前村房屋均不受溃坝尾砂流的影响，而塘棣源村内高程低于94.5 m的房屋建筑物将受到尾砂流冲击或掩埋，故需要提前制定受灾人员的疏散、撤离和避险转移方案。转移方案的制定步骤主要如下：(1)根据区域地形、交通以及事故影响范围等信息，选取合适的应急避险场所中心点；(2)根据需转移人数确定应急避险场所，可以按步骤(1)确定的中心点为圆心规划一个平面圆形区域作为应急避险场所[14]；(3)计算受灾人员到达应急避险场所的时间。

应急避险场所中心点的选取一方面应遵循安全性的原则，即此地不会被溃坝尾砂流冲击或掩埋。在能够保证本身安全的前提下，尽可能考虑能够容纳较多的人口，并远离灾害源，选择地势相对较高的地方；另一方面应考虑可达性，即受灾人员能否快速到达事先安置好的应急避险场所。在遵循上述两个原则的前提下，还应确保所选地点尽可能靠近道路，便于人员迅速撤离至应急避险场所。

经统计，塘棣源村共计15户61人将受到尾砂流威胁，根据《避险转移图编制技术要求(试行)》[15]规定，应急避险场所人均面积不小于3 m2，由此可确定应急避险场所面积应不小于183 m2。故应急避险圆形场所半径拟定为8 m，总面积约为200 m2，可以满足受灾人员安置要求。

塘棣源村内道路纵横，撤离最短路线总路程与受灾房屋到应急避险场所中心点的直线距离相差不大，故选取受灾房屋到应急避险场所中心点的直线作为转移路线，以此来估算人员撤离时间。人员撤离速度取4 m·s-1[14]。经计算，各受灾居民的转移距离在181.26～289.32 m范围内，转移时间在45.32～72.33 s内。由3.4节可知，溃坝尾砂流在漫顶发生后570 s到达塘棣源村。这表明，按照上述方法制定的避险转移方案(图11)，可在确保预警工作有效的前提下，塘棣源村可能受灾的居民有充足的时间完成疏散、撤离和避险转移。

图11 塘棣源村居民避险转移方案示意图

5 结论

(1) 从尾砂流流速、流深、淹没范围、泥沙淤积厚度等方面调查了尾砂流对下游居民、村庄公路和农田等造成的影响。燕仓尾矿库漫顶溃坝下泄尾砂共计4.82×105m3，大部分尾砂沿程淤积，其中一部分会尾砂淹没村庄公路和农田。其中，新岩前村淹没高度为110.5 m，但因新岩前村房屋建筑物高程均高于110.5 m，故该村基本没有受到溃坝尾砂流影响；相比之下，塘棣源村淹没高度为94.5 m，因该村部分房屋高程低于94.5 m，故受到了溃坝尾砂流冲击或淤埋等严重影响。

(2)在溃坝尾砂流演进模拟的基础上制定了塘棣源村居民的避险转移方案。溃坝尾砂流到达塘棣源村时间为570 s，该村居民转移时间为45.32～72.33 s。在保证预警工作有效的前提下，该村居民有充足的时间按照制定的方案完成避险转移。