2. 中国地震局地球物理研究所, 北京 100081;
3. 江苏省土木工程防震技术研究中心, 南京 210009;
4. 江苏省地震局, 南京 210014
2. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China;
3. Civil Engineering & Earthquake Disaster Prevention Center of Jiangsu Province, Nanjing 210009, China;
4. Institute of Earthquake Engineering of Jiangsu Province, Nanjing 210014, China
强震动记录和震害调查均表明,局部场地条件对地震动分布有重要影响(Bouchon M, 1973; Spudich et al., 1994; Bouchon and Barker, 1996; 王海云等, 2010; 王伟等, 2014; 迟明杰等, 2015; Meunier et al., 2008).局部场地条件包括:局部场地地形变化及局部场地介质的非均匀性.当仅考虑场地地形影响时,通常假定场地为均匀线性半无限空间模型.目前, 研究地形效应的方法主要有强震动地形效应观测、解析法和数值模拟三类方法.
地形效应观测台阵是研究场地局部地形效应对地震动影响的一个主要手段, 观测记录提供了对地震动定量分析的条件, 将地震动记录分析与震害现象相结合, 可更好地研究地震动的地形效应, 但地形效应观测依懒于台阵的分布, 成本高昂.解析法主要针对可使用简单参数描述的模型,如圆弧形(Trifunac, 1972; Yuan et al., 1992)、椭圆形(Wang et al., 1982)、U形(Gao Y et al., 2012)以及V形地形(Gao Y et al., 2013; Chang et al., 2015),这些解析方法可以从本质上解释地形对地震动的散射问题,同时也为数值模拟的精度验证提供了依据.数值方法对复杂地形及土层的场地具有良好适用性, 可以弥补解析方法在求解复杂场地时的不足.数值计算方法主要包括:有限差分法、有限元法、伪谱法和谱元法等.有限差分法是最早被用于地震波传播数值模拟研究中的.该方法算法简单, 计算速度快, 占用内存小.李小军等(1995)早年提出了适用于任何地形情况下、具有较高计算精度和技术稳定性的显式有限元-有限差分方法, 可以很好地模拟二维黏弹性场地地形对地震动的影响.这一方法避免了有限差分法只适用于相对简单的地质模型的不足, 发挥了有限元法较适合处理几何条件和物理条件复杂的问题优势, 其不足是对高频和短波长信号模拟时需提高网格分辨率, 会增加模型计算量.伪谱法的主要特点是计算速度快, 精度高, 其缺点是当在空间中任意一点发生变化时都会改变频域的所有值, 因此只能模拟横向波速变化不大的情况.谱元法融合了有限元方法和伪谱法的思想, 但其在理论上还有待进一步完善.
随着国家对西部开发的大力支持, 高速公路、大型水坝、长输管道、高铁等跨河谷和山谷的工程大量涌现, 局部场地地形的影响越来越受到重视.目前, 数值方法是处理地形效应的首选方法, 能够满足各种实际工程中复杂地形场地的地震动响应分析, 并给出定量结果.荣棉水等(2007)研究了两个相邻凸起台地地形对入射P波傅氏谱特性的影响; 赵瑞斌等(2015)针对阶梯地形地震动波动问题, 在黏弹性人工边界的基础上提出了“虚拟对称自结构体系”地震动输入技术, Chen等(2015)利用ABAQUS软件研究由河漫滩、阶地、波状起伏的丘陵地带和残丘等组成的二维复合地形场地的地震效应特征; 唐辉等(2012)通过对西山自贡公园的流动观测台阵数据分析, 结合数值模拟的研究方法, 对山体地形地震动影响作了系统的研究;蒋涵等(2015a)针对庐山地区, 采用三维谱元法和一阶透射边界模拟了地震波传播, 研究了该地区地震动的频谱特征;梁建文等(2014)采用了黏弹性边界和以等效节点力输入地震动的方法研究了三维凹陷地形地震响应.这些研究中大多采用的是二维模型(荣棉水等, 2007; 赵瑞斌等, 2015; Chen et al, 2015; 唐辉等, 2012), 少数考虑了三维的情况但也往往对模型做较大简化处理(蒋涵等, 2015b; 梁建文等, 2014), 致使三维局部地形影响反应不明显.
本文在显式有限元方法的基础上用真实的二维侧边界模型求解三维复杂地形场地的边界自由场, 将三维转化二维的解法结合黏弹性人工边界, 建立了三维复杂地形条件下地震动垂直输入的一种简化方法.并以四川省桃坪羌寨地区作为研究对象, 利用前处理软件HyperMesh进行有限元网格划分, 根据地表高程数据建立了实际地形的二维和三维场地模型.研究了真实河谷地形场地对地震波传播的影响.同时对三维模型与二维模型进行了对比, 反映了二维模型的局限性以及对河谷地形条件下建立真实三维场地模型的必要性.
2 三维复杂地形地震动输入方法及验证 2.1 三维黏弹性人工边界本文所采用的数值方法是显式有限元法结合黏弹性人工边界的时域整体分析方法.采用显式有限元法进行复杂场地地震动研究时必须考虑人工边界条件, 首先需将一无限域转化为有限域, 并在有限域上设置人工边界条件模拟外部无限域的影响(梁建文等, 2014), 在设置人工边界后, 必须保证波在人工边界处的传播特性与在原连续介质中一致, 才能实现对原连续介质的精确模拟.目前应用较多的人工边界有透射边界(廖振鹏等, 1982; 李小军等, 1996)、黏性边界(Lysmer et al., 1969)和黏弹性人工边界(梁建文等, 2014; 何建涛等, 2010; Deeks et al., 1994; DU et al., 2010)等.其中黏弹性人工边界能够较好地模拟无限地基的弹性恢复能力和辐射阻尼效应, 同时黏弹性人工边界也易于在通用有限元软件上实现.黏弹性人工边界相当于在人工边界上设置一系列由线性弹簧与黏滞阻尼器并联的弹簧-阻尼物理元件.刘晶波等(1998;2005;2006) 给出了三维黏弹性人工边界参数, 为三维无限域波动问题的求解提供了方便.三维黏弹性人工边界的弹簧-阻尼元件参数为法向:
(1) |
切向:
(2) |
其中, Kbn、Cbn为法向边界的弹簧-阻尼元件参数; Kbt、Cbt为切向边界的弹簧-阻尼元件参数; ρ为介质密度; cp=
刘晶波等(1998)通过将地震波转化为节点等效力的方式实现了波动的输入.得到了在模型边界点上所需施加等效荷载的一般计算公式:
(3) |
其中, σl(t)为边界节点l处的等效力; σ0(x, y, t)为连续介质中由原自由场产生的应力; K、C分别为黏弹性人工边界的法向与切向的弹簧和阻尼系数, 由式(1) 和式(2) 得到; u(x, y, t)、
在三维复杂地形条件下直接求解其边界各节点的自由场比较困难, 本文用二维侧边界模型求解三维复杂地形场地的边界自由场, 将其退化为二维问题求解各边界节点的等效应力.
如图 2a所示, 通过设置一个水平底边界面以及相互垂直或平行且都垂直于底边界面的四个侧边界面, 将实际无限域三维地形场地变成三维有限空间分析模型.对于侧边界外被切除的区域近似地视为自由表面向外无限延伸的半空间.如图 2b所示, 此时每一个侧边界面上的自由场可以视为该边界面形成的一个二维场地的自由场.以垂直底边界入射的平面SV波(y轴方向位移时程为u0(t))的情况为例, 考虑平面内的介质为线弹性体.则底边界面和边界面①、② 为平面内波动问题(SV波), 边界面③、④ 为出平面波动问题(SH波).根据弹性波在均匀介质中的传播特性可以得到平面①、② 的二维模型底边界线和两侧边界线上各节点等效应力的以下计算公式.
底部人工边界线处:
(4) |
左侧人工边界线处:
(5) |
右侧人工边界线处:
(6) |
式中, hl1、hr1分别为左右两个侧边界线上节点与底边界的高程差; hl、hr分别为边界面①、② 左右两个边界线的长度.利用上面公式确定边界各节点的等效应力后, 利用章小龙等(2016)给出了二维模型自由场的有限元求解方法,求解出底边界面和边界面①、② 上各节点的自由场; 平面③、④ 为出平面问题(SH波), 这两个侧边界面SH波垂直入射的自由场计算, 采用自编的二维Fortran有限元程序实现, 其平面内各节点的自由场求解方法可参见陈少林等(2014), 将其各节点的自由场代入公式(3) 便可以得到四个侧边界面上各节点的等效应力, 对于上文平面SV波输入的特殊情况, 三维空间有限模型底边界为一致输入面, 底边界面上各节点的等效应力可由式(4) 确定.至此, 实现了三维空间有限模型的地震动输入.
2.4 地震动输入方法的精度验证为了验证SV波垂直入射方法的模拟精度, 分析平面SV波垂直入射条件下三维均匀弹性半空间的动力反应问题.入射SV波选用一宽0.25 s、持时10 s的近似狄拉克脉冲, 其位移时程曲线及傅里叶谱如图 3所示.有限元截取2000 m×2000 m×2000 m的正立方体, 离散单元网格大小为25 m×25 m×25 m.波动分析介质参数如表 1所示.人工边界和人工边界处等效节点力按1.1—1.3节所述方式施加, 显式有限元方程的时间积分步长取为0.00001 s.
图 4a为平面SV波沿着z轴正方向入射时计算模型自由表面中心点的位移时程; 图 4b为自由表面中心点y轴方向输出的位移时程曲线与解析解(Achenbach J D, 1973)的对比, 可以看出本文所述方法的计算结果与理论解吻合较好.图 4c—4e三幅图为1.14~1.50 s这一时间段内自由表面中心点的空间轨迹图, 可知误差在可控范围之内, 说明了本文处理方法具有良好的模拟精度.
本文采用四川省桃坪地区一山谷地形为研究对象, 计算区域如图 5所示, 矩形区域四个端点的经纬度分别为:(103°25′E, 31°30′N)、(103°30′E, 31°30′N)、(103°30′E, 31°36′N)和(103°25′E, 31°36′N), 地表东西向和南北向的跨度分别为7.3 km和9.6 km, 地表高程差2.2 km, 从图 5可以看出, 该场地地形复杂且山体起伏剧烈, 为典型的山谷地形场地.在S-N向分别设置了四条测线SL1—SL4, 四条测线在E-W向上的坐标分别为1.7 km、2.6 km、4.8 km和5.7 km.根据四条测线SL1—SL4上各点的数据高程分别建立二维模型, 点Ai~Ei(i=1~4) 分别为四条测线上的观察点, 其中点Ai、Ei位于两侧山峰, 点Bi、Di位于两侧山坡, 点Ci坐落在谷底, 四个剖面如图 7所示.模型底部与地表高程最低点的距离为2.2 km, 模型的五个人工边界面均为平面, 设置为黏弹性人工边界.在建立三维模型时采用四面体单元对模型进行离散, 单元网格大小为25 m, 为尽可能还原真实地形, 三维大场地模型在通用有限元软件中直接进行网格划分往往比较困难, 本文利用有限元网格划分前处理软件HyperMesh根据该计算区域地表高程数据建立三维场地模型.三维模型节点和单元个数分别为9743380和17663381.对二维模型采用三角形单元进行离散, 网格最大尺寸为25 m.如图 6所示, 首先建立三维曲面模型, 在表面进行布种网格划分生成表面三角形网格部件, 再利用HyperMesh根据推进波前法生成模型内部四面体网格即可.
仅考虑局部场地地形的影响时假定模型为均匀介质, 模型介质材料参数见表 1.为了模拟真实岩体中的耗能情况, 采用瑞利阻尼, 阻尼系数根据场地1阶自振频率确定, 阻尼比取为0.002.输入波动位移时程如图 3所示, 截止频率为10 Hz, 沿着S-N(横河谷)向运动的地震波垂直底边界输入.有限元计算的时间积分步长取为0.0001 s, 所有模型计算总时程均为10 s.
4 分析与讨论 4.1 三维场地地震反应图 8为地表各节点总位移云图, 图 8a、8b两幅位移云图显示了地震波从下向上传播过程; 图 8c显示了地震波在山峰处汇聚可以观察到“鞭稍效应”, 8d可以观察到地震波向下传播时消散的现象, 表明了对该模型的计算分析满足波动要求.根据云图的特征, 可以发现在地震波传播的过程中位于山脊上的位移峰值相对于山坡上的结果相对较大, 位于山顶上的位移峰值总体相对较大.地震波在自由表面反射后, 向下传播, 当散射场通过人工边界后山体震动最后归于静止, 符合实际场地地震波的传播情况, 进一步表明了本文针对复杂场地地震动输入方法的合理性.
图 9为三维模型地表PGD的分布平面图.由图 9可知, 地表上任意一点在S-N向的PGD都较其他两个方向上的反应强烈, 说明地表地震动沿着输入地震动方向的影响比其他方向的影响更大, 高PGD区分布于山顶、山脊, 且各高PGD区分布比较离散, 主要是由于地震波在复杂陡峻地形条件下的折射和散射作用所造成的.正是由于地形作用, 并非所有的峰顶、山脊处都得到放大.
图 10a—10d分别为四个二维模型地表与三维模型的四条测线上的地震动(记录点水平间距均为118 m).每条测线上的地震动都用此测线上所有记录点的PGD的最大值归一化,图右上角注明了每条测线上所有记录点的PGD最大值, 单位cm.
从图 10可以看出, 二维与三维模型在S-N方向的输出位移整体趋势有相似性, 并且二者在SV波初至时间上也基本一致, 当SV波抵达地表后, S-N向位移逐渐增大, 各记录点波幅随记录点的位置而改变.在Z向位移反应中可以清楚地观察到P波, 由于P波波速大于SV波, P波较先抵达地表.不同测线所对应的最大位移幅值也不同.地表记录点在沿着地震动输入方向的影响比其他方向的影响更大, 同时三维模型在接近边界处的位移时程曲线振荡持时明显延长; 充分说明与已有的二维概化模型计算结果相比, 三维真实地形起伏多变, 地形对地震动的影响比二维模型更加复杂.
图 11为四条测线在地表S-N向上的峰值位移分布平面图, 可以看出二维模型与三维模型测线上的峰值位移分布都随着坡度比的增大而增大, 也基本满足在谷底处衰减, 在山峰和凸起处增强的趋势, 二维模型的地表位移幅值整体上比三维真实地形大, 在局部区域三维计算结果明显高于二维计算结果.
图 12ai—12ei分别为剖面1—4上观察点Ai~Ei(i=1~4) 的位移频谱.在四条测线中, 测线3相对其他的三个剖面的地表起伏比较复杂, 测线3地表上5个记录点的移频谱整体变化趋势相对于其他剖面的记录点较复杂; 测线1、测线2和测线4为简单的V型河谷地形, 可以看出由于地形条件的不同, 每个剖面记录的位移频谱也各异, 在山峰处记录点的位移频谱幅值整体上比其他点的反应大; 同一记录点, 二维模型相对于三维模型在高频处影响更大, 二维模型明显地放大了场地反应, 说明三维真实地形对地震动的影响不能简单地用二维模型来替代.
本文在黏弹性人工边界和显式有限元方法的基础上, 建立了三维复杂地形条件下地震波垂直输入时人工边界自由场求解的一种简化方法, 将三维场地的地震波垂直入射时边界面上自由场退化为二维问题求解, 通过数值算例验证了该方法的合理性.
同时以四川桃坪地区一山谷地形作为研究对象, 采用二维与三维模型进行对比分析.结果表明: (1) 可用二维侧边界模型求解三维复杂地形场地的边界面上自由场, 本文数值计算结果表明该方法的合理性.因此, 在求解具有不规则地形的三维场地地震反应时, 可采用二维侧边界模型来求解其边界面上的自由场; (2) 在输入脉冲波情况下, 自由面对地震动有局部放大作用, 地震反应在沿地震动输入方向比其他方向的影响更大; (3) 山谷的散射作用会影响地表位移时程曲线的波形并使其振荡持时延长, 该现象三维模型比二维模型模拟更明显; (4) 不同地表测线可能会获得不同幅值变化结果, 二维模型与三维真实地形测线的PGD分布基本满足在谷底处衰减, 在山峰和凸起处增强的趋势, 但二维模型地表整体上比三维真实地形测线的PGD大, 局部区域三维计算结果反而高于二维结果, 采用二维模型时地表的位移频谱在高频处相对于三维模型有明显的放大现象.
研究表明, 在三维复杂的地形中地形特征对于峰值位移的分布影响很大.建立真实地形起伏的三维模型对开展地形条件下地震动的影响是非常必要的.本文方法的思路可推广应用于三维复杂地形条件下地震动斜入射问题求解.
致谢感谢南京航天航空大学土木工程系陈少林教授提供的二维场地SH波垂直入射时自由场计算Fortran有限元程序, 感谢匿名审稿专家提出的宝贵建议.
Achenbach J D. 1973. Wave Propagation in Elastic Solids. Amsterdam: North-Holland Publishing Company. | |
Bouchon M. 1973. Effect of topography on surface motion. Bulletin of the Seismological Society of America, 63(2): 615-632. | |
Bouchon M, Barker J S. 1996. Seismic response of a hill:the example of Tarzana, California. Bulletin of the Seismological Society of America, 86(1A): 66-72. | |
Chang K H, Tsaur D H, Wang J H. 2015. Response of a shallow asymmetric V-shaped canyon to antiplane elastic waves.//Proceedings of the Royal Society of London A:Mathematical, Physical and Engineering Sciences. The Royal Society, 471(2174):20140215. | |
Chen G X, Jin D D, Zhu J, et al. 2015. Nonlinear analysis on seismic site response of Fuzhou basin, China. Bulletin of the Seismological Society of America, 105(2A): 928-949. DOI:10.1785/0120140085 | |
Chen S L, Zhang L L, Li S Y. 2014. Numerical analysis of the plane SH waves scattering by semi-cylindrical alluvial valley. Engineering Mechanics, 31(4): 218-224. | |
Chi M J, Li X J, Chen B, et al. 2015. The effect of topographical and soil condition on housing damage during the Yingjing earthquakes. Earthquake Engineering and Engineering Dynamics, 35(2): 94-102. | |
Deeks A J, Randolph M F. 1994. Ax symmetric time-domain transmitting boundaries. Journal of Engineering Mechanics, ASCE, 120(1): 25-42. DOI:10.1061/(ASCE)0733-9399(1994)120:1(25) | |
Du X L, Zhao M. 2010. A local time-domain transmitting boundary for simulating cylindrical ealstic wave propagation in infinite media. Soil Dynamics and Earthquake Engineering, 30(10): 937-946. DOI:10.1016/j.soildyn.2010.04.004 | |
Gao Y, Zhang N, Li D, et al. 2012. Effects of topographic amplification induced by a U-shaped canyon on seismic waves. Bulletin of the Seismological Society of America, 102(4): 1748-1763. DOI:10.1785/0120110306 | |
Gao Y F, Zhang N. 2013. Scattering of cylindrical SH waves induced by a symmetrical V-shaped canyon:near-source topographic effects. Geophysical Journal International, 193(2): 874-885. DOI:10.1093/gji/ggs119 | |
He J T, Ma H F, Zhang B Y, et al. 2010. Method and realization of seismic motion input of viscou-spring boundary. Journal of Hydraulic Engineering, 41(8): 960-969. | |
Jiang H, Zhou H, Gao M T. 2015a. The characteristics of frequency domain of ground motion in 3-D topography-A case study of Lushan area. Technology for Earthquake Disaster Prevention, 10(1): 59-67. | |
Jiang H, Zhou H, Gao M T. 2015b. A study on the correlation of ridge line and slope with Peak Ground Velocity amplification factor. Chinese Journal of Geophysics, 58(1): 229-237. DOI:10.6038/cjg20150120 | |
Li X J, Liao Z P. 1996. The drift instability of local transmitting boundary in time domain. Acta Mechanica Sinica, 28(5): 627-632. | |
Li X J, Liao Z P, Guan H M. 1995. An explicit finite element-finite difference method for analyzing the effect of visco-elastic local topography on the earthquake motion. ACTA Seismologica Sinica, 8(3): 447-456. DOI:10.1007/BF02650573 | |
Liang J W, Qi X Y, Ba Z N. 2014. Seismic respones analysis of 3D canyon based on the viscous-spring boundary. Earthquake Engineering and Engineering Dynamics, 34(4): 21-28. | |
Liao Z P, Yang B P, Yuan Y F. 1982. Artificial boundary in analysis of transient elastic wave propagation. Earthquake Engineering and Engineering Vibration, 2(1): 1-11. | |
Liu J B, Gu Y, Du Y X. 2006. Consistent viscous-spring artificial boundaries and viscous-spring boundary elements. Chinese Jounal of Geotechnical Engineering, 28(9): 1070-1075. | |
Liu J B, Li B. 2005. Three dimensional static-dynamic unified viscoelastic artificial boundary. Science in China:Series E Engineering & Material Science, 35(9): 966-980. | |
Liu J B, Lu Y D. 1998. A direct method for analysis of dynamic soil-structure interaction. China Civil Engineering Journal, 31(3): 55-64. | |
Lysmer John, Kuhlemeyer Roger L. 1969. Finite dynamic model for infinite media. Journal of Engineering Mechanics Division, ASCE, 95(4): 859-878. | |
Meunier P, Hovius N, Haines J A. 2008. Topographic site effects and the location of earthquake induced landslides. Earth and Planetary Science Letters, 275(4): 221-232. | |
Rong M, Li X J, Lu T. 2007. Effect analysis of topography on the spectrum property of incident P waves. Northwestern Seismological Journal, 29(4): 297-302. | |
Spudich P, Hellweg M, Lee W H K. 1996. Directional topographic site response at Tarzana observed in aftershocks of the 1994 Northridge, California, earthquake:implications for main shock motions. Bulletin of the Seismological Society of America, 86(1B): S193-S208. | |
Tang H, Li X J, Li Y Q. 2012. Site effect of topography on ground motions of Xishan Park of Zigong City. Journal of Vibration and Shock, 31(8): 74-79. | |
Trifunac M D. 1972. Scattering of plane SH waves by a semi-cylindrical canyon. Earthquake Engineering & Structural Dynamics, 1(3): 267-281. | |
Wang H L. 1982. Effect of surface topography on the diffraction of P, SV, and Rayleigh waves. Bulletin of the Seismological Society of America, 74(4): 1167-1183. | |
Wang H Y, Xie L L. 2010. Effects of topography on ground motion in the Xishan park, Zigong city. Chinese Journal of Geophysics, 53(7): 1631-1638. DOI:10.3969/j.issn.0001-5733.2010.07.014 | |
Wang W, Liu B D, Liu Xin, et al. 2014. Characteristics and mechanism of earthquake damage at the foot of the hill in Wenchuan Earthquake. Technology for Earthquake Disaster Prevention, 9(4): 863-871. | |
Yuan X, Men F L. 1992. Scattering of plane SH waves by a semi-cylindrical hill. Earthquake Engineering & Structural Dynamics, 21(12): 1091-1098. | |
Zhang X L, Li X J, Chen G X, et al. 2016. An improved method of the calculation of equivalent nodal forces in viscous-elastic artificial boundary. Chinese Journal of Theoretical and Applied Mechanics, 48(5): 1126-1135. | |
Zhao R B, Zhang J J, Liu Z X, et al. 2015. Analysis of seismic wave problem based on visco-elastic artificial boundary under the ladder terrain. Word Earthquake Engineering, 31(3): 220-227. | |