云南某矿山磷石膏尾矿库溃坝灾害影响分析

2023-02-22陈俊智马毓婷陈文倩

有色金属（矿山部分） 2023年1期

陈祯，陈俊智，马毓婷，安南，陈文倩

(昆明理工大学国土资源工程学院，昆明 650032)

尾矿库溃坝事故是影响较大的工程灾害事故，在矿山的安全生产中备受矿山企业关注与重视。在矿山生产中由于选矿后产生的尾砂材料需要长期积压堆存。当溃坝发生时库内的尾砂与水混合后会以类似“泥石流”的形式下涌，严重危害所涉及区域的人员以及财产安全[1]。此外，若灾害发生，尾砂中暂时不能处理的有害物质泄露会造成周围环境的污染[2]。且我国现存的上万座尾矿库大约有一成存在各种问题[3]。因此，加强对尾矿库溃坝事故影响范围的研究并建立相关数值模型，能够有效减少溃坝灾害发生时造成的生命以及财产损失，对保障尾矿库的安全运行具有重大意义。

迄今为止，就溃坝事故的产生机理以及危险系数和影响波及的范围来看，国内外学者已经做了较多成果显著的研究工作。BLIGHT[4]通过对环形尾矿坝的溃坝分析得出尾矿所处地表的湿度情况会影响尾矿浆的扩散范围，地表湿度较高的地方尾矿浆的移动扩散范围较远；PASTER等[5]以弹塑性有限元法为原理建立相关模型，该模型可以用在滑坡、泥石流、溃坝等灾害的模拟中推算流体流动随时间堆积的程度；PIRULLI等[6]通过对比分析碎屑流与溃坝泥流，探究两者共性，并在尾矿库溃坝事故分析中运用RASH 3D软件提出一个流变模型。

国内对溃坝等灾害现象的研究成果颇丰，陈青生、孙建华等[7]以建立数学模型的方式找出了溃坝现象发生后引起的泥浆流动的表现规律；周彪[8]通过分析研究不同工况对渗流场的影响，并结合FLAC3D软件分析坝体的渗流稳定性；余乐文等[9]利用FLO 2D软件对溃坝现象模拟分析，研究了溃坝后泥石流的堆积情况，并且预测了溃坝发生后的波及范围；金佳旭等[10]利用ANSYS CFX软件分析辽宁某矿山的尾矿库溃坝事件，得出了溃坝所需时间、影响范围和尾砂堆积范围以及高度等情况。

据上，国内外学者主要考虑尾矿溃坝泥浆动力特征、渗流和溃坝后的影响范围及尾砂堆积角度采用了相关理论分析和数值模拟方法分析尾矿坝稳定性[11]，并利用基于岩土工程及边坡工程成熟的理论和工程经验与相应可靠的分析软件进行数值模拟仿真。

因此，本文以云南某矿山磷石膏新建尾矿库为例，根据现场测量参数，应用相关公式设计计算三维模型大小，运用Gambit软件建立实际尺寸的三维模型并划分输入网格，运用ANSYS FLUENT软件进行溃坝数值模拟。本文着重探讨了溃坝后的影响范围，为建库安全评价和灾害预警防范提供有力依据及参考。

1 磷石膏尾矿库工程概况分析

云南某矿山磷石膏尾矿库库区位于山顶冲沟处，微地貌为山间沟谷、山脚坡地，场地东西两面为山，沟谷近东西向，东高西低。从库区到外界的交通比较方便，有建设好的公路。库区地表分水岭为近似圆弧形连绵山脉，南端最高山峰高程约为2 004 m；以沟谷为中心线，形成基本对称的两侧岸，坡度角约为10°～20°，表面有植被覆盖。但当尾矿库发生溃坝时，由于其位于山顶冲沟处，其库浆泥沙有极大的重力势能。所以，发生溃坝时库浆将在极短的时间内以极快的速度对下游造成毁灭性的破坏。

1.1 尾矿性质

尾矿以粉粒为主，塑性指数IP为7.8，小于10；按照《尾矿堆积坝岩土工程技术规范》(GB50547-2010)中粒径组类的划分标准：拟堆积的磷石膏尾矿为尾粉土。并且磷石膏尾矿渗透系数为2.52×10-5cm/s，渗透性等级属弱透水，尾矿pH值介于1.95～3.04，浸出液pH值介于1.46～6.4，属酸性。

1.2 尾矿库库容设计概况

为增加库容和方便防渗工作展开，云南某矿山企业将堆填区域内的库底及库岸进行清理平整、削坡。清理完成后场区与初期坝、库尾挡水拦渣坝一同围成一盆状：库底平缓，纵坡度不超过10%，库岸坡度不陡于 1∶2。尾矿库的计算容积见表1。

表1 尾矿库容积计算表Table 1 Tailings pond volume calculation table

