地震波斜入射下混凝土重力坝的塑性损伤响应分析

李明超，张佳文，张梦溪，闵巧玲，史博文

（1.水利工程仿真与安全国家重点实验室，天津大学，天津 300350；2.华电重工股份有限公司，北京 100070）

1 研究背景

我国西南部强震区拥有全国约80%的水能资源，随着近年来西南地区不断修建200米级、300米级的高坝，对高坝大库的抗震安全评价提出了严峻的挑战[1-2]。以往大坝安全性评价在地震输入机制方面考虑略有不足[3]，高坝坝趾地震动参数确定方法复杂[4]。用有限元法模拟结构地震响应时，多将坝体-地基体系作为封闭系统，并采用一致性激励方法输入地震，未考虑到行波效应和地基能量辐射。为此，Deeks等[5]推导了二维时域人工边界来模拟无限地基辐射阻尼；王进廷等[6]验证了考虑地基辐射阻尼对坝体非线性地震反应的必要性；刘晶波等[7-8]结合球面波动理论推导了三维黏弹性人工边界，并将波动问题转换为等效荷载的输入。波动输入方法较封闭系统改进较大，但由于地震波经过地壳中复杂介质时要进行多次折射、反射，很难确定入射方向，目前多假定为垂直向上的平面体波进行输入，这对于远场波动是合理的，但震源距场地较近时地震波通常是倾斜入射的[9]，地震动呈现出更为复杂的空间变化特性。随着大型结构和建坝地区地形的复杂性，单一方向的垂直入射无法真实反映地震输入状态，且地面运动的非一致变化对大型结构的影响难以预测。廖河山等[10]利用特征线法分析了地震SH波斜入射时层状半空间的动力响应，结果表明入射角度对地表加速度放大倍数有显著影响。Heymsfield等[11]利用边界积分方程法求解，指出倾斜基岩在SH波斜入射角度为60°时位移幅值最大。在大跨度桥梁结构[12]、复杂场地条件[13]的研究中也证明了斜入射的不利影响，因此在坝体抗震评价中有必要引入对地震入射角度的考虑。

在这方面一些学者开展了相关研究并取得了进展，苑举卫等[14]将地表地震动时程分量分解为斜入射的平面SV波和P波，证明了斜入射对重力坝结构影响明显，尤其在坝-基交界面上；岑威钧等[15]分别研究了SV波和P波入射下高面板堆石坝的地震反应特性，证明了斜入射对加速度反应谱幅值和坝体局部稳定性的影响；杜修力等[16]结合小湾拱坝研究了SV波斜入射下拱坝的地震响应，表明入射角度对低频区幅值放大系数影响显著；孙奔博等[17]研究了平面SV波斜入射下重力坝的地震响应，证明入射角度存在对坝体位移和应力有显著影响；何卫平等[18]研究了斜入射地震动确定性空间差异对重力坝动力响应和破坏模式的影响。但是目前针对高烈度区混凝土重力坝结构在斜入射地震波影响下的分析不多，且局限于传统的线弹性分析阶段，难以预知强震下坝体的真实破坏状态，在地震破坏动态分析时需从非线性的有限单元法进行突破。

目前，由Lubliner等[19]提出的混凝土塑性损伤模型（Concrete Damage Plasticity，CDP）是发展最成熟的非线性材料模型，国内外许多学者都基于此进行了研究与改进[20-23]。例如，沈怀至等[24]在弹塑性损伤力学的基础上，建立了一种坝体地震破坏评价模型并验证了它的合理性；范书立等[25]采用塑性损伤力学对混凝土重力坝进行动力分析，建立了包含能量特性的大坝整体损伤评价指标；李金友等[26]结合CDP模型，首次建立以DCR为评价指标的反应谱评价方法。还有许多学者对损伤本构和地震评价模型做出了研究与改进[27-30]，但目前还并没有广泛应用到斜入射等复杂地震输入机制下坝体的响应分析中。

