基于Geo-Studio的某钨尾矿库应力变形数值分析

许志发，龙 翠，唐永骏，高亚伟，赵怀刚，王光进

（1.昆明理工大学 国土资源工程学院，云南 昆明 650093；2.云南铝业股份有限公司安全环保部，云南 昆明 650502；3.昆明理工大学建筑工程学院，云南 昆明 650504）

0 引 言

矿山三大控制性工程之一的尾矿库是一种特殊的工业建筑物，是矿山企业在生产建设过程中的重要设施。中国是钨资源大国，储量占世界总量的68%，排名世界第一[1]。但随着经济的发展，对钨的需求量不断增长，钨尾矿库的数量和等级逐年增加[2]。因此，尾矿库的运行状况不仅与矿山能否顺利进行生产建设，还关系到尾矿库下游的周边环境安全和群众的生命财产安危[3-4]。自1910年以来，已经有数百起尾矿库事故在世界各地发生[5]。其中三分之一的尾矿库属“病险库”，因尾矿库是具有高势能的人造危险源，一旦失事，易造成非常严重的事故。例如，2008年9月8日，位于山西省的新塔矿业有限公司发生尾矿库事故造成281人死亡，直接经济损失9 619.2万元的特别重大事故；2010年9月21日，位于广东的信宜紫金矿业有限公司所属的银岩锡矿尾矿库发生尾矿库溃坝事故造成坝体下游4人死亡，直接经济损失约1 900万元，这些事故产生了极其恶劣的社会影响[6-8]。这些都是人为责任事故，如果能够规避尾矿库建设过程中的违规行为，并对筑坝过程中的应力变形进行实时有效的监控，异常变形和坝体开裂等异常现象能够被及早发现，筑坝作业中的不规范得到调整，矿山的生产建设和下游群众的生命财产安全将得到有效的保障。

尾矿库是由各种不同粒径尾矿颗粒的堆砌形成的人造土体，近年来国内外学者对其应力变形特性进行了大量的试验和探讨。Floresc[9]分别对墨西哥两种尾粉砂进行共振柱试验，应变范围为 10－5～10－2，用Davidenkov模型对试验结果进行拟合，得出了两个参数的试验数据；袁泽川[10]运用剪应变归一化和等效线性模型的方法研究围压对无黏性土、黏性土和一般性土的最大剪切模量的影响，得到尾矿在不同围压下的阻尼比和动模量衰减曲线；曹进海等[11]针对余震对尾矿坝稳定性的影响问题，以某上游式尾矿坝为例，基于完全非线性有限差分法，从不同余震烈度、网格变形、位移变化及超孔压比等方面对尾矿坝进行数值模拟分析；柳厚祥等[12]研究了尾矿颗粒在静力、动力和共振情况下的物理力学性质，模拟了在不排水有效应力法和排水有效应力法下，对尾矿库在地震作用下的应力变形进行分析；周健[13]将尾矿库看作三维两相体，研究了振动孔隙水压力增长规律，进而运用三维两相有效应力动力分析方法对尾矿库的应力形变进行分析。

研究选取四川省某钨尾矿库，利用Geo-Studio对尾矿库的应力变形运用有限元法进行模拟，分析尾矿库在原设计标高和在二期设计标高堆筑下尾矿库整体的应力应变等状态，为尾矿库的设计和施工提供指导性依据。

1 应力变形数学模型

有限元法、离散元法、拉格朗日元法及边界元法等在尾矿库、边坡等岩土工程都有较为广泛的应用。目前，尾矿库与边坡已广泛应用基于有限元法的数值模拟软件对应力应变和稳定性进行分析计算。有限元法是将复杂的物体或者系统划分为有限不重叠的单元，这些单元相互间独立又互相联结。所有作用的力在每个单元上都等效到节点上，最后将各变量或其导数的节点值和插值函数组成的线性表达式代替微分方程中的变量。有限元分析法可以用数值模拟软件来分析难以人工计算的巨量数和方程组。

1.1 有限元法基本原理

在进行分析之前需要对物体或系统进行网格划分，在二维有限元分析中可以选用三角形、矩形、四边形等多边形进行系统分析。在这几种网格划分中运用三角形进行网格划分，具有网格密度大、求解精度高等优点。常应变状和刚体位移态是单元的位移函数必须满足的条件，位移必须连续[14]。因此，函数u（x，y），v（x，y）可写为式（1）所示。

将位移值和坐标值代入式（1），得到6个求解方程，将其联立获得各个系数的求解值。将其再代入式（1）得到式（2）。

