某加高扩容尾矿库三维静力敏感性及地震反应数值分析

陶东良，王 硕，吴蒙蒙，宁 伟，徐力群，孙一清，冯亚新

(1.长沙有色冶金设计研究院有限公司，长沙 410011；2.河海大学 水利水电学院，南京 210098；3.洛阳栾川钼业集团股份有限公司，河南 洛阳 471500)

我国现存尾矿库总量约12 655 座，其中病库、险库和危库占比高达39%，正常库仅有7 745 座[1]。尾矿库一旦失事，将对其下游人民群众的生命财产安全造成不可估量的损失[2]，此外由于尾矿库筑坝材料强度存在一定的不均匀性和随机性，并且地震作用是影响尾矿坝稳定性的主要不利因素，因此对尾矿坝进行静力敏感性分析及地震作用下的动力特性分析，预测其加高扩容后的坝体稳定性及抗震性能对于确保其未来的长久稳定运行具有重要的现实意义。

国内外学者在该领域研究较为深入，如王凤江[3]通过大量采用上游法修筑的尾矿坝震害实例，对上游式尾矿坝的抗震加固研究及其存在的问题进行了详细梳理；徐志英等[4]对德兴铜矿尾矿坝在遭遇高强度地震时的安全稳定性及液化可能性进行了有限元计算，结果符合工程实际；柳厚祥等[5]开发了基于变分原理的二维稳定分析程序，得出采用对数螺旋滑动面进行尾矿坝安全稳定分析更加贴近实际；周健[6]考虑了流固耦合作用和尾矿砂石的非线性，采用三维有效应力法分析了某尾矿库在地震作用下变形时所产生的残余空隙水压力。

本文采用河海大学编写的三维静动力有限元分析程序对某尾矿库进行了静力稳定性分析、参数敏感性分析及动力特性分析，较为全面地验证了该尾矿库实施加高扩容方案的可行性，具有实际工程意义。

1 地质概况

某尾矿库地处秦岭山脉东北缘低中山区，由沟口向库内进入约700 m后分为南北两条支沟，南侧支沟向沟谷内进入约350 m后又分为南北两条支沟，将沟谷从北往南编号为(1号沟、2号沟、3号沟)，三道沟三面环山，库尾及两侧均为山坡。该尾矿库所在主沟长约1.9 km，其水文地质条件较为简单，地表水主要为库尾蓄水、尾矿水主要分布在库尾，按正常情况每天排入尾矿库尾矿水含量为1.4×104m3；库址区抗震设防烈度为Ⅵ 度，设计基本地震加速度值为0.05 g，设计地震分组为第一组；建筑场地按构造活动性、边坡稳定性和场地地基条件综合评价为对建筑物抗震不利地段，建筑场地类比为Ⅱ 等，特征周期Ts为0.45 s。初期坝为碾压堆石坝，坝顶标高1 310.00 m，坝高67 m，顶宽4 m，坝轴线161 m，上下游边坡均为1∶1.75。现状尾矿堆积坝坝顶标高约1 358.00 m，总坝高115 m，后期尾矿采取上游法堆坝，平均堆积边坡为1∶5，设计最终堆积标高为1 442.00 m，总坝高199 m，总库容2 672.8×104m3，等别为Ⅱ 等。

为验证该尾矿库加高扩容方案的可行性和合理性，本文采用非线性有限元法建立了该加高扩容完成后尾矿坝的三维有限元模型，对影响坝体变形和应力的关键性施工填筑过程进行模拟，计算并分析了加高库容完成后该尾矿库的静力和动力稳定性，结果表明加高扩容方案的设计是合理的。

2 计算模型及设计参数选取

2.1 三维有限元模型

尾矿库的实际排矿沉积过程与堆积坝的填筑过程是交叉进行的，在实际建模中很难精确模拟实际的填筑过程，因此本文为尽可能提高静动力计算结果的准确度并适当降低建模的复杂性，在充分考虑到该工程库址区的水文地质情况和现状坝高堆积坝填筑情况的基础上，采用非线性三维有限元法模拟了在后续加高扩容方案中可能影响尾矿堆积坝坝体的应力和变形的关键性施工过程。该有限元模型考虑初期坝已填筑完成且尾矿料已堆积至现状坝顶标高1 358.00 m，将尾矿堆积坝的堆积过程大致按排矿沉积的顺序进行模拟，通过分级加载的形式逐层提高坝体标高，直至最终坝顶标高1 442.00 m。

