考虑弱化效应的加高扩容尾矿坝地震可靠度分析及风险评价

李春立,李亮*,徐亮,赵民,王超,蔡凤珍

(1.青岛理工大学土木工程学院,青岛 266520; 2.青岛市勘察测绘研究院,青岛 266033)

诱发尾矿坝溃坝的因素和机理复杂繁多,究其原因,主要与尾矿坝筑坝材料的不确定性以及地震、强降雨等复杂的因素相关[1-4]。由于中国地震频发,给尾矿坝构成了极大威胁。地震对尾矿坝稳定性的影响主要表现为地震的惯性力使处于极限平衡状态的尾矿坝产生变形以及地震循环剪切作用使坝体震动液化[5-10]。鉴于此,国家有关规范中明确提出[11-12],为保证尾矿坝的安全,要求对设防烈度为Ⅶ度及以上、尾矿库等级在三级及以上的尾矿坝进行抗震性能分析。在尾矿坝地震稳定性分析方面,国内外学者开展了一系列研究,以获得尾矿坝的安全系数、位移、液化变形及大坝的应力应变状态。例如,杨安银等[13]采用有限元Quake/W软件分别研究了尾矿坝在三种地震作用下的加速度、位移、液化区域、安全系数的动力响应特性。秦晓鹏等[14]采用拟静力法和有限差分法计算地震作用下的边坡时程安全系数,并采用PL-Finn模型结合Flac3D研究了尾矿坝的液化变形。Naeini等[15]分别采用有限元Quake/W软件和有限元Sigma/W软件对大坝进行动态和应力重分布分析,并确定坝体的永久位移。朱明明[16]考虑随机地震动和尾矿材料参数的二维空间变异性,对尾矿坝进行了动力可靠度分析和地震敏感性分析。研究成果有助于合理评估尾矿坝地震安全性能,有力保障矿山安全。

近年来,由于矿山继续扩大再生产,尾矿库亟需进行加高扩容,评估尾矿坝加高扩容后的地震安全性具有指导意义。尹光志等[17]采用有限元时程分析法,分别对云南大沙河尾矿库加高扩容坝体在现状条件下、加高中增设排渗措施和直接加高后三种工况进行了动力分析和抗震性能研究。Huaman等[18]使用FLAC2D对处于秘鲁高地震活动区、采用下游法筑坝的扩容尾矿坝进行了地震灾害评估。

以往的研究中重点关注尾矿坝地震安全系数或永久位移等指标,忽略了溃坝后滑动体的大小、运移过程以及堆积状态对下游构筑物带来的风险研究,此外,地震过程中导致尾矿坝筑坝材料的弱化,进而会加剧尾矿坝溃坝风险。因此,为合理考量上述因素对尾矿坝地震安全性能影响,现采用极限平衡拟静力法与光滑粒子流体动力学方法(smoothed particle hydrodynamics SPH)联合进行尾矿坝地震可靠度分析[19-27]。考虑地震对尾矿坝筑坝材料的弱化效应,利用极限平衡拟静力法与蒙特卡罗方法寻求可能的失效样本,进而采用自主研发的SPH程序模拟失效样本下尾矿坝失稳过程以及最终堆积状态,以滑动面积评估尾矿坝的风险水平。

1 尾矿坝地震可靠度分析

1.1 极限平衡拟静力法

极限平衡拟静力法是用于评价地震作用下岩土结构物稳定性的方法之一[28-29]。该法忽略坝体的动态响应,将地震作用简化为一静力作用在坝体上,结合极限平衡方法进行稳定性分析时,将地震作用等效为水平和竖向地震力分别施加于每个土条的重心。主要考虑水平地震力的作用,如图1所示,某滑动区域被划为n个垂直土条,其中作用于土条i的水平地震力代表值Qh,i为

Ti为土条底部法向力;Zi为土条底部剪切力;Ei和Ei+1为土条间法向应力;βi为土条底部倾角,(°),i=1,2,…,n ;r为圆弧滑动半径;O为圆弧滑动圆心

