0 引言

尾矿库是具有高势能的人造泥石流危险源，尾矿库溃坝极易造成重特大事故。尾矿库溃坝的发生、发展和溃决程度等受诸多因素影响（如溃坝原因，坝体材料和库容等），模拟难度非常大。尽管如此，在过去 40 多年中，溃坝模拟还是取得了相当大的进展。大坝安全问题的关注推动了对溃坝水力学的研究，开发了许多数学模型。溃口模拟的主要有 DAMBRK 模型、BEED 模型、BREACH 模型、 Cristofano 模型、HW 模型、Lou 模型和 Nogueira 模型等。

目前，学者们通过数值模拟的方式对尾矿库溃坝进行了大量的模拟研究，卢春维（2023）基于流体颗粒理论的颗粒流模型开展尾矿库溃坝泥石流研究工作，得到了泥石流运动规律，为下游科学避灾提供指导；朱远乐等（2023）利用PFC3D模拟研究了尾矿下泄可能淹没的范围及影像区域；周薛淼 （2023）基于构建的尾矿库漫顶溃坝动态模型，模拟研究了堤防溃口处的展宽、崩塌情况、刷深和溃口周围冲刷坑发育变化过程和水流规模等信息，为溃坝堵口提供了技术支撑；邓金鹏等（2023）基于尾矿库溃坝数值模拟建立了溃坝尾砂泥石流对下游高速公路桥梁的损伤和滑动破坏评估模型，实现对高速公路损伤和滑动破坏进行定量评估；孙黎明等 （2023）以 2008 年陕西襄汾陶寺乡 980 沟尾矿库溃坝事件为背景，建立了襄汾尾矿库溃坝前后的高精度三维模型，为少资料地区尾矿库三维模型快速重建提供参考；黄泰宇等（2024）利用 Flow-3D 数值模拟软件模拟研究了溃水淹没高度和断面流量随溃坝时间的变化规律；邱子源等（2024）通过建立三维溃坝有限元模型并提出了刻画尾矿库溃坝浆体流动规律的数学模型，定量分析评价了尾矿库溃坝后对下游的影响；闫崭等（2024）构建三维数值模型，对尾矿库多级拦截坝滞留尾砂泥石流效果和对下游敏感目标防护作用开展深入研究。

综上所述，本文利用 MIKE 21 平面二维水动力模块模拟研究甘肃某新建尾矿库在极端条件下尾矿库溃坝对尾矿库周围 1 km 范围内各类设施的影响规律，分析了尾矿库溃坝演进规律及发生的前提条件，帮助企业更全面地了解尾矿库对周边环境和人类活动的潜在风险，为相关决策提供科学的依据。

1 模型的建立及前期准备

该尾矿库东北 1. 001 km 外有交建建材园，580 m 处有二冶临时搅拌站（待拆除），东侧及东南侧 1. 0 km外有嘉策铁路，尤其是东南控制区范围内存在嘉峪关交投公司采矿区。尾矿库所在地形地势为西高东低，东方向的尾矿库总坝高较大，特殊工况下的东北侧截面的坝体稳定性最弱，低于相邻截面的坝体稳定性，故该截面易发生边坡局部失稳。经现场核实尾矿库周边环境，新建尾矿库用地的库区东南角直接压覆嘉峪关市交建集团建材园的南 2#采砂场，设置禁采线 50 m，总共压覆矿产面积约 99277. 04 m²，一旦溃坝将对砂场带来影响。因此，模拟尾矿库因局部失稳而导致溃坝现象是本文研究重点内容之一。溃口选择东北侧断面所在位置。

1.1 二维模型原理及构建

（1）模型计算原理

水动力模型，描述水流受力与运动相互关系的数学模型。依据流体力学基本方程，建立数学模型，对流动水的动力过程进行数值模拟。水动力数学模型采用的二维非恒定流作为基本控制方程，方程基本假设水体为不可压缩水体且压强分布沿水深方向为静压分布，二维非恒定流方程由水流连续方程和水流运动方程组成，如公式1~6所示：

x方向水流运动方程：

y方向水流运动方程：

式（1~3）中：u、v 为垂线流速分量在 x、y 方向的分量（m/s）；h—水深（m）；ζ—水位（m）；f—科氏力系数（f=2ωsinφ）；v_e—有效粘性系数（Pa·s）；v_t —紊动粘性系数（m² /s）；τ_bxτ_by—底部切应力在 x、y 方向分量 （Pa），可表示为：

