强震区隧洞洞口段地震响应数值模拟研究

刘国庆，肖明，陈俊涛，陈世杰



强震区隧洞洞口段地震响应数值模拟研究

刘国庆1, 2，肖明1, 2，陈俊涛1, 2，陈世杰1, 2

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

假定地震波为垂直入射的弹性平面波，探讨隧洞洞口段平面波输入方法；基于拉、压条件下混凝土不同的损伤演化过程，建立适于编程的混凝土衬砌弹性动力损伤本构模型；针对地震作用下衬砌与围岩相互作用特点，建立衬砌与围岩联合承载分析模型；以处于强震区的滇中引水隧洞某出口段为实例，分析平面P波和SH波同时作用下洞口段的地震响应特性。研究结果表明：衬砌各部位位移与地震波的振动方向有关，峰值位移随着与洞口距离增大而减小，并在48 m处趋于稳定；在地震作用下，衬砌的损伤破坏以拉损为主，损伤系数随时间逐渐增大；衬砌中脱开区、滑移区的分布与其震损较严重区域的分布基本一致，主要分布于距洞口21 m之内的拱肩、拱腰和拱脚部位。

洞口段；地震波输入；混凝土衬砌；联合承载；地震响应；数值模拟

汶川大地震后，通过对震中区域公路隧道震害调查发现，洞口段在地震发生时极易破坏，是仅次于断层破碎带结构震害的抗震薄弱段[1−4]。隧洞洞口段埋深浅，所处地质条件差，岩层多为风化严重的堆积体，容易造成坡体失稳。洞口段遇强震时，常发生地表开裂、落石、坡体滑塌等，进而造成衬砌开裂、错位、垮塌等现象，从而影响隧洞的正常运行。为此，进行洞口段的抗震稳定研究是当前隧洞中迫切需要解决的问题之一。针对洞口段的地震响应分析已有一定的研究成果，研究方法主要采用模型试验和数值模拟。在模型试验方面，申玉生等[5]以雅安—泸沽高速公路山岭隧道为依托工程，对洞口段结构动力响应进行了振动台模型试验研究，分析了衬砌的破坏形态；蒋树屏 等[6]通过大型振动台模型试验，对嘎隆拉隧道洞口段衬砌的加速度及应变响应进行了研究。在数值模拟方面，蒋树屏等[7]采用动力有限元法，分析了嘎隆拉隧道洞口段衬砌的加速度、应力和位移响应特性；李育枢等[8]以国道318线某隧道为实例，基于FLAC3D分别从洞口边坡、洞门建筑及衬砌结构对洞口段的地震安全性进行了评价。由此可知，在数值计算中，洞口段抗震稳定分析的关键在于模拟混凝土衬砌的动力响应，因此，建立合理的混凝土动力本构模型及衬砌与围岩联合承载模型显得尤为重要。然而，目前的研究大多假定衬砌材料是线弹性的，且衬砌与围岩共用节点，没有考虑二者之间的动接触行为，这与实际情况不符。本文作者假定地震波为垂直入射的弹性平面波，探讨洞口段平面P波和S波输入方法；针对拉、压条件下混凝土不同的损伤演化过程，建立1种简单的混凝土衬砌动力损伤模型；基于动接触力算法，建立1种衬砌与围岩联合承载分析模型；结合滇中引水隧洞某出口段实例，对衬砌的地震响应特性进行模拟，以便为洞口段结构的抗震减震设计提供参考。

1 隧洞洞口段平面地震波输入方法

1.1 基本思路

基于波场分解原理，总波场可分解为内行波场和外行波场，由于外行波场可在人工边界处实现透射，对有限元计算没有影响，故不进行专门处理。因此，为实现地震荷载的准确输入的关键在于获取人工边界处的内行波场。考虑地震波垂直入射，模型底边界的内行波场等于入射波场。由于侧向边界的入射场和经地表自由面反射的波场均平行于该边界，鉴于局部人工边界不能模拟平行于侧向边界的波动场，此时，可将侧向边界的内行波场作为自由场计算。当入射波场已知时，将入射波场和反射波场进行叠加，即可得到侧向边界的自由场。基于黏弹性人工边界[9]，可将波的输入问题转化为求解作用于人工边界节点的等效节点力问题，以此实现地震波的输入。

图1 地震波入射模型

1.2 平面P波入射

已知位移场，相应的速度场可通过求导或差分得到。根据广义虎克定律，由人工边界处的内行波场可计算对应的应力场，进而得到人工边界处的等效节点力[10−11]。