式中：｛｝δ =[u1，v1，u2，v2，u3，v3]T为 节 点 位 移 矩 阵 ；称为形状函数矩阵，其，i＝1，2，3称为形状（插值）函数。

确定位移函数之后，求得单元的应变和应力。单元应变如式（3）所示。中

式中：[B]为单元应变矩阵。

将式（3）代入，求得单元应力如式（4）所示。

式中：[D]为弹性矩阵。

设坝体发生虚位移，网格单元虚位移为{δ*}，与之的网格单元虚应变为{ε*}，则有：

式中：A 为单元面积，m2；t为单元厚度，m；S 为单位虚位移，m；T为转置矩阵。

经变换得：

式中：{U}为总体位移列阵；[K]为总体刚度矩阵；{P}为总体荷载列阵。

总体刚度方程通过边界约束条件进行修正后，总体位移列阵{U}可以求解得到，然后计算各单元的应力和应变分量，进而对物体或系统的应力和应变进行模拟分析。

1.2 初始模量

不同初始应力下，土的初始模量Ei是变化的。根据三轴试验，若把Ei和相应的σ3用大气压力Pa归一成无量纲数，可得到关系式（7）：

式中：k、n是取决于土的性质参数，分别称为模量系数和模量指数，可根据三轴试验测出。

1.3 体积模量和泊松比

材料变形中，反应体积变化的参数一般为体积模量B或泊松比μ，两者之间的关系见式（8）所示。

体积模量和泊松比也可以表示为与应力状态有关，即：

式中：kb、m、G、F、D 皆为试验常数；c为内聚力，kPa；φ 为内摩擦角，（°）。

2 尾矿库的应力变形模拟分析

2.1 Geo-Studio概况

加拿大GEO-SLOPE公司在20世纪70年代开发了Geo-Studio，该软件是可以应用于岩土、地质、采矿等各行各业数值分析的软件之一。它具备操作方便、模块化计算、建模简便、计算分析快等优点，矿业、岩体工程等许多领域已广泛应用。SIGMA/W是Geo-Studio系列软件的重要组成部分，是应力变形分析中专业的有限元分析软件，即可实现线弹性分析，也可以完成非线性弹塑性分析等高度复杂的各种分析。Geo-Studio软件采用增量形式进行求解，可以选择线弹性、各向异性弹性、双曲线弹性、弹塑性等多种材料模型，能够很好地进行边坡、坝体等问题的模拟计算分析。

2.2 应力变形模拟

2.2.1 尾矿库概况

尾矿库在2012年完成设计和施工，由于原设计不能满足矿山扩大生产的要求，2017年进行尾矿库增高扩容工程的设计。地貌单元为山地中山，沟谷相间地形，两侧山坡地势较缓，随地形变化25°～35°，表面大部为第四系残坡积层覆盖，局部可见基岩出露，谷底为第四系冲洪积层覆盖。尾矿库堆积坝一期设计最终标高2 035 m（总坝高92 m），使用段沟底纵坡8.7%。因一期设计不能满足未来尾矿的排放需求，因此对尾矿库进行增高扩容，增高扩容后的设计坝高为147 m（即加高55 m），加高后的坝顶平台标高为2 090 m。根据《尾矿设施设计规范》（GB50863—2013）[15]对尾矿库的分类，原尾矿库为三等库，加高坝高后，按坝高考虑尾矿库为二等，故尾矿坝等尾矿库主要构筑物的级别为二级。

2.2.2 尾矿库稳定性计算

库容增加后尾矿库的稳定性应该满足规范要求，因此针对库容增高后的稳定性进行研究。由于尾矿库所在区域的地震基本烈度为Ⅷ度，所以需考虑地震作用下坝体的稳定性。本文选择三种运行工况开展尾矿库稳定性分析，一是尾矿库处于正常运行工况；二是在洪水运行工况，其洪水位按二等尾矿坝的最小安全超高（1.0 m）考虑；三是在特殊运行工况，即考虑地震基本烈度Ⅷ度时的稳定性。前两种工况运用Slide软件进行尾矿库稳定性分析。采用毕肖普法对尾矿库进行稳定性分析，计算结果如图1、表1所示。结果显示在坝高抬高至2 090 m时，尾矿库稳定性仍满足规范要求。

图1 2 090 m时两种工况毕肖普法计算结果Fig.1 Bishop method results of two conditions at 2 090 m

表1 高程2 090 m时的抗滑稳定性计算结果表Tab.1 Results of anti-sliding stability at an elevation of 2 090 m

