山西宝山矿业有限公司

小地沟尾矿库渗流与稳定研究报告

 

 

 

 

 

 

 

 

 

 

大连理工大学水利工程学院

20122

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

前言

山西省繁峙县宝山矿业有限公司位于山西省忻州市繁峙县岩头乡官地村。该公司选矿厂所在地区地形崎岖,山多川少,地质构造较为复杂,形成山高坡陡,沟谷狭窄的地形条件,找到较大的山沟作为存贮尾矿的尾矿库库址十分困难。宝山矿业在选厂以东选定小地沟为尾矿库址,条件适宜。该沟三面高山环抱,沟口窄小,沟的纵深较长,主沟长约600m,主沟的上游段有两条岔沟,岔沟窄长,地势高容积不大。沟的两岸山体较高,坡度较陡,沟的断面呈“V”型地形。小地沟做为贮存尾矿的尾矿库,在生产初期库容量小,尾矿堆积上升速度快,随着尾矿堆积高度的增高,库容量也随之增大时小地沟的地形特点,显然,利用这个沟作为贮存尾矿的尾矿库,并且要求尾矿库的服务年限10年以上,势必出现尾矿高堆才能多装方能满足尾矿库服务年限的尾矿贮存量的设计方案。在7度地震区同时采用“上游式”筑坝,基于目前我国尾矿筑坝的设计和试验研究技术水平,采用“上游式”筑坝尾矿高坝的稳定和所采用的对策都能得到解决。

若满足宝山矿业的尾矿库服务年限,依小地沟的地形条件必须尾矿高堆,高堆近200m可容纳近20年的尾矿排出量。按规范大、中型选矿厂的尾矿库的使用年限不宜少于10年,使用小地沟作为贮存尾矿的尾矿库才有其经济意义和使用价值。

尾矿库采用尾矿自身材料堆积贮存,利用尾矿自身材料的物理力学性能维持尾矿堆积坝体的稳定性,是尾矿堆积坝的构造准则。尾矿筑坝采用“上游式”筑坝,坝上小流量多管分散放矿充填。尾矿堆积坝外坡平均坡比采用1:4.2。为防止尾矿堆积坝外坡尾矿裸露面粉尘污染环境,在每年尾矿筑坝完成后,在尾矿堆积外坡面铺0.3m 厚山坡土防尘、防止雨水冲刷坡面。并在坡面种草、草皮或植紫穗槐、沙棘等灌木。

尾矿堆积没上升10m 平台,平台宽1.5m ,平台上做纵横排水明沟。由于尾矿筑坝的堆积速度初期生产上升快,中、后期减低,平均每年上升高度在14.4m 左右。对尾矿坝的稳定性影响很大,为此,设计采取RCP渗垫进行尾矿排水固结的工程措施,即在尾矿堆坝过程中,每筑坝升高10-15m 铺一层RCP排渗垫和渗水排水管,进行加快尾矿排水固结,防止尾矿堆积上升速度快而发生尾矿脱坡现象。

尾矿堆筑分为初期堆积和后期堆积,初期堆积到标高1324m ,尾矿库服务年限为5年。后期堆积到标高1400m

本报告包括尾矿坝典型断面二维渗流数模计算与整体三维渗流数模计算、静力稳定分析以及坝体抗震稳定与动力反应分析等内容。针对两种典型标高的尾矿坝,在正常水位运行、洪水水位运行以及特殊运行工况下进行了渗流计算与稳定分析,同时考虑到坝体内的水平排渗垫有效与否对浸润线与坝体稳定的影响。

1 渗流的数值模拟计算

       渗流数值模拟计算的目的是提供坝体浸润线的位置,了解坝体的渗流情况,为坝体的设计及稳定性研究提供依据。渗流模拟计算采用固定网格有限元法计算模式。采用有限元方法对宏达矿业小地沟尾矿坝进行了渗流数值模拟计算。

1.1 渗流数值模拟计算方法

1.1.1渗流数值模拟方法

Darcy定律提出以来,对多孔介质渗流问题研究已比较充分,并建立了较为完善的多孔介质渗流理论与计算模型。综合国内外的研究现状,现今渗流计算模型主要有:等效连续介质模型、离散裂隙网络模型、裂隙-孔隙双重介质模型和随机渗流模型,另外,多场耦合模型研究和应用也发展迅速。从工程应用的实际出发,等效连续介质模型有许多成熟的经验可以借鉴,计算参数的准备相对较少且较易获取,但对于许多工程中的一些关键导水设施是不可以简化和忽略的,在土石坝渗流的研究中,结合土石坝本身的特性,大多采用等效连续介质模型。因此本文采用等效连续介质模型,模拟对渗流场分布起主要控制作用的渗流控制设施、断层及软弱结构面,对土石坝渗流问题进行研究。根据饱和、非饱和渗流理论建立了以测压管水头和饱和度为变量的饱和-非饱和渗流数值模型,该模型能够处理土石坝渗流、降雨入渗等工程渗流问题。

