基于双滑面模型的混凝土重力坝深层抗滑动力稳定分析方法

龚亚琦，苏海东，陈 琴

(1.长江科学院 材料与结构研究所，武汉 430010； 2.长江科学院 国家大坝安全工程技术研究中心, 武汉 430010)

1 研究背景

大坝设计过程中，坝址岩体普遍存在软弱夹层、缓倾角裂隙等地质缺陷，对于抗剪强度低的岩体，坝体有可能在坝基内部沿软弱结构面滑动,引发大坝深层失稳。特别是在地震作用下，坝基的深层滑动稳定问题更为突出，因此，开展地震条件下坝基软弱结构面的深层抗滑动力稳定分析非常必要[1-2]。以小南海重力坝为例，前期工程勘察表明，坝基下发育有缓倾下游的剪切带，且坝后抗力体中存在连通率不大的缓倾角结构面，坝体可能沿剪切带及坝后的缓倾角结构面、岩体破坏面发生滑动，深层抗滑问题突出。

在实际工程中，深层抗滑稳定的计算主要采用刚体极限平衡方法和强度折减法[3-4]。其中，强度折减法能考虑材料的塑性破坏过程且无需设定滑裂面，但是由于涉及非线性计算，计算量较大；刚体极限平衡法的双滑面滑动模式是设计规范推荐的主要滑移模式[5]，该方法模型简单，计算方便，适合有特定滑动面的深层抗滑问题。

目前，有限元方法与刚体极限平衡法相结合的方法在坝基深层抗滑稳定计算中应用较为广泛[6-8]。然而，应用刚体极限平衡方法时，需要得到滑动面和阻滑面上的各分项作用力荷载，但是，有限元块体单元只能给出应力，无法直接给出计算需要的内力。一般做法是直接通过应力积分得到内力，而应力积分方法对单元网格形态要求较高，积分精度受有限元单元尺寸和单元类型的影响较大。同时，对于深部岩体地震作用中还存在辐射阻尼效应，包含不同岩层的非均匀岩体，为了简便，使用黏弹性等地基模型进行辐射阻尼计算时，通常先将地基等效为均质材料。然而，地震波在不同结构面内的传播特性、地震动荷载的施加都存在一定的差异，这些因素会对计算结果产生何种影响，如何评估，是计算分析过程中需要关注的问题。

本文探讨了在地震作用下应用刚体极限平衡方法和有限元法进行重力坝深层抗滑稳定分析时的相关问题，研究了计算过程中模型的建立、荷载的施加、结果的处理和稳定分析步骤等相关方法，并结合小南海水电站工程实例进行了验证。

2 理论与方法

2.1 黏弹性地基模型

抗震计算实例表明：辐射阻尼效应对地震响应的影响可达20%～40%，并随坝体体积增大和地基变形模量降低而增长[9]，因此，在进行强震作用下坝基深层抗滑稳定分析时，考虑地基辐射阻尼特性的影响将更加符合实际[8]。黏弹性地基模型通过在边界上设置单向弹簧和阻尼器吸收系统在振动过程中向外辐射的能量[10]，可以很好地模拟岩体的辐射阻尼特性，同时该方法与有限元法结合方便，在地震动分析中应用较为广泛[11-12]。

2.1.1 地震动等效荷载

杜修力等[11]通过等效节点力的方式实现了黏弹性地基模型中的地震动输入，模型中考虑了地表侧向边界反射波的吸收影响。采用相同的思路，可以进一步推导出考虑底部边界地表反射波吸收影响的地震动等效荷载公式[13]。

(1)对于S波从底部垂直入射的情况，其公式为：

u′S(t-Δt1)]+AiρcS[u′S(t)-u′S(t-Δt1)] ；

(1)

(2)

(3)

(2)对于P波从底部垂直入射的情况，其公式为：

u′P(t-Δt4)]+AiρcP[u′P(t)-u′P(t-Δt4)] ；

(4)

式中：F为等效地震荷载，其上标表示边界的外法线方向，与坐标轴方向一致，x，y，z分别为地基水平侧向和竖向，下标表示节点号和节点力分量；Ai为人工边界上节点的影响面积；K，C分别为边界弹簧刚度和阻尼，其中下标t表示切向，n表示法向；ρ为密度；cS，cP分别为S波和P波的波速；uS，u′S和uP，u′P分别为S波位移、速度和P波位移、速度；λ为拉梅常数；t为时间；Δt1为地震波从地基底部入射到达地表后反射至底部所用时间；Δt2，Δt′2分别为x向两侧边界节点处入射S波和地表反射S波的到达时间；Δt3，Δt′3分别为y向两侧边界节点处入射S波和地表反射S波的到达时间；Δt4为P波从人工边界底部入射到达地表后反射至底部所用时间；Δt5，Δt′5分别为x向两侧边界节点处入射P波和地表反射P波的到达时间；Δt6，Δt′6分别为y向两侧边界节点处入射P波和地表反射P波的到达时间。

