地震动斜入射下层状岩体隧洞接触响应分析

赵勐 肖明 陈俊涛 金汉城

摘   要：為研究层状岩体对隧洞地震响应的影响，考虑地震动斜入射特性和层状岩体层间非线性接触特性，建立一种层状岩体水工隧洞地震动力响应数值模拟方法. 首先，基于三维黏弹性人工边界条件和波场分解理论，将地震动转化为作用于人工边界上的等效节点力，建立了一种层状岩体中地震动三维空间斜入射输入方法. 其次，针对地震作用下层状岩体层间动力相互作用特点，建立了一种考虑接触面黏结滑移特性的动接触力算法. 将该模拟方法应用于巴基斯坦阿扎德帕坦水电站输水隧洞抗震稳定计算，对比分析地震动竖直入射、地震动斜入射、地震动斜入射且考虑动接触3种工况的计算结果，结果表明，地震作用下隧洞结构的应力和位移响应受地震动入射角影响明显;层间剪切、挤压破碎带的存在加剧了隧洞的地震反应，接触面附近破坏区发展较大;考虑接触作用后，衬砌腰部的应力和位移响应相比顶拱较大，首先发生开裂损伤破坏，成为水工隧洞衬砌结构抗震设计的薄弱部位，隧洞结构的损伤区主要分布于软岩穿过部位和层间接触部位.

关键词：层状岩体;水工隧洞;地震动斜入射;动接触力法;地震响应;数值模拟

中图分类号：TU45                              文献标志码：A

Abstract：In order to study the influence of layered rock mass on tunnel seismic response， considering the oblique incident angles of seismic motion and the nonlinear contact characteristics at the interface， a numerical simulation method for seismic response of hydraulic tunnel in layered rock mass was proposed. First， based on the 3D viscoelastic artificial boundary conditions and the wave field decomposition theory， the input method of an obliquely incident earthquake in layered rock mass was put forward. It can transform the seismic waves into equivalent nodal forces acting on the nodes of artificial boundaries. In view of dynamic interaction characteristics between interlayers in layered rock mass under seismic action， a dynamic contact force algorithm considering the bond-slip characteristics of the interface was presented. Then， the methods were applied to the anti-seismic stability calculation of the hydraulic tunnel of AZAD PATTAN hydropower station in Pakistan. The calculation is divided into three different working conditions， with vertically incident earthquake， with obliquely incident earthquake but no dynamic contact force， with obliquely incident earthquake and dynamic contact force. The results indicate that the stress and displacement of the tunnel structure under seismic action are greatly affected by the angle of incidence. The existence of interlayer shearing and crushing fracture zone exacerbates the seismic response of the tunnel， resulting in the fact that the failure zone near the interface develops further. After considering the contact effect， the stress and displacement response of the haunch is larger than that of the vault， so that the haunch of lining where the cracking damage first occurs is the weak part of lining structure under the action of earthquakes. The damage zone of the lining mainly distributes in the place where the soft rock passes through and interlayer contacts.

Key words：layered rock mass;hydraulic tunnel;obliquely incidence earthquake;dynamic contact force method;seismic response;numerical simulation

为缓解巴基斯坦国家电网严重缺电的局面，中国投资并帮助巴基斯坦规划建设了一大批水利水电工程，如卡洛特水电站、SK水电站和正在规划设计的阿扎德帕坦水电站，从而形成了为数众多的地下水工隧洞群. 水工隧洞往往深覆于山体中，大多具有大尺度、大埋深、洞线长等特点，不可避免地要穿越包括强震区和层状岩体区等在内的复杂地质区域，面临突出的抗震安全稳定问题[1]. 在层状岩体中进行地下工程建设，尤其对于软硬互层状岩体而言，当岩层倾角较大，结构面发育明显时，层间剪切、挤压破碎带较为常见，是水工隧洞抗震的薄弱部分. 一旦发生地震，岩体层间极易发生剪切滑移破坏，对隧洞结构造成不可逆的损伤破坏. “5.12”汶川大地震震害调查表明，深埋于地质条件较差部位的隧洞结构极易发生衬砌开裂及错位等损伤破坏[2]. 因此，研究层状岩体中水工隧洞地震响应特性和破坏机理具有重要的现实意义.