基本方程

质量守恒方程:

设流体的密度为r,土体的孔隙率为n,饱和度为S,对于一般渗流,由质量守恒定律可导出

                                            1-1

其中方向的渗流速度,对于二维问题i=1,2,对于三维问题,i=1,2,3

运动方程:

考虑各向同性渗流,其运动方程满足达西定律

                                                              1-2

K为渗透系数,H为水头函数,其定义为

                                               1-3

这里x2铅直方向的高度,p为流体压强,r为液体容重,h为压强水头。

非饱和渗流渗透系数的本构关系:

渗透系数一般是饱和度的函数,Gardner最早提出了渗透系数K与饱和度之间为幂函数关系:

                             1-4

其中,为介质饱和渗透系数,式中指数m,根据试验结果为2~3

在饱和区S=1,将运动方程带入连续性方程(1-1),对于不可压缩流体得,

                     1-5

上式即饱和渗流的Laplace方程。

在非饱和区,当空隙与大气相连通时,相对压强为零,由(1-3

可以把多孔介质的渗流域分成饱和区与非饱和区两部分,分别考虑饱和区饱和度为1和非饱和区的情况,(1-1)式可以退化为如下的分区方程:

         1-6

 

式(1-6)为描述饱和-非饱和渗流的一般微方程式。对于饱和渗流,式(1-6)右端为零。联立求解(1-6)可确定整体渗流区水头函数H及饱和度度函数S的空间与时间变换规律,进而利用达西定律(1-2)确定渗流速度场。

上述双变量述饱和-非饱和渗流模型,可以避免带自由面渗流问题中需要不断调整自由面(线)变计算域的困难。渗透自由面(线)可通过确定饱和度函数S阶跃突变而自动捕获。

方程的求解需要已知定解条件,包括初始条件和边界条件。

1)初始条件

时刻,初始条件为坐标的函数,可以写成:

                     1-7

2)边界条件

边界条件包括水头边界和逸出面边界。将饱和-非饱和区统一起来形成总的渗流研究区域,无需考虑自由面流量补给边界,也不存在自由面位置迭代问题[106]

水头边界

                         1-8

透水边界

                         1-9

逸出面边界

                         1-10

为给定边界水头,为沿边界的外法线方向,边界点坐标,逸出边位置水头。

 

若将(1-2)直接代入(1-1),得

                        1-11

 

nCs一般为渗流压强的函数,对于刚性骨架,n=Const,且,故式可改写为

                                              1-12

其中                           

为贮水率。(1-12)为地下水渗流的基本方程。目前,在国内外的渗流数值模拟研究中,很多采用的都是以水头(或体积含水量)为变量的单变量抛物型方程,而公式中的其它参数则多采用以水头为变量的函数来近似表达。单变量抛物型方程的优点在于公式变量统一求解方便,但是贮水率参数与水头函数关系不明确,基本以经验公式拟合或实验测定值描述为主;

1.1.2有限元离散与求解

将求解域划分为若干个有限单元。设单元插值函数为,则单元结点测压管水头和饱和度可表示为

                            1-13

对公式(1-6),时间采用隐式差分格式做离散,空间采用Galerkin有限元法做离散,得:

          1-14

其中为权函数,表示m单元积分后求和。经分部积分并考虑边界条件,得

 1-15

化简并整理得方程组:

                      1-16

其中

                        1-17

                         1-18

                             1-19

                          1-20

以上给出了本文饱和-非饱和渗流数值模型的有限元计算格式,问题的解可由满足时的解得到,其中表示范数,为控制精度,且 。上述计算格式采用固定网格法可避免自由面调整时网格重生成时所面临的困难,采用隐式时间推进格式求解可以在网格尺度相差较大的情况下,采用较大的时间步长仍能保证数值解的收敛性。

在饱和-非饱和渗流数值模型中,将渗流域分为饱和区和非饱和区,其中非饱和区中忽略了毛细作用的影响,使得测压管水头等于位置水头,因此非线性抛物型渗流控制方程在饱和区与非饱和区分别退化为关于测压管水头的拉普拉斯方程和关于饱和度的对流方程。其中,拉普拉斯方程采用预处理共轭梯度法(PCG)求解;对流方程采用GaussSeidel方法求解。

