地震作用下高尾矿坝三维有限元动力分析

李 瑜 宁 伟 吴蒙蒙 林 杰

（1.栾川县环境保护局；2.洛阳栾川钼业集团股份有限公司；3.河海大学水利水电学院）

金属冶炼及火力发电过程中产生的大量尾矿及其他工业废渣，往往采用筑坝拦截谷口等方式形成尾矿库进行贮存［1］。尾矿坝作为一种人造的、具有高危势能的特殊工业构筑物，不仅关系到矿山企业的经济效益，而且对其下游居民人身及建筑设施构成极大威胁［2］。地震作用下尾矿坝易发生溃决及坡地液化，诱发泥石流等地质灾害［3］，为确保尾矿坝施工及运行期的安全，开展尾矿坝动力分析具有重要的实际意义。

地震的随机性与突发性和坝体空间结构的非均匀性及变异性，致使尾矿坝抗震问题十分复杂［4］，国内外学者主要就地震液化、残余变形及坝体稳定等进行了研究。Newmark［5］认为，最小安全系数不能决定坝体在地震作用下的稳定与否；汪闻韶［6］分析了不同机理下的土体液化，认为土体极限平衡状态将先于土体液化的发生；徐志英等［7］分别采用总应力法、不排水有效应力法及排水有效应力法分析德兴铜矿4号尾矿坝的动力特性，并进行安全评价。本研究基于等效非线性黏弹性模型，针对洛阳某拟建高尾矿坝，分析其在动力状态下坝体动加速度、有效应力、动位移、地震永久变形、地震液化区域及干滩长度变化，评价其抗震安全性。

1 工程概况

拟建尾矿库位于河南洛阳栾川县境内，该库初期坝为碾压堆石坝，坝顶标高约为1 160 m，坝高为79 m，坝顶宽为4 m，上、下游边坡比为1∶2.0，坝轴线长约为165 m，初期库容为401×104m3；后期采用上游法堆坝，设计采用粗尾砂筑子坝，每级子坝高4 m，顶宽2 m，上、下游边坡比为1∶2.5，设计堆积标高为1 280 m，平均堆积边坡比为1∶5.0，尾矿堆高为120 m，总坝高为199 m，总库容为5 742.5×104m3。尾矿库等别为二等，其主要构造物级别为2级，次要建筑物为3级，临时建筑物为4级。

2 计算工况及参数

2.1 材料分区及边界条件

根据提供的工程资料，将尾矿坝堆积坝分为尾黏土、尾粉质黏土、尾粉土、尾粉砂和尾细砂5个区，加上初期坝及基岩共有7种材料分区，如图1所示。

模型边界条件：对上、下游边界，左、右岸边界和模型底边界进行位移约束，模型上边界自由。在综合分析计算区域内地形、坝体结构等特征的基础上，生成有限元网格31 521个，结点47 125个，三维模型网格图如图2所示。

2.2 静力计算参数

进行尾矿坝动力计算时，需先进行静力计算。初期坝及尾矿堆积坝区材料按非线性材料考虑，均采用邓肯模型；基岩按线性材料考虑，采用线弹性模型。部分静力计算参数如表1所示。

?

2.3 动荷载输入

“正常蓄水位＋地震”作为坝体地震反应的计算工况。以50 a超越概率5%地震波为基本设计工况，对应的基岩地震震动水平峰值加速度取52.01 cm/s2。选择类似场地地震条件的3条实测加速度记录，并取其计算结果的最危险组合。动加速度输入曲线如图3所示，地震波时长为40 s，时间步长为0.02 s，三向地震波的峰值时刻基本都在16 s左右。

2.4 动力计算说明

动力平衡方程为

式中，δ,δ̇,δ̈分别为结点位移、结点速度和结点加速度；F(t)为结点的动力荷载；M为质量矩阵，用集中质量法求得，即假定单元的质量集中在结点上；K为劲度矩阵，用常规有限元法得；C为是阻尼矩阵，C=为第一振型自振频率；λ为阻尼比；t为时间。

动力计算分析采用等效非线性黏弹性模型，即假定尾矿堆积坝为黏弹性体，相应尾矿堆积坝区域的动应力应变关系采用等效剪切模量G和等效阻尼比λ分别代表其非线性和滞后性［8］。主要计算公式如下。

式中，G、Gmax分别为等效剪切模量和初始剪切模量；λ、λmax分别为等效阻尼比和初始阻尼比；γd、γr分别为动剪应变及参考剪应变；Pa为大气压；σ'm为有效平均应力；K'、n、λmax3个参数由动三轴试验确定。

残余剪应变γp、残余体应变εpv和动剪应力Δτ的关系根据中国水利水电科学研究院进行的坝料体积变形特性的大型动三轴试验结果，采用如下公式进行计算。

式中，μ为动泊松比；Ka，na，Kv，nv为试验参数，分别是以有效固结应力σ'3、固结比Kc、等效振次N为参变数的系数和指数；残余轴应变εpa、残余体应变εpv以百分数表示，动剪应力Δτ和有效固结应力σ'3采用相同的量纲。