针对混凝土重力坝在斜入射地震波下的响应研究存在的不足，且与非线性动力分析和安全评价方面未能有机结合的现状，本文以Konya重力坝为例，建立坝体-地基三维有限元动力模型，首先将一致激励地震输入的数值结果与实测破坏效果进行对比，验证混凝土塑性损伤模型的可靠性；再基于隐式有限元结合黏弹性人工边界的方法，将地震波动输入转化为作用于人工边界的等效节点力；结合混凝土塑性损伤模型，进一步分别考虑P波和SV波多角度斜入射下坝体的非线性动力响应，针对坝体关键点位移、应力与坝体损伤重点分析，并提出地震破坏评价模型，最后根据损伤指标及损伤等级对震后坝体破坏进行安全评价，为后续研究重力坝在斜入射下的非线性响应提供一定的参考。

2 总体结构

结合实际工程分析了地震波斜入射下混凝土重力坝的非线性动力响应，并对坝体进行安全评估，总体结构如图1所示。

图1 地震波斜入射下坝体动态响应总体结构

本文主要研究内容如下：（1）考虑坝体振动特性，建立三维非线性动力有限元模型，坝体部分采用CDP模型，并对模型参数进行验证。（2）采用黏弹性人工边界模拟地基辐射阻尼，并基于此推导P波和SV波斜入射下的波动输入等效荷载公式。（3）以坝体位移、应力和塑性损伤为研究对象分析其动力响应，并重点结合损伤区长度、位置及损伤因子建立地震破坏评价模型，综合评估重力坝震后破坏状态。

3 研究方法

3.1 基于黏弹性边界的地震波动输入方法

3.1.1 三维黏弹性边界 用有限元法分析结构-地基动力相互作用时，在建模中由于计算量等限制，通常只取部分近域地基，为保证外行波能量能逸散到远域地基，通常在切取的边界上建立全局或局部人工边界来模拟连续介质的辐射阻尼，保证散射波穿过边界时不发生反射。目前常采用刘晶波等[7]提出并改进的黏弹性人工边界，其考虑到了介质的弹性恢复能力，克服了黏性边界的低频失稳问题，具有良好的稳定性。其实质是在边界节点上施加连续分布的并联弹簧-阻尼器系统，弹簧的刚度系数和阻尼器的阻尼系数计算公式如下：

式中：KBT、KBN分别为弹簧切向与法向刚度系数；CBT、CBN为阻尼器切向和法向的阻尼参数；R为散射波源至人工边界点的距离；Cs、Cp分别为S波和P波的波速；G为介质剪切模量；ρ为介质密度；αT、αN分别为黏弹性人工边界切向、法向参数，本文经过反复对比各参数下结果拟合精度将αT和αN取为0.75和1。

3.1.2 P波入射时的等效荷载 地震激励下的复杂波场中总共包含入射波、反射波和散射波三种波，其中入射波场和反射波场称为自由波场，散射波能量通过人工边界来吸收。刘晶波等[7]提出了适用于黏弹性人工边界的波动输入方法，并将地震波动转换为作用于边界节点上的等效荷载，计算公式如下：

如图2所示，在有限区域内划定高度为h、长度为l的研究范围。以一斜入射的P波为例（左下角点（x0，y0，z0）为起振点），由于斜入射的P波在经过边界时会对其产生剪切和挤压作用，因此不但会产生反射的P波，还会产生反射的SV波。

参考周晨光[31]的波动输入推导方法，假定入射的P波与Y轴的正方向夹角为θ1，入射波与反射波确定的平面与XY平面夹角为α，反射的SV波与Y轴正方向夹角为θ2。由波动理论和单元应变状态，代入式（3），得到左面边界3个方向的等效节点力：

图2 P波斜入射半无限弹性空间波场

3.1.3 SV波入射时的等效荷载 如图3所示，一斜入射的SV波在经过边界时会产生一反射的SV波和一反射的P波。假定入射的SV波与Y轴的正方向夹角为θ2、入射波与反射波确定的平面与XY平面夹角为α、反射的P波与Y轴正方向夹角为θ1。