对于−边界：

对于−边界：

对于+边界：

对于−和+边界：

1.3 平面S波入射

各人工边界处的等效节点力如下。

对于−边界：

对于−边界：

对于+边界：

对于−和+边界：

2 混凝土衬砌动力损伤模型

混凝土是一种非均质的人造石材，内部具有天然微裂缝，在地震作用下易产生损伤特性。

2.1 混凝土动力损伤本构模型

假定混凝土为各向同性介质，基于连续损伤力学基本理论，混凝土应力−应变关系可表示为

2.2 混凝土动力损伤演化模型

图2 混凝土动力损伤模型

2.3 混凝土三维动力损伤模型

3 衬砌与围岩联合承载分析模型

衬砌与围岩联合地震动响应过程中存在复杂的动接触行为，地震作用可能导致衬砌损伤开裂、脱开或滑移。工程实践表明，在不含断层破碎带的隧洞中，衬砌结构的脱开、滑移等破坏主要是局部的，大滑移破坏现象较为少见。因此，本文假定地震加载过程中衬砌与围岩接触面节点始终或近似处于点对接触状态，并假定地震加载前接触节点对处于黏结接触状态[15]。

式中：，和分别为接触节点质量、阻尼和刚度矩阵；，和分别为节点加速度、速度和位移向量；为外荷载向量；为动接触力向量。

每1时步计算完毕后，需要校核接触面的接触状态，并对动接触力进行修正，具体方法如下。

4 工程实例

4.1 计算模型与参数

滇中引水工程主要由水源工程、输水工程、配套工程等组成，是一项水资源综合利用的水利工程。输水总干线全长848.18 km，工程规模巨大，沿线地质条件复杂，抗震设防烈度较高，隧洞稳定问题突出。

本文选取沿线某出口段进行地震响应分析。该段隧洞最大埋深为35.00 m，坡角约50°，设计水深为7.51 m，围岩以Ⅳ类为主。开挖断面为马蹄形，最大开挖宽×高为8.30 m×9.87 m，采用C25钢筋混凝土衬砌结构，衬砌厚0.50 m。

建立出口段的三维有限元模型，如图4所示。模型全部采用八节点六面体单元进行离散，共剖分 110 532个单元和117 280个节点，其中衬砌单元5 760个。模型轴沿水平向与洞轴线垂直；轴与洞轴线重合，以顺水流向为正；轴与大地坐标一致。

数据单位：m

4.2 计算条件

计算程序采用自主开发的三维弹塑性损伤动力显示有限元计算平台[17]，并将本文地震波输入方法、混凝土动力损伤模型及衬砌与围岩联合承载模型嵌入其中。

表1 材料力学参数

模型底部采用黏弹性人工边界，四周采用自由场人工边界。地震波选取美国强震记录EI-Centro波，峰值加速度为3.417 m/s2，截取其中20 s作为入射波，然后经滤波、基线校正得到加速度时程曲线如图5所示。地震波由模型底部垂直入射，鉴于SV波对隧洞结构影响较小[11]，计算同时考虑P波和SH波对结构的作用。其中，SH波采用如图5所示的入射波，P波加速度取为SH波的2/3。

图5 入射波加速度时程曲线

沿隧洞轴线，距洞口每隔6 m布置1个监测断面，共计12个监测断面。在每个监测断面内衬砌结构上布置5个监测点，以监测地震加载过程中衬砌的应力、位移等特性。衬砌监测点布置如图6所示。

4.3 结果分析

本文主要从衬砌的位移、应力、损伤及脱开、滑移特性4个方面对衬砌的地震响应进行分析。以第6个监测断面为例，分析上述指标的时程特性。

4.3.1 衬砌位移分析

在地震加载过程中，衬砌监测点向和向的位移时程如图7所示。由图7可知：各监测点位移时程曲线变化规律基本一致，表明衬砌各部位处于同步振动状态。在SH波作用下，拱肩和拱腰位移较大，底拱位移最小；在P波作用下，底拱位移最大，拱顶位移最小，表明衬砌各部位位移与地震波的振动方向 有关。

图6 衬砌监测点布置

方向：(a) x向；(b) z向