式（4~5）中：c_f为底部摩擦系数：$c_f=n^{2}g/H^{\frac{1}{3}}$；n为河床糙率；τ_xτ_x—表面风应力在x，y方向分量；

式（6）中：k_s为系数，计算过程中不计风应力，故 τ_sx、τ_sy均为0。

（2）方程的离散与求解

将二维计算区域划分为非结构化或者不规则形状的小单元（即控制体积），每个单元进行水量及动量平衡计算，进而得到每个计算单位的水位及流量。用有限体积法计算各网格，网格中心为计算单位水位h，网格边界为通道，在通道中心点计算水流 Q，网格边界Q-h交错布置，如图1所示。

控制体积法控制方程为：

式（7）中：v为计算网格；n为单元边界 A的外法向单位向量；F（U）·n为方向数值通量。

将控制体取均值，转化为：

式（8）中：ΔV为各单元的体积（m³）；F_n，j为法向通量。

其中时间的离散有两种：分别是低阶显式Euler 方程和二阶Runge-Kutta方程，公式如下：

1.2 数值模型建立

（1）二维网格文件

二维模型中常用的网格类型有结构化网格和非结构网格两种，由于非结构划分法不受网格节点约束，能够合理的分布网格，适用于复杂地形区域网格划分，因此本研究中根据尾矿库实际情况，将研究区域划分为不规则非结构网格。在网格划分过程中，网格的大小是随地形地貌变化而变化的，对于高程变化小、边界规则的地形区域，网格划分面积较大；对于高程变化大，边界复杂的地形区域，网格划分面积相应减小或者进行局部加密网格，从而提高运算效率和运算精度。

（2）网格与地形图

尾矿库溃坝分析采用无结构三角形对尾矿库影响地区进行剖分，无结构三角形具有复杂区域适应性好、局部加密灵活和便于自适应的优点，能很好地模拟自然边界及复杂的水下地形，提高边界模拟精度，本次二维模型计算范围主要为位于尾矿库，计算面积7.43 km²，网格剖分采用不规则三角形网格，共剖得网格 10 万个，本文划分网格的大小对于尾矿库溃坝模拟研究是满足精度要求的。根据尾矿库坝体影响区地地形数据进行网格剖分，如图2所示。

1.3 模型参数设置与求解

运用MIKE 21平面二维水动力模块模拟尾矿库最终堆积标高时达到漫顶情况下出现溃口且溃口逐渐扩大而失稳的工况，生成二维溃坝演进模型。溃决后下泄量考虑最不利排洪构筑物失效的情况。

1.3.1 边界条件

在溃坝洪水模拟过程中，边界条件也是影响模拟结果精度的重要因素，溃坝洪水数值模拟结果的准确性很大程度上取决于该模型的边界条件。在二维模型计算中有两种边界条件，一种为开边界，一种为闭合边界，模型中闭合边界为陆地边界，开边界一般为水边界。

A. 闭合边界设定：在模型中闭合边界多为陆地边界，边界上流速变量为0；

B. 边界设定：采用一二维耦合计算溃坝洪水模拟计算，二维区域计算的上下游开边界为二维模型的边界条件。

1.3.2 模型参数选取

在MIKE 2l模型中，需设置的有：计算时间和计算步长、糙率、风场、涡粘参数等。

A. 计算时间和计算步长：计算时间及步长受模型网格大小和CFL数影响，本次计算步长设置为30 s，满足模型稳定要求；

B. 涡黏参数：采用 smagorinsky 公式计算，取值 0.28 m² /s；

C. 风场：由于风场对溃坝洪水模拟影响较小，本次计算忽略风场设置；

D. 底部阻力：本研究用曼宁系数的变化作为剪切力的函数，模拟出尾矿浆不同于水的流动特征。模型中剪切力和曼宁系数的关系如图3所示。

图3中，τ_c为尾矿浆开始流动的临界剪切力；τ_cu 为尾矿浆从屈服假塑性体变化为宾汉塑性体的临界剪切力；Ml 为曼宁系数的最小值，表征尾矿浆开始流动时受到的底摩擦，取 3.2；Mu 为曼宁系数的最大值，表征从屈服假塑性体变化为宾汉塑性体时受到的底摩擦，取32。

1.3.3 模拟工况选取