(1)

式(1)中:ξ为地震作用效应折减系数,依据《水工建筑物抗震设计标准》(GB 51247—2018)规定[30],除采用动力法计算钢筋混凝土结构外,ξ取0.25;αi为

土条i地震惯性力的动态分布系数,需根据图2(a)所示的坝体地震惯性力动态分布系数示意图确定;ah为水平地震加速度;kh为水平地震力系数,可通过查阅表1、表2确定;g为重力加速度;Wi为土条i的重力。鉴于Slope/W软件中的极限平衡拟静力法未考虑地震惯性力的动态分布,为便于计算,本文中将图2(a)所示地震惯性力动态分布图简化图2(b)所示的矩形分布图。

表1 工程场地地震烈度与地震动峰值加速度(PGA)对照表[31]

表2 工程场地地震烈度与水平地震力系数(kh)对照表[31]

H为坝高,αi为某坝高Hi时地震惯性力的动态分布系数;am为坝体顶部的分布系数

以Bishop法为例,采用极限平衡拟静力方法计算坝体安全系数FS,计算公式为

(2)

式(2)中:φi为土条底部的内摩擦角,(°);λi为土条长度,m;c为土条黏聚力,kN。在给定的尾矿坝筑坝材料参数下,可以根据式(1)和式(2)计算地震作用下坝体的安全系数,为考虑筑坝材料参数的随机性,需要结合蒙特卡罗抽样方法生成一系列样本,分别采用式(1)和式(2)计算每一样本下的坝体安全系数值,进而计算坝体的失效概率和溃坝风险。为提高计算效率,基于Win-BatchTM平台,编写批处理程序调用Slope/W软件进行坝体失效概率计算。

1.2 基于Win-BatchTM语言的失效概率计算

采用1.1节极限平衡拟静力法计算的安全系数确定坝体的极限状态函数,进而结合蒙特卡罗方法计算尾矿坝的失效概率,即

G(X)=FS(X)-1=0

(3)

式(3)中:G(X)为极限状态函数,X=(x1,x2,…,xm)为尾矿坝地震可靠度分析时考虑的岩土体参数变量,譬如黏聚力c和内摩擦角φ等,其概率分布特征通常假定为对数正态分布。依据土体强度参数的平均值、标准差及分布特征进行N次抽样,产生N组随机样本,X1,X2,…,XN。依次将每组随机样本输入Slope/W中计算尾矿坝的安全系数,若Gi(X)<0,则随机样本Xi视为失效样本,统计失效样本个数M,尾矿坝的失效概率pf≈M/N。pf精度计算式为

(4)

式(4)中:pf为尾矿坝的失效概率;Cpf为失效概率pf的变异系数;N为总抽样次数。

基于Win-BatchTM平台编写批处理程序,调用Slope/W软件批量计算尾矿坝的安全系数FS,并结合蒙特卡罗方法计算失效概率的基本流程如图3所示,主要步骤如下。

图3 基于Win-BatchTM语言失效概率计算流程图

(1)确定尾矿坝的几何模型,如坝高、土层分层,土体强度参数和工程场地地震烈度,在Slope/W软件中建立尾矿坝数值分析模型并确定kh,将尾矿坝分析模型保存为.xml源文件。

(2)根据尾矿坝土体强度参数平均值和标准差,抽取N组随机样本,X1,X2,…,XN,将N组随机样本保存为N个.xml源文件,并批量导入N个文件夹中;

(3)通过批处理程序调用Slope/W软件依次计算N组随机样本的安全系数FS1,FS2,…,FSN;

(4)基于尾矿坝的极限状态函数,统计失效样本个数M,估算失效概率pf=M/N。

2 考虑弱化效应的尾矿坝失稳风险评价

2.1 土体强度的地震弱化效应

研究表明,地震时循环往复的剪切作用会导致土体强度发生弱化,即“土体软化”[32-33]。为合理考虑土体强度的弱化效应对尾矿坝失稳风险的影响,定义土体强度的弱化系数为