地震工况下的尾矿库动力稳定性分析采用Geo-Studio系列软件的QUAKE/W模块进行模拟。其可以用来解决因地震而产生的动态载荷以及碰撞作用下的冲击载荷等作用下的动力稳定性问题。应用该模块，可以同时考虑因地震产生而作用在土体中的孔隙水压力增长和动应力，能够研究地震对土工建筑物强度稳定性的影响。在动力稳定性分析过程中，选用汶川地震波和人工合成地震波两组典型地震波对尾矿库开展地震工况的边坡稳定性分析。

采用有限元极限平衡法对地震作用各时段进行稳定性计算，得到汶川地震波和人工合成地震波两组典型地震波对坝体作用的最小安全系数。由图2可知，两种地震波作用过程中，坝体安全系数均大于1，汶川地震波作用过程中最小值为1.320，人工合成地震波作用过程中最小值为1.212，均符合规范要求。

2.2.3 几何模型与计算参数

坝体应力变形有限元分析计算剖面选择最不安全剖面（即沿谷底计算剖面）在堆高至2 043 m和2 090 m时开展其应力变形分析，而材料分区仍然采用二维稳定性的概化分区（即距子坝30 m的位置作为尾细砂与尾粉砂的分界线；距子坝180 m的位置作为尾粉砂与尾粉土的分界线；距子坝350 m的位置作为尾粉土与尾粉质黏土的分界线）。同时，计算的2 043 m和2 090 m高程时的浸润线埋深采用的是刚好满足规范规定最小安全系数时的临界浸润线埋深，即尾矿库在堆高至2 043m时其坝体的临界浸润线埋深值为14.0 m，而堆高至2 090 m时坝体的临界浸润线埋深值为25.0 m。依据室内土工试验数据分析，同时与矿山提供的勘察报告进行参照和国内外尾矿库设计中相关参数，进行全面的对比分析，模拟分析中所需的各材料的参数如表2所示。

2.2.4 边界条件与几何划分

考虑到尾矿库的实际情况，在分析中边界条件的确定如图3所示（一期：2 043 m，二期：2 090 m），由于在实际中库尾一般为山体、硬岩等自然地形，因此将左侧边界确定为水平方向的约束，即U=0，V≠0；尾矿库底部在堆坝前都需要提前将淤泥、腐殖土、杂填土等清除干净，使底部形成一个相对稳定的地基，因此将底部边界的水平和垂直约束都设定为约束，即U=0，V=0。尾矿库原堆积坝高2 043 m和最终堆积坝高2 090 m时应力变形计算边界条件和网格划分结果如图3所示。

表2 尾矿坝应力变形计算参数Tab.2 Calculation parameters for tailings dam stress deformation

图3 沿谷底计算剖面不同堆积坝高应力变形计算图Fig.3 Computational map of high stress deformation of different accumulation dams calculated along the valley floor

2.2.5 尾矿库应力变形数值计算结果

在计算中考虑到地下水对坝体应力变形的作用，分别对正常工况下坝高堆至2043m（坝高100 m）和最终坝高2 090 m时（坝高147 m）的水平、竖直方向位移，水平、竖直方向总应力和剪应力总应力进行应力和变形数值模拟，计算结果见图4、图5。

图4 坝高堆至2 043 m时正常工况下位移等值云图、总应力云图及剪应力总应力云图Fig.4 Displacement equivalent cloud image,total capacity cloud image and shear stress total stress cloud image under normal conditions of up to 2 043 m dams

3 结 论

本文采用数值模拟软件对某钨尾矿库不同堆高下的稳定性和应力变形进行模拟，分析其稳定性、沉降和应力应变，研究了其变化规律。

（1）分析库容增加下正常运行工况、洪水运行工况和地震运行工况下三种不同工况尾矿库的稳定性，分析结果表明库容和坝高增加的情况下尾矿库的稳定性符合设计规范要求。

（2）坝高堆至2 043 m（总坝高100 m）和最终坝高2 090 m（总坝高147 m）时，坝体变形规律基本一致。坝体水平向变形主要向着坝体下游，最大水平向位移出现在坝坡中上部，竖向变形自上往下逐渐降低，最大竖向变形出现在坝坡顶部。

（3）不同堆高情况下在初期坝内坡脚都出现一定的应力集中现象，但坝体内大部分区域剪应力不大。

图5 最终坝高2 090 m时正常工况下位移等值云图、总应力云图、剪应力总应力云图Fig.5 Displacement equivalent cloud image,total capacity cloud image and shear stress total stress cloud image under normal conditions of up to 2 090 m dams

