基于耐震时程法的特高拱坝极限抗震能力研究

马天骁， 张燎军， 张汉云， 史俊飞， 翟亚飞

(河海大学 水利水电学院， 南京 210098)

我国西南地区修建了大量的300 m级特高拱坝，这些高坝枢纽工程，规模和数量在世界上都是史无前例的，而处在地震动活跃地区的特高拱坝的抗震安全尤其值得关注[1-2]。在对特高拱坝进行抗震性能评价时，采用常规的时程分析法计算常常面临着结构模型复杂、自由度数目较多、计算耗时长的问题。特别是当模型存在同时考虑坝体-基岩损伤的材料非线性和各类缝面的接触非线性时，一次分析往往需要耗费大量的计算时间，效率较低。对于特高拱坝的极限抗震能力分析，则需要多次大量的时程分析计算才能最后确定，计算成本更高。为了解决这一问题，需要采用一种既能考虑时程分析中的动力效应，又能简便高效的对特高拱坝抗震性能和失效模式进行评价的方法。Estekanchi 等[3-4]提出的耐震时程法，是一种基于时域的新型弹塑性动力分析方法。它的优点是可以通过一条加速度时程曲线来代表不同强度的地震系列，只需要计算一次耐震时程分析就能够得到不同地震动强度下结构的动力响应。自该方法提出以来，在钢框架及桥梁结构动力分析中得到了较多的应用[5-8]。在水工结构动力分析领域，徐强等[9-10]采用该方法对重力坝的动力损伤指标进行研究，Hariri-Ardebili等[11]采用耐震时程法分析了强震作用下拱坝损伤模式，并对高拱坝动力损伤进行安全评价。但对于300 m级特高拱坝，目前尚缺少采用耐震时程法全面考虑坝体及基岩的损伤、横缝接触非线性等方面的研究。

特高拱坝的极限抗震能力研究是地震时确保不发生库水失控下泄的关键课题之一[12]。很多学者都对高拱坝的极限抗震能力进行研究。如Cheng等[13]提出位移突变、塑性区贯通及计算不收敛等判据来分析高拱坝的极限抗震能力；李德玉等[14]同时采用地震动超载和基岩结构面降强等方法分析了白鹤滩拱坝的极限抗震能力；张社荣等[15]从位移响应，塑性区贯通，坝体开裂等角度对重力坝的极限抗震能力进行研究；李静等[16]采用地震动超载法，从坝体损伤及横缝开度等角度研究了高拱坝的抗震安全性能。以上研究大多是采用地震动超载法，采用多次时程分析计算来确定高拱坝的极限抗震能力。但由于高拱坝计算模型较为复杂，大量的时程分析计算比较耗时，目前尚少有采用耐震时程法来研究高拱坝极限抗震能力的研究，对采用耐震时程法和常规多次时程分析法计算极限抗震能力的结果也缺少对比研究。

基于上述原因，本文采用耐震时程法研究高拱坝的动力损伤行为及极限抗震能力。以中国西南某300 m级特高拱坝为例，建立了可以同时考虑坝体混凝土及坝基岩体的损伤、横缝的动态接触非线性的有限元分析模型。基于耐震时程法，采用非线性最小二乘函数合成耐震时程曲线。分别采用耐震时程法和常规的时程分析法，对比分析高拱坝位移、损伤体积比、损伤耗能、横缝最大开度等动力响应。结果表明，耐震时程法能够有效的反映不同地震动强度下高拱坝的损伤情况，其结果与常规时程分析法的计算结果较为接近，计算效率较高，在节省计算时间和计算资源方面具有优势。

1 基于耐震时程法的极限抗震能力分析方法

1.1 耐震时程法理论

耐震时程法(ETA法)是通过对特高拱坝计算模型输入一条地震动强度随时间逐渐增大的时程曲线，以分析坝体及基岩系统从线性阶段到非线性阶段的动力响应，直至最终破坏的全过程。该方法最主要的特点是随着地震动持续时间的增加，地震动强度逐渐增大，且在某一个时间段内的目标加速度反应谱与该时段的持续时间t呈线性关系

