基于Rhino 和Fluent 耦合的尾矿库溃坝数值模拟分析

孙鸿昌 郝 喆 侯永莉 尹亮亮 张 颖

(1.中煤科工集团沈阳研究院有限公司,辽宁 抚顺 113122;2.煤矿安全技术国家重点实验室,辽宁 抚顺 113122;3.辽宁大学环境学院,辽宁 沈阳 110036;4.辽宁有色勘察研究院有限责任公司,辽宁 沈阳 110033)

尾矿库是高势能的泥石流源头[1],我国尾矿库具有体积小、数量众多的特点,加之我国人口多、分布广[2],尾矿库的建造场地往往很难避开人类居住区,导致多数尾矿库下游存在人类活动[3],形成了“头顶库”。根据相关统计,截至2018年底,我国仍然存在“头顶库”1 425 座,一旦发生溃坝,将严重威胁下游居民和重要工农业、交通设施及环境安全[4]。近年来“头顶库”溃坝事故屡见不鲜,例如2008年9月,山西省新塔矿业有限公司发生尾矿库溃坝事故,波及下游500 m 左右的矿区办公楼、集贸市场和部分民宅,房屋摧毁若干,造成276 人死亡;2015年11月,湖南郴州云锡矿业尾矿库因排泄设施损坏,导致4 人失联,造成7.9 亿的经济损失以及严重的环境污染[5]。为此,一些学者开展了尾矿库溃坝模拟研究,李全明等[6]以溃坝泥浆的动力学过程作为研究重点,得到了尾矿库溃坝水砂流动的计算方法;郝喆等[7]通过ANSYS 和FLCA3D相结合,对尾矿库加高增容过程中的坝体渗流进行了稳定性分析;陈殿强等[8]采用大坝溃决的经验公式,分析了某铁矿尾矿库溃坝时的尾砂淹没范围及深度等;王兴华等[9]通过Fluent 数值模拟研究了尾矿库溃坝洪水不同时刻的流动特性;HANSON 等[10]通过多次大规模逐渐溃坝模拟试验,得出漫顶引起的堤坝溃决的形成过程和时间可以极大地影响水库放水的速率;WANG 等[11]对溃坝后下泄洪水的演变进行了数值模拟。

根据现有成果分析可知,当前的尾矿库溃坝模拟大多采用基于数学模型的理论分析方式,或者进行实际地形简化后的二维数值模拟思路。溃坝的三维数值模拟在水库大坝中应用较多,在尾矿库中应用较少,而尾矿库溃坝后形成水砂两相流,其流动状态与水库的洪水单相流有很大区别。前人多采用单一数值模拟软件进行稳定性计算或溃坝流动状态分析,在精细地形描述和高精度建模方面存在不足。溃坝模拟中多将流动状态等效为带侧板的明渠流,缺乏三维真实地形下的流动过程模拟。因此,采用专业的3D造型软件和具有强大计算能力的数值模拟软件进行耦合分析,实现高精度复杂地形下的尾矿库水砂两相流三维溃坝模拟,具有重要意义。

作为各级应急管理部门正在推行的非煤矿山安全管理新方向,尾矿库溃坝数值模拟技术标准尚未建立,分析方法、建模精度、计算内容和成果要求等也未明确。本研究以辽宁省某“头顶库”为例,采用VOF多相流模型与k-ε湍流模型理论,利用Rhino 和Fluent 软件相耦合,考虑下游存在沟谷、农田、村庄、河流和公路等状态的复杂地形条件,开展尾矿库加高后溃坝的三维数值模拟,研究漫顶溃坝后的尾砂流动规律,获得尾砂在不同时刻对下游的淹没范围、堆积情况和流动状态,对尾矿库防灾减灾和环境保护具有一定的参考意义。所提出的耦合建模分析思路,可为类似“头顶库”溃坝分析和标准建立提供一定的借鉴。

1 工程概况

