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

赵勐，肖明†，陈俊涛，金汉城

（1.水资源与水电工程科学国家重点实验室（武汉大学），湖北 武汉 430072；2.水工岩石力学教育部重点实验室（武汉大学），湖北 武汉 430072）

为缓解巴基斯坦国家电网严重缺电的局面，中国投资并帮助巴基斯坦规划建设了一大批水利水电工程，如卡洛特水电站、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 所示），此时侧向边界区内的内行场位移（t）（i＝1，2，3）和应力（t）（i＝1，2，3）应由入射P 波、反射P 波和反射SV 波各自内行场的位移和应力叠加而成；在底部边界，内行场为P 波入射场.即：

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

图1 P 波和SV 波三维斜入射示意图Fig.1 3D diagram of oblique incidence of plane P and SV waves

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

式中：η0是一与围岩阻尼和弹性模量相关的经验常数，而与地震动的输入方式无关，张志国等[16]对各因素对η0值的影响做了大量数值试验，可参考取值；Δli为地震波i 传播至边界节点l 所经过的路程；Δt为入射P 波从零时刻波阵面传播到边界节点l 所需的时间；Δl1～Δl3、Δt1～Δt3可表示为：

侧向人工边界节点A 处的应力为：

底部人工边界节点A 处的应力为：

式中：λ、G 为介质的拉梅常数和剪切模量.由于各人工边界上y 方向的位移时程为零，故相应的y 方向应力也为零.

由式（2）～（6），可求解出P 波三维斜入射条件下人工边界处内行场的位移时程（t）、速度时程（t）及其相应的应力时程（t），以此实现三维斜入射条件下P 波荷载输入：

式中：Kli和Cli分别表示节点l（x0，y0，z0）在方向i 上的人工边界弹簧和阻尼系数，可参考文献[14]取值；Al为人工边界上节点l 的控制面积.

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 所示.

图2 三维有限元模型Fig.2 3D finite element model

图3 入射P 波位移时程Fig.3 Displacement time-history of incident P waves

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

图4 监测点A 的竖直向位移时程Fig.4 Displacement time-history of the observation point A

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

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

2.1 动接触力法的基本方程

对层状岩体层间接触系统进行有限元离散后，可得到包含动接触力的节点运动方程：

式中：M、C 和K 分别表示层间接触节点的质量、阻尼和刚度矩阵；、和u 分别表示接触节点的加速度、速度和位移向量；F 为接触节点的外荷载向量；R为接触节点的动接触力向量，且R=N+T，N 和T分别为R 的法向和切向分量.

采用中心差分法求解接触节点的运动方程，相应的时域积分格式为：

联立式（8）、式（9）和式（10），可得：

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

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

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

图5 接触模型和点对上动接触力示意图Fig.5 Contact model and dynamic contact force on node pairs

假定接触节点对在地震作用下没有相对滑动，则在t+Δt 时刻接触点对仍可认为处于黏结接触状态.此时相应的接触边界上的节点对i 和i′须满足接触边界条件，包括位移接触条件和接触力边界条件：

式中：ni为接触节点对的单位法向矢量，由节点i′指向i；ti为单位切向矢量.

联立式（11）、式（13）和式（14），可得：

式中：Mi和Mi′分别为节点i 和i′的集中质量.

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

1）若Δ1i<0 且‖‖＞cAi，表明接触面受拉，并且发生了张拉开裂破坏，接触点对i 和i′处于分离状态.此时，动接触力修正如下：

2）若Δ1i<0 且‖‖≤cAi，或者Δ1i>0 且‖‖≤μs‖‖，表明接触节点对处于黏结接触状态，此时无需对和进行修正.

3）若Δ1i>0 且‖‖≥μs‖‖，表明接触面受压，且接触点对处于滑动接触状态，即接触面发生剪切滑移破坏.此时，须校核切向动接触力，如下所示：

式中：μs和μd分别为接触面静、动摩擦因数；c 为接触面的黏聚力.若接触点对一直处于黏结接触状态，则c>0.若接触点对从黏结接触转变到滑动接触或分离状态，则c=0；Ai为接触节点i 的控制面积.由式（15）～（17）即可求得t 时刻动接触力Rt，将Rt代入式（13），计算t+Δt 时刻附加位移场Δut+Δt，更新接触节点的总位移ut+Δt.

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，与大地坐标系平行，竖直向上为正.

图6 层状岩体模型Fig.6 Interbed layered rock mass model

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

表1 材料力学参数Tab.1 Physico-mechanical parameters of materials

3.2 计算条件

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

式中：R 为材料的损伤常数；Di、分别为第i 个主应变方向上的损伤系数和塑性偏应变，为便于后处理显示材料的损伤状态，采用标量形式的损伤系数D=.

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

图7 入射波加速度时程曲线Fig.7 Time-history curve of input wave acceleration

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

图8 监测点布置Fig.8 Layout of the monitoring points

3.3 计算结果及分析

3.3.1 围岩破坏区分布

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

图9 震后围岩破坏区分布Fig.9 Distribution of surrounding rock failure zone after earthquake

当地震动竖直入射时，洞周围岩破坏区分布较少，从横向上看，洞周塑性破坏区在腰拱处有逐渐向深部扩展的趋势，且塑性区深度为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，隧洞衬砌结构相对变形量值较小.

图10 监测点合位移时程曲线Fig.10 Displacement time history curves of lining

图11 腰部与顶拱相对位移时程曲线Fig.11 Relative displacement time history curve of the haunch and top arch

工况②下，当考虑地震动的空间斜入射特性时，各监测点位移及腰、拱相对位移时程曲线相比工况①的差别主要表现在波动幅值上，各监测点处最大位移为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 所示，进一步说明了地震动斜入射对层状岩体水工隧洞接触响应的影响.

图12 层间接触面相对位移时程曲线Fig.12 Relative displacement time history curve of interface

当地震动竖直入射时，层状岩体层间相对位移在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 所示.

图13 衬砌监测点最大主应力时程曲线Fig.13 Maximum principle stress time history curves of lining

当地震动竖直入射时，在前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，且主要位于衬砌结构腰部.表明在进行层状岩体水工隧洞抗震设计时，需要采用抗断技术限制衬砌结构发生开裂破坏.

图14 衬砌结构损伤系数分布Fig.14 Damage coefficient distribution of lining structure

4 结论

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

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

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

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

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