张 宇，张树茂，崔 旋，周汉民

(北京矿冶科技集团有限公司，北京100160)

尾矿库作为重大危险源之一，在国际灾害事故排名中位列第18位[1-2]，地库水位超高引起的失稳破坏是引起尾矿库事故的主要原因之一[3-4]。上游法筑坝方法往往浸润线较高，大部分坝体均处于饱和状态，非常不利于库区的稳定性和安全性[5-6]。采取合理的方法计算上游法尾矿库浸润线位置并评估库区稳定性，对保证库区安全具有重要的理论意义和实际工程价值。

极限平衡法和强度折减法是目前使用最为广泛的两种尾矿库稳定性分析方法。极限平衡分析方法具有原理简单、物理意义明确等优点，是工程中应用最为广泛的稳定性分析方法。但是采用极限平衡法进行尾矿库稳定性分析时需预先假设滑裂面形状和初始滑裂面位置，且无法同时考虑渗流和变形等因素对尾矿库稳定性的影响，具有一定的局限性。

随着计算机技术的发展，以强度折减法为代表的有限元稳定性分析方法得到迅速发展。采用强度折减法进行土体稳定性分析时，不需要预先假设滑裂面的形状和位置，可以同时获得整个土体的应力场和变形场，弥补了极限平衡法的不足[7]。该方法的在尾矿库稳定性分析中应用的关键，是如何合理地选择失稳破坏判别准则。目前土体稳定性分析中使用最为广泛的失稳破坏准则有“数值收敛准则”、“位移突变准则”和“滑裂面贯通准则”，各准则在尾矿库失稳判别中的适用性尚有待考证。

基于此，本文以国内某上游法堆筑尾矿库为研究对象，模拟了干滩长度为100 m的尾矿库渗流场及变形场的分布特性。在此基础上，采用强度折减法分析了尾矿库的稳定性，评价了不同失稳判别准则的有效性和局限性，并提出了一种适用于上游法尾矿库的失稳判别准则，以期为对渗流作用下尾矿库的稳定性分析提供技术支撑。

1 计算软件与计算理论

1.1 ABAQUS软件简介

本文采用通用有限元计算软件ABAQUS进行了尾矿库渗流稳定性分析。该软件提供了非饱和土渗流应力耦合的分析功能[8]，并包含十分丰富的材料模型、单元模式和荷载边界条件，易于二次开发，对岩土工程有较好地适用性[9]。

1.2 饱和-非饱和渗流理论

由于库内回水及降雨的影响，非干堆型尾矿库可以看做不可压缩的固相尾矿颗粒和液相水，以及一定数量的可压缩气体的三相体系。浸润面以下，尾砂被认为是液—固两相体系，遵循以比奥固结理论为基础的饱和土渗流本构方程；浸润面以上，非饱和尾砂为固—气—液三相体系，其渗透本构由渗透系数函数和饱和度(含水量)函数进行表征[10]。

大型有限元软件ABAQUS基于固定网格，按照非饱和土力学理论进行求解，浸润面取为孔隙水压力为零处，在进行饱和-非饱和渗流时具有较好的计算精度和效率[11]。利用ABAQUS进行不饱和渗流的材料设定时，渗透系数的关系可由一个与饱和度相关的折减系数进行修正，如下式所示。

kunsat=ks·ksat

ks=(Sr)3Sr<1.0

ks=1.0Sr≥1.0

式中：kunsat为非饱和渗透系数，ks为折减系数，ksat为饱和渗透系数，Sr为土体饱和度。

在非饱和土体中，基质吸力直接受材料饱和度的影响，二者关系可由吸湿(脱水)曲线表征，通常以基质吸力(负孔压)—饱和度的关系曲线给出。

1.3 强度折减理论

强度折减法是目前岩土行业公认的可靠的、较为先进的土坡稳定分析计算方法。该方法采用抗剪强度折减系数Fr代替传统意义上的稳定安全系数Fs，又称为强度储备安全系数。计算中假定不同的强度折减系数Fr，根据折减之后的强度参数进行有限元分析，观察计算是否收敛。在整个计算过程中不断增加Fr，达到临界破坏时的强度折减系数Fr就是边坡稳定安全系数Fs，折减后的抗剪强度参数可表达为：