1.2 小地沟尾矿库二维渗流计算结果

1.2.1坝体断面材料属性图以及计算二维网格剖分图

坝体二维渗流计算分别选取尾矿堆筑典型标高1330m(初期)、1400m(后期)断面作为计算域,参照尾矿地层特征,偏于安全起见,将地层做不透水层处理并根据筑坝和尾矿堆放沉积特点,计算域共划分五种不同的材料类型见图1-1、图1-3各典型断面计算网格剖分结果见图1-2、图1-4,网格尺度为5m,能准确反映尾矿堆积体的渗流特征。

1.2.2坝体断面材料参数选取

    坝体各种材料渗透系数依据材料属性并参看往期邻近尾矿场确定,见表1-1

1-1  尾矿断面材料渗透系数

编号

土料

渗透系数Kcm/s

滤水堆石

5.0×10-2

尾细砂

3.0×10-4

尾粉砂

5.0×10-5

尾粉土

8.0×10-6

尾粉质粘土

3.0×10-6

1.2.3二维渗流计算工况表

1-2  小地沟二维渗流浸润线计算结果一览表

投产期别

堆筑标高

等级

演算工况

 

图号

初期

1330m

正常运行水位

250干滩)

渗有效

1-5

排渗失效

——

洪水位(100干滩)

渗有效

1-7

排渗失效

——

后期

1400m

正常运行水位

250干滩)

渗有效

1-9

排渗失效

1-11

洪水位(100干滩)

渗有效

1-13

排渗失效

——

 


 

1-1  小地沟尾矿库堆筑至1330高程断面材料属性图

 

 

 

1-2  小地沟尾矿库堆筑至1330高程网格剖分图

 

 

 

 

 

 

 

1-3  小地沟尾矿库堆筑至1400高程断面材料属性图

 

 

 

1-4  小地沟尾矿库堆筑至1400高程网格剖分图

 


 

1-5  小地沟尾矿库堆筑至1330高程正常水位排渗有效浸润线位置图

 

 

 

 

 

 

 

1-6  小地沟尾矿库堆筑至1330高程正常水位排渗有效等水头图

 

 

 

 

 

 

 

 

1-13  小地沟尾矿库堆筑至1400高程洪水位排渗有效浸润线位置图

 

 

 

 

1-14  小地沟尾矿库堆筑至1400高程洪水位排渗有效等水头图

 

1-7  小地沟尾矿库堆筑至1400高程洪水位排渗有效浸润线位置坐标

点号:

1

2

3

4

5

6

7

8

9

10

X:

0.00

131.60

135.10

135.72

139.54

152.77

177.21

216.50

308.73

402.33

Y:

1212.00

1218.70

1223.16

1229.00

1233.93

1243.13

1250.75

1257.72

1275.46

1293.65

点号:

11

12

13

14

15

16

 

 

 

 

X:

512.49

563.50

705.06

788.64

804.47

804.70

 

 

 

 

Y:

1315.42

1327.10

1362.11

1389.43

1397.67

1399.00

 

 

 

 

 

 

1.3 小地沟尾矿库三维渗流计算结果

小地沟尾矿库三维渗流分析采用实体建模方法,所建模型如图1-15所示,模型共划分750,062个网格单元和138,304个网格节点,如图1-161-17所示。三维渗流分析计算了坝体堆筑高程至1400 m,洪水位工况下的渗流情况,所取参数参照二维分析选取。

渗流计算结果三维渗流浸润面,如图1-18所示和三维渗流等水头图,如图1-19所示。为便于了解内部浸润面和等水头情况,采用切片处理方式,浸润面和等水头切片显示分别如图1-20和图1-21所示。

 

1-15  小地沟尾矿库三维模型示意图

 

1-16  小地沟尾矿库三维模型网格剖分示意图

 

1-17  小地沟尾矿库三维模型网格剖分示意图(局部图)

1-18  小地沟尾矿库三维渗流饱和度总体示意图

 

1-19  小地沟尾矿库三维渗流等水头总体示意图

1-20  小地沟尾矿库三维渗流饱和度切片示意图

 

 

1-21  小地沟尾矿库三维渗流等水头切片示意图

 

 

