张懂懂， 刘 洋， 熊 峰， 梅 竹， 李 爽

(1. 四川大学 建筑与环境学院，成都 610065； 2. 温州大学 建筑工程学院，浙江 温州 325035；3. 哈尔滨工业大学 土木工程学院，哈尔滨 150090)

山岭地区岩体中隧道受到围岩介质的约束，地震效应小于地上结构，一般认为不会出现严重的地震破坏情况。然而，美国北岭地震、日本神户地震、中国台湾集集地震与汶川地震等一系列强震[1-4]中，近场岩体隧道均出现严重震害，如衬砌裂损、拱起、洞室垮塌等，对隧道结构的抗震研究才逐渐活跃起来。目前关于隧道结构抗震的研究大多假定地震波垂直入射，而震害调查表明，当震源埋深较浅时，近场强震地震波往往以一定角度斜入射，进而引发隧道衬砌结构破坏，故有必要加以考虑。

近年，关于隧道等地下结构在地震波斜入射下的地震响应研究有了一些进展。刘晶波等[5]将地震动输入问题转化为波源问题，将地震动转化为等效荷载可处理波动斜输入。杜修力等[6]初步探索地震波入射角度对在建地铁的地震响应，发现斜入射时的结构地震响应规律与垂直入射时具有较大的差异。黄景琦等[7]研究了SV波斜入射对地铁车站地震响应的影响，结果表明地铁车站的地震响应随入射角度的增加而增加。周晓洁等[8]研究了SV波斜入射下矩形海底沉管隧道的地震响应规律，提出沉管隧道结构抗震设计时应当考虑斜入射的影响。Lyu等[9]研究了SV波斜入射对拱形隧道非线性地震响应的影响，结果表明隧道动力反应随入射角度增大而增大。巴振宁等[10]研究了平面P-SV波入射下山岭隧洞的动力响应，得出入射波频率及入射角度均会显著影响隧洞动力响应的结论。王飞等[11]分析了近断层脉冲型与非脉冲型地震动斜入射对水电站厂房地震响应的影响规律，认为脉冲型地震动斜输入比非脉冲型地震动斜输入对电站厂房破坏作用更大。以上研究多集中在剪切SV波斜入射的情况，考虑压缩P波作用的研究不足。尽管P波在土体中衰减很快，对远场地动力反应的影响微弱，可以只考虑SV波斜入射的情况。然而，对于容易出现隧道震害的近场地震情况下，P波产生的拉压破坏作用不容忽略[12]，有必要开展P波斜入射作用下的隧道动力响应研究。

地下结构地震响应分析一般涉及多种静力与动力荷载的共同作用，故采用数值仿真手段进行地下结构地震反应分析时，为模拟静、动力荷载下的耦合效应，通常需要对结构进行静动力联合分析。既有研究[13-14]大多通过设置人工边界条件转换可进行静动力联合分析，然而已有的边界转换方法仅进行一次静力分析获得初始应力状态，忽略了施工影响的场地初始位移状态，导致模型的位移场在人工边界转换前后不一致，从而影响隧道地震响应分析的准确性。

综上，本文采用结合黏弹性边界的有限元方法，首先基于ABAQUS软件建立地震动斜入射通用数值模型，再以一种可将静荷载下模型的应力、位移场同时作为动力分析初始状态的改进的耦合边界模拟方法，实现地下结构静动力联合分析，并经算例验证模型和方法的可靠性，最后以岩体隧道洞口段二维平面模型为切入点，对比分析P波与SV波分别斜入射对岩体隧道洞口段地震效应的影响规律。

1 地震波斜入射数值模型

1.1 动力人工边界

采用有限元法对地下结构地震反应进行动力分析时，为了兼顾数值计算成本与精度，需要从半无限域地基中截取适当范围的近场区域进行离散。近场结构产生的外行散射波传播至人工截断处会产生波动反射，造成计算结果的失真，故需要在截断边界处设置动力人工边界[15-16]以模拟半无限域地基的辐射阻尼。本文采用杜修力等[17]提出的黏弹性边界模拟远域地基。

二维黏弹性人工边界的弹簧和阻尼器系数为

法向

(1)

切向

(2)

式中：ρ为介质密度；cp,cs分别为介质的剪切波速和压缩波速；λ，G分别为拉梅常数和剪切模量；r取近场结构几何中心到节点所在边界面的垂直距离；A,B为无量纲参数，按较优建议取值为A=0.8，B=1.1。