1.3 尾矿库坝体设计概况

初期坝顶标高1 945.0 m，坝轴处自然地面标高1 922.0 m，沟谷段清基深度为4 m，最大坝高H初=23 m(不含清基深度)，坝顶宽B=4.0 m，坝体外坡比1∶0.5，内坡比1∶0.5，坝轴长104.5 m。外坡设置一条1.5 m宽的道路，具体设置在标高1 930.0 m处，内坡在紧挨坡面处设置堆石排水棱体，排水棱体顶面标高为1 935.0 m，顶部宽度为2.0 m，上游坡比1∶1.5，坡面上布置反滤层，可在一定程度上对上游流下来的水起引流作用，保护坡面。石坝体与岸坡的连接方式为斜面连接，这样可以避免坝体与岸坡的连接面产生集中渗流现象，也可以减少软弱夹层的形成。

由于堆石坝上下游坡比较缓，坝底纵向长度较大，清基工程量和筑坝石方量均较大，此外堆石坝方量较大，占用库容空间也较大，使得设计库容较小，无法满足设计库容要求。受制于尾矿库径深较短，因此不采用堆石坝设计。

2 溃坝尾矿浆流动过程

尾矿库溃坝后库浆下泄引起的浆流本质上属于库浆与水的混合物，类似泥石流，因此，本文的基本假定为库浆砂是各向同性的连续均匀的介质体，等同于高密度的水体。

2.1 尾矿库溃坝下泄总量

多日连续降雨会导致尾矿库内积水，雨水与尾砂形成浆体，结合库浆物理力学性质，为保证设计安全，考虑溃坝发生时产生的溃口相应高程处的全部库容为该尾矿库的下泄总量，具体库容如表1。假设尾矿库坝体全部容量泄出，取坝体全部容量为100万m3。

2.2 尾矿库溃口宽度

尾矿库的溃口宽度参考黄河水利委员会提出的溃口宽度计算公式。

溃口宽度计算公式见式(1)。

(1)

式中，B0为尾矿库溃口平均宽度，m；K1为坝体材料系数，尾矿库坝体为石质，取0.65；W为溃坝时的下泄总量，m3；B为溃坝发生时坝前的水面宽度，m；H为尾矿库溃坝时的水深，m。

2.3 尾矿库溃口的最大流量

假设溃坝时溃口在一瞬间完全形成，即瞬时全溃，不考虑溃口形成过程对最大流量的影响。本文采用谢任之提出的计算方法，这一方法在对尾矿库类型的溃坝中运用较为广泛，溃口最大流量计算公式见式(2)[12]。

(2)

式中，QW为溃口最大流量，m3/s；λ为流量系数，由断面形状系数、沉溺系数、堰宽比等系数计算得出，瞬时全溃时取8/27；g为重力加速度，m/s2；H0为坝前平均水深，m。

2.4 库浆下泄时间

溃坝发生时，溃口底部所在水平以上的所有库容对应的库浆会全部泄出，下泄时间受制于库容量和溃口高度。由于假设库浆全部泄出因此溃口位于尾矿库坝体底部。

库浆下泄时间计算公式见式(3)。

(3)

式中，W为下泄总量，m3；Qm为尾矿坝溃口最大流量，m3/s；K2为修正系数，取3.5。

2.5 库浆最大流速

溃坝发生前往往多降雨，尾矿库内积存大量雨水。雨水与尾矿混合后会形成饱和浆体。溃坝发生后，浆体从溃口流出，其流变性质类似于稀性泥石流[13-14]。

库浆最大流速计算公式见式(4)。

(4)

式中，Mm为河床粗糙系数；HC为尾矿浆的平均泥深，m；IC为尾矿库的纵坡降，%。α为阻力系数，按式(5)计算。

(5)

式中，φC为泥石流修正系数；γg为尾矿堆积干，t/m3。

2.6 库浆最大冲击力

库浆最大冲击力计算公式见式(6)。

(6)

式中，γ为尾矿浆的饱和容重，kg/m3；Vmax为尾矿浆最大流速，m/s；g为重力加速度，m/s2；K3为泥浆不均匀系数，取3。

2.7 尾矿库溃坝理论危害面积

尾矿库溃坝后，库浆从溃口流出并扩散，扩散后的溃坝库浆危害面积公式见式(7)。

(7)

式中，h为尾矿浆的最大堆积厚度，m。

根据以上公式计算得出的仅是理论溃坝参数，但是在实际情况中，尾矿的堆积厚度受制于地形影响，无法达到设计厚度，因此在计算时，实际值取理论值的85%，结合理论计算公式与尾矿库工程实际计算，得到所有参数的最终计算结果见表2。

表2 理论溃坝参数和实际溃坝参数Table 2 Theoretical and practical dam break parameters

3 溃坝数值模拟

