周 涛 李启航 黄宜超

(江西理工大学资源与环境工程学院，江西 赣州 314000)

尾矿库是金属或非金属矿山尾矿、工业废渣的储存场所，通常通过筑坝或者围地构成。矿山的安全生产及周边居民的生命财产安全与尾矿库坝体的稳定性密切相关。尾矿坝的安全在世界各国的矿山建设过程中都受到密切关注，各国都十分重视尾矿坝的建设与管理。

长期以来专家学者利用室内相似材料模拟实验和数值模拟技术对尾矿坝开展了广泛而深入的研究。尾矿坝溃坝的失稳模式可分为洪水漫坝失事和坝体结构失事，坝体结构失事模式又可分为渗透破坏失事模式、坝体失稳模式和地震险情失稳模式[1]。张兴凯和陶东良都采用了相似物理实验法分析了尾矿库漫顶溃坝的全过程[2,3]。李旭[4]研究了降雨诱发的尾矿坝溃坝，利用Slide软件对降雨下的边坡进行定量分析。敬小非等[5]基于云南拉拉铜矿小打鹅尾矿库工程设计资料,设计了相似材料模型实验，模拟了洪水作用下尾矿坝的垮塌和溃决机制，结果表明：尾矿堆积坝的浸润线变化存在滞后性；在水位上升过程中，坝坡中部出水平方向的总应力增量较垂直总应力增大；洪水导致尾矿坝所受的渗透力、孔隙水压力、重力增大，削弱坝体材料的抗剪强度，加大自身荷载，从而发生结构失稳。

尹光志等[6]指出尾矿坝组成材料有明显的颗粒物质特征，在数值模拟中，相对于连续介质方法，采用离散元法探索尾矿颗粒组成的结构体的破坏特征以及规律更加合适。众所周知在数值模拟中，正确的计算参数对最终结果的正确性起着决定性作用。在尾矿库溃坝数值模拟中，张千贵等[7]基于离散元软件PFC2D，采用了粘结弹性接触模型，通过不断调整细观参数并进行双轴数值模拟实验来拟合实际室内实验的曲线，从而确定最终的细观参数值。刘春等[8]利用其开发的MatDEM，进行了茂县滑坡模拟，在实测野外岩土体宏观力学参数后，输入软件进行自动材料训练，软件自动给出最符合宏观力学性质的细观参数。利用MatDEM的自动训练功能可以免去繁琐的试算测试过程，可以节省大量时间精力。王学良等[9]在连续介质软件Massfloss中根据前人对滑面摩擦参数确定了相应参数，同时在PFC3D考虑尾矿坝为饱和颗粒的极端情况，在软件中将离散元颗粒视为弹性接触模型，不考虑其粘聚力的影响。张晓宇等[10]基于离散元软件MatDEM模拟了黏性土失水开裂过程，在根据参数完成建模后，利用前人总结的公式实现了对计算参数的动态计算，数值模拟效果十分理想。此方法实际是根据实验公式决定材料的参数，依赖于前人对公式总结的正确性。目前，通过数值模拟方法物理模型的手段的主要研究对象是洪水漫顶导致的溃坝模拟，对其他溃坝模式的实验以及数值模拟研究很少。利用离散元法进行细观模拟的研究也很少。本文基于离散元软件MatDEM利用预设参数建立了尾矿坝离散元模型。在模拟中，模拟了尾矿坝逆流渐进式滑坡的动态过程。通过对比发现坝体破坏特征以及颗粒流动特征与前人室内实验结果较为一致。

1 离散元接触模型

1.1 单元间的接触模型

离散元法通过堆积和胶结具有特定力学性质的颗粒来构建岩土体结构。在此基础上通过时间步迭代算法来进行数值模拟。本文模拟的尾矿坝溃坝实验采用最基本的线弹性接触模型，两颗粒间靠法向与切向的弹簧来模拟相互作用力。通过以下公式来表征颗粒间的弹簧在法向与切向上的破坏准则。

Fn=KnXn

(1)

Fn=KsXs

(2)

Fsmax=Fs0-μpFn

(3)

其中，Kn为法向刚度；Ks为颗粒间的切向刚度；Fs0为颗粒间的初始抗剪强度；μp为单个颗粒的摩擦系数；Fn为垂直于切向弹簧上的法向压力。当式(1)中的Xn超过断裂位移Xb时两个颗粒间的法向弹簧断裂；当单元间的切向力超过式(3)中的Fsmax时两颗粒间的切向弹簧断裂，颗粒开始滑动。宏观表现就是断开的位置连成了一条弧线岩土体出现滑移面。因此式(1)以及式(3)中的参数的正确标定对模拟结果起决定性作用。