某尾矿库位于辽宁省凤城市,尾矿库平面如图1所示。该尾矿库初期坝设计为堆石坝,坝底标高222.70 m,坝顶标高232.00 m,相对高度9.30 m。主坝与副坝采用上游法筑坝,内坡比为1 ∶2,外坡比为1 ∶1.5。该尾矿库为山谷型,设计加高增容后的堆积坝顶标高为285.00 m,堆积后总坝高约48.00 m。西侧堆积坝体为主坝平均外坡比为1 ∶4.5,南侧副坝外坡比为1 ∶2.4,东侧坝体为挡水堆石坝体,外坡比约1 ∶2.0,库容约1 300 万m3。该尾矿库下游为农田、公路以及农村,副坝南侧为厂区,综上,根据《尾矿堆积坝岩土工程技术规范》(GB 50547—2010)[12]的规定,按现状坝高、库容等综合考虑,该尾矿库等级为Ⅲ等库,尾矿坝级别为3 级。

图1 辽宁某尾矿库平面(单位:m)Fig.1 Plan of a tailings pond in Liaoning Province

2 溃坝数值模拟

2.1 模拟软件评述

本研究数值模拟软件采用Rhino 建模工具[13]和Fluent 流体分析软件。

(1)Rhino 软件具备较传统建模思路更为出色的NURBS 建模方式[14],可以更好地控制模型表面的曲线度,创建出更逼真、生动的复杂模型和效果图,实现高精度的复杂地形真实模拟。

(2)Fluent 软件具有较高的计算速度、稳定性和精度,强大的后处理功能,可以模拟各种类型流体的真实流动状态,但软件自带的建模功能有很大的局限性[15]。

(3)Rhino 软件和Fluent 软件相耦合,可有效解决目前数值模拟建模能力差、地形精度低、模型太过理性化、缺乏真实性的问题,能够有效体现出下游复杂地形对溃坝水砂流动的影响。

2.2 溃坝模拟区域

在Rhino 软件中建立尾矿库加高增容扩建工程终期条件下的溃坝三维计算模型,模型内包括初期坝、堆积坝及所在区域的实际地形。研究区域长度为3.1 km,最大宽度1.8 km,研究范围约4.65 km2,研究区域内主要有青城子村、新式街、松树沟、双下等4个村庄以及溃坝下游农田。尾矿库建模范围如图2所示。

图2 尾矿库建模范围Fig.2 Modeling scope of tailings pond

2.3 初始条件假设

考虑到尾矿库生产时间不断加长,排洪设施有可能发生尾砂下陷、尾砂沉积造成堵塞,或者排水斜槽发生坍塌,预制管接口渗漏等情况[16]。在这种极端情况下排洪能力下降或消失,库内水位升高,浸润线埋深抬升过高,从而导致尾矿库发生溃坝[17]。为简化计算并进一步提高模拟结果与实际情况的契合度,本研究做出两点假设:① 假设500 a 一遇洪水过程超过10 h,库内水位持续上涨达到坝顶,洪水位上涨的时间过程在数值模拟中忽略;② 假设库内排洪系统失效,计算过程中的库内水砂主要从溃口下泄,且溃坝模拟求解中无其他外来水量叠加计算。

2.4 模型建立及边界条件

Fluent 软件的建模能力有限,而Rhino 软件中可以通过等高线建立复杂溃坝模拟的地形图。为此,本研究利用Rhino 建模和Fluent 计算进行耦合,使溃坝数值模拟结果更符合实际情况。将构建的模型导入Fluent 软件中进行计算,在Fluent 中设置计算模型和边界条件,具体如下:

(1)读入网格模型,检查网格模型质量,若有体积量为负的情况,需要重新定义网格。该模型的网格数量为350 万。地表形态及主要网格见图3。

图3 地表形态及主要网格Fig.3 Surface morphology and main grid

(2)模型设置。选取Transient 非稳定瞬态求解,计算模型选取VOF 多相流模型,相数为2,k-ε湍流模型。

(3)设置材料参数。材料为溃坝尾砂和空气两种,由于溃坝常伴随着暴雨等条件的发生,溃坝时的尾砂常处于饱和状态,所以根据实际调查和试验结果,选取尾砂的平均密度为1 750 kg/m3,黏度系数为15 Pa·s,空气选取标准大气压下的空气参数。将空气设置为主相,溃坝尾砂设置为次相。

(4)边界条件。将模型上部空气界面设置为压力入口,四周界面为压力出口,山体地形和坝体为墙体。