cm=c/Fr

φm=arctan(tanφ/Fr)

式中：c和φ为土体抗剪强度，cm和φm为维持平衡所需的或土体实际发挥的抗剪强度；Fr为强度折减系数。

2 计算模型及材料参数

某尾矿库采用初期坝+上游法堆坝，总坝高210 m。初期坝为碾压堆石坝，坝高60 m，坝顶宽为4 m，上游边坡坡比1∶1.6，下游边坡坡比1∶1.8；堆积坝平均堆积边坡为1∶5；两侧岩体为中风化白云质灰岩，尾矿库分区如图 1所示。尾矿坝二维有限元模型被离散为3 153个节点及3 020个单元，采用具有孔压自由度的四节点平面应变单元和少量退化的三角形单元进行模拟。

图1 尾矿库分区及有限元网格示意图Fig.1 Zoning and finite element mesh of tailing reservoir

计算中有限元的网格固定于土骨架，气体或液体在流过网格的过程中满足流体连续方程。土体采用基于有效应力指标的摩尔库伦弹塑性模型模拟，液体渗流采用Forchheimer渗透定律模拟，参数见表1及图2，所有材料泊松比取为0.3。需要注意的是，一般岩土工程勘察报告以及《尾矿设施设计规范》主要提供土体的压缩模量Es，而压缩模量难以直接应用于有限元数值计算中。贾堤[12]等人探讨了利用场地的岩土工程勘察报告数据获得数值分析中土体弹性模量的三种方法，建议采用公式E=αEs来估算弹性模量E，其中α建议取为8.2，Es为土体压缩模量，本文采用该方法计算尾砂材料的弹性模量。

表1 计算参数表Table 1 Calculation parameters of various materials

图2 土水特征曲线Fig.2 Characteristic curves of soil water

3 非饱和渗流条件下尾矿库流固耦合分析

浸润线的位置是决定尾矿库渗流情况及稳定性的关键因素[13-15]。本文考虑了饱和度及基质吸力对尾砂渗透系数的影响，图3给出了干滩长度分别为100 m时尾矿库浸润面位置分布情况，即总水头为0的水压力等值线分布情况。可以看出，尾矿库浸润线基本沿下游坝坡面分布，在坝中部位置略有抬高，最小埋深5 m。由于初期坝为透水堆石坝，因此浸润面均经过初期坝上游坡面流向下游坝外。

图3 干滩长度为100 m时尾矿库浸润线分布Fig.3 Distribution of phreatic line in tailings reservoir with the beach width of 100 m

图4和图5分别给出渗流条件下尾矿库变形场分布及变形矢量图。可以看出，最大水平向变形主要发生在库区1/3高度的下游坡面处，最大竖向变形主要发生在库顶滩面处。尾砂的变形矢量主要指向库底，即库区内尾砂主要发生向下的竖向变形；上游坡面处尾砂略微出现指向上游的水平变形。

图4 尾矿库变形分布图Fig.4 Distribution of deformations in tailings reservoir

图5 尾矿库变形矢量图Fig.5 Deformation vector diagram of tailings reservoir

4 尾矿库渗流稳定性分析

在采用强度折减法对尾矿库进行稳定性分析时，强度折减过程尾矿库失稳判别准则的选择至关重要，目前常用的判别准则有：1)“数值收敛准则”：以数值计算收敛与否作为评价标准，其与有限元的算法有关；2)“位移突变准则”：以特征部位的位移拐点作为评价标准；3)“滑裂面贯通准则”：以是否形成连续的贯通区作为评价标准。其中，影响数值计算收敛与否的因素较多，如数值计算方法的稳定性、计算模型网格质量等，采用“数值收敛准则”作为稳定性判据可能带来误差[16]。因此本文主要采用“位移突变准则”和“滑裂面贯通准则”作为失稳判别准则，对干滩长度为100 m的尾矿库进行了稳定性评价，计算结果如图6～7所示。

