尾矿库不同溃坝形式的数值模拟研究

林子扬, 王卫红, 杨宏珺

(1.西南科技大学环境与资源学院, 四川 绵阳 621010; 2.国家遥感中心绵阳科技城分部, 四川 绵阳 621010;3.西南科大四川天府新区创新研究院, 四川 成都 610299)

1 前言

我国境内存在1.2万余座尾矿库,其中大多数为上游式筑坝法构建的尾矿库。为了使每年尾矿排量能低成本处理,同时延长尾矿库的使用寿命,需要在现有尾矿库的基础上加高现有堆积坝[1]。从技术管理的问题来看,坝高越高,溃坝后影响范围越大,尾矿库的安全存在很大的潜在风险。

在尾矿库溃坝方面,国内外学者已开展系列研究。李全明等[2]建立尾矿库溃决泥浆运动模型,提出尾矿库溃坝范围计算的评估方法;刘嘉欣等[3]研究漫顶溃坝机理,通过沿深度取积分平均值建立溃坝的二维数学模型;崔旭阳等[4]构建了一种包含动态指标和静态指标的评价指标体系,利用动态加权贝叶斯网络和时间权重、指标权重来评估尾矿库溃坝风险;胡航等[5]强降雨洪水条件下下泄尾砂总量与溃坝后果之间的关系,并结合FLO-2D模拟软件对溃坝灾害进行了预测分析;吕松峰等[6]模拟了尾矿库在瞬时局部溃坝条件下尾砂的运动情况。Rico M 等对事故尾矿库的坝高、库容与溃坝量等资料进行回归分析,得到尾矿库溃坝量约占尾矿库总库容的1/3(包括水)。Concha Larrauri P等提出尾矿库在溃坝后不会全部泄完。

2001—2018年我国共发生95起尾矿库事故,事故中约380人丧生,其中溃坝事故占比高达56.9%[7-8]。因此,需要在加强尾矿库日常监管,降低事故隐患的同时,对尾矿库进行溃坝模拟实验,对实验结果和溃坝影响范围进行数据处理,得到应对尾矿库溃坝灾害时能提前对影响范围内的居民进行预警,保障尾矿库下游居民人身财产安全[9-10]。

综上所述,目前研究仅限于单一洪水溃坝,而实际上,在遭遇地震或其他自然灾害时,有可能发生单点或多点瞬时溃坝。因此,结合尾矿库洪水溃坝数值模拟,局部尾矿库瞬时溃坝数值模拟,提出尾矿库多点瞬时溃坝数值模拟,研究采用尾矿库安全评价项目,通过FLO-2D模拟尾矿库洪水漫顶,单点瞬时溃坝和多点瞬时溃坝,来研究尾矿库不同形式的溃坝在强度和影响范围的差异,得到单点瞬时溃坝破坏力强,影响范围略小于多点瞬时溃坝,但两种溃坝形式在强度和影响范围上都大过于洪水漫顶溃坝成果,弥补尾矿库溃坝实验形式单一的情况[11]。

2 研究区现状

尾矿库位于米易县垭口镇安宁河西岸烂坝山下的冲沟中,与选厂隔河相望,直线距离1 km,距撒莲镇政府5 km,距米易县政府20 km,具体位置如图1所示。

图1 尾矿库位置图

该尾矿库修筑方式为上游式筑坝法,2022年之前该尾矿库存在3座坝体,经过不断堆高修筑,现如今合并为一座坝体的尾矿库,设计总坝高为246 m,总库容设计为5 079.5万m3,有效库容为4 800万m3,该尾矿库由三等库升级为二等库,位于山谷中,初期坝下游是安宁河,有少量耕地,无工矿企业、水产基地及名胜古迹等重要设施。

库区属热带、亚热带森林气候,季节干湿分明,夏天炎热、多雨、潮湿,而冬天温暖、多雨、干燥。年平均气温21～22 ℃,最高38.7 ℃,雨季期间,陡峭的河谷两侧频繁发生山洪、河流暴涨和小型泥石流。

3 数据和方法