图3 SV波斜入射半无限弹性空间波场

等效荷载计算方法与P波类似，以底面荷载为例，推导3个方向的等效节点力：

式中：B1、B2和A2为入射SV波、反射SV波和反射P波的势函数幅值；Δt4、Δt5、Δt6分别为入射SV波、反射SV波和反射P波从起振点到边界上各点的延迟时间，推导过程与上节类似，只是注意P波与SV波波速的区别。

将式（4）—式（9）计算的等效节点力输入有限元模型中即可模拟出斜入射的P波和SV波。

3.2 混凝土塑性损伤模型混凝土重力坝的破坏通常由混凝土的损伤开始，进而发生开裂及裂缝扩展所致，CDP模型是由Lubliner等[19]提出，并由Lee等[32]改进和发展得来。该模型重点关注了材料拉、压性能的差异，用于模拟混凝土等准脆性材料在反复荷载下由损伤引起的材料退化，主要表现在抗拉压屈服强度的差异。混凝土单轴受拉与受压应力-应变及开裂应变关系如图4、图5所示。

图4 混凝土单轴受拉应力-应变及开裂应变关系

图5 混凝土单轴受压应力-应变及非弹性应变关系

如图4所示，在混凝土单轴受拉情况下，当应力未达到极限抗拉强度时，混凝土应力-应变为线性关系；当超过破坏应力后，进入软化阶段。图5为在单轴受压情况下，应力未达到σc0时，应力-应变也表现为线性关系，进入塑性状态后先是出现硬化，在达到极限抗压强度后进入软化阶段。两种情况下的应力应变曲线如下式[33]：

式中：E0为初始弹性模量；dt、dc分别为受拉损伤因子和受压损伤因子，取值在0（无损伤）到1（完全损伤）之间。

本文中混凝土重力坝的破坏主要体现在受拉损伤方面。此外，为合理考虑损伤开裂过程中，由于循环荷载作用下反向受压导致的裂缝闭合和弹性模量恢复现象（即“单边效应”），在设置参数时假定受拉刚度恢复参数ωt为0，受压刚度恢复参数ωc为1。

3.3 混凝土坝体地震破坏评价模型目前混凝土坝的抗震评价基于线弹性的应力标准，混凝土材料良好的抗压性能使得设计标准中对压应力的评价基本能够满足，但是强震作用下的拉伸损伤需要更好的定量评价方法。

沈怀至等[24]基于弹塑性损伤力学模型，由大坝损伤因子和损伤区沿坝体宽度方向的长度，建立了重力坝局部破坏评价模型。由于损伤破坏程度（由拉伸损伤因子做判别）与贯穿深度呈正相关性，所以本文对沈怀至的模型进行简化，以大坝拉伸损伤因子（dt≥0.8）及损伤区贯穿深度为依据，建立强震作用下大坝的局部损伤判断指标：

式中：T为权系数，损伤开裂区位于坝踵T取1.5，损伤开裂由下游面向上游发展T取1.0（贯穿上游面时取1.2），损伤开裂由上游面向下游发展T取1.2；li为局部破坏区域内的第i个水平截面上损伤因子大于等于0.8区域的长度；Li为此水平截面在坝体上的总长度。

根据Dm的大小将坝体损伤程度分为四级：

局部损伤判断指标概念清晰，便于结合有限元的结果进行计算与分析，其中dt≥0.8的情况下可以判断出此处已出现较大损伤破坏，应用简便且有一定理论依据。分级标准参考了陶能付等[34]提出的坝体震害等级分类，并符合其宏观描述。此外，许多学者也不同损伤指标对震害程度分级[35-36]，各有其优势。本文主要以混凝土重力坝的局部损伤指标为依据，实现快速合理的分级。

4 地震波斜入射下重力坝动力响应分析

