流-固耦合作用下某尾矿库稳定性数值模拟研究*

陈艳 雷波 李向阳，2 杨晨

(1.南华大学资源环境与安全工程学院 湖南衡阳 421001；2.湖南省铀尾矿库退役治理工程技术研究中心 湖南衡阳 421001)

0 引言

尾矿库是一个人造的特殊场地，用来存放金属矿山或非金属矿山分选矿石后剩余的残渣或其他废物[1]。影响尾矿库安全的一个重要因素就是水文和地下水。根据不完全统计，在尾矿库溃坝的直接原因中，50%是洪水，20%在于坝体稳定性不足，20%为渗流破坏，剩下10%为其他原因[2]。

流固耦合是研究变形固体在流场作用下的各种行为以及固体位形对流场影响这二者相互作用的一门科学[3]。流固耦合的一项重要特征就是流体与固体之间的相互作用，流体的运动受到固体的影响而发生变化，这种变化又使得固体受到流体荷载的作用而发生变形。薛世峰等[4]研究建立了非混溶饱和两相渗流与孔隙介质耦合作用的数学模型,推导出解耦的渗流系统压力方程，并对流固耦合效应进行分析。沈珠江[5]首次利用次塑形模型及比奥理论，提出了异步交叉迭代显式差分有限元数值解法。一些学者在流固耦合的基础上，结合现代科学的计算机数值模拟技术进行研究分析[6-7]。王晓鸿等[8]建立等效流固耦合的混合数学模型，用数值模拟分析研究了岩土边坡的稳定性。基于ABAQUS软件，荆周宝[9]采用有限元数值模拟研究降雨入渗对土坡稳定性的影响，得出非饱和土质边坡基质吸力、饱和度、变形和稳定性受降雨影响的变化规律。本文利用COMSOL有限元软件，研究南方某尾矿库在不同降雨强度下的稳定性，找到尾矿坝最容易发生事故的部位，为工程实际中尾矿库的治理维护提供参考依据。

1 模型建立

1.1 地下水位的确定

由于降雨入渗，地下水位会发生一定的变化，随降雨时间的持续，地下水位会持续上升。地下水对滑坡发展的各阶段都有很大影响，边坡地下水位上升会导致饱和区域增大、水头压力上升，形成渗透压，引起不饱和区域岩土强度降低，导致边坡变形[10]。

根据Rooso坡地水文模型[11]，可以计算在一定的降雨强度下尾矿库的地下水位。

式中，h为地下水高度，m；h0为降雨初始潜水深度，m；sr为土体饱和度，%；e为土体孔隙比；t为降雨历时，d；T=kz，是沿边坡垂直深度的饱和渗透系数，m2/d，其中k为土体饱和的渗透系数；z为滑坡体厚度，m;θ为坡面与水平面夹角；p为降雨强度，mm/d；a为某网格以上上游集雨面积，m2；b为网格宽度，m。

收集尾矿库所在地区的月降雨资料(见图1)[12]，作为初始数据计算地下水位。分析该地区的月降雨量，选取月降雨量最大的月份分析，即2017年6月的275.5 mm、2014年6月的262.3 mm和2002年7月份的341.4 mm，计算得到地下水位为Hw1=2.83 m，Hw2=2.70 m，Hw3=3.51 m。

图1 月平均降雨量

1.2 模型构建

选取南方某铀尾矿库东北某个坝段，长约335 m。采用渗透性良好岩石的初期坝高约21 m，坝顶宽约4 m，坝体内坡坡比1∶2.0，外坡1∶1.55。尾矿库坝高约42 m，坝坡坡度为1∶4.9，尾矿坝由初期坝和堆积坝构成，材料为土石材料。根据前文计算的地下水位的结果，代入模型，分析尾矿坝的稳定性。

尾矿坝计算模型的力学边界条件为：限定尾矿坝模型两侧的水平位移，模型坝体底部所有位移都固定。地下水渗流条件：尾矿坝底部和左右两侧为不透水边界。此次计算中，将初始应力近似地处理，将尾矿坝本身的自重应力作为初始应力，则此次模拟就是计算尾矿坝在其自身重力作用下的稳定性。

