尾矿坝地震动力响应可靠度分析

王梓鉴，杨 俊

（九江市水利电力规划设计院，江西 九江 332000）

0 引言

尾矿库是由金属或非金属等矿山选矿的废弃物所构成的，由于我国对矿产资源的需求与日俱增，尾矿库的数量也在不断增多。近年来，地震、洪水等原因导致我国尾矿库溃坝事故频繁发生，尾矿料的倾泻可能导致下游居民房屋受到冲击、农田遭受破坏，同时尾矿料中的一些重金属离子将会造成水污染，严重影响当地居民的生产生活。目前国内外学者针对尾矿坝的稳定性问题开展了大量有益的研究工作，如龙明魁[1]等对元江尾矿坝进行静力稳定性分析，研究尾矿坝的位移、应力变化规律；陈燕[2]等开展了不同降雨强度下的尾矿坝稳定性分析，尾矿坝稳定性会随着降雨强度的增加而下降；李良振[3]进行渗流-应力-边坡耦合，得出尾矿边坡稳定性系数，为尾矿坝具体施工提供了技术支持。然而，尾矿材料是性质极其复杂的地质介质，在施工过程中人为或非人为因素作用下，其物质组成与内部结构会发生一定的变化。实际上尾矿坝一般处于饱和的疏松状态，颗粒细小，并且比重较大，由于颗粒的组成、矿物的成分等因素的影响，它们看似稳定，实际上对扰动特别敏感，在地震中很容易发生液化和破坏变形[4]，因此，研究我国尾矿坝的抗震稳定性就显得尤为紧迫和必要。

由于尾矿坝的动力失稳机制复杂，目前多数的研究集中于将动力问题简化为拟静力问题来计算尾矿坝的安全系数，如李再光[5]使用拟静力法分析了南京钢铁尾矿坝在地震作用下的稳定性。当考虑简化地震作用效应时，用拟静力法简单有效，但地震作用下的动力特性却得不到真实地反映。动力有限元分析法将每一时刻的动应力施加到相应的静应力上，并对尾矿坝进行静力分布状态分析判断静力稳定性，作为坝体动力反应的基础，通过输入地震时称曲线，比较适用边界条件和岩土体结构的复杂动力问题，从而对尾矿坝进行地震过程中的变形分析和安全评价。另外使用结构可靠度分析方法可以使得稳定分析结果更具说服力。

本文结合动力有限元模型和结构可靠度理论，使用GEOSTUDIO软件，并将Seep/w、Quake/w、Slope/w三个模块相耦合，对江西省德安尾矿坝工程进行渗流分析、地震响应分析以及尾矿坝的稳定计算，最后采用蒙特卡罗法计算德安尾矿坝的失效概率，分析相关研究成果，对于德安尾矿坝的安全生产、生态环境的维护以及周边居民的健康生命安全具有十分重要的意义。

1 尾矿坝动力稳定性分析原理

在尾矿库地震失稳状态下，尾矿土的应力应变关系表现为非线性关系，本文采用等效线性本构模型[6]。等效线性模型是黏弹性模型的一种类型，特点在于将不同应变幅值下的滞回特性用阻尼比随应变变化，即将土体视为黏弹性体，反映土体滞后性和动剪应力－应变关系非线性的两个基本特征[7]。等效线性模型计算见式（1）。

式中，AL为应力应变滞回环面积；λ为等效阻尼比；AT为滞回环拐点与x轴、坐标原点O所围成的直角三角形面积。

如果震动频率很低，滞回环不包括黏性阻尼，会导致阻尼比很小，则需要修正公式进行修正，见式（2）、式（3）。

式中，λv为黏性阻尼比；δi、δi+1为相邻两个周期的位移振幅。

动力分析的步骤与静力方法大致相似，都要先把尾矿坝坝体离散为单元，然后对单元进行动力分析，动力有限元方程式为：

式中，[M ]为质量矩阵，可用集中质量法求得；[C]为单元阻尼矩阵；[K ]为单元刚度矩阵，通过常规有限元求得；a为加速度；u为速度；s为单元节点的位移；｛F（ t） ｝为单元节点上的动力荷载，可以表示为式（5）：

式中，｛Fb｝为单元结构上的体力；｛Fs｝为边界上的面力；｛Fn｝是节点上集中力；｛Fg｝为地震荷载引起的地应力。采用Rayleigh阻尼假设单元阻尼矩阵，即阻尼矩阵是刚度矩阵和质量矩阵的线性组合[8]：

式中，α、β 为比例系数。另外 α=λγω0，β=λγ/ω0，λγ为由单元结构剪应变所决定的阻尼比，ω0为坝体与尾矿库无阻尼系统的频率。

目前，通常采用Wilson线性加速度法求解式（4），然后逐步积分，求出每一时刻节点上的位移｛st｝、速度｛ut｝和加速度｛at｝。

