强震下高土石坝不同地震动输入特性评价

冯燕明，王 琪，朱 晟

(1.中国电建集团昆明勘测设计研究院有限公司，云南 昆明 650032；2.水电水利规划设计总院，北京 100120；3.河海大学水利水电学院，江苏 南京 210024)

0 引 言

现阶段的土石坝地震反应分析，大多采用单点输入法，假定到达地表面各点的地震动时程是相同的。但实际的地震波受到不同场地条件、路径条件等因素的影响，到达地表面上各点地震动都是不相同的；同时，限于目前人类对地震现象的认识水平和强震观测的技术条件，难于对地震的发生和地震波传播作出准确的预报，因此对于大跨度的土石坝有必要通过人造地震动时程来模拟真实地震过程，对不同地震动输入特性进行评价。Dibaj等[1]进行了土石坝在行进波地震作用下的动力分析，分析了一个理想均质的弹性土坝，计算中假定地震波是水平方向传播的剪切波，得到行进波计算的结果与按照刚性基岩的结果有较大差别，且偏于不安全；沈珠江等[2]深入研究了行进波在土石坝中的动力反应，认为土石坝不考虑行进波的影响偏于不安全；Hao等[3]考虑了相干效应和波的传播性，提出了基于相干函数模拟空间变化的地面运动方法；屈铁军和王前信[4-5]在Hao研究的基础上，提出了一套较完整的空间相关多点地震动合成方法，求解地震动沿管线的反应时程，为确定延伸型结构物各激励点的地震动时程提供了方法。

2008年5月12日，四川省汶川县境内发生了里氏8.0级地震，距震中仅17 km的紫坪铺混凝土面板堆石坝经受了地震烈度Ⅹ度的浅源近震考验，并且地震动输入模式有别于常规模式，因此本文根据文献[4-5]，并考虑到地震动不同方向的相关性，合成多点多向地震动，基于实测波功率谱合成多点地震动对紫坪铺大坝进行多点输入动力有限元计算，并与汶川地震大坝坝址实测波输入进行计算比较，对地震动不同输入方式进行评价。

1 单点及多点地震动合成理论

1.1 单点地震动合成

采用Scalan和Sachs提出的三角级数法[6]来合成单点地震动。非平稳的地震动加速度时程可表示为如下目标函数

(1)

式中，f(t)为强度包络函数；Ck、ωk分别为第k个频率分量的幅值和频率；φk为(0，2π)区间内均匀分布的随机相位角；t为时间。

根据公式(1)，可以构造一个平稳高斯过程

(2)

通过给定的功率谱密度函数可以求得公式(2)中三角级数各分量幅值

Ck=[4S(ωk)Δω]1/2

(3)

式中，Δω=2/T，ωk=2πk/T，T为平稳随机过程a(t)的总持续时间；S(ωk)为功率谱密度函数。

式(3)的过程是具有功率谱密度函数S(ω)的零均值平稳高斯过程。如果目标函数模拟的是加速度反应谱，因此需建立功率谱与反应谱的转化关系。采用Kaul[7]提出的功率谱与反应谱转化的经验公式，即

(4)

文中利用实测波功率谱合成地震动时，直接通过实测波时程求出地震波功率谱，并代入到公式(3)中计算。

三角级数的求和采用快速傅里叶变换(FFT)。文中由于模拟“5·12”汶川地震地震动，地震平稳段持续时间相当长，采用单峰值包络函数不能够很好地说明地震动强度随时间的变化，因此，采用Amin和Ang[8]提出三段包络线的多峰状强度包络函数，即

(5)

式中，c为衰减系数；t1、t2分别为平稳段的开始时间和结束时间。

霍俊荣等[9]考虑到包络函数受到多种因素的影响，如场地条件、震中距、震级、强震持时等，提出了公式(6)中地震动时程强度包络参数确定公式，即

logY=C1+C2M+C3log(R+R0)+ε

(6)

式中，Y为地震动参数值；M、R分别为震级、震中距；C1、C2、C3、R0、ε为包络参数衰减函数的回归系数，霍俊荣等给出了不同条件下的回归系数值。

根据ts=t2-t1，可确定地震动峰值平稳段结束时刻t2，峰值衰减系数c按振动结束时振幅为峰值的5%确定，地震加速度总持时td可以根据式(7)来确定。

(7)

由于初始时程得到反应谱只是用一些点去拟合设计反应谱，吻合程度上是概率平均的，可以按式(8)的方法通过迭代调整以提高精度，即

(8)

式中，FS(ωk)为幅值谱，通过迭代可使计算反应谱Sa(ωk)向目标反应谱Sa(ω)逼近，从而满足精度要求。

1.2 多点非平稳地震动合成

对土石坝不同地震动输入方式评价时需要合成多点地震动。基于功率谱合成多点地震动，需先生成功率谱矩阵，对于频率为ωk的地震动分量，有

(9)