4.1 坝体-地基-库水三维有限元模型Koyna混凝土重力坝位于印度西南部，1967年10月11日该地区发生了里氏6.5级地震，震中距13 km，震中位于大坝以南偏东2.4 km，重力坝坝长850 m，坝高103 m，地震发生时的坝前水位91.75 m。本文结合CDP模型对其真实震害进行模拟，并进一步分析了P波和SV波斜入射下坝体的损伤情况。取一个挡水坝段进行研究，图6为整个坝体-地基的三维有限元模型。

图6 Koyna坝三维有限元模型及关键点

坝体的尺寸为长103 m，坝基交界面处宽70 m，地基范围顺河向宽度取为40 m，顺河向自坝踵向上游侧和自坝趾向下游侧各延伸2倍坝高，竖直向自建基面向下延伸2倍坝高。坝体和地基部分均采用ABAQUS中的C3D8R单元来划分，整个模型共有5236个单元和7100个节点，其中地基边界采用弹簧和阻尼器单元来模拟黏弹性人工边界，坝体部分在Z轴方向进行法向约束。坝体及基岩各项主要力学参数如表1所示。此外，经过试算，确定合适的混凝土塑性参数为：膨胀角ψ=36.31°，流动势偏移值m=0.1，双轴极限抗拉强度和单轴极限抗压强度比fα=1.16，第二应力不变量在拉伸子午面上和压缩子午面上的比为γ=0.6667，黏性系数μ=0.0025。

表1 坝体与基岩主要力学参数

本文考虑的静荷载主要为坝体自重和静水压力，动荷载为斜入射地震波。P波和SV波分别从上游地基底部入射，P波入射角度分别为 0°（垂直）、15°、30°、45°、60°、75°、90°；SV波由于具有临界角度θcr（一般为35°），入射角度选为0°（垂直）、5°、10°、15°、20°、25°、30°，入射方式及关键点位置如图6所示。地震荷载统一选取Koyna纵向地震波记录，峰值加速度为0.312g（入射地震波在传播到地表时，要与由自由边界条件产生的反射波叠加，因此把在基岩中的入射地震波幅值取为设计地震动幅值之半[3]，文中取为0.157g），综合地震波采样频率和高阶模态振动周期的影响，选取时间步长为0.02 s，持续时间为10 s，对加速度曲线进行基线调整后积分，得到位移和速度的时程曲线，如图7所示。

图7 输入地震波的加速度、速度和位移时程曲线

4.2 模型参数验证为验证模型建立和参数取值的合理性，首先模拟真实地震条件下三维模型的损伤情况。各项材料参数和荷载与4.1节基本相同，除了动力荷载取为一致激励输入的Koyna坝的强震记录，考虑水平向和竖向的地震作用，持续时间为10 s，地震波从模型底部输入。

拉伸损伤结果输出如图8（a）所示，与图8（b）中振动台实验结果[37]作对比，破坏集中发生于折坡和坝踵处，并在折坡处沿着与下坡面呈90°的方向向坝内延伸约15 m，与试验值相近，证明文中模型和参数的选取可以较准确模拟坝体的损伤破坏。此外，如图8（c）和图8（d）所示，选用线弹性模型重复上述操作，发现时程结束时应力分布与CDP模型下的应力分布显著不同，且线弹性模型无法反映坝体的损伤累积，证明了CDP模型的优势。后续研究将从坝体关键点位移、应力和坝体塑性损伤三个方面分析P波和SV波在不同入射角度情况下重力坝的动力响应。

图8 Koyna重力坝地震响应

4.3 斜入射下坝体关键点位移图9（a）（b）分别为P波入射下坝顶A点顺河向和竖直向的位移时程曲线，其他各点变化趋势相同，只是幅值略有差异，此处不一一给出。从图9（a）（b）可以看出，顺河向最大动位移随入射角度的增加呈现先增大后减小趋势，极大值出现在入射角为60°时；竖直向最大位移随着角度的增大逐渐减小，极大值出现在入射角为0°时。图9（c）（d）为不同角度下坝体各关键点与坝踵处D点的最大相对位移。从图9（c）（d）可以看出，随着点高程增大，相对位移呈现增大的趋势；随入射角度变化规律复杂，但许多在60°时出现最大值。由此可初步看出，P波斜入射与垂直入射相比对坝体的动力响应影响较大，对坝体安全性有不利影响。