通过对小地沟尾矿库堆筑高程至1400 m洪水位工况下的三维渗流分析研究,并结合二维典型断面渗流的模拟结果来看,三向模拟总体上符合二维模拟的结论,但考虑到尾矿库安全等别较高,为控制浸润线高度,保证尾矿坝的安全运行,保证各级水平排渗设施有效畅通,严格控制干滩长度,保证初期坝堆石排渗体畅通是十分重要的。

2 坝体静力稳定分析

为了综合评价坝体在各种工况运行的静力稳定性,采用了两种方法进行静力稳定分析。第一种方法为规范规定的圆弧滑动法,稳定性校核将以此法为主。第二种方法为有限元应力应变分析,此法近年来开始应用于坝体的稳定性分析中,本报告将应用此法对坝体的应力变形状态,剪应力水平及整体稳定性进行分析。

2.1 圆弧滑动法静力稳定分析

2.1.1分析方法原理简述

土体的滑动与否决定于促使土坡的运动的滑动力与滑动面上的抗滑力这一对矛盾抗衡的结果。经典的滑弧稳定分析假设土体的应力应变关系为刚塑性的,即滑动土体为刚性体,滑动面为塑性曲面,然后取土体为脱离体,分析其在各种力作用下的稳定性。

对于外形比较复杂,特别是土坡由多层构成时,要确定滑动土体的重量及其重心位置就较复杂。滑动面上的抗剪强度又分布不均,而是与各点的法向应力有关。故在土坡的稳定性分析中,常将滑动土体分成若干垂直土条,求各土条对滑弧园心抗滑力矩和滑动力矩,分别求其和,然后求出土坡的稳定安全系数。采用瑞典圆弧法,按土石坝设计规定稳定安全系数的总应力法表达式为

         K=                  2-1

式中:第i条土滑裂面长度

       :第i条土滑裂面上的粘聚力与磨擦角

       :第i条土条实重

       :浸润线以上的条块实重

       :浸润线以下,静水位以上的饱和重

       :静水位以下的浮重

        i条土条滑裂面中心处的孔隙水压力(由浸润线算至条底)

        :第i条土条滑裂面的倾角

2.1.2体材料的力学指标

2-1给出了计算中采用的坝体土料的物理力学指标。

 

 

2-1  体材料物理力学指标

编号

土料

粘聚力

(水上)C(KN/m2)

内磨擦角(水上)j()

粘聚力

(水下)C(KN/m2)

内磨擦角(水下)j()

天然容重g(KN/m3)

容重g(KN/m3)

饱和容重gsat(KN/m3)

滤水

堆石

0

38

0

36

20

10

20

尾细砂

7.84

33

7

30

18

8.5

18.5

尾粉砂

9.8

30

8.8

28

18.5

8.5

18.5

尾粉土

9.8

28

8.8

27

19.5

10

20

尾粉质粘土

10.78

16

10

15

19

9

19

 

二维静力稳定计算工况以及计算结果见下表2-2

 

2-2  小地沟二维各工况静力稳定安全系数一览表

投产期别

堆筑标高

等级

演算工况

 

图号

K

初期

1330m

正常运行水位

250干滩)

渗有效

2-3

1.83

排渗失效

——

——

洪水位(100干滩)

渗有效

2-4

1.44

排渗失效

——

——

后期

1400m

正常运行水位

250干滩)

渗有效

2-5

1.43

排渗失效

2-6

1.40

洪水位(100干滩)

渗有效

2-7

1.33

排渗失效

——

——

 

 

2.1.3计算结果

2-2小地沟各运行工况的静力稳定安全系数表,图2-1、图2-2为稳定计算条分示例图,图2-32-7为各工况下静力计算的滑弧位置示意图,图中标明了滑弧的圆心位置,滑弧半径以及稳定安全系数。

 

 

 

 


 

 

 

2-1 小地沟尾矿库堆筑至1330高程条分示例图

 

 

 

 

 

2-2  小地沟尾矿库堆筑至1400高程条分示例图


2-3  小地沟尾矿库堆筑至1330高程正常水位排渗有效静力稳定滑弧计算结果

 

2-4  小地沟尾矿库堆筑至1330高程洪水位排渗有效静力稳定滑弧计算结果

 

 

 

 

 

2-5  小地沟尾矿库堆筑至1400高程正常水位排渗有效静力稳定滑弧计算结果

 

2.2 尾矿坝静力有限元分析