层状岩体水工隧洞地震响应分析主要包括两方面内容：1） 地震动的输入方法;2） 层间动力相互作用模拟. 在已有的针对地下隧洞群进行的动力时程分析中，多是假设地震动为从模型底部竖直向上入射. 但是根据近年来强震动观测记录的统计，发现基岩场地的地震波入射角平均为60°，从而引起结构的非一致性变形[3]. 杜修力等[4]研究了地震波斜入射条件下隧洞洞身段的地震響应特征，结果表明地震波斜入射时隧洞地震响应规律与竖直入射时明显不同. 李山有等[5]研究了地震波斜入射条件下竖直、倾斜台阶地形引起的波形转换，分析说明了研究斜入射的必要性. Heymsfield[6]分析了二维条件下斜入射SH波对倾斜基岩自由面位移的放大效应. Stamos等[7]采用一种新的频域内边界元法研究斜入射体波作用下长大隧道的地震响应. Naggar等[8]研究了地震波斜入射角度对隧道衬砌弯矩和轴力的影响. 从上述文献可看出，目前对地震动斜入射的研究已取得一些成果，但均未考虑地震动在层状岩体中传播时其幅值受岩体阻尼影响而衰减这一情况.

层状岩体层间动力相互作用属典型的动接触问题. 在有限元分析中，进行动接触迭代的数值模拟方法很多，主要有Lagrange乘子法[9]、罚函数法[10]、线性补偿法[11]和动接触力法[12]等. 其中，刘晶波等[12]提出的动接触力法以其计算效率高和稳定性好而被广泛运用，但是该方法忽略了接触介质界面的黏聚力. 本文针对地震荷载作用下层状岩体层间循环往复相互作用特点，建立一种考虑界面黏结滑移特性的动接触力算法，可以反映动力作用下层间非线性接触特性.

综上所述，本文建立了一种层状岩体水工隧洞地震动力响应数值模拟方法，该方法有效考虑了地震动在层状岩体中的斜入射特性和层状岩体层间非线性接触特性. 将本文分析方法应用于巴基斯坦阿扎德帕坦水电站输水隧洞抗震稳定计算，分析软硬互层状岩体对隧洞地震响应的影响，以期对复杂层状岩体中水工隧洞抗震设计进行有益探索.

1   层状岩体地震动斜入射方法

水工隧洞相比地表、地基和城市地下建筑物的特殊性在于其完全埋置于无限域的山体介质之中. 在对深埋水工隧洞进行动力时程分析时，由于其所在的工程区相对整个地震区域是微小的，远域地震的波动对工程区域影响较小，因此可将工程区外无限域地震区假设为弹性无限介质体，从中截取隧洞所在的有限区域进行有限化模拟. 在时域分析中首先需对计算模型各边界设置人工边界，以模拟外行波的透射、内行波的入射、边界处无限域波动场及弹性位移场，且需考虑地表自由面反射波对工程区的影响.

基于三维黏弹性人工边界[13]，可将无限域地震波动场问题转化为求解作用于人工边界节点的等效节点力问题，以实现模型内外波动场的交互. 本文基于波场分解原理，假定计算模型外无限域为均匀弹性介质体，地震波为倾斜入射的弹性平面波，将波动场分解为内行场和外行场，外行场主要由自模型内部向无限域透射的外行波动场构成，可由模型内显式有限元逐步积分计算求得，故地下隧洞地震动斜入射实现的关键是求解无限域的斜向内行场以及相应的地表反射内行场.

对于空间任意入射角度的地震波，可分解为质点振动方向与波的传播方向相一致的P波（压缩波）和质点振动方向与波的传播方向相垂直的S波（剪切波）. 下面详细阐述P波三维斜入射下地震荷载的计算方法.

1.1   P波三维斜入射

平面P波在半空间自由表面经反射后会发生波形转换，产生反射P波和反射SV波（如图1所示），此时侧向边界区内的内行场位移uRli（t）（i=1，2，3）和应力σRli（t）（i=1，2，3）应由入射P波、反射P波和反射SV波各自内行场的位移和应力叠加而成;在底部边界，内行场为P波入射场. 即：

