基于MIDAS的尾矿边坡三维流固耦合稳定性分析

李良振，张延军，秦胜伍，邓浩，倪金

1.吉林大学 建设工程学院，长春 130026；2.中国地质调查局 沈阳地质调查中心， 沈阳 110034

0 引言

矿产资源作为现代工业的一大支柱，在国民经济的发展中具有重要作用。一般矿石开采之后，提取出目前对人类有用的成分，将暂时无法提取的矿物质以及废渣选取合适位置进行集中排放，从而形成尾矿坝[1]。尾矿坝是一种特殊的工业建筑物，其安全运行对矿山的效益、周围人民的生命财产安全以及环境具有重要的影响。许多尾矿库处于非常险要的地理位置，如江河湖泊、重大交通设施、密集居民区和耕地的上游。一旦坝体发生破坏就会造成重大的财产损失和人员伤亡[2]。因此，研究尾矿坝边坡稳定性对于生命安全和企业效益具有重要意义。

在渗流-应力耦合及边坡稳定性方面，国内外大量学者形成了一套完备的计算理论。在渗流方面，Darcy通过砂桶实验，提出了一定流速条件下的达西定律，从而奠定了层流渗流的理论基础[3-4]。Perttson提出边坡稳定性计算方法-圆弧法，在 Taylor和 Fellenius的改进下形成瑞典圆弧法[5-7]。此后，相关学者根据破坏面形状以及划分条块后条间力处理方式的不同，相继建立了 Bishop 法，Janbu 法、morgenstern-price 法等[8]。在流固耦合方面，Terzaghi 提出了有效应力原理，得到了总应力、有效应力以及孔隙水压力之间的关系，从而提出了一维固结模型[9-10]。Zienkiewicz 和 shiomi 在线性固结理论的基础上考虑了非线性问题，改进了 Biot 公式，使得公式可以解决非线性问题[11-13]。李锡肈等考虑到饱和土固结效应的土体和其结构相互作用及相互影响的问题[14-15]，其研究使得土体固结在饱和土和非饱和土方面都有了相应的本构模型。为后续的直接计算和数值模拟提供了相对完备的理论基础和数学模型。

但是，国内外对于尾矿坝的渗透稳定性分析相对较少，且一般只分析渗透情况，对于坝体边坡稳定性只给出定性判定，缺乏定量分析。因此，笔者在渗流分析的基础上，结合吉林集安石墨尾矿坝实际情况和库区工程地质条件，采用有限元方法对该尾矿坝的流固耦合稳定性进行综合分析。

1 尾矿坝三维数值模型建立

本次模型建立基于谷歌大比例尺高清地形云图，通过区域等高线获取22 500个地形控制点，使用曲面拟合的方法重建区域内地面。模型高152 m，长3 250 m,宽1 410 m(图1)。初期坝采用透水堆石坝，坝底标高为505 m，坝顶标高为536 m，坝高31 m。后期堆积坝最终标高648 m，坝高120 m,坡度31.2°。总库容达1 991.3万m3。

图1 尾矿坝三维网格模型Fig.1 3D mesh model of tailings dam

本次模型网格共生成的网格节点121 859个，生成的网格单元181 467个，保证了计算精度要求。由于该地区有较多石墨尾矿坝，且尾矿颗粒成分及组成大致相同。因此通过测量其他石墨尾矿坝的尾矿颗粒及初期坝物理力学参数，得到本模型的物理力学参数，如表1所示。

表1 尾矿坝体各部分物理力学参数

2 尾矿坝渗流分析

2.1 基本假设

(1)假设在渗流过程中，固体颗粒不发生运移，渗流仅仅是水的位移。

(2)假设渗流过程中，渗流速度在一定区间内，即雷诺数(Re)上限在1～10之间，水流为层流，渗流满足达西定律。

(3)假设渗透系数为常数，忽略温度的变化

(1)

式(1)中，ρ为液体的密度，g为地球重力加速度，μ为水的动力粘滞系数，k为渗透率。

2.2 计算原理

根据Darcy定律,饱和砂体渗流方程：

(2)

式(2)中，K为渗透系数，Q为总流量，A为过水面积。H1,H2为水头，L为渗透路径长度。

在三维介质中，饱和渗流微分方程为：

(3)

式(3)中，Kx,Ky,Kz为三个方向的渗透系数，为水头。

在三维空间坐标系中，则有速度分量

(4)

Gd=i×γw

(5)

式(5)中，i为水力坡度，γw为水的容重，Gd为动水压力。

通过水力比降，可以判定渗透力分布特点以及渗透力的大小。

2.3 计算结果

尾矿坝运行过程中，坝顶水位随降水量的变化而变化，主要有干旱水位、正常水位和洪水位3种水位。水位变化影响坝顶干滩长度(坝坡顶到水位线的距离)。根据测量多个临近石墨矿尾矿坝干滩长度随水位的变化情况，得到干滩长度随水位变化情况如表2所示。

表2 干滩长度随洪水位变化

考虑到洪水位时为最危险水位，如果坝体在洪水位时是安全的，则干旱水位和正常水位时一定安全。所以，只需要计算洪水时的渗流情况，并设置安全措施以保证坝体边坡稳定性。

通过计算，得到孔隙水压力分布如图2所示。孔隙水压力的零势面即为浸润线位置。浸润线切过尾矿坝边坡，说明有水溢出。溢出高度41 m，溢出区域长度115±3 m。边坡溢出水会冲刷尾矿坝边坡坡面，对于坝坡稳定性极为不利，如果不设置排渗措施，尾矿坝边坡一定是危险的。