《生产过程危险和有害因素分类与代码》（GB/T13861—2009）明确了该尾矿库存在的主要危险因素有溃坝、洪水漫顶、渗流破坏。库区周边、库底均为第四系冲洪积层、冲洪积卵石层，厚度大，渗透系数为 35.20~36.11 m/d，具有强透水现象，具渗漏可能，但在尾矿库的建设阶段，库内已做足够防渗设施，发生渗透破坏的可能性极小。库区周边卵石呈中密—密实状态，地形与岩土体结构面组合属稳定结构，故库岸稳定性良好。

勘察区及其外围在大地构造上属中朝地台北祁连褶皱带与阿拉善台隆之间的走廊坳陷带内，从北向南依次有金塔中新生界断陷、鸳鸯池拱断束、走廊拗陷带和北祁连褶皱带，构成相间排列的隆起与凹陷，本区地质结构稳定，褶皱及断裂构造均不发育。拟建场地区域稳定性良好，无发生地质灾害的可能。

尾矿库从地形地质条件、周边环境、洪水计算及调洪演算、抗滑稳定分析方面均满足尾矿库建设要求。因此，假设尾矿库在极端条件下因排洪系统失效导致的局部冲刷溃坝和局部边坡失稳引起溃坝的可能性相对较大，本次模拟针对这两种溃坝形式的影响特征进行分析。

1.3.4 溃口宽度

溃口是尾矿库溃坝时形成的缺口，影响溃口形态包括大坝的体型及尺寸、坝体材料粒径及级配、施工方法及施工质量以及外来因素影响等。

本计算报告参考水库溃坝的溃口形状，考虑最严重后果，将新建尾矿库溃口形状设定为倒梯形。溃口位置选取对下游最不利位置尾矿库东北坝。溃口宽度参照黄河水利委员会科学研究院的经验公式，经过修订，尾矿库溃坝决口平均宽度为：

式（10）中：b 为溃口宽度（m）；W 为水库总库容 （m³）；B 为主坝长度（m）；K 为经验系数（粘土取 0.65，壤土取1.30）。W=6775.7 万m³，B=1216.3 m， H=1600 m，按照式（10）计算可得b=139.3 m。

2 尾矿库漫顶溃坝影响范围数值模拟

2.1 地形展示

根据初期坝以及堆积坝的布置方案和现状，对尾矿库的地形进行处理，得到平面图形如下图4 所示。

2.2 洪水漫顶影响范围计算结果

尾矿库东北侧坝高最大，布置的排洪管道最多，在极端条件下发生洪水漫顶的可能性也最大，因此漫顶溃口选择在坝体东北侧。漫顶溃坝工况状态下，首先在坝顶形成溃口，进而库内洪水携带泥砂冲向下游，本次分析 100 min内溃坝范围情况，溃坝范围见图5。

本次模拟考虑极端最不利情况下的溃坝，即尾矿库达到漫顶情况下出现溃口且溃口逐渐扩大的失稳破坏。该工况下模拟结果每间隔 20 min 淹没范围演进图如图5所示。

由图5可知：（1）溃坝初始阶段尾矿库内的水沙流由东北侧溃口流出坝外。20 min时尾矿砂到达下游高地并开始向东北方向蔓延，32 min 时尾矿砂包围下方高地，60 min时尾矿砂到达砂场东北角，开始向正北方向移动。92 min 时，尾矿砂在北方向和东方向的流动缓慢，开始回转向西方流动。100 min时模拟结束，淹没区域内大部分尾矿砂高程在 1536~1540 m。（2）尾矿库溃坝初期，尾矿砂大部分向东方以及北方移动，尾矿砂只向南方前进了一小部分。 （3）尾矿到达最远点至溃口位置约为 960 m，溃坝对东北侧砂场影响大于东南侧砂场。直至结束尾矿砂仍未到达建材园并与其保持一定距离，未对建材园造成影响。

2.3 洪水漫顶淹没高程计算结果

对该工况下模拟结果每间隔 20 min 淹没高程演进图如图所示（图6）：溃坝后尾矿砂在4 min时溃口处的尾矿高度突破 2 m、12 min 时东南侧尾矿高度突破 3 m、尾矿沿东南侧尾矿坝轴向方向高度最高、36 min 时尾矿砂东北角高度突破 4 m、72 min 时由于地形原因最大淹没高度出现在尾矿砂东北角、 84 min 时东北侧最大高度突破 6 m、100 min 时模拟停止，此时尾矿最大高度为6.5 m、溃坝后尾矿在前 44 min 内对其下游地形高程改变较大，44 min 之后尾矿砂虽依旧在前进，影响范围变大，但淹没高度增加不明显。