2 可靠度分析方法

可靠度是指规定的时间范围内，工程在规定的条件下所能够完成预期功能的概率，对评价工程整体具有一定的指导意义，而无法完成预期功能的概率就定义为失效概率[9]。可靠度的分析方法有很多，如JC法、一次矩阵法、蒙特卡罗法等，本文采用蒙特卡罗法（MCS）。MCS方法采用大量重复抽样的基本思想，通过模拟实验次数来获取失效概率值，当设定的次数越高，所得到的失效概率值也就越精确。失效概率Pf=m/n，式中m为失效次数，n为实验次数。使用MCS对尾矿坝土体输入参数进行随机抽样后，可将输入参数特性值分别赋给尾矿坝稳定性有限元分析模型，再进行尾矿坝稳定性有限元分析。然而大量的稳定性分析计算繁琐、效率低下，本文通过将MATLAB软件编写程序与GEOSTUDIO软件实行对接，并在DOS环境下进行相关批处理计算，无需打开GEOSTUDIO软件界面，有效地避免了重复操作有限元软件，大大提高了计算效率。随后通过批处理计算获得相应的n个临界安全系数FSmin和滑动面等信息来计算尾矿坝的失效概率。

3 工程案例

江西省九江市德安县拥有铅锌锑矿矿床一座，为了解决选矿废渣的排放堆积问题建立尾矿库。德安尾矿库采用上游法筑坝，坝顶标高64.30m，初期坝坡比为1：2，有 3 座堆积子坝，1#子坝上、下游坡比为 1：2，2#、3#子坝的上下游坡比为1：4，堆积坝目前标高75.00m。

表1 尾矿坝各土层物理力学参数

图1 德安尾矿坝模型

3.1 有限元模型及计算参数

本文通过 GEOSTUDIO中 Seep/w、Quake/w与Slope/w模块耦合建立工程地质模型来进行模拟分析。为了能更直观展现，德安尾矿坝有限元模型一共划分为2 314个节点和2 211个网格大小为1.00m的四边形和三角形混合单元，并设置了A、B、C、D四个地震响应历程点监测地震过程中的位移和加速度，工程地质模型和特征点位置见图1。边界调节为：上游水位（左侧）22.50m，下游水位（右侧）10.00m，各土层所采用的物理力学参数见表1。首先在Seep/w模块建立德安尾矿坝坝体模型，进行饱和渗流分析，随后将渗流分析结果导入Quake/w模块，输入地震波，进行动力响应分析。

3.2 地震作用下尾矿坝的动力响应分析

尾矿坝动力稳定计算分析需要在计算程序中输入地震动时程曲线（地震加速度时程曲线），且不同的地震动时程曲线对坝体动力计算结果影响较大。根据《中国地震动峰值加速度区划图》（GB 18306-2015），该尾矿所在区地震设防烈度7度，设计基本地震加速度值为0.10g，需要对加载地震波进行人工合成。根据德安尾矿坝的场地信息，以场地反应谱为目标谱生成人工加速度时程曲线，卓越周期为0.2s，时长为10s，时间步长为0.02s，如图2所示。

图2 地震加速度时程曲线

本次地震历程响应点选取A、B、C、D四点，分别代表尾矿坝初期坝的坝底、坝顶和堆积坝的坝中、坝顶四个部分。根据该尾矿坝所在的位置加入地震参数，进行位移响应分析，位移时间图见图3。

从图中可以得到4个响应点中最先产生X方向位移的是堆积坝坝顶的D点，然后是堆积坝坝中的C点，最后是初期坝的A点和B点。位于尾矿坝顶点的位移峰值大于其它各点的峰值，四个地震响应点的位移变化趋势基本相一致。可以预判，在同一地震波的影响下，尾矿坝坝顶受地震的影响最大。随着地震事件的增加，A、B、C、D四点处的峰值位移都在不断的增加。其中1.7s和6.7s时分别产生了向左和向右的最大位移，为了更加直观的看到尾矿坝的情况，分析了这两个时刻的尾矿坝相对位移图（见图4、图5）。

图3 X-位移时间图

分析A、B、C、D四点X方向速度见图6，可以看出A点的峰值速度出现在1.98s，其值为0.076m/s；B点峰值速度出现在4.50s，其值为0.100m/s；C点峰值速度出现在4.05s，其值为0.067m/s；D点峰值速度出现在4.33s，其值为0.120m/s，四个点的峰值速度有所差别，但是B、C、D点峰值速度出现的时间点大致相同，堆积坝坝顶的D点在X方向速度波动最为剧烈，各响应点最大位移从坝顶方向到坝底逐渐减小，四点总体的震动趋势与地震波行较为相似，而尾矿坝不同位置在同一地震波影响下反映出的剧烈程度不一样。图7为地震过程中尾矿坝的总应力随时间变化的分布图，震后总应力大小基本没有发生变化，堆积坝坝顶D的总应力最大，初期坝坝底A的总应力最低。