图10（a）（b）分别为SV波斜入射时坝顶A点的位移时程曲线，其他各点规律也基本相同。从图10（a）（b）可以看出，顺河向最大动位移随入射角度的增大而缓慢减小，极大值出现在入射角为0°时；竖直向最大位移随着角度的增大而增大，极大值出现在入射角为30°时。图10（c）（d）为不同角度下坝体各关键点与坝踵处的最大相对位移，图10（c）（d）可以看出，随入射角度的增大相对位移也随之增大，可以看出斜入射SV波对坝体的不利影响。

图9和图10中各点绝对和相对位移的变化规律体现了考虑入射角度的波动输入方法的必要性，同时刻相对位移差异越大，区域应力也越大，当应力超过极限抗拉强度时会产生塑性损伤破坏，这对后文中进一步探索损伤区提供了方向。

4.4 斜入射下坝体关键点应力P波斜入射时坝体各关键点第一、三主应力极值如图11所示。从图11可以看出，随着入射角度的增大，各点的主应力绝对值呈先增大后减小的趋势，基本在45°～75°这一区间达到最大值；地震动垂直入射时的最大主应力要略大于水平入射的情况，但均小于斜入射情况下的应力。其中折坡、坝踵和坝趾处为应力较大的区域，折坡处在15°～75°区间与坝踵处在60°～75°区间的第一主应力极值均有达到过极限抗拉强度（2.9 MPa），应该为后文中损伤区集中出现的区域和角度。

图9 P波不同入射角度下各测点位移曲线

图10 SV波不同入射角度下各测点位移曲线

SV波斜入射时坝体各关键点第一、三主应力极值如图12所示。从图12可以看出，随入射角度增大，各点第一主应力基本不变，第三主应力随入射角度增大呈现不断减小的趋势，可以初步判断垂直入射时坝体损伤最大。同P波斜入射时的情况类似，折坡处依然为拉应力最大的区域，且在各种角度下均达到极限抗拉强度，其他各点处均未达到，所以推测在SV波斜入射时损伤总发生在折坡处。

图11、图12中第三主应力都远大于第一主应力，这反映了混凝土材料的特性，抗压强度远大于抗拉强度；两种波型斜入射下，坝体始终未达到极限抗压强度，因此在后文中将以拉伸破坏变量为重点研究对象。针对折坡、坝踵等应力较大区域的损伤情况应重点分析，并进行安全评估。

图11 P波不同入射角度下关键点第一、三主应力极值

图12 SV波不同入射角度下关键点第一、三主应力极值

4.5 塑性损伤分析在有限元分析中，损伤因子及损伤区范围可以有效表示坝体在荷载作用下的刚度退化大小及程度，也可以展示大体积混凝土结构中裂缝萌生、扩展及失稳的全过程。图13、图14为不同角度的P波和SV波分别从上游面斜入射时坝体的损伤情况（用拉伸损伤因子来表示，图例中其变化范围为0～1）。

从图13可以看出，P波入射下坝体损伤破坏出现在入射角度为15°～75°时，在此区间损伤范围随着入射角度的增大呈先增大后减小的趋势，在折坡处损伤区贯穿深度也随之变化，在60°时达到了最大；坝体损伤破坏集中在折坡和坝踵两部分，这与上节中依据关键点应力判断出的结果吻合；在垂直入射（0°）和水平入射（90°）时，坝体保持弹性变形，均未出现损伤破坏，可见斜入射对坝体具有不利影响，特别是在强震作用下角度对坝体损伤的影响十分关键。在60°和75°两个角度下，压缩波对坝体的影响十分剧烈，在坝踵处也出现了少量破坏，且在60°时折坡处受拉应力和重力共同影响，沿水平向下弯曲方向发生了贯穿性损伤，这也是造成坝体头部左右摇摆、相对位移过大的主要因素之一。