(5)流场初始化。选取压力入口作为流场初始化条件,初始流场基本处于静止状态,将尾砂部分设置为第二相即尾砂材料,使其体积分数为1。计算在自重作用下的溃坝流动情况。

(6)迭代求解。设置时间步为0.01 s,每次迭代循环步数20 s,迭代总步数18 万步。

3 模拟结果与分析

本研究采用VOF 计算模型与k-ε湍流模型理论相结合的计算方法[18]。VOF 法在计算域内,引入了相体积分数F,即在网格内各相流体的体积与流体空间总体积的比值[19];选取适用于剪切流动强、边界层分离等复杂流动情况[20]的k-ε湍流模型。

尾矿库发生溃坝后,库区内水砂向下游倾泻而下,水砂的淹没范围能够直观地反映出不同时刻水砂对下游的淹没情况;水砂对地面的压力能够计算出水砂的淹没深度;下泄水砂端头速度大小能够反映出水砂对下游构筑物冲击力的大小。因此根据模拟结果,对溃坝后水砂的淹没范围、压力、速度进行进一步分析。

3.1 溃坝水砂淹没范围

水砂流的演进情况可表示尾矿库溃坝的全部过程。根据分析结果,提取了尾矿库发生溃坝过程以及溃决水砂流中不同时刻流动演进情况的体积分数(Volume fraction)云图,如图4所示。

图4 不同时刻溃坝水砂流演进云图Fig.4 Nephogram of water and sand flow evolution of dam break at different times

由图4 可知:当尾矿库排洪系统失效同时遭遇500 a 一遇洪水发生溃坝,坝体破坏状态表现为库内水砂漫过主坝坝顶,从坝顶倾泻而下,随后主坝顶破碎,受洪水拖曳力的作用,坝坡尾砂受冲刷下泄,库内洪水诱发尾砂流态下泄,演化运动过程位移变化明显。t=0 时,尾矿堆积坝体处于相对静止状态,随着时间不断增长,堆积坝坝顶溃口持续发展,在0～100 s 时,溃口附近尾砂由于水流拖曳力作用冲刷形成初期水砂混合流,后方尾砂随之开始滑动,临界流动面向后方扩展;在100～300 s 时,溃口后方尾矿堆积体逐渐滑出,在后方滑出水砂的流动推动下,坝坡尾砂出现冲刷沟壑;400 s 左右前端砂流淹没初期坝,此后阶段为水砂流在初期坝下游的演进阶段;在600～900 s 时,水砂流在初期坝外坡脚至所在沟谷口范围流动运移;在900～1 500 s 时,溃坝的势能转化为动能过程中伴随尾砂料之间的碰撞摩擦导致能量消耗,前端水砂流到达低势能区,在下游村庄和公路所在平地减速缓慢演进淤积;t=1 800 s 之后尾矿坝体滑动逐渐达到新的平衡状态,坝体垮塌运动趋于停止,此时水砂到达下游的1 035 m 位置。此次溃坝水砂量约297万m3,水砂覆盖下游面积约0.98 km2。

3.2 溃坝水砂压力(淹没高度)情况

尾矿库下游地形条件的复杂情况对溃坝的影响程度较大[21]。若下游高差大,相应的溃出距离及影响范围更大;若尾矿库下游区域为平原开阔区域,在纵坡度较缓的情况下,相对其他下游地形条件,影响距离较小。因此,根据分析结果,提取了1 800 s 尾矿库溃决后水砂对地面压力云图并选取距初期坝外坝脚250、500、750、1 000 m 4 处作为水砂压力监测点,如图5所示。

图5 水砂对地面的压力分布及监测点位置Fig.5 Pressure distribution of water and sand on the ground and the position of monitoring points

由图5 可知:大部分尾砂在距离初期坝外坡脚500～600 m 区域内淤积。根据压强公式,折合尾砂最大淹没深度为7.2 m,位于下游500 m 处的沟谷内,溃坝下游其他位置水砂的淹没高度为2～3 m。

本次溃坝模拟泥石流渐进式泄向下游地区。根据监测点处的压力变化规律,绘制了水砂压力随时间的变曲线,如图6所示。

图6 与坝脚不同距离处的压力—时间曲线Fig.6 Pressure-time curves at different distances from the foot of the dam

