基于PFC3D磷石膏尾矿库溃坝数值模拟研究.pdf

工业安全与环保24Industrial Safety and Environmental Protection2023年第49 卷第8 期August 2023基于PFC3D磷石膏尾矿库溃坝数值模拟研究朱远乐1,2王淇萱1,2廖文景1.2赵新勇3（1.长沙矿山研究院有限责任公司，湖南长沙41 0 0 1 2;2.金属矿山安全技术国家重点实验室，湖南长沙41 0 0 1 2；3.广西华锡矿业有限公司铜坑矿业分公司，广西河池5 47 2 0 4)摘要为探究尾矿库溃坝灾害作用对下游区域的影响，利用GOCAD/Rhino软件对某磷石膏尾矿库地形进行建模，采用颗粒流离散元方法研究设计1 40 m堆积高度尾矿库在4个不同方向上溃决砂流演进过程，分析了尾砂的位移和流速随溃坝时间的变化规律，为尾矿库下游合理规划、企业应急处置和科学防控提供依据。结果表明：1 40 m堆积高度初期坝溃决工况为溃坝影响最大、下泄体量最大的工况，下泄体最大宽度42 8 m，最大深度1 5 m左右，体积达2.0 1 7 1 0 m，尾砂最终滑移距离为1 2 2 9 m，途中淹没回水泵房及多处房屋，对下游居民造成严重影响；1 40 m堆积高度北部、西部堆积坝溃决对下游影响较小；1 40 m堆积高度东部堆积坝溃决，未影响任何下游建筑物。关键词尾矿库溃坝模拟PFC3D离散元方法Numerical simulation of phosphogypsum tailings dam failure based on PFC3DZHU Yuanlel-2 WANG Qixuan*-2 LIAO Wenjing*-ZHAO Xinyong(1.Changsha Mining Research Institute Co.,Ltd.,Changsha Hunan 410012,China)Abstract In order to explore the impact of tailings pond dam break disaster on the downstream area,the terrain of aphosphogypsum tailings pond was modeled by GOCAD/Rhino software,and the particle flow discrete element methodwas used to study the design of tailings ponds with a stacking height of 140 m in four different directions.The evolutionprocess of the burst sand flow is analyzed,and the variation law of the displacement and flow velocity of the tailingswith the dam burst time is analyzed,which provides a basis for the rational planning of the downstream of the tailingspond,emergency disposal of enterprises and scientific prevention and control.The results show that the dam failurecondition at the initial accumulation height of 140 m is the condition with the greatest impact of dam failure and thelargest discharge volume.The discharge body has a maximum width of 428 m,a maximum depth of about 15 m,a vol-ume of 2.01710 m,and a final sliding distance of tailings of 1 229 m,the backwater pump house and many houseswere submerged on the way,which had a serious impact on the downstream residents;the collapse of the northern andwestern accumulation dams with the 140 m accumulation height had little impact on the downstream;the 140 m ac-cumulation height of the eastern accumulation dam collapsed and did not affect any downstream buildings.Key words tailings pond dam break simulation PFC3D discrete element method0引言近年来，随着磷化工企业新增产能的急剧扩大，磷固废总量也随之增加。在当前磷化工企业磷石膏利用率总体仍低于总堆存2 0%的情况下，磷石膏的处理方式仍然以堆存为主，因此磷矿渣尾矿库仍是保证磷化工企业正常生产运行的重要设施2 。而尾矿库作为一个具有高势能的人造泥石流库，是一个潜在的污染源、危险源。尾矿库溃坝时泥沙流速大，从坝脚到下游1 km处往往只需几分钟，应急时间非常短3。尾矿库溃坝造成的事故和危害十分巨大，尤其在我国的矿山企业中存在着大量的头顶库，*基金项目：湖南省应急管理厅科技专项（2 0 2 0 YJ002）；金属矿山安全技术国家重点实验室开放基金（CKY/YJ202101）。尾矿库的下游通常分布着居民住宅和交通等设施，一旦发生溃坝事故，将严重威胁到下游村民的人身和财产安全4。2 0 0 7 年，辽宁省鞍山市海城某选矿厂5 号尾矿坝发生溃坝，约5 4万m尾矿下泄，造成其下游约2 km处的部分房屋被冲毁，1 6 人死亡、39人受伤5 ;2 0 0 8 年，山西省襄汾铁矿尾矿坝发生溃坝事故，约2 6.8 万m尾矿下泄，造成2 5 8 人死亡、34人受伤0 。目前，国内外学者对于尾矿库溃坝方面的研究，主要有理论计算法、模型试验法和数值模拟法等，其中数值模拟法以其高效率、低成本、条件设置灵活等25优势应用最广。WANGK等 7 利用无人机摄影测量式中,F为线性力；F为阻尼力；M为滚动阻力矩。和光滑无网格粒子SPH三维数值模拟方法，模拟溃线性力和阻尼力的更新与线性模型相同，而滚坝泥浆在下游真实地形上的演进过程，并揭示了坝动阻力矩则通过以下步骤进行更新。高、溃口宽度、泥浆黏度影响因素对于溃坝演进后果首先，滚动阻力力矩增加为：的影响规律。罗昌泰等 8 利用RAMMS碎屑流软件M:=M-k,0 b建立了数值三维模型进行数值仿真溃坝模拟，并建式中，,为相对弯曲旋转增量，k,为颗粒转动刚度，立了相应的1：1 50 物理模型，综合研究了溃坝的尾与切向接触刚度k,相关。砂流对下游的演变进程和影响范围。梁萱等 9 利用另外，滚动阻力矩没有扭曲分量。可以在接触FLOW-3D软件建立了尾矿库及其下游3km范围内平面坐标系中表示为：的实体地形及房屋和村庄三维数值模型，对尾矿库逐M=MbsS.+Mb,t渐溃坝的动态过程进行了模拟计算，分析了溃坝对下滚动阻力刚度定义为：游的影响。孙银华等 1 0 1 利用FLUENT数值模拟软件ki=k,R?建立了三维尾矿坝溃坝模型，并对尾矿坝下游3个式中,R为接触有效半径，定义为：不同位置6 个不同时间段的冲击压力、掩埋高度、流11+1动速度进行了分析。张家勇等 1 以贵州省开阳县鱼RR+RO式中,R)、R(2)分别为触点的端部1 和端部2 的半径鳅坡滑坡为研究对象，采用颗粒流离散元PFC3D对(用于球面触点）。其破坏运动过程进行数值模拟，对滑坡不同关键部然后根据阈值极限检查更新的滚动阻力力矩的位颗粒进行位移、速度监测，阐明其破坏运动特征。大小：本研究以某磷石膏尾矿库为工程背景，借助GO-CAD/Rhino软件对坝体及其下游2 km范围内的地形及重要建筑物进行建模，并利用PFC3D三维离散元模拟软件对磷石膏尾矿库溃坝后尾矿的演进过程进行模拟分析。1颗粒离散元方法颗粒流的计算方法是基于离散单元法（DEM），离散元方法由PeterCundall于1 9 7 1 年提出的 1 2 1,经过多年时间的发展以及后续学者不断的完善，离散单元法已成为研究散体介质的一种强有力工具。PFC颗粒流离散元软件因其以刚性球体为计算基本单元，能够很好的符合尾矿矿渣的力学特性，所以用PFC软件来模拟滑坡和尾矿库的溃坝能够很好地展现颗粒从坡上滑落的全过程。接触模型是PFC(ParticleFlowCode)计算理论的核心问题，其中PFC5.0中的抗滚动接触模型有以下优点：采用圆形颗粒，且避免对不规则颗粒间接触的反复判断，可提高计算效率；引入具有明确物理意义的形状参数，表示颗粒间的接触情况；引入滚动阻力力矩并与滚动刚度联系。另外，由于尾砂颗粒具有离散特性及不规则特性，考虑到颗粒之间存在的相互作用(咬合、摩擦等)对颗粒滚动产生的阻碍。因此，溃坝颗粒离散元模拟的接触模型考虑选用抗滚动接触模型。抗滚动接触模型的接触力F和力矩M通过力-位移定律更新：F=F+F,MC=M(2)(3)(4)(5)M,IMIM*M=MMMT),/MIM*极限转矩M*定义为：M=u,RF,式中，,为滚动摩擦因数；F为颗粒间的法向接触力。滚动滑移状态s，也更新为：ture,IlMII=M*S,=false,其他2工程概况某磷石膏尾矿库为典型的“头顶库”（下游1 km距离内有居民或重要设施的尾矿库），初期坝采用碾压土石坝，坝顶高程1 0 0.0 m，坝底高程7 7.0 m，坝高23.0m，顶宽5.0 m，坝轴线长度约为40 3.6 m，上、下游每级坡比均为1：2.5，堆积坝采用上游胶结法磷石膏筑坝，平均堆积坡比为1：3,设计最终堆积标高140.0m，总坝高6 4.0 m，总库容为4.2 2 51 0 m，属于三等库。磷石膏库区周边环境如图1 所示，库区下游东南侧3.1 km为长江，场区东侧约0.3km为老渣库，场区西侧约0.6 km有水库，西侧1.8 km为S225省道，均有天然山脊分隔，北部距离3.7 km为31 8 国道,场区东侧2 km为姚港三路，东偏北方向约3km为化工厂，库区下游0.1 1 km范围内有大量居民。3溃坝数值模拟3.1溃决工况(1)为分析溃坝对下游居民或重要设施的影响，本(6)(7)(8)26研究将计算在洪水情况下尾矿坝发生全溃瞬时溃坝的灾害情况，同时考虑设计最终1 40 m堆积高度初期坝溃决、1 40 m堆积高度东、北、西部堆积坝溃决4种工况。积高度初期坝溃坝范围情况，通过数值模拟可以得出不同时刻磷石膏尾矿库下游尾砂位移和流速随溃坝时间的分布规律。汇水区施区边界调节水油图2 尾矿库区模型图1 尾矿库周边环境如表1 所示为本研究数值计算所采用的物理力学参数，PFC软件中需要设置的参数与一般岩土数值计算软件不同，PFC需要通过设置颗粒与接触模型的细观参数来表征岩土体的宏观参数。基于宏细观参数之间的关系，对本研究的细观参数进行标定。其中摩擦系数因考虑洪水的原因取修正后的摩擦系数0.1。表1 尾矿物理力学微观参数参数数值球半径比Rmax/Rmin1.5球-球弹性模量E/GPa11081108球-墙弹性模量Ew/GPa球-球摩擦系数f球-墙摩擦系数f刚度比k孔隙率n阻尼D3.2尾矿库区模型构建取自初期坝坝址向下游延伸约1.7 km，东西延伸1.3km作为陆地边界。根据尾矿库库区影响区1：1 0 0 0地形图,利用GOCAD软件对尾矿库地形进行建模，模型如图2 所示。网格按照三角形划分，共有网格节点数42 2 5，划分三角形单元数量8 1 9 2 个。再将GO-CAD网格模型导出成DXF格式文件并导入到PFC3D软件中用墙（wall）生成网格作为地形，并填充球颗粒作为尾矿，颗粒数量为51 349 个（+1 40 m），如图3所示。3.3140m堆积高度初期坝溃决尾矿坝溃决后，尾矿如泥沙般倾泻而出，尾矿的势能转换为动能，但因尾矿的黏度较高，尾矿的流速会有一个升高的过程，随着地势平缓后因摩擦力流速会逐渐减小。本次分析1 1 0 0 s内尾矿库1 40 m堆PFC3D5.60T图3PFC3D软件中生成的尾矿库三维计算模型（+1 40 m)0s时，初期坝坝体发生溃决，库区尾砂泥石流开始向下游宣泄；溃坝40 s时，尾砂泥石流最大速度增0.028至9.6 9 m/s，位于尾砂最前端，尾砂最大位移达到36 8 m;0.028溃坝1 0 0 s时，最前方的尾砂泥石流到达回水泵房，2最大速度减小至6.35m/s，最大位移为40 9 m；溃坝0.7200s时，尾砂泥石流完全淹没回水泵房及其周边场0.3内建筑物，初期坝西部坝脚处少量尾砂泥石流未沿河谷向下游演进，而沿周边地势低洼处向西南方向建筑物演进，最大速度继续减小至4.55m/s，最大位移为538 m；溃坝30 0 s时，尾砂泥石流继续向下游演进，最大位移为6 7 4m；溃坝50 0 s时，从初期坝西部坝脚处漏出的少量尾砂泥石流演进至西南方向建筑物，尾砂泥石流继续向下游演进，整体演进速度为1m/s，最大演进位移为8 7 8 m；溃坝1 0 0 0 s时，尾砂泥石流局部仍有0.5m/s左右速度的流动，尾砂泥石流最大位移达到了1 1 2 4m；溃坝1 1 0 0 s时，尾砂整体停止流动最终位移1 2 2 9 m。图4为1 40 m堆积高度初期坝溃决工况下尾砂位移时程图，图中标注了尾矿库下游的村落和建筑物的具体位置，其余未标注的为公路或耕地。3.4140m堆积高度东部堆积坝溃决本次分析1 8 0 s内尾矿库1 40 m堆积高度东部堆27积坝溃坝范围情况，通过数值模拟可以得出不同时刻堆积坝溃坝范围情况，通过数值模拟可以得出不同磷石膏尾矿库下游尾砂位移及流速随溃坝时间的分时刻磷石膏尾矿库下游尾砂位移及流速随溃坝时间布规律。当东部堆积坝溃决时，由于东部局部地势的分布规律。当西部堆积坝溃决时，由于场区所处抬升程度大，尾砂泥石流几乎无法从库区内向外倾地形总体呈现西北高东南低的特征，西部堆积坝坝泻，但由于堆积高度为1 40 m，只有上部的部分尾砂脚地势抬升尤为明显，尾砂泥石流仅少量从库区内滑移1 8 0 s后便停止滑动，最大位移51 m，未影响任向外倾泻，影响下游建筑物范围很小。图7 为1 40 m何下游建筑物。图5为1 40 m堆积高度东部堆积坝堆积高度西部堆积坝溃决工况下尾砂速度时程。溃决工况下尾砂速度时程。Balldisplacement_magBalldisplacement_mag图4尾砂位移时程（初期坝溃决）Ballvelocity_magBall(49419)0.000E-00图6 尾砂位移时程（北部堆积坝溃决)BallallBallvelocity,0.0000E+00Balls(49410)9445E-Ballalls(49826500088888Balls(49826)888Ball velocity2000E-0098260000DOOOE6.2606E.3.2027EOOOOE0088831图5尾砂速度时程（东部堆积坝溃决)3.5140m堆积高度北部堆积坝溃决本次分析38 0 s内尾矿库1 40 m堆积高度北部堆积坝溃坝范围情况，通过数值模拟可以得出不同时刻磷石膏尾矿库下游尾砂位移及流速随溃坝时间的分布规律。0 s时,堆积坝坝体发生溃决；溃坝2 0 s时，由于场区所处地形总体呈现北高南低的特征，只有少量尾砂从北部堆积坝坝脚处漏出；溃坝6 0 s时，从北部堆积坝坝脚处漏出的尾砂演进至下游建筑物；溃坝1 8 0 s时，尾砂泥石流淹没部分下游建筑物；溃坝38 0 s时，尾砂泥石流基本静止，最大位移2 1 0 m。图6 为1 40 m堆积高度北部堆积坝溃决工况下尾砂位移时程。3.6140m堆积高度西部堆积坝溃坝本次分析2 6 0 s内尾矿库1 40 m堆积高度西部图7 尾砂速度时程（西部堆积坝溃决）3.77不同工况溃坝影响分析根据以上溃坝结算结果可知，在1 40 m堆积高度初期坝、北部堆积坝、西部堆积坝溃决工况下，尾矿位移较大，会影响到下游的建筑物，因此将分析这3种工况下的尾砂影响范围。由图8 溃坝影响范围图可知，在设计堆积标高+1 40m溃坝情况下，初期坝溃坝后大量尾矿开始不断向下游滑移。尾矿在滑移至距离尾矿坝最远1 2 2 9 m处停止，此工况为下泄体量最大的情况，下泄体最大宽度42 8 m、最大深度1 5m左右，体积达2.0 1 7 1 0 m，位于尾矿坝坝脚下，随着下泻距离增大，深度差逐渐趋于平缓，途中淹没多处房屋。另有一小部分尾矿沿着下游西边28的低洼地势移动，因低洼地势狭窄，尾矿量很小。在1 40 m堆积高度北部堆积坝溃决工况下，因尾矿库整体北部地势较高在北部堆积坝溃坝后有部分尾矿开始不断向下游滑移。尾矿在滑移至距离尾矿坝最远2 1 0 m处停止，下泄体最大宽度8 1 3m、体积为3.9 51 0 m，最大深度8 m左右，位于尾矿坝坝脚下，随着下泻距离增大，深度差逐渐趋于平缓，途中淹没多处临近尾矿库的房屋。在1 40 m堆积高度西部堆积坝溃决工况下，因尾矿库整体西部地势较高在西部堆积坝溃坝后有部分尾矿开始向下游滑移，但位移量较小。尾矿在滑移至距离尾矿坝最远1 1 3m处停止，下泄体最大宽度6 6 4m、体积为6.31 1 0 m，最大深度9 m左右，在滑移最远处淹没了一处建筑物。FD初期频图8 设计1 40 m堆积高度不同工况下溃决影响范围4应急预案通过分析磷石膏尾矿库设计1 40 m堆积高度在4个不同方向上、不同时刻尾砂位移、流速随溃坝时间的分布规律，可进一步开展尾矿库应急预案研究，在编制传统应急预案中增加相应内容：1）尾矿库排水设施的可靠性和排水能力、尾矿库浸润线控制是影响尾矿库溃坝的主要因素，在尾矿库的日常运行管理与检查中应重点对待。2）圈定撤离范围，明确圈定范围内的人员必须撤离。3）制定尾矿库下游群众遇险疏散转移路线，并在库区内设置逃生路线表述标牌，组织库区及尾矿库下游居民学习。4)标注尾矿库从溃坝到下游各建筑物的初始淹没时间，明确撤离时间。5结论本文通过运用PFC3D颗粒流离散元软件，建立相关磷石膏尾矿库溃坝影响分析数值模型，并在设计1 40 m堆积高度初期坝、东部堆积坝、北部堆积坝及西部堆积坝溃决情况下，开展了溃坝影响分析，得出如下结论：1）1 40 m 堆积高度初期坝溃决工况下，尾砂最终滑移距离为1 2 2 9 m，下泄体最大深度1 5m左右，途中淹没回水泵房及多处房屋，对下游居民造成严重影响。且此工况为溃坝影响最大的工况。2）1 40 m 堆积高度东部堆积坝溃决工况下，由于东部局部地势抬升程度大，尾砂泥石流几乎无法从库区内向外倾泻，只有上部的部分尾砂滑移1 8 0 s后便停止滑动，最大位移51 m，未影响任何下游建筑物。3）1 40 m 堆积高度北部堆积坝溃决工况下，尾砂最终滑移距离为2 1 0 m，下泄体最大深度8 m左右，途中淹没多处临近尾矿库的房屋。4）1 40 m 堆积高度西部堆积坝溃决工况下，尾砂最终滑移距离为1 1 3m，下泄体最大深度9 m左右，只在滑移最远处淹没了一处建筑物。5）溃坝模拟分析给出了溃坝泥石流到达下游重点建筑物的时间，为磷石膏尾矿库制定应急救援预案和应急响应提供了依据。参考文献1何秉顺，付永祥.岩溶地区磷化工固体废物堆场的设计思路：兼议干堆、湿堆方案对比 刀.金属材料与冶金工程，2012,40(S1):94-99.2王国清，顾正聪.磷石膏渣场闭库措施浅析 J.环境科学导刊,2 0 1 6,35(S1):95-97.3刘婷，廖文景，朱远乐.尾矿库溃坝模拟范围及应急预案西部堆积坝研究 J.采矿技术,2 0 2 2,2 2(1)：1 1 9-1 2 3.4鹿守山.尾矿库漫顶溃坝模型试验及数值模拟 D.南昌：南昌工程学院，2 0 1 9.5周振民，李香园.我国尾矿坝溃坝事故成因及生态环境影响评价 J.金属矿山，2 0 1 2(1 1)：1 2 1-1 2 4.6闪淳昌，张振东，钟开斌，等.襄汾“9 8”特别重大尾矿库溃坝事故处置过程回顾与总结 .中国应急管理,2 0 1 1(1 0)：13-18.7JWANG K,YANG P,YU G M,et al.3D numerical model-ling of tailings dam breach run out flow over complex ter-rain:a multidisciplinary procedure J.Water,2020,12(9):2538.8罗昌泰，李栋伟，余国平，等.模型试验及数值模拟下尾矿库溃坝尾砂流演变预测 J.振动与冲击，2 0 2 2，41（3)：307-315.9梁萱，曾智超.基于FLOW-3D的尾矿库数值模拟对下游影响研究 J.江西水利科技，2 0 2 1，47（1)：1 1-2 0.10孙银华，袁利伟，高巍.基于FLUENT的云南某尾矿库溃坝数值模拟研究 .化工矿物与加工，2 0 2 0,49(1 2)：4-8.11张家勇，邹银先，杨大山.基于PFC3D的鱼鳅坡滑坡运动过程分析 J.中国地质灾害与防治学报，2 0 2 1，32（4)：33-39.12JCUNDALL P A.A computer model for simulating pro-gressive large-scale movements in blocky rock systemsC/France:Proceedings of Symposium of Rock Mech-anics,1971:129-136.作者简介朱远乐（1 9 8 7 一)，男，硕士，高级工程师，主要从事尾矿库安全技术的设计、评价、论证等研究工作。E-mail：。(收稿日期：2 0 2 2-0 5-1 7)