尾矿坝是应用很广的一类坝型。在该类坝设计中,除了分析稳定性外,一般还须分析其应力和变形[6] 。稳定性的丧失由于局部地区应力达到了强度,变形增长,而逐步发展的。此外土体内局部出现拉应力的趋势,亦会导致裂缝的存在。通过对土体内应力变形的全面分析,可以更好地分析坝体的稳定性和估计裂缝的发展。尾矿坝静力分析的主要目的在于:(1)根据计算所得的剪应力水平,帮助判断坝的稳定性及可能的失稳面,并检查是否出现拉应力。(2)提供震前平均应力分布,为动力分析中的动模量的计算提供基本数据。(3)根据主应力比,为进一步的液化危险性分析做准备。

 

2.2.1固体力学有限元法

尾矿坝一般可作为连续介质考虑,其应力和变形过程可以通过平衡律(守恒律)与本构关系进行力学描述。这种描述在数学上可以表示为一组偏微分方程。广义的有限元法可以在数学上被理解为求解偏微分方程弱解的近似方法,各类加权剩余法及泛函形式都可作为弱解有限无格式的特例。对于固体力学中的问题,有限元离散可以等价地通过变分原理导出坝体的有限元应力——应变分析,就是把连续的坝体用网格划分为有限数目的单元体。这些单元体之间在结点相互铰接,形成离散结构,用这离散结构来代替原来的连续结构,进行分析求解。

设单元的位移,可由相应的结点位移来{u}表示。

                              2-2

在此[N]为插值函数。单元内的应变为

          2-3

这里[B]代表应变矩阵。

设应力、应变关系为

                                  2-4

D为弹性矩阵,它决定于弹性常数。

应变能为

             2-5

若不考虑集中外力,外力所作的功可表示为

                  2-6

其中为体力,为面力,S为物体上面力分布的面积。

物体的势能可表示为

    2-7

可导出

         KU=F                             2-8

其中

                              2-9

                           2-10

K为刚度矩阵,F为载矩阵。

 

2.2.2本构关系

为了描述土体受力后的行为,必须给出其本构关系,即应力与变形之间的物理规律。实际土体的应力——应变关系非常复杂,它具有非线性,非弹性,压弹性和剪胀性等等特性,并受到土体结构,应力历史,孔隙水及时间效应的影响。目前,人们提出了多种弹——塑性模型来描述土的本构关系,其中在我国应用较多的方法是——张等提出的非线性弹性模型[8]

根据肯等人的建议,三轴试验的应力与应变可近似为双曲线关系。如果用切线模量表示,有如下所述的模式E—B模式。

 

2.2.2.1 模式

在此模式中,土体的本构关系,可用切线杨氏模量和泊松比来描述

                       2-11

                    2-12

其中初始切线模量A

                              2-13

                             2-14

式中Kn为试验确定的参数,为大气压,单位与相同,以便使K成为无因次量。表示破坏比。为应力水平。根据摩尔——库仑破坏准则。得

       2-15

其中C为土的凝聚力和内摩擦角。

上述各式中共有8个参数,即KncFGD,都是由三轴试验确定的。

增量应力——应变关系在平面应变条件下按弹性理论为

    2-16

 

2.2.2.2 模式

模式中的切线泊松比与实验实测数值有时不在符合,所以有人建议取消泊松比而引入体变形模量B。本文静力计算采用E—B模型。应用E—B模型时,的定义与模型相同,体积变形模量为

                                      2-17

这样,模式中需要试验提供七个材料参数: Kncm

    对于平面应变问题,E—B模型的应力应变增量形式为

    2-18

2.2.3分析步骤

由于上述本构关系是非线性的,土体的有限元应力应变分析属于材料非线性问题。几何非线性问题较为复杂,通常在工程问题中较少考虑。

求解非线性总是的方法可分为三类,即迭代法,增量法和混合法。由于荷载是逐步级施加的,增量法可用于模拟加工过程,因而应用较广。在本报告中采用了精度较高的中点增量法,其计算步骤是,首先施加第三步荷载增量的一半进行计算,求得中点的刚度矩阵,然后再应用此中点刚度矩阵求出第三步的位移增量。中点刚度法相当于预估校正法,刚度矩阵的计算需要采用切线模量。

 

2.2.4计算方案

本次有限元应力应变分析选择尾矿坝后期(坝顶高程1400m)在洪水运行排渗有效时的工况进行计算分析。

分析是按全坝分成若干施工层进行载荷增量计算。采用四结点(可退化为三结点)等参元。底部面为固定边界,坝面为自由边界。

划分的有限元网格如图2-13。在高度向上的方向共分20个施工层,共计1043个单元,1120个结点。

