混凝土的单轴拉压行为由塑性损伤描述，分别采用拉伸损伤因子σ_t和压缩损伤因子σ_c来表示混凝土受到拉伸和压缩时的刚度退化和强度下降的特点^[8-9]。具体公式如下所示^[10]：

式（1）中：s为应力；d为损伤变量（0≤d≤1）；e为应变； 为塑性应变；D^el₀为无损弹性刚度。

损伤后的弹性模型表示为：

E=( 1 - d ) D ^{e l}₀（2）

混凝土塑性损伤模型采用非关联流动法则，其塑性势函数为^[11]：

式（3）中： ξ为混凝土塑性势函数的偏心率； σ_{t 0}为破坏时的单轴应力；ϕ为混凝土屈服面在强化过程中的膨胀角，本文取36°。

在地震荷载作用下，坝体—基岩系统的失稳主要表现为坝基滑裂面的粘结强度逐渐破坏直至丧失，处在粘结状态的区域逐渐缩小，进一步发展到滑移、张开，最终因整个滑裂面都处在滑动状态而导致基岩滑裂面的震后错动位移增大，导致整体失稳^[12-15]。由于地震荷载的随机性和往复性，地震中不同时刻的接触面状态不断变化，故很难以地震过程中某一时刻的接触面状态作为坝基失稳的判定标准^[16]。但是，震后的接触面的状态及关键点错动位移，仍然可以作为坝基稳定的参考指标^[17]。这种方法可以全面考虑基岩中滑裂面的整体状态，避免局部区域位移变化的影响。

在非限性动力接触模型中，接触面应满足法向接触条件和切向接触条件。法向接触条件是：①法向不可贯入性条件；②法向接触力为压应力条件。法向接触条件可表示为^[18-19]

g≥0; λ≥0; gλ=0（4）

式（4）中：λ为两个接触面间的法向接触力（接触力以受压为正），g为接触面法向间距。当λ=0，并且g>0时，表示接触面张开；当λ>0，并且g=0时，表示接触面闭合λ>0。

切向接触条件是：设接触面的切向摩擦满足库伦摩擦定律，则切向接触力f_cT的大小不能超过它的极限值μλ。即

引入一无量纲变量τ，τ= ||f_cT/μλ。接触面间相对滑动的速度为u̇T，定义为

在考虑震后滑裂面接触状态的同时，也可参考坝体—基岩系统震后坝基塑性区分布，综合评价大坝的动力稳定性。

某混凝土重力坝，坝高142.00 m；坝宽122.20 m，坝段厚22.00 m，坝基存在一条明显的滑裂面。为此考虑两种滑移模式建立有限元模型。

坝体与基岩均采用3D-Solid单元，坝基网格充分考虑了岩基的倾斜夹层特性和岩体材料分区，根据岩层走向分层建模，滑移面采用接触单元模拟，上、下游方向地基及其深度均取1.5倍坝高。整个模型共划分单元18 738个，节点22 573个，如图1所示。坝体采用混凝土塑性损伤材料，坝基岩体采用采用Mohr－Coulomb屈服准则。根据水利水电工程抗震设计规范^[20]要求，混凝土动态弹性模量标准值较表1中的静态弹性模量标准值提高50%。三种标号的混凝土受拉应力应变关系曲线如图2所示，损伤因子与应变关系如图3所示。坝基岩体参数如表2所示。

图1 三维有限元动力计算模型Fig.1 Three-dimensional finite element dynamic calculation model

图2 混凝土受拉应力应变曲线Fig.2 Tensile stress-strain curves of concrete

图3 混凝土损伤因子与应变关系曲线Fig.3 Relationship between damage factor and strain of concrete

表1 坝体材料参数表Table 1 Material parameters of dam body

表2 坝基岩体材料参数Table 2 Material parameters of dam foundation rock mass

计算荷载包括坝体自重、上游和下游的静水压力、坝基扬压力、库底的淤沙压力等，其中：上、下游水位分别为138 m和58.41 m；坝基面扬压力按规范要求选取；淤沙的浮容重6.0 kN·m^-3，内摩擦角12°。

坝址场地地震基本烈度为Ⅷ度，抗震类别为甲类。设计地震条件下，基岩水平向峰值加速度为0.3 g，竖向峰值加速度取水平向的2/3。以场地谱为目标谱，拟合顺河向和竖向地震波，地震动持续时间为30 s，计算时间步长0.01 s，输入的地震加速度时程曲线如图2和图3所示。对应的加速度谱与目标谱的对比如图6、图7所示。采用Westergaard附加质量法来模拟地震过程中坝体和库水的动力相互作用^[21]。坝体动力计算阻尼采用瑞利阻尼法，根据规范要求，坝体阻尼比取10%。采用地震强度超载法计算动力稳定性时，分别对设计地震的峰值加速度超载1.0倍、2.0倍、2.5倍、3.0倍、4.0倍。

