胡良才,郭大平,李哲辉,李玉雷
(中核第四研究设计工程有限公司,河北 石家庄 050021)
尾矿库是矿山的三大控制性工程之一,用于贮存水冶厂选冶矿石后剩余的大量尾矿[1],受筑坝材料特性及缺乏管理等因素影响,尾矿库发生溃坝事故的概率远远高于传统水库[2]。尾矿库一旦溃坝,产生的泥石流将快速冲泄向尾矿库下游,破坏农田以及交通设施等,给下游的居民造成严重的生命及财产损失[3]1。
国内外诸多学者采用经验公式法、模型试验、数值模拟等手段对尾矿库溃坝问题进行了研究。陈殿强等[4]针对凤城市某尾矿库,采用经验公式计算得到尾矿库溃口宽度、溃坝坝址最大流量、溃坝最大流量沿程演进、砂流传播时间等,并根据计算结果提出了尾矿库预防及应急措施。连国玺等[5]采用经验公式计算得到某铀尾矿库溃坝影响范围。尹光志等[6]采用模型试验方法,研究了不同高度下尾矿坝瞬时全部溃坝情况下泥浆的流态、演进规律以及动力特性,研究了泥浆在尾矿坝下游流动过程中的冲击力强度、淹没范围及演进规律。张力霆等[7]利用自主研发的尾矿库模型试验平台,研究了尾矿库溃坝破坏形式。刘磊等[8]采用物理模型试验方法,对尾矿坝漫顶溃坝洪水及溃口变化过程进行了研究。林江等[9]在溃坝洪水数学模型的基础上,加入非恒定流输沙方程,建立能同时模拟水流和泥沙运动的二维溃坝水流泥沙数学模型,并对云南某尾矿库溃坝后的洪水演进及矿砂在下游的淤积过程进行了二维数值模拟。金佳旭等[10]采用非牛顿流体模型中的Bingham模型,针对辽宁某尾矿库溃坝后的砂流演进过程进行了数值模拟,重点分析了溃坝后泥石流流态变化、速度矢量变化及最终堆积形态。霍紫玮[11]采用强度折减法分析某尾矿库坝体稳定性,计算得出坝体溃坝范围,利用计算流体力学软件对尾矿库溃坝进行了数值模拟,得出了溃坝泥石流淹没范围、流动速度和淹没深度等重要指标。Marina Pirulli等[12]采用数值模拟方法,研究了意大利Stava Valley尾矿库溃坝泥石流流动特性及影响范围。陈俊等[13-14]采用数值模拟方法研究了某铜矿拟建尾矿库溃坝情况,通过设置监测点,从溃坝泥石流流场、水深和尾砂淤积方面分析溃坝泥石流对下游村庄的影响。
经验公式法不能考虑地形影响,计算结果往往与实际情况存在较大误差;模型试验法费用高、试验周期长,同时存在有缩尺效应,无法准确研究诸多工程中关心的实际问题。数值模拟法费用低、计算周期较短,且可以不需考虑缩尺效应的影响,直接对工程原型进行模拟研究,可以获得大量详实的资料,还可以进行方案比较选择最优设计或试验方案,指导模型试验的进行,节省大量试验时间[3]5。随着计算机硬件及数值模拟理论的发展,数值模拟技术越来越多的应用于尾矿库的溃坝研究。笔者针对某铀尾矿库,建立包括尾矿库坝体及下游三维精细地形的数值模型,进行尾矿库溃坝三维数值模拟研究。
RNGk-ε紊流数值模型可更好地模拟流体流线弯曲程度较大的流动[15],因此采用RNGk-ε紊流数值模型来封闭基本方程组,对铀尾矿库溃坝进行三维数值模拟研究,控制方程[16]如下。
连续性方程:
(1)
动量方程:
(2)
紊动能k方程:
(3)
耗散率ε方程:
(4)
采用有限体积法来离散控制方程,VOF(Volume of Fluid)方法[18]追踪水流的自由表面。VOF方法定义一个体积函数F,若单元体内充满流体,F=1;若单元体内无流体,F=0;单元体内同时含有流体以及空气,0 相间界面的追踪方程可表示为 (5) 采用Martin等[19]的水槽试验验证溃坝数值模型。根据数值模拟软件特点,建立与试验模型尺寸相同的三维数值模型,水槽试验及数值模型示意图如图1所示。 模型计算区域长度l=12.7 cm,宽度b=5.72 cm,高度h=6.72 cm。采用结构化网格,顺水流方向为x方向,重力方向为z方向,网格的单元尺寸Δx=0.025 cm,Δy=Δz=0.1 cm,设置计算软件自动调整时间步长大小,确定初始的时间步长Δt=0.001 s。采用VOF方法追踪流体自由表面,RNGk-ε紊流数值模型封闭控制方程。数值模型上游、两侧以及底部的边界条件设定为法向速度为零、无滑移的固体边壁;模型出口为自由出流,设置为压力出口边界条件;模型顶部给定为压力入口边界条件。确定初始水体长l0=2.86 cm,宽b0=5.72 cm,高h0=5.72 cm。水槽试验与数值模拟水流波前同时间关系曲线如图2所示。 由图2可见,对于水流与时间关系的模拟,数值模拟与模型试验结果规律一致,误差较小,说明采取的数值模型真实可靠。 采用三维建模软件CATIA建立尾矿库库区及下游河道三维地形模型:1)使用Autodesk Civil 3D软件,根据实测地形图建立地形曲面,并提取库区及下游河道地形坐标及高程,导出点云数据;2)进入CATIA的Digitized Shape Editor模块,导入点云数据,过滤、修补点云数据(图3a);3)根据导入的点云数据生成mesh网格面,分析、检查、修补网格面(图3b);4)进入Quick Surface Reconstruction模块,生成实体地表面(图3c);5)进入Part模块,将地形边界投影到基准面上作为草图边界,拉伸草图至曲面,形成三维地形实体(图3d)。 某铀尾矿库主要由初期坝、尾矿堆积坝、副坝、排洪设施等组成。初期坝高25 m,坝顶标高125 m,尾矿堆积坝高27 m,总坝高52 m。初期坝上游坝坡坡度为1∶2.5,下游坝坡坡度为1∶2~1∶2.5;初期坝内库容堆满后,后期采用尾矿砂堆积筑坝,尾矿堆积坝坡度为1∶3.5[20]。基于实测的1∶1 000比例地形图,建立铀尾矿库局部溃坝三维模型,模型长4 000 m,宽3 000 m。模型计算区域包括尾矿库库区、坝体及尾矿坝下游5 km范围内河道。为了避免因计算过程中发生封顶现象而导致计算失败,设定库区最大模拟高程为160 m,高于尾矿库现状最大坝高[17]140。 假定尾矿库溃坝产生的泥石流为介于流体和散粒体之间的一种特殊的运动形式,可采用流体流动的能量方程和连续方程近似描述,数值模型采用如下基本假定:1)尾矿库库区基岩及周围山体为不透水边界;2)溃坝泥石流演进过程中,不考虑单个颗粒的体积变形;3)溃坝泥石流为各向同性连续介质流体,其流动符合宾汉流动模型的规律;4)不考虑溃坝泥石流对下游河道的冲刷。 尾矿库溃坝从形成决口到基本形成稳定的溃决断面,时间过程非常短暂,为安全考虑,尾矿库溃坝按瞬时溃坝处理。根据某铀尾矿库坝体稳定性分析计算得到的坝体滑动面位置[21],确定尾矿坝滑移破坏区,即为局部溃坝范围。对校核工况下,尾矿库瞬时局部溃坝情况进行三维数值模拟。根据黄河水利委员会水利科学研究院实际资料分析求得溃坝决口平均宽度b,计算公式为 b=k(W1/2B1/2H0)1/2, (6) 式中:b—溃坝决口平均宽度,m;W—溃坝时库容,万m3;B—坝顶长度,m;H0—坝前水深(水头),m;k—与坝体土质有关的系数,对黏土k值约为0.65,壤土k值约为1.3。计算得到尾矿库坝体溃口宽度为102.9 m。 采用结构化网格对数值模型进行网格划分,自尾矿库坝体向下游方向为x方向,重力方向为z方向。经过多次试算后确定计算网格尺寸为Δx=Δy=6.0 m,Δz=2.0 m,网格总数约1 800万。设定初始时间步长Δt=0.001 s,设置计算软件可以自动调整时间步长大小。采用VOF方法追踪流体自由表面,RNGk-ε紊流数值模型来封闭紊流基本方程组。尾矿库库区模型顶面设置为压力入口边界,压力为101.325 kPa;库底及下游河道设定为法向速度零、无滑移的固体边壁,采用标准壁面函数法处理近壁的黏性底层;下游河道出口方向为自由出流,设置为压力出口边界,压力为101.325 kPa。尾矿库下游为水田、灌木丛等,下垫面情况较为复杂,根据国内外相关研究成果,下游河道曼宁系数n取0.03。计算初始条件为库区内坝顶标高以下水的体积分数为1,即库区被水体充满。在尾矿库坝址处至坝下游5 km处,每1 km设置1个监测点,共计6个监测点,以监测溃坝泥石流流速、传播时间、泥石流深等水力要素变化情况。 图4为尾矿库坝址处流速变化过程线。由图可见,尾矿库局部溃坝情况下,坝址处溃坝泥石流流速随时间增加迅速增大至峰值,之后逐渐下降,最终趋于零。溃坝5 s后坝址处流速达到峰值,最大流速约为14.12 m/s。 尾矿库溃坝后,溃坝泥石流快速冲泄向下游,从洪水到达至最大洪峰到达时,两者时间间隔仅为数秒至数十秒,且越靠近尾矿库,间隔时间越短,从而留给下游民众安全撤离的时间越短。因此,铀尾矿库运行管理过程中,应加强监测,提前预警,为居民安全撤离争取宝贵的时间。局部溃坝情况下,尾矿库坝址及下游各特征断面洪水起涨时间、最大洪峰传播时间、洪水消落时间见图5。 尾矿库坝址及下游各特征断面监测点泥石流深变化见图6。溃坝事故发生后,坝址处及下游各特征断面溃坝泥石流深随时间增加迅速增大,之后逐渐下降。在溃坝40 s后,坝址处泥石流深达到最大值,约为12.31 m。 根据数值模拟计算得到的各特征断面最大泥石流深,结合尾矿库实测地形图,绘制局部溃坝泥石流影响范围见图7。在尾矿库下游的潭元等几个村庄设置有监测点,各监测点均未监测到溃坝泥石流,说明局部溃坝情况下,溃坝泥石流未直接冲击影响到这几个村落。 基于计算流体力学软件Flow-3D建立溃坝三维数值模型,采用水槽试验验证数值模型,结果表明建立的数值模型能很好的模拟水流运动情况。采用三维设计软件CATIA建立某铀尾矿库坝体三维模型及下游河道三维精细地形,基于VOF方法、RNGk-ε紊流模型,建立某铀尾矿库溃坝三维数值模型,对尾矿库局部溃坝进行的三维数值模拟表明:尾矿库溃坝后,溃坝泥石流流速、深度等迅速增大到峰值,之后逐渐平稳下降;局部溃坝情况下,尾矿库溃坝泥石流不会直接冲击影响下游各村庄;RNGk-ε紊流模型可用于对尾矿库溃坝泥石流水力特性进行三维数值模拟,研究成果可为工程设计、运行管理提供很好的借鉴。1.2 模型验证
2 尾矿库溃坝数值模拟
2.1 地形建模
2.2 尾矿库三维数值模型
3 结果分析
3.1 流速变化分析
3.2 传播时间分析
3.3 泥石流深结果分析
3.4 影响范围分析
4 结论