2.1.2 软弱结构面中地震波传播特性

对于深层抗滑问题，基岩内部往往被多条剪切带和软弱夹层等结构面切割而形成不同的材料区域，不同材料中地震波传播过程中波速和阻尼均不同，应用上述等效节点力公式时需要特别考虑。

本文采用分段处理的方法，即分段计算从人工边界底部垂直入射地震波进入不同材料区域的入射波和反射波到达的时间，应用式(1)—式(6)，分别计算不同材料区域内部不同节点处的等效地震荷载，如图1所示。该方法可以最大限度地考虑材料分区的影响，与等效均质材料方法相比，计算结果更接近实际地质环境。

图1 结构面中地震波传播时间示意图Fig.1 Schematic diagram of seismic wave propagationtime in structural plane

具体步骤如下：

(1)根据结构面的分布将地基沿高程、顺水流向和横河向分为不同的材料区域。每个材料区域内设置对应的波速、弹簧刚度和阻尼等参数。

(2)假设从人工边界底部垂直发射一条射线穿透各材料区域，分别记录入射波和反射波到达材料边界的时间。

(3)搜索不同边界处节点坐标，并根据记录的边界点处入射波和反射波的时间计算不同节点处Δtij值(i表示节点号，j表示基岩材料分区)。

(4)代入式(1)—式(6)，计算不同节点处等效地震力，将其作为地震动荷载施加在对应区域的节点上。

2.2 深层抗滑稳定分析方法

深层抗滑稳定问题主要采用双滑面模型和单滑面模型，基于刚体极限平衡方法，将外部荷载沿滑动面分解为滑动力和阻滑力，计算抗滑稳定安全系数。下面以双滑面模型为例，简要说明深层抗滑稳定分析方法。

2.2.1 双滑面模型

双滑面模型是目前深层抗滑稳定分析中常见的缓倾角结构面的滑移模式[5]，如图2所示。分析过程中，采用被动抗力法，首先假定抗力体BCD达到极限平衡状态，按照传统的刚体极限平衡方法，根据有限元静、动力计算成果，计算推力Q，然后再计算作用于ABD滑动体上的荷载(如水压力、自重、扬压力、地震力等)，将这些作用力沿滑动面分解成切向和法向分力，求得沿该滑动面抗滑稳定整体安全系数。

图2 双斜滑动面模式示意图Fig.2 Sketch of double oblique sliding surfaces

地震作用下的抗滑稳定分析具体计算过程如下：对于地震输入过程中的每个时刻，将有限元法计算得出的外力荷载沿指定结构面的滑裂面的切向和法向分解，得到滑动力和阻滑力，再按照抗力函数和相关公式，计算沿该滑裂面方向的抗滑稳定整体安全系数。

2.2.2 作用力的计算

应用双滑面模型时需要计算滑动体和抗力体上的作用力荷载，由于有限元计算中实体单元只能提供应力和位移成果，为了整理截面作用力，目前最常用方法是应力积分法，即直接将截面上的应力通过数值积分计算滑动力和阻滑力。由于受应力集中和应力精度的影响，滑裂面上的切向力与外力平衡的精度不高，并且积分结构受单元形态影响较大。

苏海东等[14]提出直接应用有限元节点力计算截面内力的方法：根据有限元计算的位移结果和刚度矩阵直接计算节点力。根据有限元理论，每个有限单元的平衡方程为

(7)

取某一截面一侧的单元为脱离体进行处理，则

(8)

从式(8)可以看出，内力计算的精度仅依赖于位移的结果，而对有限元计算而言，位移的精度要高于光滑后的应力结果，同时，位移结果受单元形态的影响相对较小，因此相比应力积分方法，本方法的计算精度较高。

3 工程概况

小南海重力坝最大坝高68 m。采用100 a基准期超越概率2%，相应基岩水平地震动峰值加速度的最大值为1.06g；壅水建筑物采用100 a基准期超越概率1%进行抗震校核，相应基岩水平地震动峰值加速度的最大值为1.20g。坝基内普遍存在软弱夹层、缓倾角裂隙等，且抗剪强度低。坝基以下30 m范围内主要分布有多条结构面，以其中一条剪切带(3S620)为例，如图3所示，可能的滑动模式有2种：双滑面滑动模式，沿A-B-C-O-D滑动；单滑面模式，沿A-B-C-O-E滑动。

图3 结构面分布示意图Fig.3 Schematic diagram of structural plane distribution

4 重力坝深层动力抗滑分析

采用场地谱反演的人工地震波进行时程法动力分析。计算过程中，剪切带按地质勘测的走向、产状和倾角建立模型，模型中滑动带的抗剪断综合参数是考虑剪切带连通率、层理面和混凝土置换的综合作用平均取值，抗力体抗剪断综合参数按缓倾角结构面连通率与岩桥参数平均取值。

