基于FLOW-3D的尾矿库数值模拟对下游影响研究

2021-03-13曾智超

江西水利科技 2021年1期

梁萱，曾智超

（1.江西省水利科学院，江西南昌330029；2.江西省水利规划设计研究院，江西南昌330029）

0 引言

尾矿库即是矿场提炼矿石中的精矿后将剩余废渣排入库内进行堆存和沉积的贮存场所，一般在山谷口部、河道或凹地的周边筑坝。由于尾矿库溃坝不仅会使得库区内的水流泄出，并且其裹挟尾矿等固体废渣对下游造成的危害更甚于一般的堤坝，从而引起下游地表水及地下水体的污染[1]。而且，由于尾矿坝规模及范围较庞大，且库内以包含重金属等污染物质的尾矿为主要部分，一旦发生溃坝、渗流等事故，库内尾矿砂大量冲泄，壅塞河道，冲垮房屋、农田、道路和桥梁，必将严重威胁下游群众的生命财产安全，并危害库区周边的生态环境，破坏之大，往往令人触目惊心[2]。

洪水漫顶是导致尾矿库溃坝的主要原因之一，但是由于尾矿库一般处于偏僻恶劣地区，且溃坝发展相对较快，大多数对尾矿库溃坝的研究都是溃坝后通过查阅资料、问访等方式进行调查，因此现场的监测数据尤为匮乏[3]。采用数值模拟研究可以方便直观的呈现尾矿库溃坝进程，同时，也可以克服试验不易修改和耗费时间的缺点，并且，也可利用计算机强大的运算能力，解决复杂的工程问题。姜清辉、胡利民[4]等人通过总结分析近年来国内外尾矿库溃坝研究进展，提供了预防尾矿库溃坝的可能方法；马海涛、张亦海[3]等人通过对国内外尾矿库物理模型试验研究资料的收集和分析，总结尾矿库溃坝机理及影响因素，为监测预防尾矿库溃坝提供了建议；林长强[5]基于FLUENT软件提供的自定义函数功能，对于溃口实时变化编制了动网格程序，实现了水流的三维逐渐溃坝数值模拟，计算得到溃口流量过程及下游河道洪水演进情况。

因此，本文以实际工程为例采用数值模拟的方法对尾矿库动态溃决过程进行模拟，通过采取工程措施前后对比分析，分析采取工程措施的有效性和必要性，该成果可应用于降低和规避尾矿库溃坝的风险，对其设计、施工和运行全过程都具有重要的指导意义。

1 尾矿库溃坝统计分析

通过对美国、欧洲和中国从2001年至2015年发生的尾矿库溃坝事故的统计数据[6～8]的分析可以看出，造成尾矿库溃坝的主要原因大致可以分为六类：洪水漫顶、稳定破坏、渗流破坏、结构破坏、其它原因以及不明确的原因。图1中给出了不同类型溃坝事故模式所占的百分比。其中，稳定破坏造成的溃坝事故数量最多，共有52例，占28.57%；结构破坏导致的次之，共39例，占总数的21.43%。结合图2可以看出美国和中国发生尾矿库溃坝失事的主要原因是稳定破坏和结构破坏，而欧洲发生尾矿库溃坝失事的主要原因是其他各因素导致的。当尾矿库发生洪水漫顶时往往伴随极端灾害天气，应急抢救难度大，漫顶溃口一旦启动，溃坝难以控制。因此，本次溃坝研究的主要针对洪水漫顶溃坝工况。

图1 尾矿坝溃坝失事原因分析

图2 尾矿坝溃坝事故原因分析

2 尾矿库溃坝数值模拟

2.1 尾矿库基本情况以及三维建模

本文中某尾矿库中的尾矿坝的初期坝型采用透水堆石坝[9]，后期堆积坝采用上游式放矿堆坝工艺。初期坝顶高程325.0m（黄海高程，下同），坝高25.0m，坝轴线长120m，坝顶宽5.0m，上游坡比1：2，上游坝坡布置反滤层，下游坝坡坡比1：1.75。尾矿库后期堆坝平均坡比1：3.5，尾矿库最终堆积高程取361.0m，具体剖面尺寸图如图3所示。

图3 尾矿坝剖面图

根据尾矿库实际地形图生成地形曲面，而后建立初期坝及堆积坝三维模型形成尾矿坝三维模型。如图4所示。

图4 尾矿库整体枢纽及下游村庄整体模型实体

2.2 逐渐溃坝计算模型及其边界条件

2.2.1 尾砂的基本参数

在尾矿库溃坝的数值模拟的泥沙模型中，尾砂特性控制参数如表1所示。

表1 尾砂特性控制参数

其中，临界希尔兹系数可用Soulsby-Whitehouse方程[10]计算：

其中，θcr为临界希尔兹系数，无量纲；d*为泥沙颗粒的当量粒径，mm；sn=ρn/ρ，vf为流体的运动粘度，m2/s；τcr为泥沙床面的临界剪切力，Pa；ρn为泥沙颗粒的密度，g/cm3；ρ为水流密度，g/cm3。

挟带系数、推移质系数和休止角取默认值。