监测点向和向峰值位移同时出现在2.46 s，将2.46 s时所有监测断面5个监测点的峰值位移沿洞轴线变化规律绘制成曲线，如图8所示。由图8可知：各监测点向峰值位移沿洞轴线呈逐渐减小的趋势；在距洞口0~48 m段，峰值位移减小幅度较大；当距离大于48 m时，峰值位移变化不大。这主要是因为洞口边坡对地震波有反射作用，并伴随波形转换现象，致使该段波场复杂，对隧洞结构影响较大[18]，在该时间点拱肩与底拱最大相对位移达到0.98 cm。类似 地，各监测点向峰值位移沿洞轴线变化规律与向一致，在该时间点底拱与拱顶最大相对位移达到 0.55 cm。

方向：(a) x向；(b) z向

4.3.2 衬砌应力分析

衬砌监测点最小、最大主应力时程如图9所示。由图9(a)可知：各监测点最小主应力呈波动变化，波动较剧烈的时间段为2.6~9.4 s，震后压应力较小。整体来看，拱肩、拱腰最大压应力比其他部位的大，最大压应力为9.73 MPa，没有超过混凝土抗压强度 (11.9 MPa)，表明衬砌在地震加载过程中没有出现压损破坏。由图9(b)可知：在0~2 s，各监测点最大主应力不大；2 s以后，最大主应力迅速增大至混凝土抗拉强度(1.27 MPa)，表明衬砌出现拉损破坏。在2~7 s，最大主应力波动较剧烈，可能会加剧衬砌损伤，造成衬砌疲劳破坏。地震完成后，最大主应力与初始时的相差不大。

主应力：(a)最小主应力；(b)最大主应力

4.3.3 衬砌损伤分析

由衬砌应力分析可知衬砌在地震中的损伤破坏以拉损为主，下面主要分析衬砌的拉损状况。

衬砌监测点损伤系数变化时程如图10所示。由图10可知：在0~2 s，各监测点没有出现损伤；2 s以后，损伤系数随时间逐渐增大；在2~7 s时，地震波振动较为剧烈，损伤系数增幅明显，7 s以后，损伤系数增幅较小。由此可见，衬砌的损伤演化过程与最大主应力变化规律是紧密相关的。整体来看，衬砌拱腰损伤最严重，拱肩和拱脚次之，拱顶和底拱损伤最小。

地震完成后，衬砌结构整体(取距洞口60 m范围)损伤系数分布如图11所示。由图11可知：衬砌不同部位损伤程度不同；从横向看，损伤较严重的部位分布于拱肩、拱腰和拱脚；从纵向看，洞口处损伤最严重，局部损伤系数达到1.0。若定义损伤系数大于0.2的区域为震损区，则震损区主要分布于距洞口48 m范围内，其中震损较严重的区域(＞0.5)分布长度约为21 m。

1—拱顶；2—拱肩；3—拱腰；4—拱脚；5—底拱。

图11 震后衬砌损伤系数分布

4.3.4 衬砌脱开、滑移分析

地震作用前，衬砌与围岩处于黏结接触状态；在地震过程中，衬砌与围岩接触节点的动接触力一旦突破黏聚力，极易造成接触面脱开或滑移，且该过程是不可逆的。地震完成后，衬砌结构(取距洞口30 m范围)中脱开区和滑移区分布如图12和图13所示。

由图12可知：衬砌脱开区主要分布于两侧拱腰、拱顶和底拱部位；拱腰脱开区从洞口向隧洞内部延伸较长，达到21 m；拱顶和底拱脱开区分布范围较小，延伸范围限制在9 m以内。由图13可知：衬砌滑移区主要分布于两侧拱肩、拱脚部位，这4处的滑移区从洞口向隧洞内部延伸长度为21 m。由此可知，衬砌脱开、滑移区主要分布于距洞口21 m范围内的拱肩、拱腰和拱脚部位，这与衬砌震损较严重区域的分布规律基本一致。

图12 震后衬砌脱开区分布

图13 震后衬砌滑移区分布

综上可知：距洞口48 m范围内衬砌位移响应较大、损伤较严重，因此，可将洞口段隧洞结构抗震设防长度设定为48 m；另外，衬砌脱开、滑移区主要分布于距洞口21 m范围内，该范围可作为重点设防 长度。

5 结论

1) 地震作用下衬砌各部位处于同步振动状态，各部位位移与地震波的振动方向有关；距离洞口越远，衬砌峰值位移越小，并在48 m处趋于稳定。

2) 在地震加载过程中，衬砌最大主应力极易达到抗拉强度，因此，衬砌的损伤破损以拉损为主；最大主应力的剧烈波动可能加剧衬砌损伤，造成其疲劳 破坏。