模拟全坝的材料性质,各单元赋予不同的材料号,分别模拟初期坝筑坝料,尾矿料,及堆石体等。根据浸润线的位置,某种材料的水上部分与水下部份被赋予不同的材料标号,以模拟该材料具有不同含水量时力学参数的变化及渗流力的影响。本次静力分析所选材料的参数见表2-4,这些材料参数系根据现有的实验资料与以往工程经验综合确定。

2-4  有限元静力计算参数

土料

计算容重(KN/m3)

粘聚力C(KN/m2)

内磨擦角j()

破坏比Rj

杨氏模量K

杨氏模量n

体积模量Kb

体积模量m

滤水堆石

(水下)

10

0

36

0.8

200

0.32

130

0.2

尾细砂

(水下)

8.5

7

30

0.73

106.0

0.85

20.0

0.83

尾粉砂

(水下)

8.5

8.8

28

0.45

75.0

0.80

118.0

0.77

尾粉土

(水下)

10

8.8

27

0.45

75.0

0.80

118.0

0.77

尾粉质粘土(水下)

9

10

15

0.45

75.0

0.80

118.0

0.77

滤水堆石

(水上)

20

0

38

0.8

200

0.32

130

0.2

尾细砂

(水上)

18

7.84

33

0.73

106.0

0.85

20.0

0.83

尾粉砂

(水上)

18.5

9.8

30

0.45

75.0

0.80

118.0

0.77

尾粉土

(水上)

19.5

9.8

28

0.45

75.0

0.80

118.0

0.77

尾粉质粘土(水上)

19

10.78

16

0.45

75.0

0.80

118.0

0.77

 

2.2.5分析结果与小结

2-82-13给出了小地沟尾矿坝后期有限元的计算结果。图2-8给出了x方向(水平方向)应力等值线。图2-9给出了y方向(铅直方向)应力等值线。图2-10给出了剪应力等值线。图2-11与图2-12分别给出了最大与最小主应力分布,最大与最小主应力均大于零表明无拉应力出现。剪应力水平见图2-13,剪应力水平最高不超过0.9

以上的有限元应力应变分析表明,尾矿坝内应力水平不高,坝体内未发现接近极限平衡状态仅存在拉应力的区域,因此坝体是稳定的。静力有限元分析与滑弧稳定性分析的结论相一致。

 


 