（4）堆积时向下固结沉降和尾矿沉降压缩导致侧向膨胀，引起水平方向位移。随坝体增高，坝体变形逐渐增大，且下层堆积坝会在上层堆积后出现局部“隆起”。同时，干滩面长度减小也会使坝体水平位移和垂直位移有所增大。

[1] 邬海滨，李继福，袁勤智，等.回收某硫化矿尾矿中白钨的选矿试验研究[J].中国钨业，2017，32（3）：36-41.WU Haibin，LI Jifu，YUAN Qinzhi，et al.Beneficiation of scheelite from a tailings of a sulphide ore[J].China Tungsten Industry，2017，32（3）：36-41.

[2] 李茂林，王 旭.甘肃某白钨浮选尾矿再回收白钨试验研究[J].中国钨业，2015，30（2）：17-20.LI Maolin，WANG Xu.Study on the recovery of scheelite from a scheelite flotation tailing in Gansu province[J].China Tungsten Industry，2015，30（2）：17-20.

[3] 魏作安，尹光志，沈楼燕，等.探讨尾矿库设计领域中存在的问题[J].有色金属（矿山部分），2002，54（4）：44-45，48.WEIZuoan，YINGuangzhi，SHENLouyan，etal.Discusstheproblems existing in the design of tailings ponds [J].Nonferrous Metals（Mining Section），2002，54（4）：44-45，48.

[4] MILTON A，COOKE J A，JOHNSON M S.A comparison of cadmium in ecosystems on metalliferous mine tailings in Wales and Ireland[J].Water，Air and Soil Pollution，2004，153（1）：157-172.

[5] OZCANNT，ULUSAY RESAT，ISIKNS.A study on geotechnical characterization and stability of downstream slope of a tailings dam to improve its storage capacity （Turkey）[J].Environmental Earth Sciences，2013，69（6）：1871-1890.

[6] 谢旭阳，田文旗，王云海，等.我国尾矿库安全现状分析及管理对策研究[J].中国安全生产科学技术，2009，5（2）：5-9.XIE Xuyang，TIAN Wenqi，WANG Yunhai，et al.Analysis of safety status and management countermeasures of tailings in China[J].

[7] 张兴凯，王启明，相桂生.金属非金属尾矿库安全现状及分析[J].中国安全生产科学技术，2006（2）：60-62.ZHANG Xingkai，WANG Qiming，XIANG Guisheng.Safety status and analysis of metal nonmetal tailings storage[J].Journal of Safety Science and Technology，2006（2）：60-62.

[8] 李青石，李庶林，陈际经.试论尾矿库安全监测的现状及前景[J].中国地质灾害与防治学报，2011，22（1）：99-106.LI Qingshi，LI Shulin，CHEN Jijing.Current status and prospect of tailings dump safety monitoring[J].The Chinese Journal of Geological Hazard and Control，2011，22（1）：99-106.

[9] FLORES O C，ROMO M P O.Dynamic behavior of tailings[C]//Shamsher Prakash.Geotechnical Engineering Commons.California：University of Missouri—Rolla，2001：1-6.

[10] 袁泽川.尾矿的动力应力-应变特性 [J].工业建筑，2017，47（10）：139-145.YUAN Zechuan.Dynamic stress-strain characteristics of tailings[J].Industrial Construction，2017，47（10）：139-145.

[11] 曹进海，胡 军.基于FLAC3D余震对尾矿坝稳定性影响的研究[J].中国钨业，2017，32（4）：68-73.CAO Jinhai，HU Jun.Study on the influence of FLAC3Daftershocks on the stability of tailings dam [J].China Tungsten Industry，2017，32（4）：68-73.

[12] 柳厚祥，廖 雪，李 宁，等.高尾矿坝的有效应力地震反应分析[J].振动与冲击，2008，27（1）：65-70，92.LIU Houxiang，LIAO Xue，LI Ning，et al.Effective stress seismic response analysis of high tailings dam[J].Journal of Vibration and Shock，2008，27（1）：65-70，92.

[13] 周 健.尾矿坝在地震作用下的三维两相有效应力动力分析[J].工程抗震，1995（1）：39-43.ZHOU Jian.Three-dimensional two-phase effective stress dynamic analysisoftailingsdamunderseismicaction[J].EarthquakeResistant Engineering，1995（1）：39-43.

[14] 李慧谦.尾矿库堆坝试验研究及稳定性分析[D].重庆：重庆大学，2008.LI Huiqian.Experimental study and stability analysis of tailings dam[D].Chongqing：Chongqing University，2008.

[15]中华人民共和国国家标准编写组.GB50863—2013尾矿设施设计规范[S].北京：中国计划出版社，2013：1-50.