三维有限元模型坐标系采用笛卡尔坐标系：选取西安80坐标系(X=37 551 603.932 9，Y=3 754 867.381 8)作为计算坐标系原点O；高程系统为1985黄海高程系；模型上、下游边界分别截取至现状尾矿料沉积滩向外约800 m、初期坝坝趾下游约25 m；左、右岸边界分别截至左、右岸坡向外约200 m和150 m处；垂直向地基截取至高程1 160 m，位于坝体建基面以下约150 m处，向上至加高扩容完成后最大坝顶标高1 442.00m。

依据伽辽金有限元理论[7]，采用“控制断面超单元有限元自动剖分技术”将全部计算区域离散成互不重叠或交叉的八结点六面体空间等参单元，有限元模型网格共包含79 238 个结点，78 245 个有限单元。图1为最终坝高尾矿坝三维有限元模型网格示意图，部分剖面有限元网格见图2。

图1 最终坝高尾矿坝三维有限元模型Fig.1 3D finite element model of tailings dam with final dam height

图2 部分剖面有限元网格图Fig.2 Section finite element mesh

2.2 计算参数

碾压透水堆石初期坝以及上游尾矿堆积料均按非线性材料考虑并采用邓肯张E-v模型；基岩按线性材料考虑，采用线弹性模型。根据提供的工程地质及水文地质资料，坝基岩体线性本构关系计算参数见表1，其余材料非线性本构关系的标准设计参数见表2。

表1 线性本构关系材料计算参数

表2 非线性本构关系材料计算参数

3 静力计算结果分析

计算工况为该尾矿库已堆积至最终坝高，坝体所受荷载包括自重及上游沉积滩积水压力。计算时，基岩作为一次性加载，并在对上游堆积坝坝体实施分级加载前将基岩区域全部结点的位移均初始化为0，仅保留其单元应力，从而获得地基初始应力场，因此下文所述位移均指填筑开始后坝体的位移，水平位移以顺河向为正，垂直位移以向上为正，应力以压为正。

3.1 位移

图3为该工况下坝体典型断面的顺河向位移分布图。由图3可见，由于尾矿库地形较为复杂，因此坝体及其在自重作用下的位移呈非对称分布，此外，由于静力计算时考虑了正常运行水位工况下的尾矿库上游沉积滩积水对坝顶的静水压力作用，并且坝体下基岩整体有向下游倾斜的趋势，因此坝体向下游移动的趋势较为明显。由于库区基岩向下游的平均倾斜幅度由坝中向左右两岸逐渐降低，因此最大顺河向位移发生在坝中部结点处，为-1 732.78 mm。

图3 静力作用下堆积坝坝体典型断面顺河向位移分布(单位：mm)Fig.3 Displacement distribution along the river of a typical section of accumulation dam under static action(Unit：mm)

图4为该工况下坝体典型断面的垂直位移分布图。由图4可知，坝体最大沉降发生在坝体中部，约1/2坝高处，为-2 245.39 mm，约占最大坝高的1.13 %，导致坝体发生沉降的主要原因是尾矿库在实际的排矿堆积过程中，尾矿砂持续发生排水固结，此外由于采用湿堆法堆存的尾矿砂具有较高的初始含水率，因此该工况下的坝体沉降大于一般的土石坝。

图4 静力作用下堆积坝坝体典型断面垂直向位移分布(单位：mm)Fig.4 Vertical displacement distribution of typical section of accumulation dam under static action(Unit：mm)

3.2 应力

图5、6分别为该工况下坝体典型断面的第一和第三主应力分布图。由图5、6可见，该工况下的坝体第一及第三主应力均为压应力，等值线分布较为均匀且近似于坝坡平行，数值规律表现为从坝坡向坝脚逐渐增加，最大第一和第三主应力分别为2 468.85 kPa和1 635.44 kPa，均发生在坝中部断面的最底部附近。应力水平最大值为0.827，表明坝体在目前荷载情况下是稳定的。