式中，功率谱矩阵中的对角线元素为自功率谱，非对角线上元素为互功率谱。为了使人工合成的多点地震动与实测波具有可比性，自功率谱可由实测波转换得到。实测波地震动时程首先作傅里叶变换，得到傅里叶幅值谱，由公式(10)可得到自功率谱

(10)

考虑相干效应和行波效应的任意两点j、k的互功率谱可表示为

(11)

式中，Sj(ω)、Sk(ω)为自功率谱；ρjk(ω，djk)为相干函数；va(w)为视波速；djk为j、k两点沿地震传播方向的距离；θ为传播方向与j、k两点连线的夹角。

公式(11)中相干函数参考冯启民等[10]根据1975年中国海城地震余震的观测记录分析提出的相干函数经验模型

ρij(dij，w)=e-(ρ1w+ρ2)dij

(12)

式中，ρ1、ρ2为相干性参数，本文参照海城地震参数取值。

式(11)中视波速采用屈铁军等[11]对加速度进行带通滤波得到的经验公式

(13)

式中，a，b为拟合参数。取屈铁军等根据SMART-1台阵波速拟合结果的平均值，a=3 344，b=1 095。

通过式(10)、(11)得到自功率谱及互功率谱后，对功率谱矩阵式(9)三角分解，最终可求得幅值anm和相位角θnm

(14)

(15)

式中，lnm(ωk)为功率谱矩阵式(9)三角分解后下三角矩阵第n行第m列元素。

最终合成的地震动加速度表达式可表示为

(16)

上述多点平稳随机过程乘强度包络函数式可得到非平稳随机过程。

2 地震反应分析方法

2.1 筑坝材料的动力本构模型

室内大型动三轴试验资料表明：在复杂的高应力条件下，试验粗粒料的动应力-应变关系具有硬化特性，其阻尼比Hardin假定值小，采用基于指数型动应力-应变关系模型的动剪模量与相应的阻尼比计算公式[12]。

2.2 永久变形模型

朱晟等[12]考虑初始固结围压的影响，将残余体积应变和残余剪应变表示为振次、动剪应力、应力水平以及初始固结围压的函数得到永久变形计算模型。

2.3 多点地震动输入动力平衡方程

对于平面尺寸较大的建筑物，由于地震波在结构基础面上的传播要经历一定的时间，这样，在同一时刻，结构各支承点所承受的地面运动时不同的。在这种情况下，必须考虑各支承点间相对运动所引起的结构内的准静力位移。多点地震动输入的动力方程可以表示为

(17)

本文采用Dibaj和Penzien[1]建议的方法，将绝对位移分解成准静力位移和动力位移，即

{ut(t)}={us(t)}+{u(t)}

(18)

准静力位移可以由下式确定

(19)

式中，[γ]为准静力影响矩阵。

将式(18)、(19)代入(17)，并考虑采用集中质量矩阵，可得到动力位移

(20)

3 不同地震动输入特性评价

3.1 工程简介

紫坪铺混凝土面板堆石坝位于中国四川省成都市西北60余公里的岷江上游，坝体为钢筋混凝土面板堆石坝，最大坝高156 m，坝体地震设计烈度Ⅷ度，地震具有震级大、震源浅(约14 km)、断层长(接近300 km)、持续时间长等特点[13]，地震对附近大坝造成一定影响[14-15]。2008年5月12日汶川发生里氏8级地震，大坝经受了远高于其设计水平的X度浅源近震考验，紫坪铺大坝距汶川地震震中17 km，距发震断层地表破裂带约为8 km，是距地震震中最近且工程规模最大的堆石坝工程。考虑到没有坝址基岩的实测加速度记录，选择距离最近的茂县地办地震台(051MXT)测得的基岩加速度时程(数据由国家强震动台网中心提供)作为参照对象，参考于海英等[16]给出的衰减关系，考虑上下盘效应，推求坝址基岩峰值加速度，得到NS、UD方向加速度峰值分别为0.46g、0.43g；然后根据基岩加速度记录采用比例法推求坝址基岩输入加速度曲线，如图1所示。大坝有限元计算网格剖分及材料分区如图2所示，其中坝体结点和单元数分别为427个和398个。大坝的静、动力计算参数见文献[17]。

图1 坝址基岩输入地震加速度时程

图2 大坝有限元计算网格及材料分区

3.2 地震动单点一致输入与多点输入评价

3.2.1 多点地震动合成

选取大坝与坝基接触的26个约束点进行多点地震动输入，合成各点的地震动峰值如图3所示。从图3可知，水平地震加速度各点峰值在4.6 m/s2上下波动，竖直向各约束点地震动峰值在4.3 m/s2上下波动。

图3 各约束点合成的峰值加速度

3.2.2 不同地震动输入评价

利用3.2.1中基于实测波功率谱合成地震动进行多点输入动力有限元计算，并与实测波单点输入进行比较，结合典型结点加速度反应、大坝绝对加速度极值反应及坝体变形对2种输入方式进行评价。