尾矿坝模型如图2所示，区域1为基岩，其材料为砂岩；区域2为堆积坝，材料为尾矿砂；区域3为初期坝，其材料为黏土。具体参数如表1所示。

图2 尾矿坝几何平面

表1 材料参数

2 模拟结果分析

2.1 Von Mises应力分析

Von Mises 应力图用于评价应力的分布，Von Mises应力大表示该处的应力情况不太理想，大多数情况下此处为材料变形较大的地方，也一般是材料容易破坏的地方[13]。

从图3可以看出，堆积坝与基岩接触的尾端是Von Mises应力最大的地方，基岩与尾矿坝上端的Von Mises应力最小。因此，堆积坝与基岩接触的尾端是最容易发生变形和破坏的地方。

(a)hw=2.70 m (b)hw=2.83 m (c)hw=3.51 m

图3 尾矿坝稳定模拟结果——Von Mises应力

比较3张应力图，随着地下水位的变大，尾矿库坝体受到的Von Mises应力会随之变大，应力变化主要集中在堆积坝内部。取初期坝、尾矿坝、基岩三者交汇的点分析可知，随着地下水位的增大，该点处的Von Mises应力逐渐增大。

2.2 滑动面分析

目前，很多边坡稳定性的研究方法都是先假设滑动面的位置，然后计算稳定性系数，最后进行判断[14]。此尾矿坝的滑动体主要在堆积坝部分，其中高度越大，滑动体的深度越深。

将3种情形的结果(见图4)进行比较发现，当hw=3.51时，滑动体体积时最大；当hw=2.70 m时，滑动体的体积最小。取堆积坝、初期坝、基岩3个区域的交点进行分析，得到如下结论：当hw=3.51 m时，滑动面距离基岩越近，滑动面越深，滑动体体积越大，一旦发生溃坝滑坡，危险性是最大的；当hw=2.70 m时，滑动面距离基岩最远。

(a)hw=2.70 m (b)hw=2.83 m (c)hw=3.51 m

图4 尾矿坝稳定模拟结果——尾矿坝滑动面示意

2.3 塑性应变分析

由图5可知，塑性应变主要集中在初期坝和堆积坝上，其中堆积坝与基岩交界的后端塑性变形较大，其次是初期坝的前端，基岩处没有塑性变形。

总体来看，整个尾矿库的塑性应变不超过4×10-3，将三者进行比较，可以发现当hw=3.51 m时，初期坝前端塑性应变较大；当hw=2.70 m时，整体的塑性变形都变小，尤其是在初期坝内，变化较其他区域更加明显。取堆积坝、初期坝、基岩3个区域的交点进行分析发现，hw=3.51 m时，塑性应变较大；当hw=2.70 m时，塑性应变较小，说明地下水位越深，该处塑性应变就越大。

(a)hw=2.70 m (b)hw=2.83 m (c)hw=3.51 m

图5 尾矿坝稳定模拟结果——塑性应变

2.4 位移场分析

选取位移场的xoy平面进行分析。从总位移(见图6)来看，可以发现总体位移最大处发生在坡顶部分，越垂直往下位移越小，初期坝上位移最大在与堆积坝接触的平面。但是从整体来看位移都较小。

(a)hw=2.70 m (b)hw=2.83 m (c)hw=3.51 m

图6 尾矿坝稳定模拟结果——位移

将3种情形得到的结果进行比较发现：当hw=3.51 m时，堆积坝顶端深色范围面积最大，灰色范围面积也最大，当hw=2.70 m时，堆积坝顶端深色范围面积最小，浅色范围面积是3张图中最大的。由此可以得到：随着地下水位的变深，整个尾矿坝的位移变化也变大。取堆积坝顶端的点进行分析发现，当hw=3.51 m时，位移量最大；hw=2.70 m时，位移量最小。

3 结论

(1)通过模型计算分析，Von Mises应力主要集中在尾矿坝与基岩交接处，表明该处为材料变形最大的区域，也是最可能发生破坏的地方。因此，日常的检查、维护重点要注意这些地方有没有发生裂缝或位移。

(2)尾矿库的塑性应力主要集中在初期坝前端部分，降雨量越大，将导致地下水位上升，增加了初期坝前端的塑性，增大了堆积坝部分大体积的土体滑动的可能性。