(5)

式(5)中:su为地震前尾矿坝各土层黏聚力和内摩擦角,可依据实际地勘报告或工程经验确定;s′u为土体弱化后的黏聚力和内摩擦角,土体强度弱化系数η可由现场试验或实验室试验进行测定。本文中假定不同的土体强度弱化系数η分别为0.75、0.5和0.15,结合工程实例探讨其对尾矿坝失稳风险的影响。

2.2 光滑粒子流体动力学方法确定滑动面积

尾矿坝失稳后的运动过程需通过大变形方法来模拟。光滑粒子流体动力学方法作为一种无网格方法,克服了传统有限元方法网格畸变对计算结果的不利影响,通过一系列携带质量、位移、速度等变量的粒子模拟计算区域,各粒子之间的相互作用通过核函数来模拟,基于质量守恒和动量守恒,引入不同宏观本构,该方法可模拟坝体滑动过程。限于篇幅,关于SPH方法的基本原理在此不再赘述,读者可参阅文献[34-37]。本团队自主研发了SPH程序并成功地进行了香港翡翠道滑坡模拟,并与现场监测结果对比验证了所研发程序。除此之外,应用所研发的SPH程序进行了一系列的滑坡过程模拟。

采用SPH方法计算坝体滑动面积的大体思路如下:针对M组失效样本,分别将每组样本值输入自主研发的SPH程序,保存每次迭代计算的粒子滑动位移uj,对粒子滑动位移进行降序排列,u1,u2,…,uq,q为粒子总数。设置粒子的临界位移值θ,当某粒子滑动位移uj超过临界位移值θ时,该粒子滑动面积计入尾矿坝的失稳滑动面积Si。表达式为

(6)

式(6)中:Si为尾矿坝的失稳滑动面积;dj为第j个粒子面积;q为粒子总数;uj为粒子的滑动位移;θ为粒子的临界位移值,可由工程经验给出或者取尾矿坝坝高H的某一分数值,本文中暂取θ=5%H、10%H、15%H。

2.3 尾矿坝失稳风险计算

尾矿坝失稳风险R可通过其失效概率pf和失稳后果C的乘积来表征。采用蒙特卡罗方法计算尾矿坝M种失效模式的失效概率pf,i,由第1节可知,尾矿坝失效概率pf=M/N,其中每种失效模式的失效概率为pf,i=1/N(i=1,2,…,M)。用式(6)计算所得Si量化其失稳后果Ci(i=1,2,…,M),尾矿坝的失稳风险R为

(7)

本文所提尾矿坝失稳风险的计算流程,具体如图4所示。

图4 尾矿坝失稳风险的计算流程图

3 案例分析

3.1 工程概况

云南省斑毛沟扩容尾矿库位于哀牢山脉西侧,属山谷型尾矿库。在原有旧尾矿坝(简称“旧坝”)的下游选址建设高达198 m的新尾矿坝(简称“新坝”),使得原有尾矿库由三级库扩容为新二级尾矿库,总坝高H=230 m,总扩容达1 087.3万m3。根据地勘报告显示,尾矿坝土层分别为坝基、初期坝、子坝、尾粉土和尾粉质黏土,场地地震动峰值加速度PGA为0.1g,地震动反应谱特征周期为0.45 s。扩容尾矿坝施工期间主要分为三个阶段,分别为新库区扩容至原库区初期坝坝顶、新库区扩容至原库区子坝中点、新库区扩容至原库区坝顶。本文分别对扩容尾矿坝的三个施工阶段进行稳定性分析。

3.2 扩容尾矿坝拟静力稳定性分析