由图6 可知:各位置处水砂压力在不同时刻出现峰值,反映了溃坝水砂在下游逐步淹没的进程。水砂流到监测点位置后,压力以约1.1 kPa/s(合0.064 m/s)的高梯度增加到最大值,然后迅速减低,随着时间增加,监测点压力衰减速度越来越缓慢,最后趋于稳定不再变化。下游750 m 处的农田地势较低,压力略高于250 m 和1 000 m 处,下游500 m 处地势低洼且有大量尾砂在此堆积,所以水砂压力最大。随后由于流距增长以及地面阻力的共同作用,水砂压力开始减弱,最后趋于稳定。

3.3 水砂速度演进情况

溃坝后不同时刻的水砂速度云图如图7所示。由图7 可知:t=300 s 时,溃坝水砂到达初期坝底部,由于溃口与坝底存在48 m 的高差,溃口水砂在重力势能和堆积坝体摩擦的共同作用下从0 迅速增大;随后由于溃口增大,水砂流量加大;t=600 s 时,水砂流经下游500 m 处,水砂流速进一步增加,快速向下游流动;随着时间不断增长,堆积坝坝顶溃口持续发展,受下游地势平缓和地表面粗糙程度的影响,水砂的整体流速降低,尾砂在500 m 处堆积;此后由于尾矿库水砂重力势能降低;t=1 800 s 时,溃口与库区内水位逐步持平,水砂的流速非常小,不会向下游产生较大区域的流动。

图7 不同时刻溃坝水砂速度云图Fig.7 Water and sand velocity nephogram of dam break at different times

溃坝水砂端头速度变化曲线如图8所示。分析图8 可知:0～300 s 时,速度迅速增加,300 s 时溃坝水砂端头速度达到2.9 m/s;300～600 s 时,速度持续增加,在600 s 时端头速度达到最大值3.7 m/s,此时水砂的冲击力最大;之后水砂流速逐渐降低,到1 800 s 时流速趋于0 停止不动。

图8 溃坝水砂端头速度—时间曲线Fig.8 Velocity-time cures of water and sand at the end of dam break

4 结 论

针对精细真实三维地形下尾矿库溃坝模拟分析存在的不足,通过Rhino 建立高精度地形图,与Flunet 模拟软件相耦合,分析了溃坝水砂受真实地形影响的流动状态,得出以下结论:

(1)采用VOF 多相流计算模型与k-ε湍流模型理论,利用Rhino 建模和Fluent 计算相耦合,可以有效实现高精度地形下的水砂两相流溃坝三维数值模拟,获得尾砂在不同时刻对下游的淹没范围、堆积特征和流动状态。

(2)“头顶库”发生溃坝后,在下游地势较平缓和谷口较宽阔的情况下,发生溃坝后水砂的淹没范围广泛,其下游水砂淹没形状近似为三角形:距尾矿坝较近位置淹没范围宽,水砂前端受沟谷地形影响淹没范围窄。尾矿库溃坝后的溃决持续时间长,水砂的淹没范围达到下游的1 km 左右。

(3)溃坝水砂渐进式泄向下游地区,溃坝水砂在下游各处压力都具有先增大后减小,最后趋于稳定的特点。大部分尾砂会淤积在溃坝下游沟谷内,导致压强高于其他地区。

(4)发生溃坝后,溃口水砂在重力势能和堆积坝体摩擦的共同作用下从0 迅速增大,在某一时刻,流速达到峰值,此时水砂冲击力最大,对下游破坏力最强;随后,溃口水位不断降低,受下游地势平缓和地表面粗糙程度的影响,水砂流速缓慢降低不会向下游产生较大区域的流动。

(5)在汛期来临前,建议拓宽加深下游河道,增强沟谷的输水性和堆积能力,降低水砂流通性,减缓溃坝水砂流速和淹没范围,防止对下游产生较大的冲击力和危害;在下游建立拦挡导流设施,降低水砂的淹没深度,延长水砂到达下游村庄的时间,保障居民生命财产安全。根据本研究尾矿库溃坝水砂淹没范围的模拟结果,建议在下游淹没范围的1.5 倍区域内禁止建设厂矿和村落。