假设入射P波在零时刻的位移时程为u0（t），入射波零时刻波阵面与水平面夹角为入射角α，节点 l（x0，y0，z0）为模型人工边界上某一节点，L为模型底部至自由面高度. 需注意本文研究基于入射波波前平行于隧洞轴线，故根据入射P波波前和人工边界节点l的几何空间位置关系，可求得模型各人工边界处内行场位移时程：

式中：β1为反射SV波在半空间自由面的反射角，β1 = arcsin（cs sin α/cp），cp和cs分别为P波和SV波的波速;A1和A2分别为反射P波和反射SV波的幅值放大系数，其值参考文献[14]取得;η（·）为考虑岩体阻尼情况下，地震动的幅值沿传播距离的衰减系数，胡进军等[15]研究表明地震动在复杂岩层中传播时可认为其幅值受岩体阻尼的影响而呈线性衰减，张志国等[16]进一步将其表述为：

SV波和SH波三维斜入射的计算公式可根据上述P波斜入射条件下地震荷载的计算公式进行类似推导，不再赘述. 由此本文建立了一种与三维黏弹性人工边界相适应的层状岩体中地震动三维斜入射的输入方法.

1.2   算例验证

为验证本文层状岩体中地震动斜入射输入方法的合理性，建立了一个有限元模型来分析斜入射条件下半无限域三维弹性介质体的动力响应问题. 三维有限元模型的尺寸为800 m×800 m×800 m，共计46 656个六面体单元，如图2所示. 模型弹性模量为10 GPa，泊松比为0.3，密度为2 000 kg/m3. 计算模型底部和側向边界施加三维黏弹性人工边界，顶部为自由面，取顶部自由面中点A（400，400，800）为监测点. 入射P波的位移时程曲线如图3所示.

图4所示为当平面P波入射角分别为15°和30°时监测点A的竖直向（z向）位移时程的理论解与数值解. 需要说明的是本文入射角α是指入射波零时刻波阵面与水平面夹角，故P波入射角分别为15°和30°时波阵面的法向向量分别为（0.259，0，0.966）和（0.5，0，0.866）. 从图中可看出不同入射角下监测点A的竖直向位移时程的计算结果与理论值比较符合.

2   层状岩体层间动接触系统分析模型

在构造应力和地震荷载联合作用下层状岩体层间易发生相对错动，从而对隧洞结构造成严重的损伤破坏. 根据巴基斯坦阿扎德帕坦水电站现场观测资料，工程区层间剪切、挤压破碎带较为常见. 本文针对层状岩体层间循环往复作用特点，建立一种考虑界面黏结与滑移特性的层状岩体接触系统动力响应分析模型，用来模拟层间非线性滑移破坏.

2.1   动接触力法的基本方程

由式（11）～（13）可看出，由t时刻接触节点的运动状态和动接触力，可求解出t + Δt时刻接触节点的运动状态. t时刻接触节点的运动状态是已知的，而动接触力Rt 是未知量，需根据t ～ t + Δt时刻的接触状态基于相应的接触条件计算.

2.2   考虑界面黏结滑移特性的动接触力算法

假定在地震作用前层状岩体层间接触良好，考虑界面黏聚力，则接触点对处于黏结接触状态，如图5所示. 在强震过程中，层间可能会发生相对滑动，导致接触节点与其对应的单元某一面发生接触，接触面进入滑动接触状态，此时不再考虑黏聚力.

显然，式（15）中法向和切向动接触力是在层状岩体层间处于黏结状态下求得的，而层间接触面是水工隧洞中薄弱部分，其在地震循环作用下的损伤是不可忽略的. 实际上，在动力作用下，层状岩体层间接触面存在黏结接触、滑动接触和分离等多种接触状态，因此在每一时步计算完毕后，需要对接触节点对的接触状态进行判别，并对动接触力进行修正[17].  层间接触面的破坏形式主要包括沿切向的剪切滑移和沿法向的张拉开裂，具体方法如下：

3 工程实例

3.1   工程概况和计算模型