为便于数值模拟,简化其坝基部分进行建模(图5)。考虑到尾矿库内水位对尾矿坝稳定性的影响,首先利用有限元Seep/W软件对各施工阶段进行渗流稳定性分析,之后再耦合Slope/W软件进行安全系数的计算。采用极限平衡拟静力法模拟扩容尾矿坝的地震性能,场地地震动峰值加速度为0.1g,根据中国地震动参数区划图(GB 18306—2015)[31]的地震动峰值加速度(peak ground acceleration,PGA)与工程场地地震烈度对照表(表1)确定尾矿坝的地震烈度为Ⅶ度,由《水工建筑物抗震设计标准》(GB 51247—2018)[30]可知,坝顶动态分布系数am=3.0。再依据表2确定水平地震力系数kh为0.1,地震作用效应折减系数ξ=0.25。本算例坝高H=230 m,0.6H高程的动态分布系数取1.0+(am-1)/3,坝底动态分布系数取1.0,等效后的地震惯性力动态分布系数为1.73,为安全起见,最终确定地震惯性力的动态分布系数为2.0。

图5 扩容尾矿坝简化模型

以扩容尾矿坝的三个施工阶段为分析工况,工况一为新库区扩容至原库区初期坝坝顶、工况二为新库区扩容至原库区子坝中点、工况三为新库区扩容至原库区坝顶。根据地勘报告的各土层参数(表3),分别计算三种工况条件下的最小安全系数FS。

表3 扩容尾矿坝土层物理力学参数

图6汇总了三种工况条件下扩容尾矿坝的最小安全系数及其滑动面。其中,工况一下最小安全系数为1.56,其对应的滑动面(称之为临界滑动面)位于旧坝;工况二下最小安全系数为1.60,其临界滑动面位于新坝;工况三下最小安全系数为1.22,其临界滑动面位于新坝。结果表明:随尾矿坝扩容高度的上升,旧坝的稳定性基本保持不变,而新坝的稳定性持续降低,自2.03降低至1.22。该尾矿坝总库容1 087.30万m3,坝高230 m,由构筑物抗震设计规范(GB 50191—2012)[12]尾矿坝抗震等级可知(表4),其抗震等级为二级,由表5可知,规范要求的最小安全系数为1.15。三种工况条件坝体最小安全系数均大于1.15,满足规范要求。但需要注意的是,工况三下,其最小安全系数与规范要求值接近,需要重点关注其稳定性,下节重点针对工况三进行可靠度分析与风险评价。

表4 尾矿坝抗震等级[12]

表5 尾矿坝地震稳定性最小安全系数值[12]

图6 扩容尾矿坝三种工况下尾矿坝地震稳定安全系数

3.3 扩容尾矿坝工况三失效概率计算

考虑尾矿坝筑坝材料的不确定性,采用蒙特卡罗方法对工况三下扩容尾矿坝进行可靠度分析。假定子坝、尾粉土、尾粉质黏土的黏聚力c、内摩擦角φ均为对数正态分布,以表2中的各土层参数为平均值,变异系数均假定为0.1,随机生成1 000组样本。采用图3所示流程计算1 000组随机样本下坝体的最小安全系数并计算其失效概率。图7以最小安全系数为横轴,以该安全系数的频率为纵轴,绘制最小安全系数直方图。由图7可见,最小安全系数的平均值为1.218,标准差为0.104,失效样本个数为10个,其失效概率为1%,由式(4)可知,计算所得失效概率变异系数Cpf≈0.3。

图7 工况三下安全系数的直方图

3.4 扩容尾矿坝失稳风险评价

3.4.1 尾矿坝失稳过程模拟

尾矿坝一旦溃坝,将会导致尾矿库内大量的有害物质流向下游,对下游的居民、环境、水源等造成不可估量的损失。基于此,以第10个失效样本为例,取弱化系数η=0.15,图8分别给出了t=1,3,5,10 s时的滑动状态。t=1 s时,尾矿坝的最大滑动位移为2.17 m,小于尾矿坝的临界位移值θ,此时整体处于稳定状态;t=3 s时,最大滑动位移为18.99 m,主要位于新坝处的子坝和尾粉质黏土部分,旧坝部分的滑动位移为4～6 m;t=5 s时,最大滑动位移为37.79 m,旧坝部分的滑动位移约为8～10 m,此时部分坝体已经越过初期坝,有向下游滑动的趋势;t=10 s时,最大滑动位移为241.56 m,部分坝体已经滑落至初期坝的边缘,扩容尾矿坝整体开始变形。