图4为地震动多点输入、实测波单点输入下坝顶典型结点地震加速度反应的傅里叶谱与反应谱曲线。从图4可知：①无论水平顺河向还是竖直向，两种输入方式下坝顶结点加速度反应的Fourier谱幅值主要集中在0.5～1.0 s范围内，坝顶加速度反应周期明显延长；②两种输入方式坝顶结点加速度反应Fourier谱幅值在大坝前两阶振型自振周期附近水平顺河向和铅直向分量得到显著放大，说明输入地震波在接近大坝自振周期位置的分量由于较强的共振激励作用产生了显著放大；③地震动多点输入、实测波单点输入计算的水平向地震加速度频谱特性与实测值比较接近，而竖直向地震动单点输入反应幅值比多点输入幅值大，多点输入计算值更接近实测值。

图4 坝顶典型测点地震反应

地震动多点输入、实测波单点输入动力反应极值的比较如图5所示。从图5可知，水平向加速度极值分布趋势是一致的，都是由坝基向坝顶方向逐渐增大，最大位置出现在坝顶，并且量值在坝轴线上游均一致，在下游侧实测波单点输入量值较多点输入大。竖直方向加速度反应极值实测波单点输入比多点输入大，多点输入极值的减小主要是由于各点的相干性和相位的不同步性引起的。多点输入时，加速度反应极值大致为对称分布，这是由于地震波传播时地震波速度有限，使得到达各点的时间不同，各约束点的相位不同步，坝体反应在极值分布就会为大致对称分布。

动剪应力极值等值线对比如图6所示。从图6可以看出，不同输入方式下大坝动剪应力极值分布规律一致，自坝基向坝顶逐渐减小，实测波单点输

入下的动剪应力极值比地震动多点输入大，分别为0.42、0.38 MPa；多点输入与实测波单点输入极值出现位置相同，均在坝轴线偏上游处，多点输入的最大动剪应力极值较小是由于各点之间的相干性。

实测波单点输入及多点地震动输入永久变形计算结果比较如图7所示。从图7可以看出，沉降量随着坝体高程的增加逐渐增大，两种输入分布规律一致。实测波单点输入计算得到850.0 m高程坝轴线处沉降量为102.0 cm，地震动多点输入计算沉降量为85.0 cm，多点输入总体计算量值比实测波单点输入偏小。计算得到紫坪铺大坝震后永久变形较大，主要原因为，一方面“5·12”汶川地震持时与强度均较大，等效振次达到64次，大于seed建议的30次；另一方面震中距较小，坝基竖向地震加速度极值与水平方向较为接近，使得大坝竖直方向的地震惯性力较大，坝料可能发生破碎和颗粒重组产生了塑性体积变形。

图7 大坝竖向永久变形

由于地震惯性力为瞬时作用的循环荷载，即使坝坡潜在滑块的抗滑安全系数在短时间内小于1.0，引起塑性滑移，但是当加速度减小或反向时，这种位移的趋势又将停滞，这样一系列数值大、时间短的惯性力的全面影响的结果将是坝坡的累积滑移，因此可以通过计算土石坝坝坡在地震中发生的永久滑移量，来评价抗震稳定性。

利用地震过程中静、动应力叠加结果，分别采用实测波单点输入和多点输入计算坝坡滑块的抗滑安全系数。潜在危险滑弧条数实测波单点输入、多点输入分别为130条和67条，取安全系数小于1.0滑弧条数出现频次最多的滑块作为潜在最危险滑块，在地震过程中，计算所得到的下游坝坡潜在最危险滑块逸出深度约24 m。针对潜在最危险滑块，利用NEWMARK方法计算下游坝坡的累计滑移量，结果如图8所示。从图8可知，实测波单点输入和多点输入最大滑移量分别为15.7、11.7 cm。为了与多点合成人造波具有相同的持时，实测波计算时选取地震前60 s过程，通过地震加速度时程曲线可以看到在50 s时出现了强度较大的震动，故累积滑移量在50 s时突然增大。

图8 下游坝坡最危险滑块累积滑移量

4 结 论

结合汶川地震基于实测波功率谱合成地震动对紫坪铺大坝进行多点输入动力有限元计算。与地震动单点输入相比，多点输入时地震反应有明显减小，这主要由于多点输入时考虑了地震波的相干效应和行波效应。对不同的地震动输入进行频谱特性分析，基岩输入地震波的高频、短周期分量基本上被坝体滤波，坝顶加速度反应周期明显延长，输入地震波在接近大坝自振周期位置的分量由于较强的共振激励作用产生了显著放大。实测波单点输入、多点输入计算大坝竖向永久变形均位于坝顶，在量值上多点输入较小；最大动剪应力出现在坝基面坝轴线附近，多点输入较单点一致输入略小；大坝下游坝坡表面潜在危险滑块出现的位置较高，多点输入与实测波计算得到典型滑块在地震过程中的最大滑移量一致。