2-8  尾矿坝后期水平方向应力分布(单位KPa

 

2-9  尾矿坝后期铅直方向应力分布(单位KPa

 

1-10  尾矿坝后期剪应力分布(单位KPa

 

2-11  尾矿坝后期最大主应力分布(单位KPa

 

2-12  尾矿坝后期最小主应力分布(单位KPa

 

2-13  尾矿坝后期剪应力水平


3 坝体抗震稳定与动力反应分析

       抗震分析按七度烈度设防。抗震分析包括两部分:拟静力法抗震稳定分析和有限元地震动力反应分析。抗震稳定分析主要给出在设计地震裂度时滑弧稳定的安全系数。地震动力反应分析主要计算在输入地震波作用下,坝体的应力,变形及加速度放大因子。

 

3.1 抗震稳定性分析

3.1.1拟静力法抗震稳定分析原理

拟静力法抗震稳定分析就是在坝体静力滑弧稳定分析的基础上,再考虑到地震惯性力对稳定的作用。根据“水工建筑物抗震设计规范SDJ10-78[7],抗震稳定安全因数K可写为

K=              3-1

上式中分子表示抗滑力(力矩),分母表示滑动力(力矩)。重力与渗透力的影响是通过“替代容重法”考虑进去的。各项的具体含义是:

    就条土滑裂面长度,

    条土第条滑裂面上的粘聚力与摩擦角,

    i条土条实重

    i条土条实重

    :浸润线以上的条块实重

    :浸润线以下,静水位以上的饱和重

    :静水位以下的浮重

    :浸润线以下,静水位以上的浮重

    i条土条滑裂面的倾角

    :作用在条块重心处的水平地震惯性力,即条块实重乘条块重心处的为水平地震系数,为地震水平最大加速度的统计平均值与地震重力加速度的比值,对于设计地震裂度为七度,=0.1是综合影响系数取1/4为地震加速度分布系数,其形式见SDJ10-783,其值按坝高向由1.02.5变化。

    :作用在条块重心处的竖向地震惯性力,即条块实重乘条块重心地的,其作用方向可向上(—)或向下(+),以不利于稳定的方向为准;

Mc:水平向地震惯性力Q对圆心的力矩。

 

3.1.2计算结果

动力稳定计算的材料参数按表2-1取,地震烈度为7°,表3-1列出小地沟各运行工况的动力稳定系数表3-1至图3-5给出各工况下静力计算的滑弧示意图,图中标明了滑弧的圆心位置,滑弧半径以及稳定系数。

 

3-1  小地沟二维各工况动力稳定安全系数一览表

投产期别

堆筑标高

等级

演算工况

 

图号

K

初期

1330m

正常运行水位

250干滩)

渗有效

3-1

1.51

排渗失效

——

——

洪水位(100干滩)

渗有效

3-2

1.16

排渗失效

——

——

后期

1400m

正常运行水位

250干滩)

渗有效

3-3

1.18

排渗失效

3-4

1.15

洪水位(100干滩)

渗有效

3-5

1.07

排渗失效

——

——

 


3-1  小地沟尾矿库堆筑至1330高程正常水位排渗有效动力稳定滑弧计算结果

 

3-2  小地沟尾矿库堆筑至1330高程洪水位排渗有效动力稳定滑弧计算结果

 

 

 

 

3-3  小地沟尾矿库堆筑至1400高程正常水位排渗有效动力稳定滑弧计算结果

 

 

 

 

 

 

 

 

 

 

 

3-4  小地沟尾矿库堆筑至1400高程正常水位排渗失效动力稳定滑弧计算结果

 

3-5  小地沟尾矿库堆筑至1400高程洪水位排渗有效动力稳定滑弧计算结果


3.1.3计算结果分析

根据《选矿厂尾矿设施设计规范》和《尾矿库安全技术规程》,表3-2给出了规范要求的抗滑稳定安全系数。

3-2  坝坡抗滑稳定最小安全系数

 

初期

后期

坝顶标高m

1330

1400

坝高(m

116

186

等级

正常运行

k

1.25

1.25

洪水运行

k

1.15

1.15

特殊运行

k

1.05

1.05

 

7度烈度对宝山小地沟尾矿库的稳定性进行验证,在堆筑高程1330m 特殊运行工况(洪水位运行+7度地震)下,动力稳定最小安全系数K=1.16>1.05,满足规范要求,进一步验证坝体物理参数合理性,进而将其应用于分析堆筑高程1400m 正常水位+7度地震和洪水位工况+7度地震的工况中,最小稳定安全系数分别是K=1.18K=1.07均大于1.05。划弧位置均位于后期堆积体上,需要强调的是,特殊工况运行时,同样需要严格保证各级排渗设施有效,严格控制干滩长度,密切监测初级坝体透水能力,有效控制浸润线高度,保证坝体稳定。

 


3.2 有限元地震动力反应分析

动力分析是研究地震条件下坝体的反应,主要是坝体各部位的反应加速和动剪应力,同时亦可能给出在地震动时程内坝体各部位应力与变形的反应时程曲线。

3.2.1动力分析的基本方法

经有限元离散后,运动状态中各结点的动力平衡方程如下

                  3-2

上式左端各项仪次代表惯性力,阻尼力和弹性力,右端是动力荷载向量。代表位移,若仅存在地震动力载荷,经变换后得

                                  3-3

其中a代入输入的地震加速度时程。

                          3-4

                            3-5

[B][D]为应变位移矩阵与应力应变矩阵,通过[B][D]

                 3-6

                                       3-7

分别代表应力与应变。

    阻尼矩阵[C]与材料的内摩擦特性有关,现在较常用的Rayleigh理论假设阻尼由两部分组成,一部分与应变率成正比,另一部分与速度成正比,从而阻尼矩阵可以写为

                                 3-8

为比例系数。

    对于像土这样的不均匀材料,且应力应变关系是非线性的,只能按下式算出单元阻尼矩阵后再叠加

                                  3-9

 

3.2.2动力分析的本构关系

本方法采用大量线性粘弹性模型描述土体振动时的本构关系。材料模量包括剪切模量G泊桑比,和阻尼比。在中(G),G与应变的关系是非线性的,通常表示为曲线。其中最大动剪比模量可表示为

                                 3-10

Kn为试验参数,见表3-3为大气压力,为土体平均有效应力(可由静力有限元分析确定)。利用Rayleigh阻尼系数可由下式求出

                          3-11

式中为单元阻尼比其值见表3-5为坝体的基本频率,可通过解特征值问题求出。

 

3.2.3计算步骤

由于本构关系中(G~是非线性的,有限元方程的求解工作量很大,工程应用是常采用等价线性方法来模拟,即假定(G)可由平均动剪应变确定,而经经验认为是最大剪应变的某一分数(通常根据Seed提议采用0.65),为此,需要采用迭代法计算,即先假定的初始分布,通过时间积分求出新的分布,如此迭代,直到误差小于容许值为止。动力方程的时间分采用Wilson-法求解,基频通过斯托多拉迭代法解特征值问题求出[9]

 

3.2.4材料的动力性质

考虑到地震时动荷动历时很短,孔隙水来不及排出,浸润线以下土体采用饱和容重计算。在本次计算中,筑坝材料的容重及最大动剪应力中的Kn值见表3-3。筑坝材料的及阻尼比随动应变幅的关系曲线值见表3-4其值参照已建工程经类比确定。

3.2.5有限元网格,边界条件与设计地震

动力有限元分析的网格与静力分析完全相同,由四边形与三角形网格组成。

地震动输入于基岩中,坝底固结于基岩,坝面为自由边界条件。

输入地震动为 ELCETRO (E-W) 记录地震动时程线,该地区性震时程持时T0=32秒,加速度峰值=210

库区位于忻州市繁峙县岩头乡官地村,设计按地震烈度为7度计算,设计基本地震加速度为0.1g

 

3-3  容重及参数值

材料名称

水下容重(KN/m3)

水上容重(KN/m3)

Gmax系数K

泊桑比

Gmax指数n

滤水堆石

20

20

1105

0.45

0.5

尾细砂

18.5

18

742

0.43

0.5

尾粉砂

18.5

18.5

515

0.43

0.5

尾粉土

20

19.5

515

0.43

0.5

尾粉质粘土

19

19

515

0.43

0.5

 

3-4  关系

ed

3.00E-06

1.00E-05

3.00E-05

1.00E-04

3.00E-04

1.00E-03

3.00E-03

1.00E-02

滤水堆石G/Gmax

0.99

0.99

0.94

0.87

0.75

0.56

0.36

0.28

尾细砂G/Gmax

0.99

0.93

0.83

0.66

0.52

0.35

0.20

0.15

尾粉砂G/Gmax

0.99

0.93

0.83

0.66

0.52

0.35

0.20

0.15

尾粉土G/Gmax

0.99

0.93

0.83

0.66

0.52

0.35

0.20

0.15

尾粉质粘土G/Gmax

0.99

0.93

0.83

0.66

0.52

0.35

0.20

0.15

滤水堆石 l %

2

4.23

5.47

6

10.62

15.68

20.29

21.5

尾细砂 l %

1.3

2.4

6.0

5.5

9.5

15.0

21.0

26.0

尾粉砂 l %

1.3

2.4

6.0

5.5

9.5

15.0

21.0

26.0

尾粉土 l %

1.3

2.4

6.0

5.5

9.5

15.0

21.0

26.0

尾粉质粘土 l %

1.3

2.4

6.0

5.5

9.5

15.0

21.0

26.0

 

3.2.6计算结果

计算工况为尾矿坝后期(坝顶标高1400m)洪水运行排渗有效情况下的地震动力反应。动力计算时输入的平均有效应力由静力分析的结果给出,其分布见图3-6。尾矿坝后期计算得出的动应力分布分别见图3-7~图3-93-12给出尾矿坝水平方向的加速放大因子,从中可以看出,水平向加速度基本没有放大,加速度放大因子较高的区域出现在坝顶体附近,最大值不超过为1.30

尾矿坝后期计算得到的动剪应力及动应力比见图3-103-11。最大动剪应力比出现在坝体顶部,但是数值不大,因此,对于七度烈度地震,尾矿坝基本不会产生液化。

3-14至图3-17给出了尾矿坝初期坝顶部附近的动应力与加速度反应时程。从结点反应加速度时程线可以看出,动应力反应同加速度反应是相对应的,且反应频率较输入的地震频率低。

3-6  尾矿坝后期平均有效应力分布(单位KPa

3-7  尾矿坝后期 X方向动应力分布(单位KPa

3-8  尾矿坝后期Y方向动应力分布(单位KPa

3-9  尾矿坝后期 X-Y方向动应力分布(单位KPa

 

3-10  尾矿坝后期最大动剪应力等值线(单位KPa

 

3-11  尾矿坝后期动应力比等值线

 

3-12  尾矿坝后期水平方向加速放大因子等值线

   

 

3-13  输入的ELCETRO (E-W) 地震动时程线(计算时归化到0.1g

 

3-14  尾矿坝后期坝顶附近动应力反应时程