(1)

式中：tt为目标时间点；SaC(T)为事先确定的反应谱；T为结构的自振周期；t为任意时间点；SaT(T,t)为t时刻的目标加速度反应谱。从式(1)可知，在tt和SaT(T,t)确定的条件下，采用耐震时程法得到的加速度时程曲线中，从起始时间点到t时刻的加速度反应谱值与时间t成正比。因为位移反应与加速度反应存在相关性，故耐震时程法加速度时程的目标位移反应谱可以表示为

(2)

式中，SuT(T,t)为t时刻时的目标位移反应谱。在某一固定时刻t，同时满足式(1)和式(2)的加速度时程在一定的精度许可范围下是可以获得的；为了使在加速度时程上的任何点均满足式(1)和式(2)则需要求解满足式(3)的无约束变量的优化问题[17]

α[Su(T,t)-SuT(T,t)]2}dtdT

(3)

1.2 基于耐震时程法的极限抗震能力分析

在对特高拱坝进行非线性耐震时程分析后，为了方便对比分析常规时程分析法和耐震时程法的计算结果，需要将天然地震动的强度换算成耐震时程法对应的时间，地震动的等效耐震时程换算关系可以采用式(4)

(4)

同理，耐震时程曲线上各个时间点与天然地震动强度换算关系可采用式(5)

(5)

式中：tET为不同强度地震动的等效耐震时间；tTarget为耐震时程曲线总时间；S1为地震动的调幅参数；PGAC为规范中时程分析所用地震加速度时程的最大值。

通过式(4)和(5)可以从耐震时程分析计算结果中的到在各级地震动强度下坝体及基岩的动力响应。

损伤体积比(DVR)能够定量的反映大坝的损伤程度,损伤体积比采用下式计算[18]

(6)

式中：n为坝体单元数目；i为坝体单元编号；di(t)为i号单元在t时刻对应的损伤因子；Vi为i号单元的体积。损伤体积比的范围为0～1，损伤体积比越大，表示坝体损伤情况越严重。当损伤体积比为0时，表示坝体无损伤，当损伤体积比为1时，表示坝体所有单元已完全损伤。

高拱坝在地震损伤的过程中常常伴随着损伤耗能的增加，因此损伤耗能指标也可以在一定程度上反映高拱坝动力损伤情况，损伤耗能采用下式计算

(7)

式中：MD为坝体损伤耗能；σ为坝体单元应力；εck为开裂应变；V为单元体积。

本文综合考虑在各级地震强度下拱冠梁的最大相对位移、横缝最大开度、坝体损伤程度、损伤体积比、损伤耗能等因素，进而确定大坝的极限抗震能力。

2 特高拱坝有限元模型及耐震时程曲线的生成

2.1 工程概况及有限元模型

以中国西南某300 m级特高拱坝为例，建立能够同时考虑高拱坝—基岩相互作用的非线性有限元分析模型(图1)。坝址区地震基本烈度为Ⅷ度，特征周期Tg=0.35 s。该特高拱坝最大坝高285.5 m，拱冠梁坝底处厚60 m，坝顶处厚14 m。坝顶中心线弦长602.2 m，弦高比为2.17。正常蓄水位时水深为275.5 m，相应的下游水深为53.5 m。坝体共设置横缝30条。模型共划分单元31 453个，自由度数为115 149，横缝30条。由于特高拱坝计算模型规模较大，全部采用精细网格划分会大量提高计算时间，但网格过大会造成计算精度下降。为了平衡这个问题，本文在坝体和近坝部分基岩采用精细网格划分，远端基岩采用较为稀疏的网格划分。

图1 特高拱坝非线性有限元模型Fig.1 Nonlinear finite element model of ultra-high arch dam

2.2 材料参数及荷载

坝体混凝土及基岩具体材料参数如表1所示，按照规范规定，混凝土的动态弹性模量较静态弹性模量提高50%，强度参数较静态提高20%。基岩材料弹性模量不提高，其损伤本构关系，目前尚缺乏充分研究。本文采用文献[19]提供的方法，将基岩的材料参数折算为混凝土的强度参数，采用下式计算