对于双斜滑面模式，抗滑稳定计算首先需要确定“第2滑动面的倾角β”和“起点距坝趾距离L0”2个关键参数(如图3所示)，该参数值可以通过敏感性计算分析获得。

4.1 敏感性分析

首先采用平面模型进行滑裂面参数的敏感性分析，确定双斜滑面模式参数β和L0，平面模型如图4所示。

4.1.1 第2滑动面的倾角分析

结合地质条件，分别选取第2滑动面倾角β=30°，35°，40°，45°，50°，其中第2滑动面从坝趾平面滑出，计算地震作用下的最小深层抗滑安全系数K。不同β角条件下最小深层抗滑安全系数见表1。

表1 不同β角条件下最小深层抗滑安全系数Table 1 Minimum values of deep anti-sliding safetyfactor in the presence of varying inclined angle β ofsliding surface

从表1可以看出，随着第2滑动面倾角的不断增加，无论采用极限状态公式还是抗剪断公式，结构面的最小安全系数均呈现先减小后增大的趋势，其中β=40°时，安全系数最小。

4.1.2 第2滑动面起点分析

选定第2滑动面倾角β=40°，比较不同起点与坝趾距离L0=0，5，10，20 m时地震作用下的最小深层抗滑安全系数，计算结果如表2所示。其中，滑动体和抗力体的抗剪断系数根据滑动面岩体性质平均取值。

从表2可以看出，在第2滑动面起点从坝趾平面向下游移动过程中，安全系数最小值变化不大，同时，结合其他工程实例，假定沿坝趾平面滑出。

表2 不同起点距离L0条件下最小深层抗滑安全系数Table 2 Minimum values of deep anti-sliding safetyfactor in the presence of varying distance L0 fromsliding starting point to dam site

图5 三维有限元模型Fig.5 Three-dimensionalfinite element model

4.2 最小安全系数

根据敏感性分析结果，第2滑动面按倾角40°从坝趾平面滑出考虑，进行三维有限元分析，有限元模型如图5所示。

双滑面模式和单滑面模式深层抗滑动力安全系数曲线如图6所示。计算结果表明：校核地震作用下，3S620剪切带按照双滑面模式深层抗滑稳定最小安全系数为3.036，按照单滑面模式最小安全系数仅仅在3.88 s时刻为2.286，其他时刻均>2.3。

图6 校核地震作用下3S620深层抗滑稳定安全系数Fig.6 Safety factor of deep anti-sliding stability ofstructural plane 3S620 under the action ofchecking earthquake

4.3 地基均匀性的影响

下面分别从结构应力变形和深层抗滑稳定安全系数几个方面，对比地基均匀性的影响，结果见图7。

图7 不同地基模型条件下结构应力、变形和深层抗滑稳定安全系数的比较Fig.7 Comparison of stress, deformation, and deepanti-sliding stability between different foundationmodels

由于采用的黏弹性地基模型，相比无质量的地基应力能减小30%～50%，同时坝址地震加速度峰值较小，所以，地震作用下，动应力水平不高。相比均匀地基，非均匀地基条件下，齿槽部位动应力增幅约20%～30%，而坝顶相对动位移则变化不大，如图7(a)和图7(b)所示。可见，对于黏弹性地基模型，采用均匀地基和非均匀地基参数，对基础附近应力有一定的影响，而对结构整体位移则影响相对较小。从图7(c)和图7(d)可以看出，无论单滑面模式还是双滑面模式，地基均匀性对稳定安全系数的影响均不大，这主要是由于规范公式的推导是基于刚体极限平衡方法，其出发点从整体考虑结构的稳定，因此，基岩的均匀性影响不大。

5 结 论

根据刚体极限平衡方法和有限元基本原理，研究了应用双滑面模型进行混凝土重力坝深层抗滑动力稳定分析的方法和流程：在黏弹性地基模型中根据剪切带和软弱夹层中不同材料确定地震波的传播时间，通过等效荷载实现地震动输入;根据有限元计算的位移结果和刚度矩阵直接计算节点力，确定每个地震作用时刻滑动面和阻滑面上的作用力荷载;应用规范给出的双滑面和单滑面模型计算最小安全系数，其中，对于双滑面模型通过敏感性分析确定倾角和滑动面起点，最终实现地震作用下重力坝深层抗滑动力稳定计算。计算结果表明：

(1)采用本文提出的计算流程进行地震荷载作用下混凝土重力坝的深层抗滑稳定分析是可行的。

(2)由于刚体极限平衡方法本身是材料力学的方法，从整体考虑结构的稳定，因此，对于动力深层抗滑稳定分析，基岩的均匀性影响不大，采用均匀基岩进行黏弹性分析也是可行的。