阿扎德帕坦水电站位于巴基斯坦的Jhelum河上，为该河段水电开发中的一级，以发电为主. 引水发电系统位于河流左岸，导流隧洞布置在右岸. 导流隧洞工程区地质条件复杂，隧洞所穿越的基岩岩性为砂岩与非砂岩类呈互层状分布，地层主要为单斜构造，岩层倾角∠71°～∠81°，层间剪切、挤压破碎带较为常见.

阿扎德帕坦工程区属于地震活动区，主要受印度板块持续向欧亚板块俯冲运动影响. 巴基斯坦北部和阿扎德地区为地震强烈活动地区，受多个地震板块构造影响，地震活动多由本地区活动断层运动引起，阿扎德帕坦水电站工程区就在这一区域. 时间最近的大地震为2005年10月8日的7.6级地震. 根据中国地震局地质研究所研究成果，阿扎德帕坦水电站工程区50年超越概率10%（DBE）的峰值加速度为0.315g，对应的地震基本烈度为Ⅷ度. 2#导流隧洞洞长667 m，洞身段埋深在140 m左右. 采用圆形断面结构，开挖洞径10.0 m，围岩以Ⅲ、Ⅳ类为主，衬砌采用C25钢筋混凝土结构，厚度为60 cm.

选取如图6所示洞身段含泥、砂岩互层的隧洞区域建立水工隧洞三维有限元模型，考虑到若模型建至地表，则单元数量过多，动力计算耗时将呈指数级增长. 为了提高计算效率，隧洞顶部取50 m. 层状岩层走向与洞轴线垂直，倾角为70°. 模型共剖分了48 608个八节点六面体单元和52 041个节点，其中混凝土衬砌单元3 888个. 模型范围及坐标系：x方向从-60.0 m到60.0 m，与洞轴线垂直;y方向从 -85.0 m到85.0 m，与洞轴线重合，顺水流为正;z方向从400.0 m到520.0 m，与大地坐标系平行，竖直向上为正.

三维初始地应力场根据设计院提供的实测地应力反演分析得到，侧压力系数取kx = 1.1，ky = 0.85，kz = 1.0. 层状岩体砂岩、泥岩、接触面和衬砌的材料力学参数取值见表1. 围岩临界阻尼比取5%，则参考文献[16]砂岩层的地震动幅值衰减系数取0.04%，非砂岩层的地震动幅值衰减系数取0.07%. 动力计算之前，采用三维弹塑性损伤有限元法进行隧洞的静力开挖与支护计算，其相应的计算结果作为动力计算的初始条件.

3.2   计算条件

计算程序采用课题组自主开发的大型地下洞群抗震稳定动力时程分析平台[18]，并将本文地震动输入方法和动接触力算法嵌入其中. 接触面模拟的部位分别位于非砂岩层（即软岩）与砂岩层（即硬岩）相交界面处. 动力加载前，首先基于节点分离技术[19]，通过增加砂岩层与非砂岩层接触面两侧的共用节点，并设置一一对应的接触节点对，以完成软岩与硬岩单元的节点分离. 围岩和衬砌采用基于M-C屈服准则的动力弹塑性损伤本构模型[18]，三维损伤演化方程如下所示：

模型的底部、四周和顶部均采用黏弹性人工边界，以吸收斜入射条件下的地震波及其在地表自由面的反射波. 地震波采用美国强震记录的El-Centro波，并根据阿扎德帕坦水电站工程区抗震设防烈度，将峰值加速度调整为3.15 m/s2，截取其中变化剧烈、幅值较大的20 s时段作为入射波，经滤波和基线校正处理后加速度时程曲线如图7所示. 计算时考虑斜入射（∠30°入射角）和竖直入射的地震动对围岩和衬砌的作用. 动力计算同时考虑P波和SV波对水工隧洞的作用. 其中，SV波采用如图7所示的入射波，P波加速度时程取为SV波的2/3[20].

选取软岩中间断面为监测断面，布置如图8所示的监测方案，监测点A、B、C分别位于监测断面上衬砌的顶拱、左腰部和底部几个关键部位，用以监测地震加载过程中衬砌的位移和应力等指标特性，监测点D、E分别位于层状岩体硬岩与软岩接触面两侧，用以监测层状岩体层间相对运动特征. 动力计算分3种工况：①地震动竖直入射，不考虑动接触;②地震动斜入射，不考虑动接触;③地震动斜入射，考虑动接触. 需要注意的是，工况①地震动竖直入射是指地震动自模型底部垂直入射，水平面内沿垂直水流向（x向）振动，竖直向振动取水平向振动的2/3. 地震动斜入射的入射方向矢量为（0.5，0，0.866）.