3) 衬砌各部位损伤系数随时间逐渐增大；衬砌震损区主要分布于距洞口48 m范围内的拱肩、拱腰和拱脚部位，并且距洞口越近，震损越严重。

4) 衬砌脱开、滑移区主要分布于距洞口21 m范围内的拱肩、拱腰和拱脚部位，这与衬砌震损较严重的区域分布规律基本一致。由此可见拱肩、拱腰和拱脚是衬砌结构抗震的薄弱部位。

[1] 李天斌. 汶川特大地震中山岭隧道变形破坏特征及影响因素分析[J]. 工程地质学报, 2008, 16(6): 742−750.LI Tianbin. Failure characteristics and influence factor analysis of mountain tunnels at epicenter zones of great Wenchuan earthquake[J]. Journal of Engineering Geology, 2008, 16(6): 742−750.

[2] 高波, 王峥峥, 袁松, 等. 汶川地震公路隧道震害启示[J]. 西南交通大学学报, 2009, 44(3): 336−341. GAO Bo, WANG Zhengzheng, YUAN Song, et al. Lessons learnt from damage of highway tunnels in Wenchuan earthquake[J]. Journal of Southwest Jiaotong University, 2009, 44(3): 336−341.

[3] 崔光耀, 王明年, 于丽, 等. 汶川地震公路隧道洞口结构震害分析及震害机理研究[J]. 岩土工程学报, 2013, 35(6): 1084−1091. CUI Guangyao, WANG Mingnian, YU Li, et al. Seismic damage and mechanism of portal structure of highway tunnels in Wenchuan Earthquake[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(6): 1084−1091.

[4] WANG Z Z, ZHANG Z. Seismic damage classification and risk assessment of mountain tunnels with a validation for the 2008 Wenchuan earthquake[J]. Soil Dynamics and Earthquake Engineering, 2013, 45: 45−55.

[5] 申玉生, 高波, 王英学. 强震区山岭隧道洞口段结构动力特性分析[J]. 岩石力学与工程学报, 2009, 28(增1): 3131−3136. SHEN Yusheng, GAO Bo, WANG Yingxue. Structural dynamic properties analysis for portal part of mountain tunnel in strong earthquake area[J]. Chinese Journal of Rock Mechanics and Engineering, 2009, 28(Suppl 1): 3131−3136.

[6] 蒋树屏, 文栋良, 郑升宝. 嘎隆拉隧道洞口段地震响应大型振动台模型试验研究[J]. 岩石力学与工程学报, 2011, 30(4) :649−656. JIANG Shuping, WEN Dongliang, ZHENG Shengbao. Large-scale shaking table test for seismic response in portal section of Galongla tunnel[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(4): 649−656.

[7] 蒋树屏, 蒋华, 王晓雯, 等. 高烈度地震区公路隧道洞口段地震响应分析[J]. 现代隧道技术, 2008, 45(6): 5−10. JIANG Shuping, JIANG Hua, WANG Xiaowen, et al. Seismic response analysis for a highway tunnel portal section in high-intensity seismic region[J]. Modern Tunnelling Technology, 2008, 45(6): 5−10.

[8] 李育枢, 李天斌. 高烈度地震区隧道洞口段地震安全性分析及评价[J]. 现代隧道技术, 2011, 48(3): 46−52. LI Yushu, LI Tianbin. Seismic analysis and safety evaluation of the tunnel portal section in regions of high seismic intensity[J]. Modern Tunnelling Technology, 2011, 48(3): 46−52.

[9] 刘晶波, 王振宇, 杜修力, 等. 波动问题中的三维时域粘弹性人工边界[J]. 工程力学, 2005, 22(6): 46−51. LIU Jingbo, WANG Zhenyu, DU Xiuli, et al. Three-dimensional visco-elastic artificial boundaries in time domain for wave motion problems[J]. Engineering Mechanics, 2005, 22(6): 46−51.

[10] 杜修力, 赵密. 基于黏弹性边界的拱坝地震反应分析方法[J]. 水利学报, 2006, 37(9): 1063−1069. DU Xiuli, ZHAO Mi. Analysis method for seismic response of arch dams in time domain based on viscous-spring artificial boundary condition[J]. Journal of Hydraulic Engineering, 2006, 37(9): 1063−1069.

[11] 汪精河. 基于波动理论的高陡边坡山岭隧道洞口地震响应分析[D]. 成都: 西南交通大学土木工程学院, 2016: 79−84. WANG Jinghe. Seismic response analysis for mountain tunnel portals in high steep slope based on wave theory[D]. Chengdu: Southwest Jiaotong University. School of Civil Engineering, 2016: 79−84.

[12] 李杰, 曾莎洁, 任晓丹. 混凝土动力随机损伤本构关系[J]. 同济大学学报(自然科学版), 2014, 42(12): 1783−1789. LI Jie, ZENG Shajie, REN Xiaodan. A stochastic rate-dependent damage model for concrete[J]. Journal of Tongji University: Natural Science, 2014, 42(12): 1783−1789.

[13] CEN Weijun, WEN Langsheng, ZHANG Ziqi, et al. Numerical simulation of seismic damage and cracking of concrete slabs of high concrete face rockfill dams[J]. Water Science and Engineering, 2016, 9(3): 205−211.

[14] ZHU W C, TANG C A. Micromechanical model for simulating the fracture process of rock[J]. Rock Mechanics and Rock Engineering, 2004, 37(1): 25−56.

[15] 邓建, 肖明, 陈俊涛. 高压引水隧洞运行期复杂承载过程数值分析[J]. 中南大学学报(自然科学版), 2017, 48(5): 1261−1267. DENG Jian, XIAO Ming, CHEN Juntao. Numerical analysis for complex bearing process of high pressure diversion tunnel at runtime[J]. Journal of Central South University(Science and Technology), 2017, 48(5): 1261−1267.

[16] 何川, 李林, 张景, 等. 隧道穿越断层破碎带震害机理研究[J]. 岩土工程学报, 2014, 36(3): 427−434. HE Chuan, LI Lin, ZHANG Jing, et al. Seismic damage mechanism of tunnels through fault zones[J]. Chinese Journal of Geotechnical Engineering, 2014, 36(3): 427−434.

[17] 张志国, 肖明, 陈俊涛. 大型地下洞室地震灾变过程三维动力有限元模拟[J]. 岩石力学与工程学报, 2011, 30(3): 509−523. ZHANG Zhiguo, XIAO Ming, CHEN Juntao. 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.

[18] 侯森, 陶连金, 赵旭, 等. SH波作用下山岭隧道洞口段结构动力响应研究[J]. 岩石力学与工程学报, 2015, 34(2): 340−348. HOU Sen, TAO Lianjin, ZHAO Xu, et al. Dynamic response of portal section of mountain tunnels under seismic SH wave action[J]. Chinese Journal of Rock Mechanics and Engineering, 2015, 34(2): 340−348.

(编辑 伍锦花)

Numerical simulation of seismic response for tunnel portal section in strong earthquake area

LIU Guoqing1, 2, XIAO Ming1, 2, CHEN Juntao1, 2, CHEN Shijie1, 2

(1. State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China; 2. Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering of Ministry of Education, Wuhan University, Wuhan 430072, China)

Assuming that seismic wave is vertical incident elastic plane wave, the plane wave input method for tunnel portal section was studied. Based on different damage evolution processes of concrete under tension and compression conditions, an elastic dynamic damaged constitutive model of concrete lining suitable for programming was constructed. According to the interaction characteristics between lining and rock under seismic load, an analytical model of joint loading between lining and rock was established. Taking an exit section of Dianzhong diversion tunnel in strong earthquake area as an example, the seismic response characteristics of portal section under synchronous actions of plane P wave and SH wave were analyzed. The results show that the displacements of different lining parts are related to the vibration direction of seismic wave, and the peak displacements decrease gradually from the portal to the interior, and then stay stable when the distance is 48 m. The damage failure of lining is mainly caused by tension under seismic load, and the damage coefficient increases gradually with time. The separation and slip zone distributions of the lining are basically consistent with its severe seismic damage area, which are mainly at spandrel, haunch and arch foot within 21 m of distance from the portal.

portal section; seismic input; concrete lining; joint loading; seismic response; numerical simulation

10.11817/j.issn.1672-7207.2018.11.022

TU45

A

1672−7207(2018)11−2804−09

2017−12−17；

2018−02−27

国家重点基础研究发展计划(973计划)项目(2015CB057904)；国家自然科学基金资助项目(51579191) (Project (2015CB057904) supported by the National Key Basic Research and Development Program (973 Program) of China; Project(51579191) supported by the National Natural Science Foundation of China)

肖明，博士，教授，从事地下结构稳定数值分析研究；E-mail: mxiao@whu.edu.cn