2.2.2 网格划分及监测点、截面设置

根据所建立的整体计算模型，模拟区域为库区及其下游3km范围，为保证计算精度并提高计算效率，模拟区域共划分为5个网格块，总有效网格共计310万左右。

为分析溃坝尾砂下泄过程泥沙淤积情况，共选取9个监测点进行分析。监测点位置如图5所示。为分析溃坝尾砂下泄过程中流量及速度情况，共截取12个河道断面进行分析，断面位置如图6所示。

图5 监测点位置图

图6 各流量监测断面位置

2.2.3 采取工程措施后的计算模型

采取工程措施后的计算模型如图7所示，此模型在尾矿坝坝脚下270m处设置刚性拦挡坝，坝高15m，坡度1:1，坝顶高程290m。计算模型的前处理设置、计算求解初始和边界条件、监测断面和监测点的设置相同。

图7 采取工程措施后模型图

2.2.4 边界条件及来流条件设置

（1）来流条件设置。为了提高计算效率保证计算效果，本工程数值模拟计算中对洪水过程线做如下等效处理：

（a）溃坝的起始时刻为尾矿库内洪水漫顶的开始时刻，洪水填满库容的时间过程在本次数值模拟中忽略。

（b）以库内尾水高程至尾矿库坝顶高程之间的库容为依据，在设计洪水过程线中找到设计洪水刚好填满尾矿库存水库容（即漫顶）的时刻Ts，在Ts时刻以前的洪水过程线中的洪水总量在数值模型计算中以洪水满库库水的形式考虑。

（c）针对在Ts时刻以后的洪水过程线，计算Ts时刻后设计洪水的洪水总量，以该洪水总量除以数值模型计算总时间（本文中数值模型计算时间为1 000s）即可得到换算的流量值，该流量值即为数值模型中的流量边界条件输入值。

（2）边界条件设定。网格块相连接部分均设为对称（Symmetry）边界；有山体阻挡的部分设为壁面（Wall）边界；底面均设置为壁面（Wall）边界；顶面均设置为压力（Pressure）边界，也即自由液面；尾矿库库尾上游边界考虑洪水来流条件，设为流量边界；下游出口处为模拟的流体流出边界，故设为自由出流（Outflow）边界。

（3）初始溃口的选择。由于在数值模型中无法模拟尾矿坝在洪水漫顶过程中的薄弱部位，本模型在模拟时通过人为设定初始溃口来模拟尾矿坝渐溃的薄弱部位，考虑渐溃的最不利情况，初始溃口设置在最大坝高对应的坝顶位置。该位置溃口发展至最终结束时，下泄尾砂量最大，对下游的影响最严重。

2.3 计算结果与分析

2.3.1 各断面流量对比

各断面流量最大值统计如表2所示。从表中数据可以看出，增加拦挡坝之后断面最大流量值为2 093.24m3/s，发生在溃坝后360s的2-2断面处。因为拦挡坝建在断面2-2的下游，图8给出了3-3断面之后采取工程措施前后流量对比图，从图8中数据可知，断面3-3的最大流量明显减少，结合水流流场图可以看出拦挡坝中拦住了一部分水流，下游溃坝水流减少，水流到达断面时间明显延迟，从断面7-7的数据可以发现，流量随时间变化的趋势是一致的。10-10断面之后无水流经过。由此可见，增加拦挡坝可以有效挡住水流，降低流量，延缓水流的进程。

图8 采取工程措施前后各断面流量对比图

表2 采取工程措施前后各断面流量对比值

2.3.2 各断面平均流速对比

采取工程措施前后各断面平均流速对比如图9所示，各断面流量最大值统计如表3所示。由于拦挡坝建在断面2-2的下游，对比采取工程措施前后监测断面的平均流速随时间变化图可以看出，断面3-3的平均流速明显降低，结合表3的对比数据，断面3-3的平均流速减小了近一半，发生时刻也明显延迟。采取工程措施后，10-10断面之后无水流经过，流速为0。

图9 采取工程措施前后各断面平均流速图

表3 采取工程措施前后各断面最大平均流速值对比

2.3.3 监测点水位对比

结合监测点位置图5，对比分析采取工程措施前后各监测点的水深情况，不难看出，采取工程措施后可以适当减少水流的淹没范围。图10给出了各监测点采取工程措施前后水深对比图。

图10 采取工程措施前后各监测点水深对比图

监测点1～3位于下游沟口左侧偏高的位置，增加拦挡坝后经过三个点的水流明显减少，对监测点附近房屋的影响可以忽略不计。结合整体水位场图可以看出，采取工程措施后，水流到达监测点4～6的时间明显延迟，但是由于监测点4、5、6处于山谷的最低处并且是水流的必经之处，加上拦挡坝坝内库容有限，因此对于监测点4～6来说，水深影响不大。

结合各断面流量流速图可以看出，增加拦挡坝后，水流的流速降低，淹没范围减小，大部分水流聚集于监测点6与监测点7之间的洼地上，监测点7～9的水流影响减小。采取工程措施前后各监测点水深数据对比见表4。