(8)

表1 坝体及基岩材料参数表Tab.1 Material parameter of dam body and bedrock

式中：ft为折算后基岩的抗拉强度；c和φ分别为地基岩体的黏聚力和内摩擦角。本工程坝基岩体c取2.5 MPa，φ取53.47°，计算得到折算后的抗拉强度ft为1.65 MPa。

本文采用Abaqus中的混凝土塑性损伤模型(CDP模型)来模拟坝体混凝土及基岩的损伤行为。混凝土受拉应力应变关系曲线是基于GB50010—2010《混凝土结构设计规范》推荐的混凝土应力应变关系得出的。损伤因子的计算是基于能量损失理论，将脆性固体材料的损伤因子D定义为[20]

(9)

式中：W0为无损材料的应变能密度；Wε为损伤材料的应变能密度；E0和E分别为无损材料和损伤材料的四阶弹性系数张量；ε为相应的二阶应变张量。

混凝土及基岩的受拉应力应变关系曲线如图2所示，相对应的损伤因子与应变关系曲线如图3所示。

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

图3 混凝土及基岩受拉损伤因子与应变关系曲线Fig.3 Relationship between tensile damage factor and strain of concrete and bedrock

施加的静力荷载主要有坝体自重、上下游的静水压力，淤沙压力以及温度荷载。正常蓄水位工况下，上游水位高程为600 m，淤沙高程为490 m。

对于高拱坝而言，采用修正的Westergaard法模拟坝体与库水的动力相互作用[21]

(10)

式中：p为水深h处的动水压力；a为水深h处的加速度；H0为库水深度；h为计算点的水深；ρ为水体密度。

坝体采用阻尼比为5%的瑞利阻尼计算。坝基四周设置黏弹性边界来模拟地基辐射阻尼，通过在边界施加弹簧和阻尼单元来实现。采用波动输入法，按文献[22]提出的计算公式，通过编制将底边界的入射波场和侧边界的自由波场转化为对应的解析应力场的程序，对Abaqus进行二次开发，通过在边界施加等效节点力的形式实现地震动输入。

2.3 耐震时程曲线的生成

基于2.1节的式(3)，本文采用NB 35047—2015《水电工程水工建筑物抗震设计规范》[23]中的标准反应谱为目标谱，取特征周期Tg=0.35 s，阻尼比为5%，反应谱最大代表值βmax=2.50。参照规范第5.5.8条的要求：“采用时程分析法计算地震作用效应时，应以5%的设计反应谱为目标谱，生成至少3套人工模拟地震加速度时程作为基岩的输入地震动加速度时程”。为了避免偶然性，本文采用非线性最小二乘法生成3条(ETA-1～ETA-3)时长为40 s的耐震时程曲线。如图4(a)为第一条耐震时程曲线ETA-1的加速度时程曲线，从图中可以看出，随着时间的增加，地震动强度逐渐增大。图4(b)给出了40 s的耐震时程曲线在不同时间段的反应谱，分别为：0～10 s，0～20 s，0～30 s，0～40 s。可以看出4个时间段的反应谱能与目标谱能够很好的吻合，类似的，图5(a)、(b)以及图6(a)、(b)分别表示第二条(ETA-2)和第三条(ETA-3)耐震时程曲线及对应的反应谱。经过计算，三套耐震时程曲线的平均误差值在20%以内，拟合误差较小。

(a) 耐震时程曲线(ETA-1)

(b) 加速度反应谱(ETA-1)图4 耐震时程曲线及对应的加速度反应谱(ETA-1)Fig.4 Endurance time accelerogram and corresponding acceleration response spectrum(ETA-1)

(a) 耐震时程曲线(ETA-2)

(b) 加速度反应谱(ETA-2)图5 耐震时程曲线及对应的加速度反应谱(ETA-2)Fig.5 Endurance time accelerogram and corresponding acceleration response spectrum(ETA-2)