尾矿库初期设计为沟口处修建3座坝体,随着坝体升高,将3座坝体合并为1座坝体,因此选取该尾矿库呈现出单点和多点溃坝模拟的成果会有明显。

选取的不同溃坝条件分别为洪水漫顶溃坝数值模拟;瞬时单点溃坝数值模拟;瞬时多点溃坝数值模拟。其中洪水漫顶溃坝数值是以某断面流体流量随时间变化的过程来客观反映泥石流灾害整体运动过程。瞬时溃坝数值模拟相较于洪水漫顶,主要在于其突发性和不可预见性,为了保障下游地区人民人身财产安全,需提前做好此类事故预测工作,为此选取瞬时单点溃坝。在此基础上,由于尾矿库坝顶宽度过长,在突发情况下不止一处会发生崩坏等现象,因此在实验中选取多点瞬时溃坝来模拟该情景下,下游受灾情况。

在实验软件上,采用FLO-2D软件实现尾矿库溃坝数值模拟实验,其具有对泥石流有较好的重现性,能表示泥深、流速,以及反映流体流变时空特征等优点[12]。该软件通过分析地形数字高程模型为实验蓝本,将尾矿库区域的地形制作成二维单元网格形式,计算得出每个网格单元在实验中对应的流速,泥沙淹没深度以及流量数据。

3.1 模拟流体参数设定

本次模拟是基于体积守恒和水力学原理构建流体质量守恒和动量守恒方程的微分形式,采用显式中心差分法对守恒方程进行求解运算,可以有效计算单元网格内流体流深、流速和影响范围的时间变化。

鉴于模拟模型原理,模拟流体以一种恒定的固液混合流体形式进行运动堆积,其需要设定以下流体参数:体积浓度折中取值0.5[13]。泥石流重度取值1.25 kg/m3。层流阻力系数本次模拟参照Woolhiser 实地测定的泥石流层流阻力系数设定该值为2 200[14]。屈服应力和黏滞系数:流体的屈服应力选择詹前登实验参数。

3.2 地表环境参数设定

在本次模拟中,根据高分辨率航空图像,将数值模拟区域分类为地物识别,并根据地物类别估算渗透性和阻碍能力。通过曼宁系数描述地表阻碍能力,使用Green Ampt透射模型描述表面渗透能力[15]。

4 数值模拟

4.1 洪水漫顶溃坝数值模拟流量曲线

参照防排洪验算和尾矿库建库标准,以1000年一遇的洪水频率为尾矿库溃决计算频率,依照水文计算洪水流量和模拟流体参数计算库区断面的泥石流流量过程线如图2所示。

图2 泥石流流量过程线示意图

图2中1#,2#,3#代表3座坝体,以最北方向为1#坝体依次向南方排列。洪水流量过程线表示尾矿库由于承受不住洪水容量溃坝后,洪水向坝体下游流动的过程;泥石流流量过程线表示尾矿库由于洪水漫顶后产生溃坝,洪水夹带着尾砂向坝体下游流动的过程。

4.2 洪水漫顶溃坝数值模拟成果

基于相关资料与参数,开展尾矿库洪水漫顶溃坝数值模拟,其计算结果如图3所示。

图3 尾矿库洪水漫顶数值模拟结果示意

根据数据模拟成果分析,在1000年一遇的洪水频率下,尾矿库一旦发生溃决对尾矿库下游地表环境造成严重破坏。尾矿库坝体涉及范围广,下游共有3处流域通往安宁河。依据本次数值模拟,模拟灾害流体分别从坝体上的3处位置运动至下游沟口,模拟流体平均速度约为0.3 m/s,模拟体过流的最大厚度集中在2 m以下,占灾害影响面积的82.6%,最大沉积厚度在4 m以上的区域占总沉积面积的3%,泥深最高沉积厚度为7.3 m。

4.3 瞬时单点溃坝数值模拟泥石流流量计算

1)溃口宽度计算

溃口宽度是尾矿库溃坝灾害破坏程度的主要影响因子之一,同时也是计算下泄流量的重要参数。通过水库溃决资料得到计算溃口宽度的计算公式[16]。

(1)

式中:b——溃口平均宽度,m;

