考虑矿砂淤积的尾矿坝溃坝洪水数值模拟.pdf

1994-2009 China Academic Journal Electronic Publishing House. All rights reserved. 文章编号10072228420090520150204 考虑矿砂淤积的尾矿坝溃坝洪水数值模拟 林 江1,张绪进1,张小峰2,董炳江2,李 俊3 1. 重庆交通大学西南水运工程科学研究所,重庆400016 ; 2.武汉大学水资源与水电工程科学国家重点实验室,武汉430072 ; 3.长江上游水文水资源勘测局,重庆400014 摘 要为模拟尾矿库溃坝稀性泥石流的演进,在溃坝洪水数学模型的基础上,加入非恒定流输沙方程,建立能同 时模拟水流和泥沙运动的平面二维溃坝水流泥沙数学模型。最后将该模型应用于云南某尾矿坝溃坝后的洪水演进及矿 砂在下游的淤积过程,模型计算结果可为灾后影响评估及预防措施制定提供依据。 关键词溃坝;数值模拟;尾矿坝;稀性泥石流 中图分类号TV 143. 4 文献标识码A Numerical Simulation of Mine Tailing Dam2Break Considering Sediment Deposition LIN Jiang1,ZHANG Xu2jin1,ZHANG Xiao2feng2,DONG Bing2jiang2,LI Jun3 1. Southwest Research Institute of Water Transport Engineering , Chongqing Jiaotong University , Chongqing 400016 , China ; 2. State Key Laboratory of Water Resources and Hydropower Engineering Science , Wuhan 430072 , China ; 3. Bureau of Hydrology and Water Resources Survey of the Upper Yangtze River , Chongqing 400014 , China Abstract To simulate the program of dilute debris from the break of mine tailings dams , a horizontal 22D numerical model is estab2 lished to simulate both the dam2break flow and unsteady sediment transport processes. Finally a case study is made. The dam is loca2 ted in Yunnan Province. The model is used to simulate the dam2break flow and sediment deposition processes , after overtopping of the dam result from heavy rainfall. Simulated results can be used for the planning of the countermeasures and risk assessment. Key words dike burst ;numerical simulation;mine tailings;dilute debris flow 收稿日期2008201208 基金项目“十一五” 三峡工程泥沙研究项目 115 20323 ;“十一五” 国家科技支撑项目 “三峡工程运用后泥沙与防洪关键技 术研究”2006BAB05B02。 作者简介林 江19832 , 女,硕士,主要从事水力学及河流动力学 研究。 0 引 言 汛期尾矿坝一旦溃决,水流将携带库内堆积的大量尾矿砂 冲向下游,对水库下游的村庄、 农田以及交通设施等造成巨大 影响,严重威胁坝下人员安全和当地经济发展,同时尾矿砂的 外流还会给当地环境带来污染。如1962年9月26日云南锡 业公司火谷都尾矿库溃坝,涌出尾矿330万m3,水38万m3,冲 毁耕地36. 05 km2,11座村寨和一个农场,共13 970人受灾,其 中死亡171人、 伤92人,直接经济损失2 000万元[1];2000年 10月18日,广西南丹县大厂镇铜坑矿区选矿厂发生尾矿坝坍 塌的特大安全事故,造成28人死亡,56人受伤,位于尾矿库下 游的70间居民住房不同程度地被冲毁、 冲坏,直接经济损失 340万元[2];2006年4月30日陕西商洛市镇安县金矿尾矿库 坝发生溃坝事故,造成17名群众下落不明,5名群众受伤,米粮 河上游遭受污染[3]。 为了评价尾矿坝溃坝后的洪灾风险,积极探求切实可行的 防洪减灾措施,这需要充分把握洪水淹没和下游矿砂淤积随时 间的变化。研究开发能模拟这种情况的数学模型具有十分重 要的意义。 不同的尾矿坝溃坝会产生不同流态的泥石流,当为稀性泥 石流时,不仅需要精确模拟非恒定含间断波的洪水演进,同时 还需要考虑溃坝水流作用下的矿砂外泄,及在下游河道、 农田、 村庄和工厂的淤积过程,洪水波的演进计算应与矿砂冲淤计算 同时进行。国内外对于尾矿坝溃坝数值计算研究的文献比较 少,国内陈青生[3]采用MAC法对尾矿库溃坝砂流进行计算模 拟,取得了一定的效果,但其只考虑了砂流单一流情况,如若为 稀性泥石流,则需从水砂两相流来考虑尾矿坝溃坝问题。 有限体积法求解的是浅水方程积分形式的守恒型方程组, 离散后仍能处处保持水量与动量的平衡,具有激波捕捉性,因 而能处理溃坝水流这类的间断水流问题。本文在基于有限体 积法的平面二维溃坝水流精确模拟的基础上,耦合非恒定流输 051中国农村水利水电 2009年第5期 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved. 沙的计算,考虑了悬移质和推移质作用下的矿砂淤积。在建立 精确溃坝非恒定水沙模型的基础上,准确地预测和模拟尾矿坝 溃坝后下游洪水行进过程及矿砂淤积过程,对灾后的生命损 失、 经济损失和社会环境影响的评估,以及防洪抢险措施的决 策实施等具有重要的实际意义。 1 溃坝水沙数学模型 1. 1 基本方程 在笛卡儿坐标系下,平面二维洪水运动的连续方程和动量 方程可表达如下 9Z 9t 9M 9x 9N 9y 01 9M 9t 9uM 9x 9vM 9y -gh 9Z 9x - gn2uu2 v2 h1/ 3 γt9 2 M 9x2 92M 9y2 9N 9t 9uN 9x 9vN 9y -gh 9Z 9y - gn2vu2 v2 h1/ 3 γt9 2 N 9x2 92B 9y2 2 平面二维泥沙连续方程形式如下 9hs 9t 9Ms 9x 9Ns 9y - ε9 2 hs 92x 92hs 92y -ρ ′ 9Z0 9t 3 悬移质冲淤引起的地形变形方程 ρs 9Zbs 9t α ωs - s34 推移质冲淤引起的地形变形方程 ρb 9Zbg 9t 9gbx 9x 9gby 9y 5 悬移质水流挟沙力公式 s3 Ks U ghω m 6 推移质输沙率公式 gbx Kb ρ ρb-ρ ρbu1- Uc u u3 ωC2 u2 v2 Uc gby Kb ρ ρb-ρ ρbv1- Uc v v3 ωC2 7 gbx gby0 u2 v2≤Uc8 上述各式中h为水深; Z为水位; u和v为x和y方向的流速; Muh; Nvh; n为Manning糙率系数;γt为紊动粘性系数; t 为时间; g为重力加速度; Zbs、Zbg、Z0分别为悬移质和推移质引 起的地形变形和总地形冲淤厚度; x、s3为垂线平均含沙量及水 流挟沙力;ω为泥沙沉速;ε为泥沙扩散系数, 取ε 0.15hU3, U3为摩阻流速;α为恢复饱和系数;ρs、ρb分别为悬移质和推移 质淤积物的干密度;ρ 、ρ ′ 分别为清水密度、 泥沙干密度; gbs、gby 分别为x、y方向的单宽推移质输沙率;U为x和y方向的合流 速;ks、m为水流挟沙力公式的系数和指数; kb为推移质输沙率 公式的经验系数,取kb 0.01;Uc为推移质起动流速, Ch1/ 6 。 1. 2 方程离散 本文在方程离散时空间上采用有限体积法,运用守恒格式 对水沙连续性方程进行离散,保证计算域内水量和沙量的守 恒,时间上采用蛙跳法可满足二阶精度的要求,对流项采用随 来流方向进行切换迎风方向的差分格式,扩散项采用中心差分 法,使用交错网格计算物理量的分布,如图1。模型详细离散过 程见式 9 - 式13。 图1 变量的空间布置图 水流连续方程离散后为 Zn3i1/2, j1/2-Zn1i1/2, j1/2 2Δt Mn2i1, j1/2- Mn2i, j1/2 Δx N n2 i1/2, j1-N n2 i1/2, j Δy 09 x方向水流运动方程离散后为 Mn2i, j1/2- Mni, j1/2 2Δt uni, j1/2Mni, j1/2- uni-1, j1/2Mni-1, j1/2 Δx uni, j1/2≥0, uni-1, j1/20 uni, j1/2Mni, j1/2 Δx uni, j1/20, uni-1, j1/2≤0 uni1, j1/2Mni1, j1/2- uni, j1/2Mni, j1/2 Δx uni, j1/2≤0, uni-1, j1/20, v ～n i-1, j1/2≤0 v ～n i, j3/2Mni, j3/2-v ～n i, j1/2Mni, j1/2 Δy v ～n i, j1/2≤0, v ～n i-1, j3/20 -v ～n i, j1/2Mni, j1/2 Δy v ～n i, j1/20, v ～n i-1, j3/2≥0 -g hn1i-1/2, j1/2 hn1i1/2, j1/2 2 Z n1 i1/2, j1/2-Zn1i-1/2, j1/2 Δx - g ni-1/2, j1/2 ni1/2, j1/2 2 2 2Mni, j1/2 hn1i-1/2, j1/2 hn1i1/2、j1/2 u n i, j1/2 2 v ～n i, j1/2 2 hn1i-1/2, j1/2 hn1i1/2, j1/2 2 1/3 10 151考虑矿砂淤积的尾矿坝溃坝洪水数值模拟 林 江 张绪进 张小峰 等 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved. 上式中 uni, j1/22Mni, j1/2/ hn1i-1/2, j1/2 hn1i1/2, j1/2 v ～n i, j1/2 vni-1/2, j vni-1/2, j1 vni1/2, j vni1/2, j1 / 4 y方向水流运动方程的离散方法与x方向的相同,离散结 果略。 泥沙连续方程离散后为 hs n3 i1/2, j1/2- hs n1 i1/2, j1/2 2Δt Mn2i1/2, j1/2sn1i1/2, j1/2- Mn2i, j1/2sn1i-1/2, j1/2 Δx N n2 i1/2, j1sn1i1/2, j1/2-N n2 i1/2, js n1 i1/2, j-1/2 Δy -α ω s n1 i1/2, j1/2- sn13i1/2, j1/2 CFFx- CFBx Δx CFFy- CFBy Δy 11 CFFxεh n1 i3/2, j1/2 hn1i1/2, j1/2 2 sn1i3/2, j1/2- sn1i1/2, j1/2 Δx CFBx、CFFy、CFBy相应交换x、y后给出。 悬移质冲淤引起的地形变形方程离散后为 ΔZbsi, jΔtα ωi, j s i, j- s3i,t / ρs12 推移质冲淤引起的地形变形方程离散后为 ΔZbgi, jΔt gbxi1/2, j1/2-gbxi-1/2, j1/2 Δx gbyi1/2, j1/2-gbyi1/2, j-1/2 Δy ρb13 2 模型中几个问题的处理 2. 1 计算区域 本文模型计算区域的选取为坝址及其下游地形区域。 2. 2 初始条件 模型中要求给定的初始条件包括水位、 流速。文中模型 采用下游为干河床的演进,溃口初始水位为未溃时水库水位 值,坝址下游水位等于地面高程,初始流速u 0, v 0。 2. 3 边界条件 计算区域边界分两种类型,即固体边界和开边界。其中 固体边界指周围山体、 高地等水流不能通过的边界,采用不可 穿入条件,即固体边界的法向流速为零,相邻单元交界面处无 流量通过;开边界指水流进入和流出的计算区域的边界,进口 处给定坝址溃口的流量过程和出沙过程,计算区域出口处,在 流动方向上各参数梯度变化为零。 渐变的坝体溃口其变形机理非常复杂,有推移、 悬移、 塌 落、 滚动等现象,目前还难以准确地描述其过程。本文仅以矩 形逐渐线性变化来模拟溃口渐溃的现象,矩形渐变,直到库容 基本泄空为止。 坝体逐渐溃决的流量计算采用宽顶堰溢流的流量计算 公式 Q σcmb2g H3/ 2 0 14 H0 H v20 2g 15 式中H0为行近流速水头的堰前水头;V0为行进流速;σc为侧 收缩系数;b为溃口宽度; m为流量系数,它与堰型、 堰高等边界 条件有关;σc、m等参数由水力学中相应表格和公式确定。 2. 4 计算参数选取 水流挟沙力公式的系数和指数ks取0.025, m取0.92。 泥沙冲淤恢复饱和系数冲刷时取α 1.0,淤积时α 0.25。 糙率根据下游地形条件而定,房屋、 村庄取0. 08 ,树丛取 0. 075 ,农田取0. 045。 空间步长对于山区性的尾矿坝,其地形起伏比较大,故 空间步长不宜过大。 时间步长溃坝洪水变化剧烈,时间步长过长会出现不符 合实际情况的数据结果。因以库朗特CFL条件为控制条件 Δt CtminΔtx,Δty Δtxmin i Δx | ui1/2, j1/2| ghi1/2, j1/2 , Δtymin j Δy | vi1/2, j1/2| ghi1/2, j1/2 式中Ct为Courant数,一般取小于1. 0的数;Δx ,Δy分别为 x , y方向的空间步长。 3 模型应用计算 云南某磷石膏堆场目前正在建设中位于云南滇池下游, 其下游2 km有海口河通过。场区为一狭长谷地,整个场区呈 东西走势,西高东低,为三面环山。场区内植被疏松,冲沟内农 田较多,谷口处为三汁箐水库。三汁箐水库下游地势平坦,地 面高程一般在1 900 m左右。场区内磷石膏库堆积坝坝高约 130 m ,库容约3 900万m3,下游有村庄和铁路线等。 采用上述模型对此尾矿库溃坝后的洪水淹没及矿砂淤积 影响进行模拟,结合场区附近的地形情况,模型模拟范围为坝 址下游约3 750 m3 080 m ,网格大小为ΔxΔy 10 m ,时间 步长Δt 0. 05 s ,矿砂为粒径为0. 1 mm的均匀沙 ,尾矿坝位置 为计算区域的入流口。 模型总共考虑了4种工况进行模拟计算,具体见表1。 表1 计算工况一览表 工况号 溃时水位/ m 泻砂总量/ 万m3 洪水总量/ 万m3 溃口最终尺寸 宽 高 / m 溃坝历时/ hour 12 0301 500500170605 22 070300200140305 32 070500300220605 42 0703 50010003001005 由于篇幅限制,本文取泻砂量最大值的工况4来叙述,工 况4通过渐溃的宽顶堰溢流公式计算坝址溃口处的流量过程 如图2。 251考虑矿砂淤积的尾矿坝溃坝洪水数值模拟 林 江 张绪进 张小峰 等 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved. 图2 坝址溃口流量过程图 图3为工况4尾矿坝溃坝后下游最终的矿砂淤积图,计算 结果显示当柳树菁磷石膏库溃坝后,洪水急速冲向下游的三 汁箐水库,由于洪水在水库内的滞留,导致大量的矿砂在三汁 箐水库内淤积,当入库水量超过三汁箐水库的库容后,洪水溢 过坝顶,冲入下游的平原地区,流速减小,矿砂淤积,进入海口 河的矿砂和洪水顺着河道流出计算区域。 图3 最终矿砂淤积图 图4、 图5分别为工况4代表点水深变化过程线和矿砂淤 图4 代表点水深变化过程线 图5 代表点矿砂淤积过程线 积过程线,可见,柳树菁磷石膏库一旦溃坝,其将严重危害下游 的铁路,沙锅村和达子小村。 4 结 语 本文通过对二维浅水方程和泥沙悬移质方程进行离散,建 立非恒定、 均匀不平衡悬移质泥沙数学模型,用于模拟尾矿坝 溃坝后稀性泥石流的水流矿砂运动及下游地形变形。实例计 算证明该模型可以正确模拟出尾矿坝溃坝后的洪水演进和矿 砂运动的基本规律,在实际应用中是可行且可靠的,为今后尾 矿坝的安全风险评估及防洪抢险决策的实施奠定了基础。□ 参考文献 [1] 丁言伟,于小川.尾矿库利用擅自动工万万不行[J ].当代矿工, 2002 ,9 13 - 15. [2] 梁雅丽. 10. 18南丹尾矿坝大坍塌[J ].沿海环境, 2000 ,12 7 - 7. [3] 张红平,高雪相.黄金尾矿溃坝之后[J ].中国气象报, 2006 ,5 1 1 - 2. [4] 陈春生,孙建华.矿山尾矿库溃坝砂流的计算模拟[J ].河海大学 学报,1995 ,235 99 - 105. [5] Zhang Xiao2feng , Yu Ming2hui. Numerical simulation of bed de2 ation dike burst [J ]. Journal of Hydrodynamics B ,2001 , 4 60 - 64. [6] 白玉川,许 栋,王玉琦,等.二维溃坝波遇障碍物的水流泥沙数 值模拟[J ].水利学报,2005 ,129 538 - 543. [7] 谢任之.溃坝水力学[M].济南山东科学技术出版社,1989. [8] 李义天,赵明登,曹志芳.河道平面二维水沙数学模型[M].北京 中国水利水电出版社,2002. [9] [英]迈克 剑桥.尾矿坝失事的回顾[J ].水利水电快报,2002 ,23 5 6 - 7. [10] 张新华,隆文非,谢和平等.二维浅水波模型在洪水淹没过程中 的模拟研究[J ].四川大学学报工程科学版 ,2006 ,381 20 - 25. 351考虑矿砂淤积的尾矿坝溃坝洪水数值模拟 林 江 张绪进 张小峰 等