由于各类型尾矿库的地理环境差异巨大且复杂多变，库容物质也不尽相同，因此，想建立统一模型进行理论分析的可能性微乎其微。现如今，随着科技的发展，运用ANSYS FLUENT数值模拟手段能够广泛而可靠的适用不同库型、不同条件、各种工况下的溃坝泥沙演进过程，进而为理论研究提供更加行之有效的参考数据。

3.1 三维建模

由于堆积子坝的库容最大，因此研究堆积子坝溃坝灾害对坝体下游空间造成的影响最具有代表性。使用Space Claim软件建立500 m×500 m×200 m的尾矿库下游空间简化三维模型和简化溃口三维模型105 m×100 m×35 m(其中100 m为坝体宽度)。将溃口放置在下游空间模型左侧，坝体右侧地势较左侧高，溃坝时，泥浆由于地势的阻挡，会向左下侧流动。假设溃坝事故发生时溃口在一瞬间完全溃决，溃口高度为35 m，为整个堆积子坝高度；溃口长度为整个坝轴线长，即105 m；计算区域内共有908 170 m3下泄库浆。溃坝计算模型及尺寸示意图见图1。

图1 溃坝计算模型及尺寸示意图Fig.1 Schematic diagram of dam break calculation model and dimensions

3.2 网格划分

使用ICEM CFD软件对三维模型划分网格。设置简化溃口模型一侧为入口，下游空间简化模型一侧为出口，其余模型边界设置为壁面，如图2所示。对整个模型进行结构化网格划分，网格为正六面体，网格大小为：Δx×Δy×Δz=1 m×1 m×1 m，将网格导入ANSYS FLUENT后，通过模拟计算直观得出库浆体积量在不同时段对网格的渲染程度来判断溃坝后的大致淹没区域。划分网格数量为50 367 500个，如图2所示。

图2 ICEM CFD网格划分细节Fig.2 Details of ICEM CFD meshing

3.3 模拟方法及边界条件设定

模拟求解器设置：

1)根据模型，设定重力作用的方向为z轴负方向，重力加速度大小取值为-10 N/kg；

2)库浆在运动过程中为不可压，因此在求解中选择Pressure-Based求解器；

3)库浆是非定常流的流体，因此设定瞬态计算，时间步长为0.01 s，模拟10 000步，模拟时长100 s，大于下泄时间78 s，可以观察整个溃坝事故的库浆完整运动情况；

4)库浆流态为湍流流动模型设置为标准k-ε湍流模型，该模型适合模拟流体的扩散情况，符合库浆流出后扩散的实际情况；

5)根据现场采集的磷石膏尾砂配置的库浆，其密度在1 674.7 kg/m3，黏度为0.43 Pa·s，因此设定材料参数时，以液态水为基础修改密度和黏度进行模拟；

6)求解方法设置为PISO，该求解方法适用于非定常计算；

7)设置残差为1×105，已知溃口处的速度条件，因此选择标准初始化中的入口初始化。

边界条件设置：

1)设置简化溃口模型一侧为速度入口，速度取计算值；

2)下游空间简化模型一侧为唯一自由流出口，权重为1；

3)其余模型边界设置为静止壁面边界，无滑移情况；

4)壁面边界中，壁面的粗糙高度和粗糙系数根据壁面材料设定，由于模拟的壁面为粗糙的岩土质地面，因此取粗糙高度为0.1，取粗糙高度为0.5；

5)下游计算区域左右两侧均为封闭边界。

边界条件示意图如图3所示。

图3 模型边界条件示意图Fig.3 Schematic diagram of model boundary conditions

3.4 模拟结果及分析

模拟结束后，通过生成流动动画，观测尾矿坝溃决后不同时间的库浆淹没范围图(图4)。

图4 尾矿坝溃决后不同时间的库浆淹没范围图Fig.4 Map of slurry inundation range at different times after tailings dam break

由模拟结果的不同时刻截图可知：1)在溃坝初期，库浆主要呈柱状涌出并开始向两边扩散，此时溃口库浆流速较快。2)随着淹没面积的形状和趋势可以看出，由于溃口靠近下游区域左侧，尾矿浆从坝址处决堤后先接触左侧的封闭边界，而后慢慢向下方和右侧扩散，但主要沿着左侧边界朝着溃口正对方向逐渐进入下游区域，大约在50 s内到达模拟下边界，与此同时计算区域右半部分的淹没面积显著增大。3)结合模拟可以发现，自50 s开始，涌入右半部分的泥浆开始加速蔓延扩散，然后泥浆从整个下游边界出口流出。由于库内泥浆下泄量巨大，溃坝发生后50 s时，库浆就能到达溃口下流500 m左右，这表明仅50 s库浆就能纵穿整个下游区域。4)到80 s时，整个下游区域几乎已经完全被淹没。