云南某尾矿坝渗流三维计算分析

安晓宇，程 鹏，乐 陶

（1.天津大学，天津300072；2.交通运输部天津水运工程科学研究所港口水工建筑技术国家工程实验室，天津300456；3.中交第一公路勘察设计研究院有限公司，西安710075；4.中钢集团马鞍山矿山研究院有限公司，马鞍山243000）

云南某尾矿坝渗流三维计算分析

安晓宇1，2，程 鹏3，乐 陶4

（1.天津大学，天津300072；2.交通运输部天津水运工程科学研究所港口水工建筑技术国家工程实验室，天津300456；3.中交第一公路勘察设计研究院有限公司，西安710075；4.中钢集团马鞍山矿山研究院有限公司，马鞍山243000）

针对云南某尾矿库蓄水后可能存在的绕坝渗漏问题，在充分分析坝址区工程地质和水文地质特征的基础上，建立坝址区三维有限元计算模型，计算分析二期满库洪水和正常运行工况下地下水的渗流情况。计算结果表明，下游坝坡坡脚处地基渗流溢出处的水力梯度较大，建议在下游坡脚3～5 m范围内地基土层表面布设厚度为0.8～1.0 cm的透水反滤碎石层；同时，由于坝基渗漏量较大，建议在坝体底部铺设土工膜，减少两岸绕渗，降低坝体的浸润线。

尾矿库；尾矿坝；渗流；三维有限元模型

尾矿库是堆存尾矿和澄清尾矿水的专用场地，通常由围护构筑物即尾矿坝（初期主、副坝及后期尾矿堆坝等）和排水构筑物（排水井、管、斜槽及隧洞等）组成。尾矿坝的安全不仅关系到矿山库区自身的安全，而且还关系到下游及周边人民群众的生命财产安全。在尾矿坝建设之前，需要对坝体的渗流稳定性开展研究，分析坝体浸润线的埋深状态等特征，这样可为坝基的抗滑和防渗设计提供参考。使用三维渗流分析的方法，可以更准确的得到库区坝址处的渗流情况。

本文针对云南某尾矿坝工程，为了解其自身及其地基的防渗、抗渗性能，以及坝体的总渗漏量，采用三维有限元方法，对该尾矿坝及地基进行数值计算，从而定量分析评价坝体的总渗漏量和渗透稳定性，并对防止渗透变形提出渗控措施，对尾矿坝的设计、施工提供了参考依据。

1 工程概况

云南某尾矿库位于昆明市西山区海口镇石马哨村南西约1km处石马哨沟中游，距450万t/a选矿厂北东约7.0 km，库区中心距海口镇约7 km，距昆明市35 km。石马哨沟是螳螂川中游右侧的一条支沟，坝址上游有两条分支沟，主沟总长7.0 km，沟谷纵坡约为i=0.024，流向由西南至东北，坝址上游汇水面积4.0 km2；尾矿坝址距上游分水岭2.5 km，坝址距出口螳螂川约4.5 km，距昆晋铁路4.0 km。尾矿坝位于库区下游沟口处，该库在坝顶标高+2 055 m位置时，总库容1 762.02万m3，有效库容1 497.72万m3。

1.1 地质条件

根据勘察报告，坝址处的各岩层渗透性较大。拟建的尾矿库库区内的地层主要有：植物层（Qpd）、第四系冲洪积层（Qel+pl）、第四系残坡积层（Qel+al）、上震旦统东龙潭组（Z2dl）及上震旦统陡山沱组（Z2d）岩土层。场地内分布的地层从上至下主要为：植物层（Qpd）、粉质粘土（Qel+pl）、粉砂（Qel+pl）、粘土（Qel+pl）、粘土（Qel+dl）、全风化泥岩、强风化藻白云岩、强风化石英砂岩、中风化藻白云岩、中风化石英砂岩、强风化白云岩、中风化白云岩、溶洞土等。尾矿坝典型地质横剖面设计图如图1所示。1.2 坝体设计方案

坝体分两期堆筑，二期坝体采用下游式堆筑法，从排水棱体顶部用碎石土沿一期坝贴坡加高，最终坝高80 m，坝顶宽度6.0 m，坝顶标高为+2 055 m，在坝体外坡+1 995 m、+2 010 m、+2 025 m、+2 040 m标高位置设置4条宽2.5 m的马道，内坡坡比为1：2.25，内坡透水层沿一期透水层往上堆筑，外坡坡比为1：2.5，河床处清基面高程1 970 m，坝顶长度约664 m。尾矿坝纵剖面设计图见图2所示。