1.2 地震波输入模拟

隧道地震反应研究属于土-结构动力相互作用问题，现有对于土-结构动力相互作用系统的地震动输入方法主要分两种，振动法与波动法。振动输入方法假定下卧基岩刚性，一般通过在计算模型底部以加速度时程、速度时程或位移时程的形式施加惯性力作用，来近似实现地震波均匀入射，多应用于地面结构的抗震分析，由于该方法不处理侧边界节点，忽略了地震波在模型侧边界的振动传播过程，难以考虑地震波斜入射引起的非一致输入效应。波动输入方法是一种理论上更为严密的输入方法，利用波场分解技术[18]将总波场分为自由场和散射场，再根据力的叠加原理将入射地震波转化为作用在各边界面上的等效荷载，使人工边界处位移、应力都与自由场一致来实现地震波动输入，较好地模拟了地震波在人工边界面上的自由场效应，结合吸收外行散射波的动力人工边界即可用于处理地震波斜入射问题。

基于黏弹性边界的地震波动输入方法中，人工边界处i节点力j方向的等效地震荷载计算公式如下

Fij=(Kijuij+Cijvij+σij)Ai

(3)

式中：Kij,Cij为边界i节点j方向的弹簧刚度和阻尼器系数；uij,vij为i结点j方向的自由场位移和速度向量；σij为i结点j方向的自由场应力向量；Ai为人工边界i节点的控制面积。

根据式(3)可知，等效荷载计算的计算需先求得已知地震波下各边界处的自由场反应(位移、速度和应力场)。根据波动理论，斜入射P波与SV波传播至表面时会发生波形转换，产生反射P波和反射SV波，如图1所示，故自由波场可分解为入射波场与反射波场。在不同边界面上构建自由场时选择不同的波场组合会求解出不同的等效荷载，影响最终的计算结果。已有研究[19]表明，构建边界自由场时，考虑的波动叠加越完整，计算精度越高。故各边界处自由波场反应采用一条入射波和两条反射波叠加的结果，可利用解析或者数值方法求解得到。

图1 P波与SV波地表波形转换示意Fig.1 Schematic diagram of P wave and SV wave waveform conversion on the surface

需要注意的是，由图1(b)可知，对于SV波斜入射的情形，当反射P波的反射角βrp=π/2时，入射角αis达到了临界值，根据斯奈尔定理可知

(4)

故对于SV波入射，本文仅考虑其临界入射角范围以内的情况。由于人工边界处节点通常数量众多，手动操作定义黏弹性边界、等效荷载时程的过程复杂费时，故借助ABAQUS前处理二次开发平台编制黏弹性边界自动施加及等效节点力计算辅助程序，实现批量建立弹簧-阻尼器及地震荷载的快速添加。

1.3 自由场验证

为了验证1.2节平面波斜入射方法的模拟精度及开发程序的可靠性。在均匀半空间中截取1 000 m×500 m的矩形区域，建立二维自由场有限元模型如图2所示，并设置模型顶部中点A和模型中点B为监测点。假设半空间为线弹性介质，材料参数为：密度1 800 kg/m3，弹性模量1 GPa，泊松比为0.25，剪切波速527.05 m/s，压缩波速860.66 m/s。采用满足波动有限元模拟精度要求的边长5 m的平面应变单元(CPE4)进行离散。

图2 自由场有限元模型Fig.2 Finite element model

入射P波和SV波都为狄克拉单脉冲波，位移时程如图3所示，脉冲时间跨度0.5 s，峰值1 m。入射角度都为α=20°，记录监测点的位移时程曲线，与理论值对比如图4、图5所示。

图3 入射波位移曲线Fig.3 Time history curve of incident wave displacement

图4 P波斜入射下监测点处位移时程曲线Fig.4 Time history curve of displacement at monitoring points

图5 SV波斜入射下监测点处位移时程曲线Fig.5 Time history curve of displacement at monitoring points

从图4、图5可以看出，P波和SV波下的计算结果都与基于波动理论计算的理论解波形吻合，虽幅值略有不同，但误差不超过5%，满足计算精度要求，表明了本文采用的波动输入方法和程序的准确性。

2 静-动力边界转换