1.2 两个不同单元的连接

每个单元都有自身的刚度kn，两颗粒接触时实际上是两根弹簧的串联，式(1)中的Kn是两个单元间的连接刚度，对于法向刚度为kn1和kn2的两个单元，其等效刚度Kn为：

(4)

对于切向刚度为ks1和ks2的两个单元，等效切向刚度Ks为：

(5)

MatDEM在完成材料赋值操作后，不同材料的颗粒间的连接刚度会根据公式自动调整。同一种材料不同粒径的颗粒也会赋予不同的颗粒刚度，因此颗粒间的连接刚度也不同。

2 尾矿坝模型建立

尾矿坝模型的建立步骤主要有堆积建模、重力沉积、切割地层分组与材料赋值等。边坡主要包含初期坝、尾矿坝主体等，分别给这两个赋予表1中材料砂土、水质以完成建模。

表1 建模中采用的两种土体材料的力学参数

2.1 随机堆积模型

堆积建模目的是模拟真实世界的重力沉积过程。如图1a)所示,建立5.5 m×1.28 m的长方形模型箱，设置分散系数，在模型箱内生成0.14 cm～0.78 cm之间的随机单元。堆积模型中生成的不同粒径大小自身有着不同的切向刚度与法向刚度，不同粒径颗粒间的连接的刚度也不同。重力沉积的过程中，规则的排列的单元被赋予了随机的初速度，单元运动一段时间后到达随机位置上；单元颗粒在重力的作用下向下沉积以模拟颗粒自然沉积的过程；利用上压力板产生的荷载来压实模型以模拟颗粒在重力沉积后经历的压实作用。

2.2 数字高程建模与材料赋值

MatDEM可以利用导入的坐标点生成曲线或者曲面并以此为基础划分地层，对于任何复杂的几何模型在确定了各层面的坐标后便可以很方便的构建。如图1b)所示，曲线描述了不同地层界面的坐标位置。在堆积模型建成的基础上，将这些曲线的坐标位置导入MatDEM后便可以将堆积模型中多余的随机颗粒切割去除，效果如图1c)所示。在完成了地层切割过后，对各地层赋予相应的材料参数，如图1c)中，深灰色地层被赋予了水的材料性质以模拟库中水，黑色地层被赋予了砂土材料。切割完地层并且赋予不同地层材料后，颗粒间的刚度发生了改变，地层的受力因此发生改变。再次对模型进行受力平衡，为了防止模型破坏颗粒飞出，需要加强颗粒间连接使其不可断裂，在重力平衡后再将颗粒间的连接设置为可断裂的状态。至此边坡模型的建立完成。MatDEM材料的细观参数通过宏观参数按照公式自动转化而成，详细公式见文献[11]。

3 尾矿坝饱和—非饱和模拟

3.1 渗透力与重力

实际工程中很难对每个孔隙的渗流情况进行研究，也是没有必要的，设定浸润线如图1d)所示,浸润线以上为非饱和区域，浸润线以下为饱和区域。研究坝体渗流，不考虑渗流在土壤孔隙中流动路径的迂回曲折，只考虑渗流的主要流向，并认为全部的渗流空间均被水所充满。在模拟饱和砂渗流主要作以下假设：1)坝体内部渗流为层流，认为坝内符合达西定律。2)坝体内部渗流为渐变流(杜平假定)，认为渗流中任意过水断面的水平流速和比降都是相等的。

基于以上假设通过自编程序为每个浸润线以下的颗粒赋予渗透力。根据达西定律，在水平与垂直方向赋予颗粒体力以合成与浸润线方程切线方向一致、大小符合达西定律的渗透力。根据每个颗粒的坐标，实现给每个颗粒竖直方向赋予垂直向上的力以模拟浸润线下的浮力。

3.2 饱和砂性能的模拟

浸润线下经过库水的浸泡，浸润线以下的坝体变为了饱和砂。砂土饱和后颗粒间的连接破坏要比非饱和容易，这是因为颗粒间充满了液体，颗粒间的距离增大，使得颗粒间的连接更加容易断裂，砂土的抗张拉能力大幅度下降。王海东等[12]研究了不同饱和度下砂土的抗剪强度，发现内摩擦角随含水率的增加呈线性下降；粘聚力随着饱和度增加经过上升达到峰值后呈线性减小。如式(3)，内摩擦系数与初始抗剪力是影响切向弹簧可承受最大力的重要指标。