3.3   计算结果及分析

3.3.1   围岩破坏区分布

在地震循环荷载作用下，水工隧洞洞周围岩循环加卸载，围岩应力一直处于波动状态，且塑性变形不断累积，使得围岩总应变逐渐增加，导致围岩总破坏区体积相应增大. 3种工况下震后洞周围岩破坏区分布如图9所示.

当地震动竖直入射时，洞周围岩破坏区分布较少，从横向上看，洞周塑性破壞区在腰拱处有逐渐向深部扩展的趋势，且塑性区深度为2.52 m，开裂区分布较少，深度为0.97 m;从纵向上看，开裂区主要分布在软岩穿过的区域，其余处较少，仅在腰拱处出现. 考虑地震动斜入射时，软岩及层间破碎带处的塑性破坏区显著变大，而开裂区增大不明显，塑性区深度达4.96 m，开裂区深度为1.65 m，表明考虑地震动在层状岩体中的传播特性后，斜入射地震动极大削弱了洞周围岩的稳定性，有可能导致隧洞结构的局部失稳和损伤破坏. 考虑地震动的斜入射和层间动接触相互作用时，层间破坏区进一步扩展，主要表现为塑性区和开裂区均明显增大. 与仅考虑地震动斜入射工况相比，塑性区深度增加了2.15 m，开裂区深度增加了0.86 m，且软岩穿越的洞周几乎被开裂区包围. 从图9中还可看出震后开裂区主要分布在隧洞的腰拱及其上部，这与Wang等[21]的研究结果相符.

3.3.2   衬砌结构位移时程分析

3种工况下衬砌结构不同监测点处合位移时程如图10所示. 由图10可看出3种工况下，1）顶拱、腰部和底部位移时程曲线的波形和波动规律基本一致，均出现了多个明显的波峰. 监测点位移同时出现波峰和波谷，表明水工隧洞衬砌结构各部位处于同步震动状态;2）在0～5 s时间段内，衬砌结构各部位位移时程曲线呈现大幅度波动，腰部最大位移要比顶拱和底部大. 本文采用腰部与顶拱的合位移差值来表征相对位移，相对位移可以更好地表征隧洞衬砌结构的变形特征，图11所示为3种工况下腰部及顶拱位移动力响应差值时程曲线.

工况①下，衬砌结构各监测点处最大位移为7.4 cm，腰部与顶拱最大相对位移为0.67 cm，发生在4.95 s，但震后腰、拱相对位移为0.21 cm，隧洞衬砌结构相对变形量值较小.

工况②下，当考虑地震动的空间斜入射特性时，各监测点位移及腰、拱相对位移时程曲线相比工况①的差别主要表现在波动幅值上，各监测点处最大位移为8.7 cm，腰、拱最大相对位移为0.85 cm，发生在5.5 s，震后相对位移为0.52 cm，表明衬砌变形受地震动入射角影响较大，地震动的三维斜入射特性和入射边界的非一致特性对隧洞衬砌结构位移响应影响较大. 这主要是因为考虑地震动斜入射时，地震波场与竖直入射时具有明显差异. 地震动斜入射时，入射P波和SV波会在自由面发生波形转换，各自分别形成反射P波和反射SV波，故模型人工边界处波场是由不同入射波和反射波叠加而成的，使得人工边界上各节点具有不同的振动波形，产生了放大效应. 而地震动从模型底部竖直入射时，经地表自由面反射后地震波仍具有相同的反射方向和振动幅值，模型人工边界上节点也具有相同的振动方向和振动波形. 因此，斜入射时衬砌结构各部位位移响应相比竖直入射时要大.