W——尾矿库溃坝总量,万m3;

B——尾矿库水面宽度;

Hw——尾矿库水面高,m;

K1——筑坝材料相关系数,如黏土材质取0.65,土壤材质取1.3。

2)溃坝总量计算

尾矿库灾害的关键取决于溃坝总量的分析,对于计算尾矿库溃坝总量在世界上没有统一的标准,在数值模拟中将尾矿库全库容作为溃坝总量数据,但在实际溃坝案例中,溃坝总量总是小于尾矿库库容,只存在部分库容通过溃坝泄出。不同尾矿库溃坝方式存在诸多不同的影响因素,因此尾矿库溃坝总量的计算存在差异性和不唯一性。根据Rico M等的研究,本次实验选取溃坝量为总库容的1/3。对此次溃坝模拟采用的溃坝总量为5 034万m3,选取有效库容为1 676.4万m3,占尾矿库总库容的33.3%。

3)坝址处溃坝流量过程线计算

通过对水库溃坝计算方法比对,主要是采用圣维南法,经验公式等方法。对于尾矿库溃坝采用借鉴水库溃坝计算方法[17],其公式为

(2)

式中:QM——溃坝流沙的最大流量,m3/s;

m——不同形式的沟谷断面参数;

H0——尾矿库上有水深。

对于不同形式的沟谷断面,可选取相对应或组合模型进行取值,具体参数见表1。由于尾矿库坝体形状介于二次抛物线和组合抛物线之间,m选取为5/2。

表1 不同沟谷断面形状的特征参数

瞬时溃坝条件下,在坝体发生突然破坏后,坝址处下泄尾砂流量在开始时达到最大值,随后快速下降,当到达临界值时,缓慢下降,对溃坝流量采用4次抛物线来进行概化。此次溃坝流量4次抛物线成果见表2。

表2 溃坝流量4次抛物线概化参数

表2中,T表示为溃坝砂流下泄时长,t表示为溃坝下泄时间刻度;QM表示为溃坝流砂的最大流量;Qt表示为t时刻砂流流量。

绘制溃坝尾砂流流量过程曲线,主要参考溃坝中水力学经验公式,采用以下计算参数:

(1)溃坝总量W的计算。采用总库容溃坝1/3,溃坝总量取值为1 676.4万m3。

(3)尾矿库溃坝砂流最大流量QM。

4.4 瞬时单点溃坝数值模拟流量曲线

参照水文电力手册和尾矿库建库标准,以尾矿库单点溃坝计算频率[18],依照模拟流体参数计算库区断面的尾砂流量过程线如图4所示。

图4 瞬时单点溃坝流量过程曲线图

4.5 瞬时单点溃坝数值模拟成果

基于相关资料与参数,开展尾矿库瞬时单点数值模拟,其计算结果如图5所示。

图5 瞬时单点溃坝模拟成果图

依据本次数值模拟,模拟灾害流体分别从坝体上的3处位置运动至下游沟口,模拟流体平均速度约为17.7 m/s,由于本次模拟有大量尾砂冲入安宁河内,选取模拟最大值位于地面部分,模拟体过流的最大厚度集中在6.5 m以下,占灾害影响面积的76.8%,最大沉积厚度在6.5 m以上的区域占总沉积面积的23.2%,尾砂最高沉积厚度为23.21 m。

4.6 瞬时多点溃坝数值模拟泥石流流量过程线

瞬时溃坝条件下,在坝体发生突然破坏后,坝址处下泄尾砂流量在开始时达到最大值,开始快速下降,当到达临界值时,缓慢下降。

参照水文电力手册和尾矿库建库标准,以尾矿库3条沟口坝体坝长比作为权重参数,本次实验中权重比为1#∶2#∶3#= 0.351 314 85∶0.326 197 28∶0.322 487 871,将溃口设置在三座坝体上计算频率,依照模拟流体参数计算库区断面的尾砂流量过程线如图6所示。

图6 瞬时多点溃坝流量过程曲线图

4.7 瞬时多点溃坝数值模拟成果

基于相关资料与参数,开展尾矿库多点瞬时溃坝数值模拟,其计算结果如图7所示。