部分动力计算参数如表2所示。

?

2.5 孔隙水压力说明

采用徐志英等［9-10］，根据试验资料整理出的可用于三维应力状态的振动孔隙水压力增长模式，在每一时段末增加了残余孔隙水压力或残余变形的计算。对于孔隙水压力，其初始平均有效应力取为出于安全考虑，短时地震条件下孔隙水压力的扩散与消散很少，本计算中对此部分进行忽略；采用最大孔隙比为最大孔隙水压力，为平均有效应力）作为判断初始液化的标准，即最大孔隙水压力比达到1时粉土质砂液化［11］。

3 动力计算结果与分析

表3分别给出了三维有限元计算成果下的坝体最大加速度反应、应力反应、最大位移反应、地震永久变形的最大值。

3.1 加速度反应

坝体动加速度反应如图4所示。坝体顺河向的最大绝对加速度值为1.64 m/s2，放大倍数为3.15；垂直向的最大绝对加速度值为1.11 m/s2，放大倍数为2.13。顺河向、垂直向加速度反应最大值等值线与坝坡近似呈平行分布，随坝体高程的增加而不断增大，最大值出现在初期坝顶和坝体中部高程坝坡处。此外，在堆积坝坝体内部存在局部最大加速度先减小再增大的现象，主要是由于地形、堆积坝坝坡过缓、材料分区和浸润线的影响。

3.2 动应力反应

坝体动应力反应如图5所示。堆积坝坝体最大第一动主应力为2 422.15 kPa，最大第三动主应力为1 823.01 kPa。从坝体典型断面上动主应力分布可知，堆积坝动主应力在坝体内部分布较为均匀且全为压应力，最大动主应力从坝体表面向内部不断增大，极大值发生在坝基附近，并受地形影响，出现局部应力集中现象。

3.3 动位移反应

坝体动位移反应如图6所示。堆积坝顺河向最大动位移为50.57 mm，垂直向最大动位移为15.13 mm。堆积坝顺河向、垂直向最大动位移值均发生在坝顶附近。坝顶各点数值接近，从坝顶向下动位移反应减小。从典型断面的动位移反应分布来看，其动位移反应不显著，其中垂直向的动位移反应较小，顺河向动位移相对垂直向动位移较大。由于坝基良好，顺河向、垂直向的最大动位移值较小且较为接近，坝体各部分变形均匀，故地震情况下堆积坝表面基本不会出现裂缝。

3.4 地震永久变形

坝体地震永久变形如图7所示。堆积坝坝体顺河向地震永久变形近似与坝坡呈平行分布，最大值为125.12 mm，位于顶部坝坡处；垂直向地震永久变形分布均匀且随高程增加而不断增大，最大值为-124.24 mm，位于堆积坝顶，按堆积坝最大坝高199.0 m计算，地震永久沉降约为坝高的0.06%。

正常运用条件下，尾矿坝沉积滩滩顶标高为1 280.0 m，超过尾矿库沉积池正常运用水位的1 275.5 m。在地震作用下，堆积坝坝体（含沉积池）的最大沉陷（永久变形）约0.124 m，地震涌浪高约1 m。考虑最大沉降后，尾矿坝沉积滩滩顶标高仍超过正常运用水位3.376 m，大于规范的最小安全超高为1.0 m；干滩长度为225.1 m，大于规范规定的100 m［12］。因此，不会发生地震作用导致沉积池水漫溢尾矿坝的情况。

3.5 地震液化分析

地震下坝坡特征单元孔压增长如图8所示，液化区域如图9所示。随着地震的发生，坝坡特征单元的孔隙水压力不断增大（忽略孔隙水压力的消散与扩展）；堆积坝在上游坝顶、浸润线以下附近出现小范围的地震液化区域，沉积坝约75 m宽度范围，深度为上表面2层网格，主要处在尾黏土和尾粉质黏土材料分区处。此处材料位于浸润线以下，含水率较高。由具体液化区域分布可看出，堆积坝位于整个库盆内，小范围的地震液化区域不会对尾矿库的整体稳定性造成很大影响。

4 结论

（1）坝体动加速度反应、有效应力反应及动位移等值线近似与坝坡平行且受地形影响，最大值位于初期坝顶或堆积坝顶部；堆积坝顺河向动加速度反应、动位移反应皆大于垂直向反应。

（2）地震永久变形总体随高度的增加而增大，水平向永久变形最大值位于顶部坝坡处，垂直向永久变形最大值位于堆积坝顶；其中最大沉降永久变形为0.124 m，约为坝高的0.06%；考虑地震涌浪后，坝顶安全超高及最小干滩长度仍满足规范要求。

（3）地震条件下，坝体单元孔隙水压力不断增大（忽略孔隙水压力的消散和扩展），堆积坝上游浸润线以下部分出现小范围液化区，对坝体稳定性基本无影响。