从图14可以看出，在各个角度下，坝体折坡处均出现了水平向下弯曲的损伤破坏，损伤范围和损伤区贯穿深度随着角度的增大不断减小，这与剪切波在水平方向的分量减小有关，坝踵处均未出现破坏；垂直入射时坝体损伤最为严重，在折坡处下方约13 m处还发生了水平向裂缝；在小角度剪切波的作用下，坝体左右晃动剧烈，易在折坡处等应力集中区产生裂缝，这与实际情况相吻合。相比于P波，坝体损伤破坏对SV波入射角度变化敏感性更强，5°的角度变化可以产生很大的差距，因此在分析中入射波的角度应该重点考虑。

图14 SV波从左侧斜入射时不同入射角度下坝体损伤情况

由图13和图14的变化规律可知，采用不同波型和不同输入角度并结合CDP模型，才能更真实全面地反映坝体的破坏情况。根据图中所示的损伤开裂破坏长度，可由式（12）计算损伤指标，进行定量分析与评价。由于两种波型斜入射下坝体的破坏主要集中在折坡处，此区域在强地震作用下易产生近水平向的贯穿性损伤区。所以在损伤评价中以折坡处损伤为重点来分析坝体局部损伤的变化情况，结果见表2。

表2 地震破坏指标

由表2可见，坝体损伤指标直观地表示出坝体随地震波入射角度的损伤变化情况，将图13、图14中由损伤范围得到的规律进一步深入，重点关注裂缝的相对位置与扩展深度和方向，可以看出P波入射下指标随入射角度增加呈先增大后减小趋势，最大值达到0.60；SV波入射下指标不断减小，最大值为0.63。其中P波入射角度在60°时与SV波入射角度在0°、5°、10°时破坏指数都达到了0.60以上，说明这些情况下坝体损伤较为严重，应该重点关注。各情况中指标值可以对坝体损伤进行等级划分，验证了模型的合理性。

结合所有情况可知，坝体折坡处和坝踵是损伤集中发生的部位，应该针对此薄弱区域加强防护，例如采用高性能混凝土等，提高抗震设计。此外，在地震波斜入射情况下，本文中的地震破坏评价模型可以较好的对坝体破坏程度定量评价，可在其他类似工程中推广使用。

5 结论

本文通过建立混凝土重力坝坝体-地基三维非线性动力分析模型，分别分析了坝体在地震P波和SV波斜入射下的动力反应并进行了地震破坏评价，得出如下结论：（1）地震动斜入射下的动力反应与水平和垂直入射相比有显著不同，P波入射下位移应力和损伤在60°时达到最大响应，SV波入射下在0°时达到最大，证明了考虑入射角度的必要性；地震入射角度和不同地震波型的影响也应当综合考虑。（2）坝体附加应力和损伤破坏产生的主要原因是坝体各点振动的不协调性，并且显著受到关键点高程和坝体形状的影响；强震作用下应力较大的区域易达到混凝土抗拉强度值，对坝体造成极大的破坏。与传统的线弹性分析相比，CDP模型更能真实反映坝体在强震作用下的动力响应和破坏程度。（3）采用不同波型和不同输入角度的地震波并结合CDP模型定量计算损伤指标，能全面且直观地反应坝体破坏程度，在同类型的坝体中可以推广本文中的地震破坏评价模型，并基于此加强薄弱区域的处理与抗震设计。

文中对地震波斜入射下坝体非线性动力反应及损伤破坏进行了分析，但未同时考虑地基非线性，地基材料的复杂本构关系对斜入射波动输入时边界节点等效荷载的确定提出了新的考验，将在后续研究中结合非线性的分层地基和复杂地质条件下断层等情况深入讨论坝体-地基体系的动力响应。