图7 瞬时多点溃坝模拟成果图

依据本次数值模拟,模拟灾害流体分别从坝体上的3处位置运动至下游沟口,模拟流体1#坝平均速度约为2.18 m/s,模拟流体2#坝平均速度约为2.17 m/s,模拟流体3#坝平均速度约为2.14 m/s,由于本次模拟有大量尾砂冲入安宁河内,选取模拟最大值位于地面部分,模拟尾砂的最大厚度集中在5.4 m以下,占灾害影响面积的80.6%,6.5 m以上最大冲积厚度区域占整个灾害影响面积的19.4%,尾砂最高沉积厚度为23.56 m。

5 下游威胁分析对比

依据洪水漫顶溃坝数值模拟结果,模拟灾害流体沿着尾矿库下游沟道向下运动,并在流域沟口处进行了扩展堆积,形成扇形的破坏区,从而对沟口的居民建筑、农田耕地和交通干道造成了大面积的淤埋,甚至有流体向下运动堆积至安宁河流域造成河道堵塞,引发次生灾害。现选取3处下游流域沟口区域进行模拟破坏分析,下游影响区:该沟口区域距离坝顶约1.2 km,沟口区域内模拟流体横向波及长度达220 m,约15 984 m2耕地被模拟流体淤埋,1栋居民建筑边缘被流体波及,1栋疑似临时建筑的房屋被流体淤埋破坏,但灾害流体并未继续向下堆积至安宁河道之中。

依据尾矿库瞬时溃坝数值模拟结果,溃坝的尾砂冲过安宁河达到对岸,瞬时溃坝最远处距离坝顶约为2.3 km,沟口区域内模拟流体横向波及长度达1 506 m,约有3.35×105m2农田耕地被尾砂覆盖,其中23栋居民建筑边缘被流体波及,4栋建筑的房屋被流体淤埋破坏。其中单点瞬时溃坝对溃坝点处的坝体造成严重影响,对另外两座坝体影响较小,而多点瞬时溃坝分别对坝体造成的影响相差不大,因此单点瞬时溃坝对2#坝体下的房屋耕地受灾最为严重。

6 结论

本文利用对四川省攀枝花市米易县某尾矿库DSM影像进行了处理,再通过FLO-2D软件对尾矿库进行多种溃坝数值模拟实验,得到以下结论:

(1)不同溃坝条件下的溃坝流量不同,瞬时单点溃坝的溃口流量速度为17.7 m/s,短时间内对下游地区的冲刷破坏性较强,瞬时多点溃坝的溃口流量相对有所降低,速度为2.17 m/s,具有持续性破坏作用,洪水漫顶造成的溃坝灾害流量为0.3 m/s最小,造成的破坏性最小,但由于伴随着大量降雨,溃坝灾害的同时伴随着洪水等次生灾害。

(2)不同溃坝条件下的溃坝破坏强度不同,瞬时溃坝模拟结果中,单点瞬时溃坝对对应坝体下游受灾地区破坏力最强,其余坝体下游受灾相较于多点溃坝较弱,多点溃坝灾害受灾程度较均匀,受灾面积扩大,但两种溃坝形式下受灾区域大致相同。

(3)不同溃坝条件下的尾砂沉积方式也不同,洪水漫顶溃坝中尾砂堆积主要在初期坝下游地区;单点瞬时溃坝尾砂主要堆积于主体受灾坝中后段,堵塞初期坝;多点瞬时溃坝尾砂堆积于溃坝口和初期坝位置,在河道处沉积的尾砂也是最多的。

(4)在三种溃坝模拟中,洪水漫顶溃坝破坏力最小,尾砂沉积量最少;单点溃坝模拟,破坏力最强,尾砂主要集中于主坝体口;多点溃坝模拟,溃坝破坏力较分散,破坏力由两边向中间集中,尾砂同时也主要集中于坝体中部,对受灾区域冲击面积最广。

(5)尾矿库溃坝后,尾砂大多流入安宁河中,使其破坏程度得到有效削减,减少了溃坝灾害对城镇的直接危害。