由于地震等动力荷载作用时，地下结构已在结构原有自体质量、土压力等静荷载作用下具有初始应力，为计算趋于真实情况下的地震响应，模拟复杂荷载下结构的静、动力耦合效应，通常需要进行地下结构的静动力联合分析。目前在线弹性小变形范围内，一般是利用叠加原理将地下结构单独在静力荷载或动力荷载下的各项计算结果进行线性叠加，得到结构的最终反应值；对于涉及非线性效应的分析问题，叠加原理不宜使用，需要一定的处理方法进行静动力联合分析。

针对地下结构-围岩系统的静、动力耦合效应模拟技术主要有3种：第1种是静、动力统一人工边界，如黏弹性静-动力统一人工边界、无限元静-动力统一人工边界；第2种是动力松弛法，基于静力解是外荷载下瞬态响应的稳态部分这一事实，采用一个有阻尼衰减的动力分析过程求解结构系统在重力等外载下的静力问题，当模型趋于平衡后再进行地震荷载的动力分析；第3种是静动力人工边界转换，先采用静力边界求解静力问题，然后在动力分析中将静力边界换成动力人工边界，同时将静力分析结束时边界节点的约束反力以集中力形式施加在对应的单元节点上。

第3种概念清晰且易于在通用数值软件中实施，然而已有的转换模拟方法只能使模型的应力场在静动力边界转换前后基本一致，模型的位移场在边界转换前后的差异十分明显。本文就此基于ABAQUS提出一种使分析模型的应力场及位移场在静动力边界转换前后都一致的转换模拟方法及操作流程。

已有转换方法(方法①)：先分步进行静力分析求得开挖支护后地下结构-围岩系统在重力等外载的静应力场及约束反力，然后另建一个尺寸、网格等与原模型相同的空模型，将静应力场及外荷载输入该空模型，再进行一次静力分析，得到地下结构-围岩系统初始应力状态后，最后在动力分析步中将边界条件转换为黏弹性边界并输入地震荷载进行动力分析。

新的转换方法(方法②)：分两步进行。第一步，首先计算重力作用下围岩的天然应力场，进而按照施工步骤，依次建立多个静力分析步计算并导出洞室开挖支护后地下模型边界处的约束反力，然后直接在静力计算模型后添加动力分析步，导出inp文件。第二步，修改导出的inp文件：利用关键字*include在动力分析步中添加黏弹性边界、地震荷载及在第一步中导出的约束反力，运行新inp文件即可。需要注意的是，黏弹性边界中的弹簧阻尼器元件会在所有分析步中产生作用，故需要利用生死单元技术(关键字*remove与*add)在静力步中抑制其作用，以免引起计算收敛性问题。

分别采用方法①与方法②对某隧道结构施加零动力荷载进行计算。静力分析时仅考虑自体质量，模型底部固定，侧边界法向约束，其静、动力计算结束时的应力及位移结果分别如图6、图7所示。

图6 静力开挖及零动力荷载下土-结构系统应力云图Fig.6 Stress nephogram of soil-structure system under static excavation and zero dynamic load

图7 静力开挖及零动力荷载下土-结构系统位移云图Fig.7 Displacement nephogram of soil-structure system under static excavation and zero dynamic load

由图6、图7可知，方法①应用下，人工边界转换后位移场与洞室开挖后土体沉降规律差异较大，这是因为该转换方法仅用一次静力计算获取静力平衡状态，忽略了施工过程对隧道工程动力分析初始位移状态的影响。而方法②在动力分析前进行了多步静力计算模拟施工流程，故静动转换前后模型的应力、位移空间分布及大小基本一致，相对误差低于0.7%。误差来源主要是静力开挖后侧边界会产生竖向位移，使侧边界的弹簧元件在动力分析初始时刻产生了微弱的初始应力，从而导致模型整体受到轻微的不平衡力作用，若静力计算时采用固定边界，该转换方法的误差会更小。

3 山岭隧道斜入射地震响应算例

3.1 计算模型

以西藏自治区某山岭隧道洞口段为研究对象，针对典型马蹄形隧道的横截面，建立二维围岩-隧道整体分析模型，如图8所示。模型水平宽度150 m，竖向高度100 m，在模型底部及两侧施加黏弹性边界。隧道的最大洞径13.6 m，高度10.9 m，埋深20 m，初期支护厚度0.2 m，二次衬砌0.5 m，洞周设置有锚杆加固围岩。围岩及衬砌均采用平面单元，衬砌结构采用ABAQUS内嵌的混凝土塑性损伤模型(CDP模型)模拟混凝土的非线性力学行为，该模型中的拉、压损伤因子利用GB 50010―2010《混凝土结构设计规范》规定的C30混凝土应力-应变关系，根据Sidoroff的能量等价法计算得到；Ⅴ级围岩模拟采用Mohr-Coulomb本构模型，锚杆加固采用锚固区模拟，模型材料参数如表1和表2所示。监测点布置方案如图9所示，θ为监测点和横截面形心O点连线与竖直方向之间的夹角，用以表征衬砌位置，假设倾斜地震波从模型左侧输入计算区域。