图4 1.2s左侧最大位移

图5 8.8s右侧最大位移

图6 地震历程中速度时间变化图

图7 地震历程中总应力时间变化图

图8 最危险滑面平均加速度时程图

随后对尾矿坝进行稳定性分析，将Seep/w、Quake/w和Slope/w模块进行耦合，在地震动力分析的基础上，使用进入进出法分析、计算尾矿坝的安全系数，不同时刻尾矿坝最危险滑面平均加速度和坝体安全系数时程曲线如图8和图9所示。在初始状态下该尾矿坝的安全系数为1.36，一开始安全系数趋于稳定，之后便有较大的起伏，在4.7s时，安全系数最低为0.92，根据首次破坏准则，此时尾矿坝可能发生溃坝危险。且安全系数时程图跟地震波型图亦十分相似，安全系数的峰值与地震波的峰值存在一定的滞后性。

3.3 可靠度的计算

对尾矿坝进行可靠度分析，结合德安尾矿坝工程资料，并参考郑欣[10]、张扬等[11]，确定土体黏聚力及摩擦角的均值及标准差，结果如表2所示。

图9 最危险滑面安全系数时程图

表2 黏聚力及摩擦角参数统计表

尾矿坝破坏可靠度分析方法使用了MCS法，在Seep/w模块中以参数均值建立尾矿坝有限元分析模型、划分有限元网格、设置边界条件，在Quake/w模块中建立动力模型，在Slope/w模块建立稳定性分析模型，在Seep/w、Quake/w和Slope/w模块耦合计算的基础上，通过设置滑移面进入进出口的位置来进行边坡稳定分析，求出安全系数。将尾矿坝稳定性有限元分析模型另存为名为“安全系数.xml”的计算源文件。当尾矿坝土体输入参数表示为标准正态随机变量之后，可将输入参数特性值分别赋给尾矿坝稳定性有限元分析模型，生成新的尾矿坝稳定性分析“安全系数.xml”计算文件。最后使用MATLAB软件在DOS环境下直接调用GEOSTUDIO软件内部程序对新生成的“安全系数.xml”文件进行一体化边坡稳定性有限元批处理分析，采用摩根斯坦-普莱斯法计算最危险滑动面，无需打开GEOSTUDIO软件界面，大大地提高了计算效率。最后通过可靠度理论来计算失稳风险概率，实验次数为10 000次，倘若人为的对有限元软件进行重复性的操作来计算10 000次的话，可能需要几个星期，这里使用批处理方法只需普通电脑自动计算几天，无需人为干涉。经计算，德安尾矿坝在7度地震作用下的失稳风险概率为0.37%。

4 结论

（1）本文结合动力有限元模型和结构可靠度理论，分析了德安尾矿坝在7度地震作用下的动力稳定性，探究了德安尾矿坝地震过程中的位移、速度、应力、安全系数等变化情况，比传统的拟静力方法更能反映尾矿坝在地震作用下的动力响应。在地震作用中，不同时刻坝体稳定性随地震加速度的变化而变化，随着临界滑面平均加速度的增加，尾矿坝安全系数会逐渐降低。

（2）运用蒙特卡罗随机抽样方法，并编写MATLAB程序与GEOSTUDIO软件相结合，将商业有限元软件GEOSTUDIO当做黑匣子直接调用，进行一体化尾矿坝边坡稳定性有限元批处理计算，可以快速的进行了10 000次尾矿坝的稳定性分析，大大提高了尾矿坝的计算效率。

（3）经可靠度分析，德安尾矿坝的失稳概率为0.37%。尽管德安尾矿坝在7度地震工况下的失稳概率较低，与尾矿坝设计规范给出的允许安全系数和破坏概率相比，该尾矿坝发生的失稳破坏可能性较小，但是仍需要设计一些尾矿坝安全防护抗震措施。

（4）本文通过理论分析与数值模拟相结合的方法对地震作用下尾矿坝的稳定性可靠度进行了计算分析，研究结果也具有一定的实际意义，但是分析过程中仍做了一些相应的简化，因此，作者认为针对该方面的研究还有待进一步深入：1）影响尾矿坝稳定性的不确定因素较多，本文仅考虑了尾矿材料参数粘聚力和内摩擦角的不确定性，未考虑到地震以及降雨作用的不确定性，因此，在今后的研究中可进一步考虑更多参数的不确定性，以使计算结果更加贴合实际；2）本文采用的GEOSTUDIO软件为2D软件，不能反映尾矿坝最真实的受力变形情况。因此，在今后的研究中可借助其他3D有限元软件进行计算分析，结果将会更加接近实际。