工况③下，当地震动斜入射且考虑层间动力相互作用后，各监测点处最大位移为9.8 cm，腰、拱最大相对位移达到1.41 cm，发生在5.15 s，在前期地震动波动较为剧烈的0～6 s时间段内，其腰、拱相对位移时程曲线波动幅度相比前2个工况较大，波动范围为-0.55～1.41 cm，震后相对位移为1.04 cm. 从理论上说，衬砌结构的相对位移在一次地震结束后应当回到0[3]. 然而，工况③中衬砌腰部和顶拱的相对位移值在震后为1.04 cm，表明隧洞结构在地震荷载作用下发生永久变形. 可见在构造应力和地震荷载联合作用下层状岩体层间易发生剪切滑移破坏，且斜入射地震动加剧了衬砌腰部的相对变形，致使衬砌腰部抗剪段安全问题突出.

3.3.3   层状岩体层间相对运动分析

在层状岩体层间接触系统地震动响应过程中，硬岩与软岩层间循环往复作用，发生了复杂的动接触行为，如：黏结接触、滑动接触和分离等多种接触状态，进而产生了层间错动位移，对隧洞结构造成严重的损伤破坏. 3种工况下层状岩体层间相对位移时程曲线如图12所示，进一步说明了地震动斜入射对层状岩体水工隧洞接触响应的影响.

当地震动竖直入射时，层状岩体层间相对位移在0线附近波动变化，层间最大相对位移为1.77 cm，在0～6 s时间段内，层间相对位移波动幅度较大，主要在-1.60～1.77 cm范围内上下波动，后期波动幅度逐渐减小，震后基本为0. 当考虑地震动斜入射后，在0～6 s时间段内，层间相对位移在0线上下剧烈波动，波动范围为-1.96～2.16 cm，层间最大相对位移为2.16 cm，发生在5.75 s，受斜入射下地震动输入非一致性的影响，后期波动幅度虽较0～6 s时间段有所减小，但相比工况①仍较大，震后逐渐减小到0线附近. 当考虑层间动接触力后，在前期0～1.5 s时间段内，层间相对位移在0线附近波动极小，表明此时间段内层状岩体层间接触面接触良好，处于黏结接触状态，或从黏结接触状态向滑动接触或分离接触状态转变的过渡阶段;在1.5～6 s时间段内，层状岩体层间相对位移波动较为明显，且出现了明显的错动位移，最大错动位移为-3.91 cm，发生在5.2 s;后期层间相对位移主要在-3.00～-4.00 cm范围内上下波动，震后为-3.94 cm，表明考虑层间动接触力后，层状岩体层间发生了明显的剪切滑移破坏，且受地震动斜入射的影响，层间相对位移在-3.50 cm上下波动明显.

3.3.4   衬砌结构应力时程分析

由于混凝土的抗压强度远大于其抗拉强度，地震作用下水工隧洞衬砌结构的损伤破坏主要是拉裂破坏. 因此本文主要分析地震作用下衬砌结构的最大主应力变化规律. 3种工况下衬砌结构不同监测点处最大主应力时程曲线如图13所示.

当地震动竖直入射时，在前5 s内，顶拱、腰部和底部的最大主應力变化剧烈，波动范围主要为-0.50 ～ 1.47 MPa，在峰值处腰部、顶拱和底部的拉应力量值分别达到1.47 MPa、1.20 MPa、0.96 MPa，腰部的拉应力量值超过了混凝土的抗拉强度，且腰部的应力水平明显大于顶拱和底部的应力水平.

当地震动斜入射时，受自由面反射波场叠加效应的影响，整体上看隧洞结构应力水平的地震反应大于竖直入射时的地震反应. 斜入射时衬砌腰部最大主应力相比顶拱和底部一直波动较为剧烈，在峰值处衬砌腰部、顶拱和底部的拉应力量值分别为1.61 MPa、1.27 MPa、1.23 MPa，其中腰部和顶拱最大拉应力达到了混凝土的抗拉强度，表明衬砌的腰部为衬砌结构受力的不利部位.

当地震动斜入射且考虑层间动接触力后，斜入射的地震动加剧了层状岩体硬岩和软岩间循环往复相互作用，腰部最大主应力波动范围增大，主要在0.5～1.81 MPa，在峰值处衬砌腰部、顶拱和底部的拉应力量值分别为1.90 MPa、1.69 MPa、1.45 MPa，均超过了混凝土的抗拉强度. 表明在考虑层间动接触力后，腰部的应力响应更为明显，层状岩体层间相对滑移对衬砌结构受力具有重要影响，加剧了隧洞腰部拉应力的变化，使得衬砌结构发生拉裂损伤破坏.