图1 尾矿坝沿坝轴线地质剖面图Fig.1 Geological profile along dam axis of tailing dam

2 渗流有限元计算基本原理

假定渗流从围岩及尾矿入渗表面到出逸断面的整个渗流区域内为连续水流，且渗流水流充满整个坝体、坝基渗流孔隙区的连续介质，采用部分概化的方法，整个渗流区域内各不同渗透性的岩体尽可能按实际地质分布情况模拟［1-2］。各部分坝体材料视为渗流各向同性体；各部分坝基均视为渗流非各向同性体，水平向与垂直向渗透系数有一定差异（一般水平向渗透系数约为垂直向的3～4倍）。根据坝基地层的渗透性特征，建立计算区域内非均质各向异性渗流模型。

对于研究区内的某一单元，在该单元中的渗流满足下列控制方程［3-4］

式中：θw为体积含水量；kwx、kwy、kwz为土体三个方向的渗透系数，随含水量而变化；H为总水头，是孔隙水压力与位置水头之和（H=uw+ρwgz）；ρw、g分别为水的密度与重力加速度，当忽略水的可压缩性时，ρw为常量。

当岩土体不饱和时，体积含水量θw随孔隙压力变化的函数关系可通过土水特征曲线（SWCC）来确定。Van Genuchten提出了下列表达式（VG模型）

式中：a、n为试验拟合参数；Se为有效饱和度；θs、θr分别为土体的饱和与残余体积含水量；uw为孔隙水压力，为负值。

当岩土体不饱和时，其渗透系数随含水率（饱和度）变化，可由式（3）计算

式中：Kr为相对透水率，KS为饱和渗透系数。

VG模型中相对透水率Kr与有效饱和度Se的关系如式（4）所示

对应控制方程式（1）的定解条件可以是

图2 尾矿坝纵剖面设计图Fig.2 Longitudinal design of tailing dam

式中：H(x,y,z)表示区域内任一点(x,y,z)总水头；φ0(x,y,z)和q0(x,y,z)为已知函数；n为边界Γ2外法线方向。

由控制方程（1）和辅助方程（2）～（4）以及定解条件（5）、（6），共同构成描述饱和-非饱和岩土体稳定-非稳定渗流场的数学模型。该模型可统一考虑饱和-非饱和区，计算得到孔隙水压力为零的点对应在浸润面上。该模型既可用于计非稳定渗流场，也可用于计算稳定渗流场：当式（1）的左端项为0时，该方程退化描述稳定渗流控制方程。

3 计算模型与参数

3.1 计算程序简介

渗流计算采用河海大学岩土工程研究所研制的US3D三维有限元程序。该程序有如下主要特点：程序既可用于饱和岩土体渗流计算也可用于非饱和岩土体渗流计算；既适用于稳定渗流计算也适用于非稳定渗流计算；无需变动网格即可计算得到浸润线和逸出点；可得到计算区域内的渗流水头分布、水力梯度和渗流力；可模拟土体饱和度变化和渗透系数的非线性；可根据计算要求计算通过任意给定截面的渗流量。

3.2 三维有限元计算模型

计算中X轴为沿坝轴线方向从左岸指向右岸方向，零点定在左坝肩处；Y轴为垂直坝轴线从上游指向下游方向，零点在坝轴线处；Z轴竖直向上，零点与高程零点相同。模型计算区域确定为X轴方向从左右坝肩各向外侧60 m范围；Y轴方向从坝轴线向上下游分别各取约2倍坝底宽；Z轴方向坝底清基面以下约2倍坝高，即计算模型的长宽高尺寸为：1 200×784×255 m。模型沿直坝轴线划分41个剖面，剖面间距约为20 m，局部区域进行加密。为了能尽可能模拟地层的变化情况，计算模型共设置16个材料分区。模型中主要采用8结点六面体等参单元，局部用6结点五面体或4结点四面体单元过渡，共划分23 554个单元，23 305个结点。整体三维有限元网格见图3。

3.3 渗透系数的确定

图3 三维有限元网格整体效果图Fig.3 Three⁃dimensional numerical calculation model and meshing diagram

某尾矿坝渗控分析中，必须采用真实反映实际地质条件和坝体填筑情况的坝体坝基各介质的渗透系数。在参照云南化工院和武汉勘察研究院出具的《某尾矿库岩土工程勘察报告》和长沙冶金设计研究总院的《某尾矿库坝体稳定性研究——尾矿坝固结渗流分析模型研究报告》中尾矿渗透系数数量级基础上，结合类似工程的基本特点，拟定了头石山尾矿坝体坝基的渗流分析参数，汇总于表1中。