图2 尾矿坝浸润线分布Fig.2 Distribution of infiltration line of tailings dam

图3 尾矿坝应力坡降分布Fig.3 Distribution of stress gradient of tailings dam

应力坡降范围是0～1.247 4 m2·s-2。水力坡降等值区呈现带状分布，应力坡降最大值集中在水头边界处和堆积坝坡脚处，坡脚处为最大值位置(图3)。根据渗透力，可以判定，可能发生破坏的部位在水头边界处和堆积坝坡脚处，其中堆积坝坡脚处破坏将最严重。由于上游水头边界处破坏为尾细砂破坏，而且破坏位置距离坝顶达到125 m，为安全地带。而坡脚处由于渗透力产生破坏则会对尾矿坝边坡产生非常不利影响。

为保证坝体安全，设置排渗体用以改变渗流路径和渗流量，从而提高坝体稳定性。排渗体为砂垫层，渗透性比尾细砂大。每个砂垫层呈水平状，降低了施工的难度(图4)。砂垫层厚1 m，宽约100 m，根据规范要求，结合模型基岩形状和尾矿坝形态，一共设置11块排渗体。排渗体物理力学参数如表3所示。

图4 坝体排渗体布置Fig.4 Arrangement of seepage body of dam body

名称天然含水率/%天然密度/g·cm-3孔隙比比重渗透系数/cm·s-1内摩擦角/°内聚力/kPa承载力特征值/kPa排渗体131.830.5462.491.0×10-22520110

设置排渗体之后，浸润线埋深明显加大，横向埋深增加至101±1.5 m，纵向埋深增加至35±0.5 m，尾矿坝安全性极大提高(图5)。此时，水由初期坝溢出，溢出高度5±0.5 m，符合规范要求。初步判断尾矿坝边坡安全。

应力坡降范围在0～1.100。水力坡降等值区呈现出带状分布(图6)。应力坡降最大值集中在上游水头边界、中间砂垫层以及堆积坝坡脚处，其中，上游水头边界坡降最大，为1.100 m2·s-2。但是由于距离坝顶位置129 m，距离较远，即使发生渗透破坏，对于尾矿坝边坡的稳定性影响不大。

图5 尾矿坝浸润线分布Fig.5 Distribution of infiltration line of tailings dam

图6 尾矿坝应力坡降分布Fig.6 Distribution of stress gradient of tailings dam

3 渗流-应力-边坡耦合

3.1 本构方程

忽略温度场和化学场对应力场的影响，只考虑渗流和应力之间的耦合，方程为：

(6)

式(6)中，p为孔隙水压力，t为时间，s为饱和度，v为流体容量，α为比奥系数，ε为体积应变。

3.2 强度折减法

考虑自重条件下的边坡稳定性，使用强度折减法进行计算。即将坝体材料的强度参数，按照一定比例系数，进行不断折减，直至坝体达到极限破坏状态。与此同时，坝体破坏部位会形成等效塑性应变区，塑性应变贯通的位置即为滑动面(带)位置。

(7)

式(7)中，c′为折减后的内聚力，φ′为折减后的内摩擦角，c为材料初始的内聚力，φ为材料初始的内摩擦角。R为强度折减系数，即最终的安全系数是临界破坏时的强度折减系数。

由于排渗体强度参数与尾细砂强度参数基本一致。为了计算简便，将渗流计算结果直接附到原坝体进行计算。

3.3 计算结果

图7 坝体最小主应力Fig.7 Minimum principal stress of dam body

最小主应力和最大主应力均为负值，说明均为压应力，没有拉应力(图7，图8)。最大剪应力在深度最大处达到最大值，往上呈现层状逐级递减，判断剪应力主要受自重影响比较大(图9)。

图8 坝体最大主应力Fig.8 Maximum principal stress of dam body

图9 坝体最大剪应力Fig.9 Maximum shear stress of dam body

根据塑性应变区分布，可以判定，坝顶处是最大值处，且坝坡上缘与坡脚处塑性应变区连通，形成闭合。三维边坡滑动面与二维边坡滑动面不同，为槽形滑动面。堆积区后部的塑性变形与周围山体边界大致呈平行条带状分布(图10)。判断后缘的变形是由于在相同的重力条件下，山体的强度高变

图10 坝体塑性应变区分布Fig.10 Distribution of plastic strain zone in dam

形量小，堆积砂体较为软弱变形量大，所以引起变形不协调，从而产生与山体接触边缘大致平行分布的塑性变形带。综上所述，塑性变形带由初期坝前缘延伸至堆积区后缘，一旦发生破坏，滑体规模是十分巨大的。由于初期坝强度很高，所以塑性变形小，塑性变形带被初期坝截断。现设计的尾矿坝安全系数为1.536 4，远远大于规范要求，不会发生破坏。

4 结论

(1)尾矿坝边坡中设置的排渗体主要对坝体浸润线有直接影响，设置排渗体之后，浸润线埋深显著增加，达到稳定状态。

(2)排渗体可以改变渗流途径，从而降低坡面以及坡脚的水力坡降，降低了渗透破坏的可能性。

(3)设置排渗体之后，渗流情况发生了变化。计算得到滑动带是与山谷近乎平行的槽状，由初期坝前缘延伸至堆积区后缘。尾矿坝边坡稳定性系数达到1.536 4，符合规范安全要求。