3.3.5   衬砌结构损伤分析

图14所示为考虑地震动斜入射和层间动接触相互作用后，衬砌结构震后损伤系数分布图. 由图可知，在地震动输入完成后，衬砌损伤区主要分布在层状岩体层间接触部位和软岩穿过部位，且向两侧延伸约5 m范围内，最大损伤系数接近于1，且主要位于衬砌结构腰部. 表明在进行层状岩体水工隧洞抗震设计时，需要采用抗断技术限制衬砌结构发生开裂破坏.

4   结   论

基于地震动三维空间斜入射输入方法和考虑层状岩体层间黏结滑移特性的动接触力算法，建立了一种层状岩体水工隧洞地震动力响应分析模型. 结合阿扎德帕坦水电站水工隧洞工程实例，对地震作用下衬砌结构动力响应及损伤破坏进行分析，得到如下结论：

1）层状岩体中软岩穿过部位的破碎带加剧了隧洞结构的地震响应，表现为围岩塑性破坏区及开裂区的显著增大. 斜入射地震动削弱了水工隧洞衬砌结构的稳定性，在考虑层间动力相互作用后接触面附近破坏区进一步发展且软岩穿越的洞周几乎被开裂区包围.

2）地震作用下隧洞衬砌结构不同部位监测点的位移时程曲线与输入地震动位移时程曲线相类似，表现为一种同步震动趋势. 对比地震动斜入射和竖直入射，隧洞衬砌结构的变形响应受地震动入射角影响较大. 当同时考虑地震动斜入射和层间动接触力时，层状岩体层间易表现为剪切滑移破坏.

3）与地震动竖直入射相比，斜入射时隧洞衬砌结构的应力响应更大. 当考虑层间动接触作用后，接触面附近处衬砌应力进一步加大，衬砌腰部应力响应相比顶拱较大，因而会首先发生开裂破坏，故腰部是水工隧洞衬砌结构抗震设计的薄弱部位.

4）地震作用下层状岩体层间震动不同步，极易发生相互错动，进而使衬砌损伤破坏. 衬砌损伤区主要分布于层间接触部位和软岩穿过部位，且向两侧延伸约5 m.

参考文献

[1]    HUANG J Q，ZHAO M，DU X L. Non-linear seismic responses of tunnels within normal fault ground under obliquely incident P waves[J]. Tunnelling and Underground Space Technology，2017，61：26—39.

[2]    SHEN Y S，GAO B，YANG X M，et al. Seismic damage mechanism and dynamic deformation characteristic analysis of mountain tunnel after Wenchuan earthquake[J]. Engineering Geology，2014，180：85—98.

[3]    JIN X，LIAO Z P. Statistical research on S-wave incident angle[J]. Earthquake Research in China，1994，8（1）：124—134.

[4]    杜修力，黄景琦，赵密，等. SV波斜入射对岩体隧道洞身段地震响应影响研究[J]. 岩土工程学报，2014，36（8）：1400—1406.

DU X L，HUANG J Q，ZHAO M，et al. Effect of oblique incidence of SV waves on seismic response of portal sections of rock tunnels[J]. Chinese Journal of Geotechnical Engineering，2014，36（8）：1400—1406. （In Chinese）

[5]    李山有，廖振鹏. 地震体波斜入射情形下台阶地形引起的波型转换[J]. 地震工程与工程振动，2002，22（4）：9—15.

LI S Y，LIAO Z P. Wave-type conversion caused by a step topography subjected to inclined seismic body wave[J]. Earthquake Engineering and Engineering Vibration，2002，22（4）：9—15. （In Chinese）

[6]    HEYMSFIELD E. Two-dimensional scattering of SH waves in a soil layer underlain with a sloping bedrock[J]. Soil Dynamics and Earthquake Engineering，2000，19（7）：489—500.

[7]    STAMOS A A，BESKOS D E. 3-D seismic response analysis of long lined tunnels in half-space[J]. Soil Dynamics and Earthquake Engineering，1996，15（2）：111—118.