表1 某尾矿坝渗透系数汇总表Tab.1 Permeability coefficient of tailing dam

3.4 容许水力梯度的确定

库区地表出逸处第四系土体的临界水力梯度可根据《土工原理与计算》［5-6］中所给出下式进行计算

式中：n为土体孔隙率；Gs为土粒比重；ξ为侧压力系数；φ为土的摩擦角；c为单位土体所受的凝聚力；γw为水的容重。

根据地勘报告资料，地表土体n孔隙率为0.4～0.5，土粒比重为2.72～2.74，侧压力系数取0.5，土体摩擦角为19～25°，凝聚力为30～40 kPa。从而可推算得地表土体临界水力梯度为11～16。若取容许安全系数为2～3，则地表出逸处土体容许水力梯度为3.7～8。

3.5 计算方案

经过讨论，确定共对8种方案进行了计算。计算结果表明，由于两岸基岩的渗透性很大，如果岸坡段坝体底部不设置土工膜，则坝体内浸润线较高，对坝体稳定不利。因此，建议岸坡段坝体底部也设置土工膜。所以，将岸坡段设置土工膜作为基本方案进行计算分析。由于篇幅有限，在本文中列出其中3种方案的计算结果（表2）。

表2 某尾矿坝渗流计算分析方案Tab.2 Seepage calculation scheme of tailing dam

4 渗流计算结果及分析

4.1 二期满库洪水运行工况计算结果

图4给出了二期满库洪水运行工况计算结果（A1方案）河谷处最大横剖面（坝0+320）总水头等值线分布情况。可以看出，受土工膜隔渗作用［7-8］，浸润线在坝体迎水面处转折，浸润线从上游棱体与土工膜间沿水平排水层达到坝体下游侧。表明受坝体迎水面及坝底土工膜包裹隔渗，加上坝体底部水平排水层及下游排水棱体作用，坝体内浸润线保持在坝基附近。在土工膜和下游排水棱体附近，总水头等值线分布密集，表明受土工膜隔渗和排水棱体排水作用，相邻区域水头损失较大。

图5为A1方案河谷最大横剖面（坝0+320）水力梯度等值线图。从坝体内各点水力梯度分布情况来看，水力梯度最大值发生在河谷处下游坡脚排水棱体附近的地基土中，最大水力梯度达1.12。尽管水力梯度大于0.8的单元不多，但还是应引起注意。

由于河谷段2 020 m以下的下游坡脚处的水力梯度较大，建议在下游坡脚3～5 m范围内地基土层表面布设透水反滤碎石层，厚度0.8～1.0 cm。排水棱体中的水力梯度较小，最大值近0.31。尾矿砂中的水力梯度大都小于0.3。水力梯度均不大，不会发生渗透破坏。此方案的渗漏量为95 500 m3/d。渗漏量较大，主要是两边岸坡渗透系数较大的缘故。

4.2 二期满库正常运行工况计算结果

图6给出了二期满库正常运行工况计算结果（A2方案）河谷处最大横剖面（坝0+320）总水头等值线分布情况。

与A1方案计算结果（图4）比较，坝体坝基内渗流规律基本相同。A1方案的上下游最大水头差约84 m，A2方案的水头比A1方案水头减小幅度为1.8%。

图4 A1方案河谷处最大横剖面总水头等值线（单位：m）Fig.4 Total head contour map of the valley maximum horizontal section in A1 scheme

图5 A1方案河谷最大横剖面水力梯度等值线Fig.5 Hydraulic gradient contour map of the valley maximum horizontal section in A1 scheme

图6 A2方案河谷处最大横剖面总水头等值线（单位：m）Fig.6 Total head contour map of the valley maximum horizontal section in A2 scheme

图7 A2方案河谷最大横剖面水力梯度等值线Fig.7 Hydraulic gradient contour map of the valley maximum horizontal section in A2 scheme

图8 A2方案河谷处最大横剖面总水头等值线（单位：m）Fig.8 Total head contour map of the valley maximum horizontal section in A5 scheme

图9 A5方案河谷最大横剖面水力梯度等值线Fig.9 Hydraulic gradient contour map of the valley maximum horizontal section in A5 scheme

图7为A2方案河谷最大横剖面（坝0+320）水力梯度等值线图。河谷断面下游出口处的水力梯度相比A1方案减小为1.07，幅度约4.5%。和A1方案相比，A2方案的坝体与岸坡内水力梯度略有减小。

