为了充分开发我国水能资源,西南地区的高坝水电枢纽数量不断增加,而且高坝数量不断攀升,规模史无前例(周建平等,2007).由于该地区为高烈度地震区,地震活动强烈,因此高坝工程的抗震安全性受到了极大的关注.国内学者围绕高坝抗震分析与安全评价问题开展了大量的研究工作(罗奇峰等,2003;徐艳杰等,2012;姚振兴和郑天愉,1985),但是关于高坝坝址区地震动问题的研究却相对比较滞后.地震动参数(峰值加速度、反应谱和地震历时)的确定已经成为决定高坝枢纽抗震安全评价准确性的关键课题.此外,由于大型高坝枢纽通常处于高山峡谷地区,地震动随着空间位置的变化, 亦十分显著,因此坝址区地震动动参数的三维分布特性也需要特别关注.
坝址区地震动的形成经历了震源断层破裂、地震波传播和局部场地反应等物理过程,是一个非常复杂的自然现象.对于重大高坝工程,需要通过专门的地震危险性分析确定地震动参数.地震资料丰富的地区可以采用大量实测记录回归分析得到衰减模型,如美国西部的新一代衰减关系(Next Generation Attenuation) (Power et al., 2008),但对于我国西南部等地震资料相对缺乏的地区只能借用其他地区的模型(胡聿贤等,1996).在传统的地震危险性分析中,基于经验衰减关系来估计坝址区的地震动参数,忽略了地震波产生和传播的物理机制.由于不同坝址区地质构造和地形条件的特殊性,统一套用其他地区获得的经验性规律很可能无法正确地描述该坝址区的地震动特征.而且,衰减关系反映的是某一地区的地震宏观规律,不能代表单一地震事件的特性,无法考虑每次地震由于发震机制、传播过程和局部场地效应差异而导致的场址地震动特异性.
近年来,将基于地震学理论的地震动模拟方法引入工程场址地震危险性分析中,考虑场址地震动生成的物理机制受到了研究者们的关注(Douglas and Aochi, 2008; Graves et al., 2011; Cui et al., 2013; Baker et al., 2014),正在成为发展趋势.事实上在1994年北岭地震发生后,Heaton以及Hall等就开始融合地震学和工程学,尝试“断层破裂-结构响应”一体化的计算思路(Heaton et al., 1995; Hall et al., 1995; Hall, 1998).随着数值方法和计算机技术的进步,采用谱元法(Komatitsch and Vilotte, 1998; Komatitsch et al., 1999)对地震过程进行大规模计算并为工程结构提供输入条件成为发展方向,Krishnan等(2006)首先基于谱元法完成了从断层破裂直到地震波作用于结构的“end-to-end simulation”.需要指出的是Heaton, Hall和Krishnan等研究的均为高层框架结构,与之相比较,我们所关注的大坝结构具有以下两个明显的特点:(1)参与振动的振型多,频率范围在0~10 Hz甚至更宽;(2)空间跨度大且位于高山峡谷,地震动的非均匀分布问题十分突出.He等(2015)基于这两点,利用高精度的震源和传播介质模型,对美国南加州Pacoima坝址区进行了模拟,获得了较好的结果.
然而目前我国西南地区,地震资料相对匮乏,各个大型水利工程都按照传统地震危险性分析进行设计.本文则针对现行地震安全评估方法中的一些不足,结合各个工程自身需要,通过建立起包含发震断层,波速结构,起伏地形的谱元法模型,应用超级计算机大规模模拟坝址区波场,给出坝址区地震动特性及其三维分布规律,与传统方法互为补充,进一步完善高坝抗震分析与安全评价的地震荷载输入环节.
2 重大工程场址地震动参数的确定 2.1 传统的地震危险性分析地震危险性分析方法分为概率分析法和确定性分析法.前者是美国学者Cornell (1968)提出的一种分析方法,认为某一特定区域内未来将要发生的地震活动的时空分布、强度以及所产生的地震动,都是随机的,最终以概率形式给出地震动参数(PGA,反应谱等).概率分析法的优势在于能够描述场地受地震活动影响的随机性,可以对工程场址在使用年限内遭受不同程度地震的概率进行定量估计,符合工程设计的理念.因此自提出以来,在地震工程中得到了广泛的应用,是工程场址地震动安全评估的主要方法.然而概率分析方法是对区域内所有地震影响的包络,得到的坝址区地震动不是某一次具体地震环境下的真实样本实现,失去了震源位置,震中距,震级等基本地震物理概念,不能代表任一可能发生地震所产生的地震动.
确定性分析法是在概率分析出现之前被普遍采用的一种地震危险性分析方法,不考虑地震发生的时空概率.在工程应用上,专门针对某一场或几场特定地震,研究特定场地的地震动参数(US Army Corps of Engineers, 1995).确定性分析方法,可以认为是概率分析法的一个特例.确定性分析方法不含有地震发生概率的概念而是针对最不利的情况进行设计,符合地震灾害是由某一个具有明确物理特性参数的地震事件所引发的事实.对于一些重大工程,我们不仅关心在其使用年限内遭遇某一地震烈度的超越概率,同时也需要考虑最极端情况的地震,应对任何可能的不测事件,因此确定性分析方法对于大型关键工程地震风险评估同样是必不可少(US Army Corps of Engineers, 1995; Krinitzsky, 1995).
无论是概率分析法还是确定性分析法,都是依靠地震动衰减的经验关系式建立震源和场地地震动之间的关系,例如基岩水平峰值加速度和反应谱衰减关系的一种表示方式为(霍俊荣,1989)
(1) |
式中,Y代表地震动参数,M表示震级,R表示震中距,ε是均值为零的随机误差项,ci(i=1, 2, …, 6)为回归系数.这种基于大量实测记录的经验估计方法比较适合于宏观层面的地震安全性评估.然而对于具体的工程项目,即使采用确定性分析,虽然具有确定的震级,震中距等概念,但是在预测地震动参数时,地震波从激发到传播的物理过程被概化,仍然不能包含地震个性化信息.
2.2 基于谱元法模拟的地震危险性分析谱元法是一种数值方法,由于其结合了有限元法的灵活性和伪谱法的精度,且便于大规模并行计算,在地震模拟中得到了越来越广泛的应用.近年来,一些学者开始尝试在地震危险性分析中运用地震学的数值模拟方法,以替代经验衰减关系.对于重大高坝枢纽工程,特别是附近有孕震大断层的情况,考虑由此引发的某一特定地震甚至是极端情况下可能的破坏就显得十分必要.
2.2.1 谱元法简介谱元法与常规有限元法不同之处仅在于插值点位置和积分方法的选择.谱元法以正交多项式极值点为插值点,正交多项式可以选择为Chebyshev多项式(Chebyshev谱元法)或者Legendre多项式(Legendre谱元法),本文采用Legendre多项式,其递推式为
(2) |
式中,下标n表示多项式阶数,其中L0=1,L1=x.数值积分采用Gauss-Lobatto积分,积分点位置同插值点位置,对应的积分权值为
(3) |
其中n为Legendre多项式阶数,ξj为局部坐标值.图 1为一个四阶谱单元的插值点位置和插值函数示意图.
考虑震源-介质-河谷场址的综合模型由各个子模型组合而成.基于此模型的分析结果不仅针对特定地震事件,而且得到的地震动还包含了震源机制,非均匀介质和局部场地地形效应的影响.
对于能够简化成点源的震源而言,震源时间函数和矩张量能够分别描述位错震源在时间域和空间域中的性质.其中震源时间函数简化成光滑的高斯函数(震源速度时程)(Komatitsch and Tromp, 2002),矩张量的表达式则可以写成
(4) |
中大地震的断层破裂尺度较大,场址离断层较近的情况下,震源不能简化为点源,因此震源模型需采用有限断层(Irikura, 1986).有限断层模型将断层面离散成一系列子源,每个子源可视作一个点源,这些子源被激发的相互时空关系为核心和难点所在.一些影响较大的地震(如1994北岭地震)经过学者们的大量研究,其断层破裂的具体过程的反演已经获得一些成果(Wald et al., 1996),但是对于一些潜在的孕震断层,破裂模式只能基于假定.
波速结构是地震数值模拟中最核心的部分,也是数值方法区别于经验方法的本质所在.通过网格划分和材料赋值,波动方程在整个时空域被确定性地求解,不含有任何的人为假定及经验规律的套用,但是其计算量非常巨大.已有研究表明,基于高精度的波速结构,谱元法计算可以获得在较宽频带上与实测记录匹配良好的地震动(He et al., 2015).目前美国加州地区已经建立了十分精细的速度结构模型(Plesch et al., 2011).但是我国相关的资料还比较分散,缺乏系统性的集成,这方面工作还亟待进一步的研究.
不同于民用建筑,高坝跨度都在几百米数量级而且位于高山峡谷,地震动空间分布不均匀现象十分突出(Alves and Hall, 2006),是模型中必须考虑的内容.得益于谱元法本身所具备的优势,可直接整合数字高程模型(DEM)并进行精细化的地形模拟(Lee et al., 2009),考虑坝址区高山峡谷地形的放大效应.
3 算例与对比分析 3.1 模拟方法验证采用Graves (1996)的的算例,计算均质半无限空间中具有一定埋深的双力偶点源激发的波场.如图 2所示,模型尺寸沿x,y,z方向分别为40 km×40 km×20 km,均匀划分为48×48×20谱单元.传播介质设为均质体:P波波速4000 m·s-1,S波波速2300 m·s-1,密度1800 kg·m-3.震源埋深2.5 km,走向、倾角和滑动角均为90°,总地震矩M0为1.0×1016 N·m,观测点取地表面距震中10 km,方位角φ=90°时的切向速度分量以及φ=0°时的径向和竖向速度分量.计算结果如图 3所示,虚线为谱元法计算结果,实线为频率-波数积分法(Wang and Herrmann, 1980)结果.从图中可以看出谱元法和频率-波数积分法的计算结果基本重合,谱元法的双力偶震源模型和弹性波动过程模拟具有良好的计算精度,验证了本文模拟方法的合理性.
以我国四川大岗山水电站高拱坝场址为背景,假设发震断层距离坝址最小距离为5.5 km,地震震级Mw=6.0,震源深度10 km,断层尺寸为10 km (走向)×9 km (倾向)、倾角70°、正南北走向,总地震矩1.259×1016 N·m,Rake角0°(纯走滑型断层)(四川省地震局工程地震研究院等,2004).采用第2节介绍的谱元法模型,分析坝基面的三维地震动分布.
以场址为中心,选取26.25 km×40.00 km×20.00 km区域为模型计算范围,如图 4a所示,均匀划分地为240×480×50个谱单元,总自由度数约为9.96亿.波速结构采用王椿墉等(Wang et al., 2007)的成果.震源采用有限断层模型,整个断层面划分为90个子断层,每个子断层尺寸为1 km×1 km,断层平面划分及破裂方式如图 4b所示.在场址河谷处沿横河向设置9个观测点以获得空间非均匀地面运动,如图 4c所示.
图 5所示分别为考虑实际地形(图 4a)和不考虑实际地形(假定地表为平坦面)情况下的波场快照图(东西方向速度分量),时间节点分别为0.4、1.2 s、2.0、2.8、3.6 s.从图 5的波场对比可以看出:(1)从整体来看,地形的起伏对整体波场的分布规律影响不显著,较大尺度的有限断层破裂所激发的近场波场呈现出椭圆形(图 5中的条纹分布),考虑实际地形的地震波场与平坦地形的分布仍然大致相同;(2)从局部来看,考虑实际地形时的波阵面不再是平坦地形条件下的光滑曲面,而是随着地形的起伏产生局部的扭曲,相邻条纹间相互渗透,甚至出现局部区域与周围反方向的运动(如图 5时刻3.6 s).地形起伏对局部地震波场的影响机制主要是因为山脊(山谷)放大(削弱)作用,以及地形高低导致的行波效应.对高坝工程的场址河谷而言这种影响不可忽略,需要予以考虑.
沿河谷横向布置的9个观测点分别提取加速度时程最大值,如图 6和表 1所示.由此可以看出,局部地面运动的空间分布和地形有紧密的联系,沿河谷两岸明显呈现不均匀的特性.竖直方向上分布相比于横河向较为均匀,峰值在289~336 Gal之间,河谷和两岸差别不大.然而水平方向的峰值差别很大,横河向运动的峰值分布大致和地形轮廓一致,呈中间小两边大的近似“V”字形分布,分布范围从河谷处最小的141.3 Gal到左岸最大的347.1 Gal,相差近2.5倍.顺河向峰值则是呈右岸大左岸小的分布趋势,峰值范围(136.3~321.3 Gal)则和横河向接近.这表明地形效应对不同波型的影响结果不尽相同,对主要由S波控制的水平向运动的影响在宏观上仍然表现为河谷缩小,岸坡放大的规律,但是对由P波控制的竖直向运动的影响作用不是特别明显.
为了与传统基于衰减关系的分析方法进行对比,这里进一步采用确定性地震危险性分析方法,利用公式(2)计算相应的场址地震动参数,回归系数分别采用美国西部和研究区域长轴和短轴衰减系数(四川省地震局工程地震研究院等,2004),如表 2所示.取断层破裂起始点到5号测站(位于河谷中心处)距离作为震中距(即R=7.2 km),则基于三种衰减关系得到的6.0级地震场址水平峰值加速度分别为232.4、318.9和177.3 Gal (表 2).相应地,谱元法计算得到的9个观测点水平峰值加速度均值为343.3 Gal,与本研究区域长轴参数计算所得结果(318.9 Gal)相比略微偏大.由此可见使用谱元法数值计算不仅能够在平均意义上获得与现行危险性分析具有可比性的结果,更能够在局部体现出波场的空间不均匀性,具有一定的优势.
这里需要说明的是,场地河谷地面运动由震源、传播路径和局部场地共同决定,每一个因素都对结果都会产生一定的影响,本工程实例所呈现的物理现象不一定是普遍规律,但对于任何一个类似的模型模拟都具有一定的参考意义.
4 结论与展望高坝的抗震安全分析与评价包括三个环节,即场址地震动参数(峰值加速度、反应谱和地震历时)的确定、高坝地震反应的分析以及抗震安全性的评价.坝址地震动参数的确定是高坝枢纽抗震分析的先决条件,是抗震设计和安全评估的基础.基于地震学和工程学之间的学科交叉思想,本文尝试了一种新的坝址地震动参数确定思路:针对设定地震,采用谱元法模拟震源破裂的时空过程、地震波传播的物理机制以及坝址峡谷的局部场地效应,生成坝址区的三维地震动场.通过实例分析表明,场址峰值加速度计算结果与传统基于衰减关系的计算结果具有一定的可比性.与传统的坝址地震危险性分析方法相比较,该模拟方法具有以下几方面的特点:
(1)针对坝址区的特定潜在震源,采用数值方法模拟地震波发生和传播的物理过程,可以生成更符合高坝枢纽实际地震地质构造、区域岩体动力特性以及坝址河谷地形地质条件的地震动荷载,体现不同地震激发的坝址峡谷地震响应的特异性.这为进行最大可信地震评估提供了一种新方法,以合理预测极端地震作用下坝址区的地震动特性.
(2)数值模拟可以生成整个坝址枢纽区的三维地震动场,尤其是对高坝地震响应特性有重要影响的坝基面非均匀地震动荷载,克服了通过假想平坦地区基岩出露点确定坝址地震动参数的不足.
(3)根据模拟得到的坝址区三维地震动场,可以将整个高坝枢纽作为一个整体进行抗震设计和安全性评价,使各个组成部分到达相似的抗震安全度,避免薄弱环节的出现.如汶川地震中,太平驿、福堂等水电站尽管坝体结构完好,但是因泄洪系统失效造成大坝漫顶等震害.同时,也可以针对整个枢纽不同组成部分的抗震设防要求,生成不同安全等级标准的地震荷载,分别用于大坝、高边坡、地下厂房结构等的抗震分析与评价.
(4)针对一些特殊重大工程,基于这一研究思路,可以通过库区和坝体的弱震记录反演,不断对模拟完善,提高模拟精度.一旦地震发生,就可以根据地震震级和震源破裂机制,以及一些测站的地震记录,仿真模拟坝体及附属结构、地下厂房及机电设备、库区高边坡的震害情况,为震后采取应急措施提供科学的决策依据.
当然,本文只是初步提出了震源破裂-传播介质-坝址峡谷地震响应模拟的基本思路,还需要不断完善,仍然面临着诸多挑战性的问题.如坝址到潜在震源区域介质参数的确定,超大规模计算的能力限制,模型尺度耦合,坝址峡谷地质地形非线性效应的影响等,有待于进一步深入的研究.
Alves S W, Hall J F. 2006. Generation of spatially nonuniform ground motion for nonlinear analysis of a concrete arch dam. Earthquake Engineering and Structural Dynamics, 35(11): 1339-1357. DOI:10.1002/(ISSN)1096-9845 | |
Baker J W, Luco N, Abrahamson N A, et al. 2014. Engineering uses of physics-based ground motion simulations.//Proceedings of the 10th National Conference in Earthquake Engineering. Anchorage, Alaska:Earthquake Engineering Research Institute. | |
Cornell C A. 1968. Engineering seismic risk analysis. Bulletin of the Seismological Society of America, 58(5): 1583-1606. | |
Cui Y, Poyraz E, Olsen K B, et al. 2013. Physics-based seismic hazard analysis on petascale heterogeneous supercomputers.//Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis. Denver, Colorado:ACM, Article No. 70. | |
Douglas J, Aochi H. 2008. A survey of techniques for predicting earthquake ground motions for engineering purposes. Surveys in Geophysics, 29(3): 187-220. DOI:10.1007/s10712-008-9046-y | |
Engineering earthquake institute of Sichuan earthquake administration. 2004. Feasibility study report for the Dagangshan hydropower station:Attachment 3-1:Seismic safety evaluation report for the engineering site (in Chinese). | |
Graves R W. 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences. Bulletin of the Seismological Society of America, 86(4): 1091-1106. | |
Graves R, Jordan T H, Callaghan S, et al. 2011. CyberShake:A physics-based seismic hazard model for southern California. Pure and Applied Geophysics, 168(3-4): 367-381. DOI:10.1007/s00024-010-0161-6 | |
Hall J F, Heaton T H, Halling M W, et al. 1995. Near-source ground motion and its effects on flexible buildings. Earthquake Spectra, 11(4): 569-605. DOI:10.1193/1.1585828 | |
Hall J F. 1998. Seismic response of steel frame buildings to near-source ground motions. Earthquake Engineering and Structural Dynamics, 27(12): 1445-1464. DOI:10.1002/(ISSN)1096-9845 | |
He C H, Wang J T, Zhang C H, et al. 2015. Simulation of broadband seismic ground motions at dam canyons by using a deterministic numerical approach. Soil Dynamics and Earthquake Engineering, 76: 136-144. DOI:10.1016/j.soildyn.2014.12.004 | |
Heaton T H, Hall J F, Wald D J, et al. 1995. Response of high-rise and base-isolated buildings to a hypothetical Mw7.0 blind thrust earthquake. Science, 267(5195): 206-211. DOI:10.1126/science.267.5195.206 | |
Hu X Y, Zhou K S, Yan X J. 1996. A method for evaluation of ground motion in regions with few acceleration observation data. Earthquake Engineering and Engineering Vibration (in Chinese), 16(3): 1-10. | |
Huo J R. 1989. Studies on ground attenuation rule of the near fault[Ph.D.thesis] (in Chinese). Harbin:Institute of Engineering Mechanics, China Earthquake Administration. | |
Irikura K. 1986. Prediction of strong acceleration motions using empirical Green's function.//Proc. 7th Japan Conf. Earthq. Eng.. Tokyo. | |
Komatitsch D, Vilotte J P. 1998. The spectral element method:An efficient tool to simulate the seismic response of 2D and 3D geological structures. Bulletin of the Seismological Society of America, 88(2): 368-392. | |
Komatitsch D, Tromp J. 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation. Geophysical Journal International, 139(3): 806-822. DOI:10.1046/j.1365-246x.1999.00967.x | |
Komatitsch D, Tromp J. 2002. Spectral-element simulations of global seismic wave propagation-I. Validation. Geophysical Journal International, 149(2): 390-421. DOI:10.1046/j.1365-246X.2002.01653.x | |
Krinitzsky E L. 1995. Deterministic versus probabilistic seismic hazard analysis for critical structures. Engineering Geology, 40(1-2): 1-7. DOI:10.1016/0013-7952(95)00031-3 | |
Krishnan S, Chen J, Komatitsch D, et al. 2006. Case studies of damage to tall steel moment-frame buildings in southern California during large San Andreas earthquakes. Bulletin of the Seismological Society of America, 96(4A): 1523-1537. DOI:10.1785/0120050145 | |
Lee S J, Chan Y C, Komatitsch D, et al. 2009. Effects of realistic surface topography on seismic ground motion in the Yangminshan region of Taiwan based upon the spectral-element method and LiDAR DTM. Bulletin of the Seismological Society of America, 99(2A): 681-693. DOI:10.1785/0120080264 | |
Luo Q F, Shi C X, Cao B Z, et al. 2003. Introduction of structural aseismic research in China in recent four years. Acta Seismologica Sinica (in Chinese), 25(6): 637-652. | |
Plesch A, Tape C, Graves R, et al. 2011. Updates for the CVM-H including new representations of the offshore Santa Maria and San Bernardino basin and a new Moho surface.//Southern California Earthquake Center Annual Meeting, Proceedings and Abstracts, vol. 21. | |
Power M, Chiou B, Abrahamson N, et al. 2008. An overview of the NGA project. Earthquake Spectra, 24(1): 3-21. DOI:10.1193/1.2894833 | |
U.S. Army Corps of Engineers. 1995. Engineering and design:earthquake design and evaluation for civil works projects. Regulation No. 1110-2-1806, Washington. | |
Wald D J, Heaton T H, Hudnut K W. 1996. The slip history of the 1994 Northridge, California, earthquake determined from strong-motion, teleseismic, GPS, and leveling data. Bulletin of the Seismological Society of America, 86(1B): S49-S70. | |
Wang C Y, Herrmann R B. 1980. A numerical study of P-, SV-, and SH-wave generation in a plane layered medium. Bulletin of the Seismological Society of America, 70(4): 1015-1036. | |
Wang C Y, Han W B, Wu J P, et al. 2007. Crustal structure beneath the eastern margin of the Tibetan Plateau and its tectonic implications. Journal of Geophysical Research-Solid Earth, 112: B07307. | |
Xu Y J, Mu H L, Zhang C H, et al. 2012. 3D finite element modeling of seismic responses of Baozhusi gravity dam in Ms8.0 Wenchuan Earthquake.. Chinese J. Geophys. (in Chinese), 55(1): 293-303. DOI:10.6038/j.issn.0001-5733.2012.01.029 | |
Yao Z X, Zheng T Y. 1985. Estimation of strong ground motion for the Ertan dam region. Acta Geophysica Sinica (in Chinese), 28(2): 170-180. | |
Zhou J P, Chen G F, Dang L C. 2007. Review on anti-seismic safety evaluation for high dams in China. Journal of Hydraulic Engineering (in Chinese)(S1): 54-59. | |
胡聿贤, 周克森, 阎秀杰. 1996. 缺乏地震动加速度记录地区地震动估计的映射法. 地震工程与工程振动, 16(3): 1–10. | |
霍俊荣. 1989.近场强地面运动衰减规律的研究[博士论文].哈尔滨:国家地震局工程力学研究所. http://www.oalib.com/references/18564404 | |
罗奇峰, 石春香, 曹炳政, 等. 2003. 中国近4年结构抗震进展介绍. 地震学报, 25(6): 637–652. | |
四川省地震局工程地震研究院, 中国地震局地质研究所等. 2004.四川省大渡河大岗山水电站预可行性研究报告附件-工程场地地震安全性评价报告. | |
徐艳杰, 牟海磊, 张楚汉, 等. 2012. 汶川地震中宝珠寺重力坝地震响应的三维有限元模拟. 地球物理学报, 55(1): 293–303. DOI:10.6038/j.issn.0001-5733.2012.01.029 | |
姚振兴, 郑天愉. 1985. 对二滩水电站坝区场地地面运动的估计. 地球物理学报, 28(2): 170–180. | |
周建平, 陈观福, 党林才. 2007. 我国高坝抗震安全评价的现状与挑战. 水利学报(S1): 54–59. | |