图8 围岩-隧道整体分析模型(m)Fig.8 Tunnel calculation model diagram(m)

图9 监测点示意图Fig.9 Layout of monitoring nodes

表1 模型中的材料参数Tab.1 Material parameters in the model

表2 二衬混凝土物理力学参数Tab.2 Physical and mechanical parameters of secondary lining concrete

廖振鹏等[20]研究表明，计算模型能准确模拟地震波的传播过程，网格尺寸应小于输入波最高频率成分所对应波长的1/10～1/8，即要求

Δl≤λ/10

(5)

为了平衡计算精度与效率，将隧道结构及附近围岩用较细的网格划分，距离结构较远的围岩部位以稍大的网格离散。经过试算，模型网格的最大尺寸取1.8 m。

本文采用瑞雷阻尼，计算方法如下

C=αM+βK

(6)

式中：M为质量矩阵；K为刚度矩阵；α，β分别为质量阻尼系数与刚度阻尼系数，可按式(7)计算

(7)

式中：ϖi,ϖj分别为计算模型的第i、第j阶自振频率；ξ为阻尼比，结构与岩土体均取0.05；提取模型的第1阶与第3阶自振频率代入式(7)计算。模态分析时，模型的侧边界设置为竖向约束，水平自由，底边界设置为固结。

由式(4)结合地层参数可知，SV波临界入射角为28.01°，故其入射角度αis取0°(垂直入射),10°，15°，20°，25°,28°；P波入射角度αip取0°，15°，30°，45°，60°，75°。

图10 Kobe波加速度时程曲线Fig.10 Acceleration time history curve of Kobe wave

选用近场地震动Kobe波作为输入波，为节省时间，截取包含峰值加速度前后10 s的数据段，经基线校准和低通滤波后得到一条20 s的加速度时程曲线，根据工程场址要求将峰值统一调幅为0.3g。入射波的加速度时程曲线如图11所示。

图11 P波斜入射衬砌加速度峰值分布图Fig.11 Distribution of acceleration peak value of lining with oblique incidence of P wave

3.2 结果分析

3.2.1 加速度

图11、图12显示了地震P波与SV波不同入射角度下，衬砌沿着横截面方向(-180°≤θ≤180°)各监测点峰值加速度变化规律。

由图11、图12可知，P波和SV波的入射角度对衬砌结构加速度幅值的大小有着显著的影响，而对其空间分布的影响不大。P波入射时，在水平方向上，加速度峰值首先随着入射角度的增加而明显增大，60°入射时达到最大，然后随着入射角度的增加而减小；在竖直方向上，衬砌加速度峰值随着入射角增加而显著减小，如0°入射时，拱顶的加速度峰值为4.16 m/s2，75°入射时为1.37 m/s2，相比垂直入射，降低幅度为67%，这表明垂直入射是P波下垂直地震反应最不利的情况。SV波入射时，随着入射角度增加，衬砌水平加速度和竖向加速度均呈增大趋势，在28°入射时达到最大。但水平加速度仅在入射角度较大时，增幅才较为明显；而竖直加速度随入射角度增加的增大则较为均匀。表明P波下隧道衬砌的加速度对入射角度更为敏感。

图12 SV波斜入射衬砌加速度峰值分布图Fig.12 Distribution of acceleration peak value of lining with oblique incidence of SV wave

3.2.2 应 力

提取衬砌内表面的第一应力(主拉应力)峰值评估混凝土衬砌结构的受力状态。