新建尾矿库周围可能受尾矿库溃坝影响设施只有位于东北侧交建建材园以及二冶临时搅拌站 （待拆除），故在考虑尾矿库溃坝影响时只考虑建材园。在尾矿库溃坝后尾矿流动区域内选取8个断面考察尾矿库溃坝对周围影响，断面位置见图7。

2.4 各断面泥砂流量变化分析

各断面距新建尾矿库东侧交点距离分别为： DM1-溃口处，DM2-150 m，DM3-310 m，DM4-500 m， DM5-676 m，DM6-820 m，DM7-960 m 断面处的溃坝数值模型计算流量（m³ /s）过程线如图8所示。

a—0 min；b—20 min；c—40 min；d—60 min；e—80 min；f—100 min

根据图8 所示，断面 1~7 处的洪峰流量分别为 780 m³ /s、679 m³ /s、420 m³ /s、206 m³ /s、182 m³ /s、 134 m³ /s、122 m³ /s，通过上图和上述数据可以看出：离尾矿库越远的断面水砂混合体到达时间和洪峰到达时间越长，DM1~DM7 各断面处的水砂混合体到达时间分别为：1.5 min、12 min、15 min、34 min、 48 min、56 min 和 73 min；DM1~DM4 处均呈现出先迅速增加到洪峰流量后缓慢下降，说明这几处断面地势差较大，沿程阻力较小；DM4 呈现缓慢增加至一定值后缓慢降低，说明此地地势落差小流速较慢，DM6、DM7 先流量迅速增加到一定值然后迅速下降直至 0 流量；从各断面流量曲线可以看出除 DM3还依然具有一定流量外，其余断面基本趋近于 0，说明尾矿砂的流动即将停止。

2.5 下游各断面泥沙淹没高程变化分析

各断面处的溃坝数值模型计算淹没高程线见图9 所示。由图可知：各断面最大淹没高程分别为 1553.3 m、1547.8 m、1544.5 m、1541.4 m、1539.4 m、 1539.2 m、1539.2 m；各断面处溃坝后水砂混合体淹没高程随溃坝后历时变化曲线形状可以看出：断面1~4 处由于距离溃坝处较近，所以淹没高程迅速增加后趋于一稳定值；而断面 5、6、7 处于持续上涨阶段，且上涨速度较快，若超出其最大容量后极有可能淹没其周边其他设施。将淹没范围回溯至地形图（图10）。

a—0 min；b—20 min；c—40 min；d—60 min；e—80 min；f—100 min

3 尾矿库边坡局部失稳溃坝影响范围数值模拟

3.1 边坡局部失稳影响范围计算结果

本次模拟为边坡局部失稳造成的溃坝，对该工况下模拟结果每间隔8 min淹没范围演进图如图11所示。由图可知：尾矿砂于尾矿库东北侧溃口涌出，并迅速向东移动，8 min 时尾矿砂包围高地，12 min 时尾矿砂向北方蔓延；16 min 时向西蔓延， 24 min 时尾矿砂接近建材园；40 min 时尾矿流动结束，尾矿到达最远处距离溃口 1500 m；边坡局部失稳溃坝无论是淹没范围还是流速都远远大于漫顶溃坝。对东北侧砂场有较大影响；尾矿砂并未到达、接触建材园，边坡局部失稳溃坝不会对建材园造成影响。

a—DM1；b—DM2；c—DM3；d—DM4；e—DM5；f—DM6；g—DM7

3.2 边坡局部失稳淹没高程计算结果

对该工况下模拟结果每间隔 8 min的淹没高程演进图如图12 所示。由图12 可知：溃坝后的尾矿砂在初始时溃口处的尾矿高度为 4 m。8 min 时东南侧尾矿高度突破 5 m，尾矿沿东南侧尾矿坝轴向方向高度最高，12 min 时由于地形原因最大淹没高度出现在尾矿砂东北角最大深度突破 6 m。40 min 时尾矿砂流动停止，此时尾矿最大深度为9 m；边坡局部失稳引发溃坝后，尾矿在前 20 min内对其下游地形高程改变较大，20 min 之后尾矿砂淹没高度增加不明显。最大高度为洪水漫顶引起的溃坝的 1.5倍。