[8]  EL NAGGAR H，HINCHBERGER S D，EL NAGGAR M H.Simplified analysis of seismic in-plane stresses in composite and jointed tunnel linings[J].Soil Dynamics and Earthquake Engineering，2008，28（12）：1063—1077.

[9]    BATHE K J，BOUZINOV P A. On the constraint function method for contact problems[J].Computers & Structures，1997，64（5/6）：1069—1085.

[10]  PERI？譎 D，OWEN D R J. Computational model for 3-D contact problems with friction based on the penalty method[J]. International Journal for Numerical Methods in Engineering，1992，35（6）：1289—1309.

[11]  KWAK B M，LEE S S. A complementarity problem formulation for two-dimensional frictional contact problems[J]. Computers & Structures，1988，28（4）：469—480.

[12]  LIU J B，SHARAN S K. Analysis of dynamic contact of cracks in viscoelastic media[J].Computer Methods in Applied Mechanics and Engineering，1995，121（1/2/3/4）：187—200.

[13]  劉晶波，王振宇，杜修力，等. 波动问题中的三维时域黏弹性人工边界[J]. 工程力学，2005，22（6）：46—51.

LIU J B，WANG Z Y，DU X L，et al. Three-dimensional visco-elastic artificial boundaries in time domain for wave motion problems[J]. Engineering Mechanics，2005，22（6）：46—51. （In Chinese）

[14]  杜修力，赵密，王进廷. 近场波动模拟的人工应力边界条件[J].力学学报，2006，38（1）：49—56.

DU X L，ZHAO M，WANG J T. A stress artificial boundary in fea for near-field wave problem[J]. Chinese Journal of Theoretical and Applied Mechanics，2006，38（1）：49—56. （In Chinese）

[15]  胡进军，谢礼立. 地下地震动频谱特点研究[J]. 地震工程与工程振动，2004，24（6）：1—8.

HU J J，XIE L L. Spectral characteristics of earthquake sub-ground motions [J]. Earthquake Engineering and Engineering Vibration，2004，24（6）：1—8. （In Chinese）

[16]  张志国. 地下洞室群地震响应数值分析方法研究 [D]. 武汉：武汉大学，2012：86—87.

ZHANG Z G. Study on numerical simulation methods for seismic response of underground cavern complexes[D]. Wuhan：Wuhan University，2012：86—87. （In Chinese）

[17]  陳世杰，肖明，陈俊涛. 隧洞块体破坏过程及稳定评价的数值方法研究[J]. 湖南大学学报（自然科学版），2020，47（5）：31—38.

CHEN S J，XIAO M，CHEN J T. Study on numerical method for failure process and stability evaluation of rock block in tunnel[J]. Journal of Hunan University （Natural Sciences），2020，47（5）：31—38. （In Chinese）

[18]  张志国，肖明，陈俊涛. 大型地下洞室地震灾变过程三维动力有限元模拟[J]. 岩石力学与工程学报，2011，30（3）：509—523.

ZHANG Z G，XIAO M，CHEN J T. Simulation of earthquake disaster process of large-scale underground caverns using three-dimensional dynamic finite element method[J]. Chinese Journal of Rock Mechanics and Engineering，2011，30（3）：509—523. （In Chinese）

[19]  赵健，肖明，陈俊涛，等. 基于单元重构与节点分离的大型地下洞室软弱结构面模拟方法[J]. 湖南大学学报（自然科学版），2017，44（3）：134—142.

ZHAO J，XIAO M，CHEN J T，et al. Simulation methodology of weak structural planes in large underground chamber based on element reconstruction and node separation[J]. Journal of Hunan University （Natural Sciences），2017，44（3）：134—142. （In Chinese）

[20]  水工建筑物抗震设计标准：GB 51247—2018[S]. 北京：中国计划出版社，2018：16—17.

Standard for seismic design of hydraulic structures：GB 51247—2018[S]. Beijing：China Planning Press，2018：16—17. （In Chinese）

[21]  WANG Z Z，GAO B，JIANG Y J，et al. Investigation and assessment on mountain tunnels and geotechnical damage after the Wenchuan earthquake [J]. Science in China Series E：Technological Sciences，2009，52（2）：546—558.