图5 静力作用下堆积坝坝体典型断面第一主应力分布(单位：kPa)Fig.5 First principal stress distribution of typical section of accumulation dam under static action(Unit：kPa)

图6 静力作用下堆积坝坝体典型断面第三主应力分布(单位：kPa)Fig.6 Third principal stress distribution of typical section of accumulation dam under static action(Unit：kPa)

4 静力参数敏感性分析

由于尾矿库在实际的排矿堆积过程中存在诸多不确定性因素，可能导致尾矿堆积料的抗变形性能下降，因此有必要开展尾矿砂主要变形参数敏感性分析，本文考虑将坝体上游堆积料的变形模量K分别降低和提高10 %及20 %后参与计算，敏感性分析参数见表3，位移及应力计算结果见图7、8。

表3 各敏感性分析工况下的尾矿料变形模量K

图7 不同敏感性分析工况下的坝体最大水平及垂直位移计算结果Fig.7 Calculation results of maximum horizontal and vertical displacement of dam body under different sensitivity analysis conditions

图8 不同敏感性分析工况下的坝体应力及应力水平计算结果Fig.8 Calculation results of dam stress and stress level under different sensitivity analysis conditions

由图7可见，坝体变形对尾矿砂主要变形模量K变化的敏感性较高，K降低，坝体的整体力学性能变差，弹性模量系数偏低，压缩变形增大，从而导致坝体位移增加，反之则减小，K降低10%和20%工况下，坝体最大顺河向位移相较设计工况分别增大6.86%和7.72%，最大沉降分别增大7.24 %和8.35 %，K提高10%和20%工况下，坝体最大顺河向位移相较设计工况分别减小5.49%和7.34%，最大沉降分别减小6.33%和9.32%。由图8可见，坝体的应力特性对K变化的敏感性较低，相较设计工况，K降低10%和20%工况下，坝体最大第一主应力增大率仅分别为0.33%和0.41%，最大第三主应力增大率仅分别为0.12%和1.36%，应力水平增长率分别为1.09%和2.42%；K提高10%和20%工况下，坝体最大第一主应力减小率仅分别为0.97%和1.32%，最大第三主应力减小率仅分别为0.35%和1.67%，应力水平减小率分别为1.81%和3.39%。

5 动力稳定分析

5.1 动力稳定性计算力学参数

动力稳定性计算采用等效非线性粘弹性模型，土体动应力-应变的非线性和滞后性基本特征借助等效剪切模量和等效阻尼比来实现。依据Hardin和Drnevich公式，并结合工程实际确定坝料各分区及坝基岩体的动力计算参数，见表4、5。

表4 坝料各分区动力计算参数

表5 坝基岩体动力计算参数

5.2 输入地震动加速度曲线

坝址区地震烈度及尾矿坝抗震设防烈度均为VI度，50年超越概率5%的基岩地震动加速度曲线时长为40 s，水平向输入基岩地震动加速度曲线见图9，垂直向取其2/3。基岩地震动水平峰值加速度依据规范调整为52.01 cm/s2。

图9 水平向输入基岩地震动加速度曲线Fig.9 Horizontal input bedrock ground motion acceleration curve

5.3 地震反应分析

1)加速度响应

图10为在地震作用下坝体中部断面的顺河向及垂直向绝对加速度分布情况。从高程上来看，顺河向及垂直向绝对加速度反映大小由坝基至坝顶逐渐增大；从距离来看，绝对加速度数值存在由坝体内部向坝坡方向逐渐增大的趋势。此外，由于初期坝与堆积坝材料的最大动剪切模量存在一定差异，导致绝对加速度反应在初期坝坝顶附近出现局部集中现象，因此在抗震加固设计与施工阶段应重点考虑在该区域采取适当的抗震加固措施。计算得出坝体在顺河向及垂直向的绝对加速度最大值分别为1.67 m/s2和1.03 m/s2。

图10 坝体中部断面顺河向及垂直向绝对加速度分布(单位：m/s2)Fig.10 Absolute acceleration distribution along the river and vertically in the middle section of the dam(Unit：m/s2)