a—DM1；b—DM2；c—DM3；d—DM4；e—DM5；f—DM6；g—DM7

3.3 下游各断面泥砂流量变化分析

各断面距新建尾矿库东侧交点距离分别为： DM1-溃口处，DM2-150 m，DM3-310 m，DM4-500 m，DM5-676 m，DM6-820 m，DM7-960 m 断面处的溃坝数值模型计算流量（m³ /s）过程线见图13所示。

a—0 min；b—8 min；c—16 min；d—24 min；e—32 min；f—40 min

根据图13所示，断面 1~7处的洪峰流量分别为 13000 m³ /s、9567 m³ /s、5498 m³ /s、1867 m³ /s、544 m³ / s、167 m³ /s、89 m³ /s，通过上图和上述数据可以看出：离尾矿库越远的断面水砂混合体到达时间和洪峰到达时间越长，DM1~DM7 各断面处的水砂混合体到达时间分别为：1.5 min、8 min、12 min、15 min、24 min、28 min和37 min；DM1~DM3断面处均呈现出先迅速增加到洪峰流量迅速降低，说明这几处断面地势差较大，沿程阻力较小；DM4~DM6三个断面均呈现迅速增加至一定值后缓慢降低，中途某时流量有上升趋势，最后趋近于0，DM7处先流量迅速增加到一定值然后迅速下降后又逐渐增加到峰值流量；从各断面流量曲线可以看出 DM1~DM6 基本已经停止，DM7还在流动，尾矿砂流还在向外扩展，只有当各断面流量逐渐衰减趋于0时，尾矿库溃坝停止。

3.4 下游各断面泥沙淹没高程变化分析

各断面距新建尾矿库东侧交点距离分别为： DM1-溃口处，DM2-150 m，DM3-310 m，DM4-500 m，DM5-676 m，DM6-820 m，DM7-960 m 处的溃坝数值模型计算淹没高程（m）线见图14所示，将淹没范围回溯至地形图（图15）。

a—0 min；b—8 min；c—16 min；d—24 min；e—32 min；f—40 min

由各断面处溃坝后水砂混合体淹没高程随溃坝后历时变化曲线形状可以看出：各断面最大淹没高程分别为 1550.5 m、1548.2 m、1546 m、1542.5 m、1542 m、1542 m、1540.9 m；断面 1~3 处淹没高程迅速增加后趋于一稳定值；断面 4 迅速增加后趋于稳定，在后期又有所抬升，而断面5、6处于持续上涨阶段，且上涨速度较快，断面 7 虽在持续上升，但上升高度微乎其微。

4 结论

（1）通过溃坝可能性分析，认为本尾矿库防洪标准按500年一遇设计，在正常标高下运行时，在良好的管理和安全的防洪系统下，数据显示不会发生溃坝。当水库达到最高洪水位，排洪系统全部失效，且同时遇上远超 1000 年一遇的洪水，24 h 入库洪水总量超出安全运行水位至坝顶之间库容时，可能产生漫顶溃坝风险。

（2）发生洪水漫顶溃坝时，尾矿砂淹没范围最远处距离溃口位置 960 m，最大淹没高度 6.2 m，尾矿下泄量为 582.65 万 m³，占总库容的 2.58%，溃口处最大流量为780 m³ /s。

a—DM1；b—DM2；c—DM3；d—DM4；e—DM5；f—DM6；g—DM7

（3）发生边坡局部失稳溃坝时，尾矿砂淹没范围最远处距离溃口位置1500 m，最大淹没高度9 m，尾矿下泄量为1220.34 万m³ ，占总库容的13%，溃口处最大流量为1300 m³ /s。

（4）边坡局部失稳造成的溃坝，在影响范围、淹没高度、最大流量、下泄总量方面都大于洪水漫顶溃坝。

（5）两种溃坝后下泄的尾矿砂对东北侧砂场的影响均大于东南侧，洪水漫顶溃坝的影响范围未超过 1 km。边坡局部失稳溃坝工况的东北侧的溃坝影响范围超过了 1 km，但此时溃坝已经停止，泥石流不再向前移动。两种工况淹没范围距建材园最近的点依然与建材园保持一定距离，不会造成人员伤亡以及大量经济损失，风险可控。

a—DM1；b—DM2；c—DM3；d—DM4；e—DM5；f—DM6；g—DM7

参考文献