(a) 耐震时程曲线(ETA-3)

(b) 加速度反应谱(ETA-3)图6 耐震时程曲线及对应的加速度反应谱(ETA-3)Fig.6 Endurance time accelerogram and corresponding acceleration response spectrum(ETA-3)

3 计算结果分析

本小节同时采用耐震时程法和常规时程分析法对高拱坝进行非线性动力计算。首先，基于标准反应谱生成3组人工地震波(GM-1～GM-3)，以峰值加速度(PGA)作为强度指标，将每条地震动在0.1g～1.0g之间进行调幅计算，间隔为0.1g，共进行30次非线性动力分析。其次，采用2.3节生成的三条耐震时程曲线，进行三次耐震时程分析。为了全面反映高拱坝在地震过程中的动力响应，本小节采用坝体最大相对位移、坝体动态损伤分布、坝体损伤耗能、横缝最大开度等指标来对高拱坝的极限抗震能力进行评价。

3.1 坝体最大相对位移及横缝开度分析

拱冠梁最大位移和横缝最大开度是评价特高拱坝地震易损性的重要动力响应参数。为了方便分析对比耐震时程法和常规时程分析法的计算结果，对两种方法得到的最大相对位移及横缝最大开度计算结果进行汇总，如图7和8所示。

图7 拱冠梁最大相对位移Fig.7 Maximum relative displacement of crown cantilever

图8 横缝最大开度Fig.8 Maximum opening of contraction joint

从图7结果可以看出，由于选取的是从起始时刻至该时刻间最大的动力响应。因此得到的数值曲线随时间呈逐渐增大的趋势。从整体上来看，耐震时程法计算结果与常规时程分析法结果较为接近，三条ETA曲线的位移值集中分布在常规时程分析法结果的均值附近，结果吻合较好。

从图8的结果可以看出，对于坝体横缝最大开度而言，一般而言超过50 mm就很可能会造成止水撕裂，威胁大坝的安全。采用ETA法计算得到的坝体横缝最大开度整体上处于时程分析法的均值水平，耐震时程法和常规时程分析法的计算结果具有很好的一致性。当ETA时间在0～30 s时，对应最大PGA约为0.75g，三条ETA曲线值均小于50 mm，说明在罕遇地震作用下，横缝不会发生较大的张开而导致止水被破坏。

3.2 坝体—基岩损伤分析

为了能够对大坝的损伤情况进行定量的分析，本小节通过采用损伤体积比(DVR)和损伤耗能这两个指标来反映大坝的损伤程度。绘制它们随时间变化曲线。从图9、10中可以看出，两种方法得到的坝体损伤体积比很接近，结果离散性较小。从损伤耗能曲线上来看(图10)，也可以得到相似的规律。

图9 坝体受拉损伤体积比Fig.9 Tensile damage volume ratio of dam

图10 坝体损伤耗散能Fig.10 Damage dissipation energy of dam

由于三条耐震时程计算得到的坝体损伤分布规律非常相近，本小节以其中一条耐震时程曲线的损伤计算结果为例，分析特高拱坝的动力损伤规律。

当ETA曲线处在0～10 s时，对应的最大峰值加速度约为0.25g。从图11(a)～(c)可以发现，此时坝体主要处于线弹性状态，在坝体上下游面及拱冠梁均未出现明显的损伤区域。由于考虑了基岩的损伤，使得坝基附近坝体的应力得到释放，因此在拱冠梁的坝踵坝趾处未出现受拉损伤。

(a) 拱冠梁(t=10 s)

(b) 下游坝面(t=10 s)

(c) 上游坝面(t=10 s)

(d) 拱冠梁(t=20 s)

(e) 下游坝面(t=20 s)

(f) 上游坝面(t=20 s)

(g) 拱冠梁(t=30 s)

(h) 下游坝面(t=30 s)

(i) 上游坝面(t=30 s)

(j) 拱冠梁(t=40 s)

(k) 下游坝面(t=40 s)