由图13可知，地震波斜入射下衬砌内的主拉应力随入射角度增加呈整体增大趋势，说明垂直入射假定将明显低估地震作用下马蹄形断面隧道结构的应力状态。随着入射角度的不同，P波与SV波下衬砌应力的空间分布特征也不同，表现为峰值应力所在的位置发生改变；P波下最大值出现在拱顶、拱底处，尤其是拱顶，几乎从零增长到2.48 MPa，其余部位的应力较小且起伏变化不大，表明P波斜入射下拱顶和拱底更易出现受拉破坏；SV波入射角较小时，应力最大值出现在拱顶、拱底处，入射角较大时，最大值出现在拱脚，且远远大于其他部位的应力，表明除了拱顶和拱底，SV波斜入射条件更易引起拱脚处的应力集中而出现拉裂损伤；此外，相比P波斜入射，SV波斜入射下拱肩处应力水平较高。P波下大多应力峰值在60°入射时产生，SV波则大多在28°入射时产生，对比两者结果可知，地震强度相同时，斜入射SV波下隧道衬砌的应力水平高于斜入射P波的情形。

图13 地震波斜入射衬砌各监测点主拉应力峰值Fig.13 Peak value of principal tensile stress at each monitoring point of lining with oblique incidence of seismic wave

3.2.3 混凝土损伤发展

本文仅展示P波60°、SV波28°斜入射的损伤结果。如图14～图16所示，地震动垂直与斜入射作用下隧道衬砌损伤演变差异明显：垂直入射P波与SV波作用下，隧道衬砌均未发生损伤；斜入射P波与SV波作用下，隧道衬砌均发生了不同程度的损伤。

图14 地震波垂直入射衬砌拉伸损伤Fig.14 Tensile damage of lining under vertical incidence of seismic wave

图15 P波斜入射衬砌拉伸损伤发展过程Fig.15 Development process of tensile damage of lining under oblique incidence of P wave

图16 SV波斜入射衬砌拉伸损伤发展过程Fig.16 Development process of tensile damage of lining under oblique incidence of SV wave

斜入射P波作用下隧道衬砌的损伤发展过程：5.1 s左右，仰拱右侧内表面开始出现损伤；20 s时，仰拱左右侧均有小范围累计损伤。斜入射SV波作用下隧道衬砌的损伤发展过程：3.8 s左右，仰拱左右侧内表面开始出现损伤；5.2 s左右，仰拱左右侧损伤加剧，贯穿衬砌，同时拱顶附近开始出现较大损伤；8.6 s左右，左右拱脚处外表面先开始出现损伤，随后向周围扩散，逐渐形成贯穿损伤；20 s时，拱肩外表面开始出现轻微损伤，仰拱和拱部亦形成多处贯穿损伤，最大累计损伤发生在拱顶附近处。斜入射SV波下衬砌损伤程度更大，表明峰值地面加速度(peak ground acceleration，PGA)大小相同时，斜入射SV波对隧道衬砌具有更大的破坏作用，相比斜入射P波，SV波作用使衬砌横向上出现大面积的损伤，产生更大的塑性应变，从而降低隧道整体的稳定性。

为了解P波斜入射条件下衬砌结构的损伤发展过程，另外计算了PGA=0.6g时P波60°入射情况下衬砌结构的动力反应。如图17所示，P波斜入射下的损伤发展过程与SV波斜入射下的损伤发展过程类似，先后在仰拱、拱顶和拱脚处表面出现损伤，再逐渐形成贯穿损伤，最后拱肩附近出现轻微损伤。

图17 P波斜入射衬砌拉伸损伤发展过程(PGA=0.6g)Fig.17 Development process of tensile damage of lining under oblique incidence of P wave (PGA=0.6g)

3 结 论

基于黏弹性人工边界与等效节点力代替地震动的波动输入方法，在通用ABAQUS平台上建立了P波与SV波倾斜输入数值模型，提出了一种改进的静动力耦合边界模拟方法应用于地下结构的静动力联合分析，并研究了斜入射P波与SV波下马蹄形断面隧道结构的动力响应，结论如下：

(1)本文针对静动力耦合边界模拟提出一种改进的模拟方法，经算例验证，可使分析模型的应力场与位移场保持一致，结果表明，相对误差低于0.7%。

(2)随入射角度的增加，P波下隧道结构的水平加速度先增加后减小，竖向加速度则呈递减趋势，而SV波下，隧道结构的水平与竖向加速度均呈增大趋势； P波与SV波斜入射下拱顶、拱底处应力分布情况相似，而入射角度较大时，SV波斜入射下拱肩、拱脚处的应力水平明显更高于P波斜入射工况。

(3)P波与SV波斜入射作用下隧道衬砌的损伤位置类似，主要集中在上、下拱及其附近位置；通过应力和拉伸损伤结果表明，在相同地震强度下，SV波斜入射下对隧道结构的破坏程度高于P波。