4 数值模拟结果

通过建模完成了如图1a)所示的二维尾矿坝模型。利用该软件模拟了浸润线抬升后并在保持稳定的情况下，饱和砂颗粒间有效应力减小，抗剪强度减小，最后在重力作用下发生剪切破坏与滑坡。记录了尾矿坝坝体失稳的位移，总体沉降以及边坡破坏特征等结果。在不考虑降雨的情况下，边坡最终破坏形式呈现典型的逆坡牵引式滑坡(见图2c))，破坏形式与前人的室内实验结果基本一致[13]。

如图2a),图2b),图2c)所示分别表示了数值模拟过程中第5时步、第15时步以及第25时步的位移场，每一个步时在软件中迭代计算了23 250次，模拟现实时间0.186 s，由于水位已经上升至破坏时的水位且没有考虑砂土力学性质逐渐减弱的过程，即主要研究坝体破坏的过程，因此破坏发生较快。坡体中不同颜色代表了不同的位移值，具体位移值见各图的颜色条。如图2a)所示，随着坝体开始崩溃，坡角处位移相比坡体后缘与中部明显颜色更深，这意味着坡脚的位移有着明显的滑动趋势，最大水平位移达到1.6 cm左右，坝体的后缘以及中部几乎没有移动；而在竖直方向的位移上，坝体后缘部分相比坝体中部与坡脚有着明显的沉降迹象，最大沉降达到了0.8 cm左右，这是因为后缘部分浸润线上方非饱和土占比相对于坝体中部和坡脚大，浸润线以下的砂体在渗流饱和后其有效应力大大降低，在上部非饱和砂的重力作用下有非饱和砂土更容易被向下压实。

当溃坝模拟进行到第15步时，坡脚的水平移动距离已经达到了4 cm以上，坝体坡脚至坝体中部部位都有较大的水平位移，从坡脚开始向中部逐渐减小，而中部以上至坝体后缘部分水平位移相比中部以下不是很大，但坡顶最大水平位移也达到了0.5 cm左右；在垂直位移方向，如图2b)所示坡顶的沉降位移已经达到了2.5 cm左右，沉降从坡顶向坡脚处开始逐渐减小，在坝体中部沉降为1 cm左右，而坡脚附近的沉降几乎为零。此时坝体已经出现了明显的裂缝，这是因为水平方向互动位移不均导致颗粒间法向更容易断开，再加上坝体沉降不均颗粒间的切向连接也遭到了破坏。

当溃坝模拟进行到了第25步时，如图2c)坡脚水平位移超过5.5 cm，坝体中部的位移也增长到了3 cm左右，随着坝体中部水平位移的增大，原有的两条裂缝间逐渐出现了一条新的张拉裂缝，而坡顶的水平位移达到了1 cm左右；如图2c)所示，坡顶处的沉降位移已经达到了3.5 cm左右。原有裂缝的块体沉降位移差别不大为2 cm左右，但该区域仍然出现了新的裂缝，该裂缝的产生主要是由于水平位移变化不均引起的张拉裂缝。

5 结语

本模拟根据前人总结的实验公式为不同饱和程度的砂土体自动实现了参数的更新。利用离散元软件Matdem建立了尾矿模型，模拟了尾矿库沉降滑坡的细观过程，分析了产生裂隙的成因，得出以下结论:尾矿坝牵引坡脚处水平位移大，逆坡向上水平位移呈递减趋势。尾矿坝沉降位移在边坡顶部最明显，顺坡向下沉降位移逐渐减小。坝顶与坝体中部的沉降差大导致裂缝产生，在沉降差别不大的区域，裂缝产生的主要原因是水平位移不均所导致的。从整个溃坝的过程来看，坝体溃坝有着明显的牵引式滑坡特征，坝脚以及坝体中部水平位移明显大于坡顶位移，牵引了坝顶的移动，沉降不均加速了裂缝的坝体颗粒的断裂，裂缝的出现且变宽主要是由于水平位移不均导致的，但沉降不均同样能导致裂缝的产生。从合位移图可知浸润线上部整体有下滑趋势，坡脚的位移变形最大如若增加初期黏土坝将会大大阻止坝体破坏的发生。MatDEM具有建模简便，计算快速的特点。配合高性能计算机，基于其进一步利用野外边坡的实际数字高程数据、岩土体强度测试数据建立实际工程的离散元模型以研究和预测尾矿坝溃坝灾害的影响范围。同时可以通过本文的方法研究降雨带来的泥石流灾害的危害范围。