图8 η=0.15的扩容尾矿坝失稳滑动过程

根据上述尾矿坝的失稳滑动过程,对尾矿坝最大滑动速度进行量化。其中,1 s内坝体尚未滑动,其滑动速度为0,3～5 s内的最大滑动速度为9.40 m/s,5～10 s的最大滑动速度为40.95 m/s。在5 s内,尚未滑至初期坝,滑动速度相对较小,整体处于可控状态。5～10 s内部分尾矿坝逐渐滑落至初期坝的底部边缘,其滑动速度迅速增加,导致滑动位移急剧增大。10 s后滑动速度将持续增大,坝体及库内有害物质将迅速殃及下游部分区域。

3.4.2η=1时尾矿坝失稳风险评价

图9 不同θ下扩容尾矿坝失效样本对应的风险Ri

在本例中,取θ=15%H=34.5 m时,只有当失稳后尾矿坝的滑动位移超过34.5 m时,才会产生风险,如图10所示,10个失效样本下的尾矿坝的最大位移仅为25.69 m,未超过34.5 m,因此风险为0。图10还表明,随着失效样本安全系数的降低,尾矿坝滑动位移值呈非线性增长趋势,当考虑地震弱化效应后,失效样本安全系数会降低,因此需要重点关注考虑弱化效应后的尾矿坝失稳风险量化。

图10 扩容尾矿坝失效样本的最大滑动位移

3.4.3η对尾矿坝失稳风险的影响

为研究η对计算结果的影响,取η=0.75、0.5、0.15,分别采用SPH模拟不同弱化系数下,尾矿坝失效样本的失稳滑动面积,并计算在不同临界位移值下的风险。以土体强度弱化系数为横轴,尾矿坝的风险为纵轴,绘制不同临界位移值下,其风险随弱化系数的变化曲线,如图11所示。临界位移值θ=5%H时,η=1.0、0.75、0.5、0.15时,尾矿坝的失稳风险R分别为309.17、380.31、437.72、842.95 m2;临界位移值θ=10%H时,η=1.0、0.75、0.5、0.15时,尾矿坝的失稳风险R分别为36.14、300.97、418.23、824.74 m2;临界位移值θ为15%H时,η=1.0、0.75、0.5、0.15时,尾矿坝的失稳风险R分别为0、205.12、398.84、807.16 m2。

图11 不同θ下扩容尾矿坝风险随η的变化曲线图

图12 不同η下θ取值对扩容尾矿坝风险的影响

4 结论

结合某扩容尾矿坝工程实例,考虑地震对尾矿坝筑坝材料的弱化效应,采用极限平衡拟静力法与光滑粒子流体动力学方法结合进行尾矿坝地震可靠度分析及风险评价。利用极限平衡拟静力法与蒙特卡罗方法寻求可能的失效样本,进而采用自主研发的SPH程序模拟失效样本下尾矿坝失稳过程以及最终堆积状态,以滑动面积评估尾矿坝的风险水平,并得出以下结论。

(1)随尾矿坝扩容高度的上升,旧坝的稳定性基本保持不变,而新坝的稳定性持续降低,拟静力法所得安全系数自2.03降低至1.22。

(2)在工况三下,即扩容高度完成后,尾矿坝最小安全系数与规范要求值接近,需要重点关注其稳定性,经可靠度分析,当尾矿坝筑坝材料变异系数为0.1时,其相应的失效概率为1%,该失效概率值的变异系数为0.3左右。

(3)不考虑地震作用对土体的弱化效应,扩容尾矿坝的失稳风险R=309.17 m2,尾矿坝的失稳风险R随临界位移值θ的增大而减小。

(4)考虑弱化效应后,随弱化系数减小,扩容尾矿坝的失稳风险增大;弱化系数大于0.5时,临界位移值θ对扩容尾矿坝失稳风险影响显著。