在研究地震荷载作用下的坝基动力稳定时，考虑基岩辐射阻尼特性将更接近实际情况。本文坝基模型的边界上建立粘弹性人工边界，即在底边界和侧边界施加弹簧阻尼单元来吸收系统在地震过程中向外辐射的能量，进而模拟地基辐射阻尼的影响。黏弹性边界上的弹簧和阻尼单元参数可以采用下式表示：

式（7）中：K_N、K_T分别为弹簧的法向和切向刚度；C_N、C_T分别为阻尼器在法向和切向上的阻尼系数；r为波源至人工边界的距离；α_N、α_T分别为法向与切向黏弹性人工边界修正系数。

采用波动输入法，在边界以等效节点力的形式实现地震动输入。

图4 水平向加速度Fig.4 Horizontal acceleration

图5 竖向加速度Fig.5 Vertical acceleration

图6 水平向加速度反应谱Fig.6 Horizontal acceleration response spectrum

图7 竖向加速度反应谱Fig.7 Vertical acceleration response spectrum

在设计地震加速度下，坝体几乎无损伤，整体均处于完好状态（图8a）；当峰值加速度为2倍设计地震时，下游坝面首先出现损伤，损伤初步向上游坝面发展，在重力坝上部逐渐形成一条开裂通道（图8b）；当峰值加速度提高到2.5倍设计地震时，下游坝面的损伤区域较2倍设计地震时进一步扩大，开裂通道继续向上游面发展，但未完全贯穿坝体（图8c）；当峰值加速度提高到3倍设计地震时，开裂通道贯通了坝体，并且下游坝面的损伤区域进一步扩大，并且形成了多条开裂通道（图8d），此时表明大坝已经处于极限破坏状态，因此可以判定大坝可以承受3倍设计地震强度而不发生溃坝。

在坝基非线性动力接触模型中，五种接触面状态，分别是粘结状态、滑动状态、闭合状态、张开状态、失效状态。

在设计地震加速度下及对设计地震进行超载2.0倍、2.5倍、3.0倍下的震后滑裂面状态分布如图4所示。从图9（a）可以看出，在设计地震时，沿着坝踵垂直向下的接触面首先张开，处于张开状态。坝踵以下的滑裂面大部分处于粘结状态，只有局部区域的滑裂面仍处于稳定状态。当超载3.0倍时，整个滑裂面的粘结强度均已丧失，相当一部分滑裂面处于张开状态，还有一部分滑裂面处在滑移状态，此时坝体处于失稳的临界状态。因此，从接触面状态来判断，坝体—基岩系统的超载安全系数为3.0。

随着超载倍数的增加，基岩塑性区分布如10所示。当承受设计地震时，只有在软弱夹层和坝趾局部区域发生了塑性屈服，滑裂面处的岩体处于弹性状态。当承受2.5倍设计地震时，塑性区沿滑裂面从上游到下游发展，同时坝趾处的塑性区范围扩大。当超载倍数达到3.0倍设计地震时，沿滑裂面塑性区基本贯通，坝趾部位出现较大的塑性区。当超载倍数达到4.0倍时，塑性区范围进一步扩大。综合塑性区的动态变化，可知系统的超载安全系数为3.0。

随着地震强度的超载，坝踵和坝趾位移如图 11所示。从图中可以看出，出当超载倍数在2.0~3.0之间时，曲线开始出现拐点，随后位移迅速增长。但在超载2.0倍和2.5倍时，虽然位移变化速度加快，从图9（b）-（c）可知，此时仍有相当一部分接触面处在粘结状态，表示滑移面的粘聚力尚未完全破坏，坝体仍然能保持稳定而不会出现失稳破坏。当超载3.0倍时，位移变化量进一步增大，从图9（d）可知，此时整个滑裂面装态均处于张开或滑动状态，此时表征坝体已经处在临界状态。当超载倍数到达4.0~4.5时，震后位移已经很大，表征系统已经失稳。故综合位移变化和关键点位移特征，超载安全系数为3.0。

图8 不同强度地震作用下坝体损伤分布图Fig.8 Damage distribution of dam under different intensity earthquake

图9 不同强度地震作用下接触面状态分布图Fig.9 State distribution of contact surface under different intensity earthquake

图 10 不同地震强度下坝基塑性区分布Fig.10 Distribution of plastic zone of dam foundation under different earthquake intensity

图 11 关键点震后错动位移与超载系数关系曲线Fig.11 The relationship curve between dislocation displacement and overload coefficient of key points after earthquake