2)位移响应

图11为在地震作用下坝体中部断面的顺河向及垂直向位移响应分布情况。结果表明，坝体顺河向及垂直向的最大动位移值均发生在相对坝高最大的坝顶附近，顺河向及垂直向最大动位移值分别为72.38 mm和29.67 mm。此外，由于坝基处岩体分布较为均匀，坝体各部分变形均匀，因此地震时将不会出现平行或垂直于坝轴线方向的裂缝。

图11 坝体中部断面顺河向及垂直向位移响应分布(单位：mm)Fig.11 Displacement response distribution along the river and vertically in the middle section of the dam(Unit：mm)

3)应力反应及安全系数m

计算表明，坝体最大第一、第二和第三主应力反应大小(叠加静力)分别为2 533.72、1 845.72 和1 795.64 kPa，将该有限元模型上，某一单元的局部抗滑稳定安全系数Fs定义为该单元所在潜在滑移面上的抗剪强度与剪应力的比值，计算得出坝坡单元抗滑稳定最小安全系数为1.244，出现在9.38 s，平均安全系数为1.415，该典型单元在地震期间的安全系数变化情况如图12所示。

图12 坝坡某一典型单元的安全系数时程变化曲线Fig.12 Time history variation curve of safety factor of a typical unit on dam slope

4)地震永久变形

图13为在地震作用下的坝体中部断面顺河向及垂直向永久变形分布情况。依据规范，要求地震水平加速度不小于0.05g地震区的尾矿库，堆积坝滩顶与正常生产水位的高差应大于或等于对应的最小安全超高和地震沉降值以及地震壅浪高度之和。计算表明，地震后坝体的最大永久顺河向及垂直向位移分别为-500.95 mm和-847.62 mm，地震永久沉降约占最大坝高199 m的0.43%。在地震作用下，堆积坝坝体(含沉积池)的最大垂直向永久变形(沉降)约0.85 m，地震壅浪高约1 m。尾矿坝沉积滩滩顶高程为1 442.00 m，正常运用水位为1 439.00 m，考虑最大沉降后，尾矿沉积滩滩顶标高超过正常运用水位2.15 m>2.0 m(地震壅浪高1.0 m+最小安全超高1.0 m)，同时干滩长度大于规范规定的100 m。因此，不会发生因地震作用而导致的沉积池水漫溢尾矿坝的情况。

图13 坝体中部断面顺河向及垂直向永久变形分布(单位：mm)Fig.13 Permanent deformation distribution along the river and vertically in the middle section of the dam(Unit：mm)

6 结论

本文首先分析了该尾矿坝在各料区均采用标准设计参数下的静力稳定性，并对尾矿砂的主要变形参数进行敏感性分析，最后在静力计算结果的基础上计算并分析了坝体在地震作用下的动力稳定性，主要结论如下：

1)各工况下的坝体应力变形及位移分布符合一般规律，数值上与类似工程较为相近，表明该计算程序的可靠性。设计工况静力计算结果表明，堆积坝顺河向及垂直向位移呈非对称分布，坝体最大顺河向及垂直向位移分别为-1 732.78 mm和-2 245.39 mm，坝体第一及第三主应力均为压应力，应力水平为0.827，表明坝体不会在自重及上游水压力荷载下发生破坏。

2)静力参数敏感性分析表明，堆积坝坝体变形对尾矿砂的主要变形模量K的变化敏感性较强，应力及应力水平对K变化的敏感性均较弱，建议在堆积子坝时尽可能使用粒径较大的尾矿料，并提高坝体的密实度，以减小坝体变形并提高其稳定性。

3)动力稳定性分析结果表明，坝体在地震作用下的绝对加速度反应较为强烈，初期坝坝顶附近呈现局部集中现象，应考虑在类似区域采取适当抗震加固措施；位移响应分布较为规律，不会因地震导致坝体出现水平或垂直向裂缝；由于堆积坝最大坝高可达199 m，因此主应力反应也较大，基于应力计算结果计算出的坝坡典型单元的最小安全系数为1.244，平均安全系数为1.415，满足规范要求；考虑发生最大沉降后，不会发生因地震作用而导致的沉积池水漫顶等情况。