A2方案为正常运行工况，其上游水位比A1方案洪水运行工况的水位低1.5 m。计算结果显示：当上游水位略有下降时，坝体渗漏量略有下降；从A1方案的95 500 m3/d减小到92 400 m3/d，减小幅度3.2%。

4.3 二期满库洪水+降雨运行工况计算结果

图8给出了二期满库洪水+降雨运行工况计算结果（A5方案）河谷处最大横剖面（坝0+320）总水头等值线分布情况。图9为A5方案河谷最大横剖面（坝0+320）水力梯度等值线图。与二期满库洪水运行工况（A1方案）计算结果（图4）比较，坝体坝基内渗流规律基本相近。水力梯度比A1方案的计算结果略大，其中，河谷段地基土层中的最大水力梯度增大8.0%。

比较A5和A1方案，结果显示，总体上降雨作用对坝体渗流场影响不大，坝坡及岸坡水力梯度均不大，大部分区域小于0.5。

5 结论

根据前述计算结果，得出如下结论与建议：（1）采用拟定的基本计算参数和各种工况条件下，某尾矿库渗流场符合一般规律。坝体和地基中大部分区域水力梯度小于0.5，远小于其容许水力梯度；但在河谷段高程2 020 m以下的下游坝坡坡脚处地基渗流溢出处的水力梯度较大，建议在下游坡脚3～5 m范围内地基土层表面布设厚度为0.8～1.0 cm的透水反滤碎石层。（2）尾矿库总渗漏量主要取决于地基（风化石英砂岩及白云岩）的渗透性大小，坝基渗漏量较大，为了减小总渗漏量，需在坝体底部铺设土工膜，这样可以有效的减少两岸的绕渗，降低坝体的浸润线。

［1］Van Genuchren.A closed⁃form equation for predicting the hydraulic conductivity of unsaturated soil［J］.Soil Sci.Am.J.，1980，44 （5）：892-898.

［2］齐清兰，张力霆.尾矿库渗流场的数值模拟与工程应用［M］.北京：中国水利水电出版社，2011.

［3］胡云进.裂隙岩体非饱和渗流分析及其工程应用［M］.杭州：浙江大学出版社，2009.

［4］骆祖江，张弘，李会中，等.乌东德水电站坝址区地下水渗流三维非稳定流数值模拟［J］.岩石力学与工程学报，2011，30（2）：341-347. LUO Z J，ZHANG H，LI H Z，et al.Three⁃dimensional Transient Flow Numerical Simulation of Ground Water Seepage in Dam Site of Wudongde Hydropower Station［J］.Chinese Journal of Rock Mechanics and Engineering，2011，30（2）：341-347.

［5］殷宗泽.土工原理［M］.北京：中国水利水电出版社，2007.

［6］钱家欢，殷宗泽.土工原理与计算：二版［M］.北京：中国水利水电出版社，1996.

［7］SL/T231-98，聚乙烯（PE）土工膜防渗工程技术规范［S］.

［8］GB/T17643-2011，国家标准：土工合成材料聚乙烯土工膜［S］.

Analysis of 3D seepage of a tailing reservoir in Yunnan

AN Xiao⁃yu1,2,CHENG Peng3，YUE Tao4
（1.Tianjin University,Tianjin 300072,China;2.Tianjin Research Institute for Water Transport Engineering,National Engineering Laboratory for Port Hydraulic Construction Technology,Tianjin 300456,China;3.CCCC First Highway Consultants Co.,LTD.,Xi′an 710075,China;4.Sinosteel Maanshan Institute of Mining Research Co.,LTD., Maanshan 243000,China）

To investigate the leakage problem below a tailing dam in Yunnan province,a three dimensional fi⁃nite element model was developed according to the hydrogeological and engineering geological characteristics.The effects of groundwater seepage under full⁃reservoir and normal condition were studied and analyzed respectively. The calculated results show that the foundation′s seepage overflow point at downstream slope toe presents higher hy⁃draulic gradient.Thus the pavement of permeable filter gravel layer(0.8-1.0 cm)above 3-5m range of foundation soil at downstream slope toe was proposed.It is suggested that the geotechnical membrane should be laid at the bot⁃tom of the dam to reduce the seepage around both sides and lower the soaking line of the dam.

tailing reservoir;tailing dam;seepage;three dimensional finite element model

TU 46;O 242.21

A

1005-8443（2016）06-0641-05

2016-06-01；

：2016-07-11

安晓宇（1988-），男，天津市人，工程师，主要从事环境岩土工程研究。

Biography：AN Xiao⁃yu（1988-），male，engineer.