(l) 上游坝面(t=40 s)图11 各时刻坝体受拉损伤分布图Fig.11 Tensile damage distribution of dam at each time

当ETA曲线处在0～20 s时，对应的最大峰值加速度约为0.5g。此时坝体下游面首先出现损伤开裂，从拱冠梁的损伤分布图中可以看出，损伤区域主要集中出现在坝体的中上部位，并从下游面逐渐向上游面扩展。坝体上游面整体上未出现明显的损伤区域，两侧岸坡坝段的坝踵处开始出现轻微的损伤，河床坝段的坝踵处完好，未出现明显的损伤区域，但坝基损伤区域进一步扩大，并向横向发展。

当ETA曲线处在0～30 s时，对应的最大峰值加速度约为0.75g。此时的地震动强度已经达到了罕遇地震以上，从损伤分布图来看，坝体下游面的损伤区域进一步扩大，特别是中上部区域损伤较为严重。从拱冠梁及上游面的损伤分布来看，此时损伤已经从下游贯穿到上游，并在上游面的中上部位形成了明显的损伤区域。两岸岸坡坝段坝踵处的损伤区域进一步扩大，河床坝段坝踵处未出现明显的损伤。但坝基出现了较大范围的损伤，可能会危及帷幕的安全，会给坝体的稳定性造成安全隐患。

当ETA曲线处在0～40 s时，对应的最大峰值加速度约为1.0g。从图中可以看出，此时大部分坝体均处于严重损伤状态，损伤范围已经贯穿坝体，但损伤未必一定会造成坝体开裂，也有学者提出损伤因子超过0.7时可以认为混凝土发生宏观裂缝，该判据也可作为判断混凝土开裂范围的参考。考虑到此时大坝的损伤程度，可以判定大坝已失去承载能力。

需要说明的是，单元尺寸的大小对于损伤结果具有较大的影响。本研究同时采用耐震时程法和常规方法来对比分析特高拱坝的极限抗震能力，由于两种方法采用的是同一套网格，因此网格尺寸不影响两者的对比分析。

3.3 大坝极限抗震能力分析

为了全面的评价大坝极限抗震能力，需要综合考虑坝体的最大位移、横缝的最大开度、坝体及基岩的损伤程度等因素。从3.1～3.2节分析可知，采用耐震时程法得到的坝体—基岩动力响应与采用常规的时程分析法得到的结果非常相近。耐震时程法能够全面的反映出大坝和基岩从完好到逐渐破坏直至最终失效的过程。在耐震时程分析的前10 s时间段，仅基岩出现轻微损伤，坝体完好，此时大坝处于安全状态。当处在第20 s时刻时，尽管坝体的局部区域发生了损伤，但仍未贯穿坝体，横缝最大开度仍达不到破坏止水的程度，此时大坝仍能够承载。当处在30 s的时刻时，坝体的损伤区域进一步扩大，并且贯穿了拱冠梁的中上部，基岩损伤区域进一步横向发展，可能会破坏帷幕进而威胁大坝的稳定性。此时大坝承受的最大峰值加速度PGA为0.75g，综合上面的分析，可以判定此时大坝处于极限状态。

4 结 论

本文基于耐震时程法来研究特高拱坝—坝基的动力损伤行为，并进行极限抗震能力分析。计算中综合考虑了坝体混凝土及坝基岩体的损伤等材料非线性问题以及横缝开合等接触非线性问题，得到以下结论：

(1) 耐震时程法具有较高的计算效率，计算得到的坝体位移、损伤、横缝最大开度整体上与常规时程分析法接近，且计算结果离散性小。

(2) 耐震时程法能够全面的反映出大坝和基岩从完好到逐渐破坏直至最终失效的过程。通过一次耐震时程计算，即可大致判断出大坝的极限抗震能力，可以大大简化极限抗震分析的计算工作量。

(3) 采用耐震时程法对某300 m级特高拱坝进行极限抗震能力分析，综合位移、横缝开度、坝体及基岩的损伤情况等方面，确定其极限抗震能力为0.75g，为相关工程的极限抗震能力分析提供一种新思路。