图6给出尾矿库外坡面坡顶和坡中特征点处位移与强度折减系数的关系。可以看出，坡顶特征点1的强度折减系数Fr由1.41增大到1.43的过程中，位移迅速增加而发生突变，该点对应的安全系数取Fs为1.41。同理，坡中特征点2的安全系数Fs应取为1.43。可以看出，采用“位移突变准则”不同的特征点对应的安全系数不同，难以选取某个单一的值作为整个尾矿库的安全系数。另外，由于位移是一种标量指标，难以反映位移方向对边坡稳定性的影响。

图6 特征点位移与强度折减系数的关系Fig.6 Relationship of displacement and stability factor of characteristic points

图7 尾矿库滑裂面塑性应变分布Fig.7 Plastical strain distribution of slipping surface of tailings reservoir

图7给出了强度折减系数Fr分别为1.41和1.43时尾矿库内塑性应变的分布情况。这两种工况下库区内均出现了大范围滑裂带，滑裂面位置大致相同，贴近基岩面的深层滑动。这是由于尾矿库的分区导致各层尾砂的抗剪强度逐层减弱，强度最低的尾粉质黏土主要沉积在库底，形成抗剪薄弱面。需要注意的是，由于构成初期坝的堆石体具有透水性好、强度高的特点，因此塑性应变区无法通过初期坝而形成贯通面。即由于初期坝的存在，在采用强度折减法进行上游法尾矿库稳定性分析的时候，采用“滑裂面贯通准则”难以对尾矿库是否失稳进行准确判别。

本文提出一种新的尾矿库失稳评价准则，即“变形矢量图准则”，该方法以初期坝附近尾砂出现指向坝外、及沿着初期坝上游坡向上的变形矢量为整个尾矿库失稳的判别依据。图8给出Fr=1.41、Fr=1.43、Fr=1.53工况下尾矿库整体变形矢量图。Fr=1.41时，尾矿靠近外坡附近堆积尾砂主要产生指向库底及库外的变形，但变形矢量方向大部分未指出库外，仅在初期坝附近呈现平行于堆积外坡面的变形；Fr=1.43时，库区中下部尾砂整体出现指向库外的变形，由于初期坝的阻碍作用，靠近初期坝的尾砂产生沿着初期坝坡向上的变形，尾矿库发生失稳破坏。随着Fr增大为1.53，失稳破坏继续发展。因此，该尾矿库的安全系数Fs取为1.41。通过上述分析可知，“变形矢量图准则”既可以直观地呈现整个尾矿库的变形情况，又可以给出滑动体的分布情况，从而得到整体尾矿库的稳定性综合安全系数，克服了单独采用“位移突变准则”和“滑裂面贯通准则”的不足。

图8 不同强度折减系数Fr下尾矿库变形矢量图Fig.8 Deformation vectors of tailings reservoir with various Fr

5 结论

本文采用强度折减法对某上游法尾矿库饱和—非饱和渗流条件下稳定性进行了系统的分析，对比了不同失稳判别准则在尾矿库稳定性分析中的适用性和局限性，在此基础上提出了一种新的判别准则。主要结论如下：

1)本文基于饱和—不饱和渗流理论，对尾矿库渗流条件下的浸润线分布和变形特性进行了分析。研究结果表明，渗流作用下尾矿库竖直方向主要发生竖向沉降，随着高度的增加竖向沉降增大；水平向主要发生指向库内的变形，且在外坡面1/3处达到最大。

2)尾砂在库区内呈现强度较高的粗颗粒先沉积，强度较低的细颗粒后沉积的分区特性，因此越靠近库底及库内岸坡尾砂的抗剪强度指标越低，尾矿库一旦失稳，将发生沿着库底岸坡的深层失稳破坏。

3)本文建议采用“变形矢量图准则”作为上游法尾矿库的失稳判别准则。该方法以初期坝附近尾砂出现指向坝外、及沿着初期坝上游坡向上的变形矢量为标志，可综合地判断尾矿库的整体稳定安全系数。