引言

随着我国社会经济的快速发展和城市化进程的推进，土地资源紧缺问题日益突出，合理利用和开发地下空间资源是解决城市可持续化发展的重要手段。城市地下综合管廊作为容纳2类及以上城市工程管线的构筑物及附属设施，与直埋管线相比具有一定优势，是城市现代化发展的必然趋势。我国自1958年在北京天安门广场改造工程中建设了第一条综合管廊以来，已在多个城市开始建造地下综合管廊（钱七虎等，2007）。地下综合管廊工程延伸范围宽广，由于场地土特性复杂多变，地下综合管廊工程施工过程中易遇砂土液化、软土震陷等诸多问题，一旦发生破坏，将对社会经济造成不可估量的损失。因此，地下综合管廊防震减灾问题成为地震工程领域研究热点（黄雨等，2013，袁勇等，2022）。

研究地下综合管廊结构地震响应规律的主要途径有物理模型试验和数值模拟分析，如史晓军等（2008）较早地开展了地下综合管廊大型振动台模型试验研究，分析了模型结构加速度响应、应变响应、接触面土压及引起结构内力的影响因素，研究发现地下综合管廊壁板与土接触面的作用力直接导致结构产生了内力；蒋录珍等（2015）以横向有接头的综合管廊非一致激励振动台试验为基础，采用非线性弹簧单元模拟地下综合管廊接头建立了有限元模型，研究发现地下综合管廊响应分析的数值模拟结果与试验结果较一致，非线性弹簧单元可作为综合管廊结构接头进行数值模拟；施有志等（2018）建立了双仓综合管廊三维动力有限元数值模型，研究了Rayleigh波与底部地震加速度共同作用下综合管廊的动力响应特征；Zhang等（2019）依托特高压综合管廊工程，以隧道-竖井变刚度节点为对象，开展了大规模多工况振动台模型试验研究，分析了软土中竖井-隧道连接在横向激励作用下的差异响应；Duan等（2019）以西安新建预制综合管廊为原型，进行了系列振动台试验，研究黄土场地综合管廊抗震性能，并分析了含管道的综合管廊地震响应；仉文岗等（2020）利用振动台试验及数值模拟对比研究了不同地震波作用下管廊结构及周围土体动力响应特征，发现地震过程中综合管廊结构及周围土体基本保持一致的运动模式，结构拐角处会产生较大的弯矩响应；韩俊艳等（2022）基于北京工业大学多台阵振动台系统，研究了一致、非一致激励振动台试验中不同箱体内土体、管道地震响应的异同，研究发现非一致激励使管道和土体产生的相对位移较一致激励大，会放大管道应变响应；汤鹏等（2022）以越江气体绝缘输电线路（GIL）综合管廊为工程背景，考虑内部支架结构与管廊结构相互作用，采用有限元方法分析了不同剖面工况下结构抗震性能。

综上所述，已有不少学者对综合管廊地震反应进行了研究，但对考虑越江综合管廊结构上覆水影响的精细化三维有限元模型分析研究较少。为此，本文基于苏通GIL综合管廊工程，建立精细化有限元模型，采用修正Davidenkov黏弹性动力本构和Byrne孔压增量耦合模型研究综合管廊非线性地震反应特性，以期为类似工程抗震设计提供参考。

1.   有限元模型建立

1.1   GIL综合管廊工程概况

以苏通GIL综合管廊工程为依托，本工程所在区域地震活动属中强地震活动水平，有必要开展综合管廊地震反应特性研究（汤鹏等，2022）。选取苏通GIL综合管廊穿越长江段的江南至江北某段截面，该段水深约40 m，管廊土体以下埋深为22.6 m，主要穿越粉细砂和细砂层，土质较软，且有局部液化的可能。管廊内径10.5 m，外径11.6 m，整环为“7+1”分块错缝拼装，管片厚0.55 m，每段管片长2 m，管片采用C60混凝土，主筋为HRB400钢筋，接头螺栓等级为10.9级，管廊内部结构采用C30混凝土。

1.2   三维有限元数值模型

采用ABAQUS通用有限元软件进行精细化建模，如图1所示，由于分析断面处土层几何几乎呈均匀分布，将土简化为水平成层，从上至下依次为粉质黏土混粉砂层、粉土层、粉细砂层、细砂层和中粗砂层。考虑模型边界效应，横断面长度取200 m，即地基宽度为地铁车站结构宽度的10倍，高80 m，纵向选取10个管片（共20 m）进行分析（Wang等，2017）。模型土体、管廊及附属结构采用8结点六面体线性减缩积分单元（C3D8R）模拟，螺栓与钢筋采用三维线性梁单元（B31）模拟，模型共划分为154 618个单元。网格划分时考虑了截止频率对应的波长沿土层深度变化（Ruan等，2019）。本文截止频率取25 Hz，据此，竖向单元尺寸＜2 m，横向单元尺寸＜1.2 m。钢筋与螺栓截面面积参考GB 50010—2010《混凝土结构设计规范》（中华人民共和国住房和城乡建设部，2011）及GB/T 51438—2021《盾构隧道工程设计标准》（中华人民共和国住房和城乡建设部等，2021）中的相关规定，钢筋与螺栓采用嵌入约束埋入管片中，管片与管片、管片与土之间采用通用接触，力学行为表现为法向硬接触，切向摩擦系数为0.6，考虑管廊内部结构采用混凝土浇筑成型，内部结构与管廊之间采用绑定约束。

图 1  有限元模型

Figure 1.  3D finite element model

下载: 全尺寸图片 幻灯片

模型分析步骤主要为：平衡地应力→浇筑附属结构→地震动输入。为得到地震前管廊与土体应力分布，首先对土体和管廊进行地应力平衡，以消除重力作用下的位移，然后激活管廊内附属结构，以得到动力分析前的初始状态，该步骤中的边界条件为侧面法向约束、底部3个方向位移约束，由于土体上部水深40 m，不能忽略静水压力的影响，在模型顶部施加392 kPa的静水压力。通过ABAQUS重启动，实现隐式与显式算法间的数据传递，可将所有场变量作为初始条件。在step模块中勾选重启动选项，选择保存最后一帧数据，对整个模型施加initial state的预定义场。

1.3   材料模型与参数

为了描述土体在地震动作用下由于孔压积累导致的非线性行为，土体本构模型采用Chen等（2019，2021）提出的修正Davidenkov黏弹性动力本构和Byrne孔压增量耦合模型，如图2所示。土体参数取值如表1所示。

图 2  修正的Davidenkov模型剪应力-剪应变关系曲线

Figure 2.  Shear stress-shear strain curve for the modified Davidenkov model

下载: 全尺寸图片 幻灯片

表 1  土体参数

Table 1.  Soil parameters

土层名称	最大动剪切 模量Gmax/kPa	泊松比	参数A	参数B	参考剪应变γr	参数β	参考围压/kPa	参数C1	参数C2	门槛剪应变γth	参数m	参数n
④1粉质黏土混粉土层	58 644	0.4	1	0.4595	0.000 377	0.003	1 000	1.051	0.143	0.000 2	0.345	668.9
④2粉土层	123 032	0.4	1	0.4595	0.000 508	0.003	1 000	0.921	0.163	0.000 2	0.345	668.9
⑤1粉细砂层	176 400	0.4	1	0.4595	0.000 478	0.003	1 000	0.813	0.185	0.000 2	0.345	668.9
⑤2细砂层	202 752	0.4	1	0.4595	0.000 541	0.003	1 000	0.813	0.185	0.000 2	0.345	668.9
⑥1中粗砂层	219 978	0.4	1	0.4595	0.000 617	0.003	1 000	0.813	0.185	0.000 2	0.345	668.9

下载: 导出CSV 

| 显示表格

钢筋和螺栓采用理想弹塑性模型模拟，通过ABAQUS软件中嵌入的混凝土塑性损伤本构（CDP）可模拟混凝土在拉压过程中由于塑性行为导致的塑性应变积累、刚度退化等问题，模型中混凝土本构及钢筋、螺栓参数如表2、表3所示。

表 2  混凝土参数

Table 2.  Concrete parameters

混凝土强度 等级	密度/ （kg·m−3）	泊松比	弹性模量/GPa	膨胀角/（°）	偏心率	双轴与单轴抗压 强度之比fb0/fc0	屈服形态影响 参数K	黏性系数
C30	2 300	0.2	30	30	0.1	1.16	0.667	0.000 5
C60	2 500	0.2	36	38	0.1	1.16	0.667	0.0005

下载: 导出CSV 

| 显示表格

表 3  钢筋与螺栓参数

Table 3.  Rebar and bolt parameters

材料	密度/（kg·m−3）	弹性模量/GPa	泊松比	屈服强度/MPa
HRB400钢筋	7 800	200	0.3	400
10.9级螺栓	7 800	200	0.3	900

下载: 导出CSV 

| 显示表格

1.4   输入地震动

根据地质勘察资料，50年超越概率10%的地表地震动峰值加速度为0.118 g，100年超越概率10%的地表地震动峰值加速度为0.206 g。选取同次地震同个台站水平向（地震动1）和竖向（地震动2）地震动记录，如图3所示，地震动峰值加速度均调整为0.1 g，傅里叶谱中灰线为利用0.4 Hz的Parzen窗对所选加速度时程傅里叶谱进行平滑处理的结果，计算时将地震动峰值调幅为0.1 g、0.2 g和0.3 g进行输入。考虑仅沿管廊横向、仅沿管廊纵向及沿管廊横向、纵向作为输入模式。管廊横向为x轴方向，纵向为z轴方向，当输入横向、纵向地震动时，打开对应方向的位移约束，同时在边界网格结点上施加初始状态时的约束反力，防止模型动力计算初始时刻不平衡应力导致土的膨胀效应。单独分析管廊横向或纵向输入地震动时，采用地震动1；施加双向地震动时，横向采用地震动1，竖向采用地震动2。

图 3  输入地震动时程和傅里叶谱

Figure 3.  Input time history of ground motion and Fourier spectrum

下载: 全尺寸图片 幻灯片

2.   结果与分析

2.1   综合管廊应力

管片内力分布与地震损伤分析对研究管廊抗震性能尤为重要，不同形式的内力会导致管片在地震动作用下表现拉压、弯曲、扭转等变形，过大的变形会导致管片环间张开及开裂破坏等问题，因此在管片四周和45°方向布置监测点（图4），监测管片应力、轴力、弯矩及剪力等内力响应。为确保在指定位置得到应力结果，对管片网格划分进行设计，以管片中心为原点，从拱顶开始每隔5°进行划分，确保网格足够精细。

图 4  管廊监测点布置

Figure 4.  Utility tunnel monitoring points

下载: 全尺寸图片 幻灯片

管片监测点在各工况下的应力分布包络曲线如图5所示。由图5可知，管片在地震动作用下，拱腰与拱肩受力较大，管片应力包络线整体呈对称状，随着地震强度的增加，管片在45°方向的应力增长更明显，而在拱顶与拱底未出现明显的应力增长。在纵向地震动作用下，管廊应力增长并不明显，这主要由管廊与土的纵向摩擦作用引起，对管廊横断面应力影响较小；在横向地震动作用下，管廊应力增加明显，双向地震动作用下增长幅度更大，应力峰值为16.4 MPa，出现在右拱腰，这说明横向地震动作用使管廊应力增长较明显，管廊在45°共轭方向与拱腰部分更易出现张开与损伤。

图 5  管廊横断面应力分布（Mises应力）

Figure 5.  Utility tunnel cross-sectional stress distribution （Mises stress）

下载: 全尺寸图片 幻灯片

弯矩、剪力和轴力响应最大值出现时刻约在第10 s，且10 s附近响应值相差较小，故取10 s作为特定时刻进行分析。峰值加速度为0.3 g时，管片在10 s时的内力如图6所示。由图6可知，仅考虑单向地震动时，弯矩分布较均匀，管片整体呈负弯矩受力情况。3种工况条件下管片轴力与剪力分布区别较小，管片轴力最大值位于右拱脚附近，且管片横截面45°共轭方向出现拉力区，这部分管片易出现横向张开与拉裂，剪力最大值位于右拱腰附近。由此可知，管片结构在地震过程中发生剪切拉压变形时，在拱腰与45°共轭方向易出现内力增大情况，且可能出现拉力区，造成管片开裂。

图 6  管廊横断面内力分布

Figure 6.  Internal force distribution in the cross section of the utility tunnel

下载: 全尺寸图片 幻灯片

2.2   管片张开量

管片张开量是评估综合管廊损伤的重要指标，当环间张开量过大时，易造成管廊漏水、进砂等灾害。选取管片易张开的拱顶与拱底作为监测点，9种工况下管片在10 s时的纵向环间张开量如图7所示。由图7可知，随着输入地震动峰值加速度的增大，综合管廊纵向环间张开量逐渐增加，当输入同一地震动时，管廊拱顶和拱底均在5、6环间的纵向张开量最小，拱顶最大张开量基本出现在7、8环间，而拱底最大张开量分布无明显规律；在相同的峰值加速度下，拱顶监测点中出现最大张开量的工况为沿隧道纵向输入地震动的工况，而拱底监测点中出现最大张开量的工况为双向输入地震动的工况，且拱顶张开量大于拱底。通过比较同一峰值加速度水平下地震动输入方向对管片张开量的影响可知，沿纵向输入的地震动对管片张开量的影响更显著，在仅输入横向地震动工况中，各监测点环间张开量区别并不显著，由此可知，沿管廊纵向的地震动更易造成管廊纵向张开。

图 7  管片纵向环间张开量

Figure 7.  The amount of opening between the longitudinal rings of the utility tunnel

下载: 全尺寸图片 幻灯片

受限于纵向长度的影响，除研究管片环间张开量最大值外，还需分析纵向大尺度模型环间张开量在地震动作用下的变化情况。以沿管廊纵向和双向输入峰值加速度为0.3 g的地震动为例，管片拱顶与拱底3、4环之间张开量时程曲线如图8所示。由图8可知，地震动未达峰值前，环间张开量无明显变化；地震动达到峰值附近时，环间张开量增长明显，并在一定范围内波动；地震动达峰值后，随着地震动的减弱，环间张开量未明显减小，甚至在峰值加速度为0.3 g的双向输入工况中发现拱顶环间张开量有所增大。在整个地震动输入过程中，拱顶环间张开量一般大于拱底，这与图7所示结果保持一致。

图 8  管廊环间张开量时程曲线

Figure 8.  Time history curve of opening amount between utility tunnel rings

下载: 全尺寸图片 幻灯片

2.3   综合管廊内部结构损伤

混凝土塑性损伤模型中弹性模量E的降低与刚度损伤因子d相关，刚度损伤因子是应力状态和单轴损伤变量的函数。对于单轴循环状态，ABAQUS中定义刚度损伤因子与拉、压损伤因子的关系如下：

式中，d_c和d_t分别为压损伤因子和拉损伤因子，s_c和s_t分别为压、拉应力状态的函数，以模拟与应力反转相关的刚度恢复效应。因此，本文主要将刚度损伤因子作为分析管廊结构损伤的主要指标。不同地震动作用下20 s时刻管廊结构刚度损伤因子云图如图9所示。由图9可知，当峰值加速度为0.1 g时，横向和纵向地震动作用下管廊结构整体未发生受拉损伤，而管廊结构底部框架结构连接位置略受拉损伤，但仍处于弹性工作状态；当峰值加速度为0.2 g时，横向和双向地震动作用下管廊内部结构受拉破坏且贯穿整个截面，而施加纵向地震动时，管廊结构底部框架结构连接位置出现损伤积累；当峰值加速度为0.3 g时，横向和双向地震动作用下管廊内部结构底部框架结构连接位置受拉损伤严重，导致内部结构在两端连接处及中间支撑处连接点出现拉压变形，使重要连接点出现损伤积累，进一步导致整体结构倒塌破坏，而施加纵向地震动时，管廊结构底部框架结构受拉损伤范围和程度不断扩大，但并未贯穿整个横截面。综上所述，沿管廊横向的地震动作用对内部结构的影响较大，内部结构出现了挤压破坏，而沿管廊纵向的地震动作用引起了内部结构与管廊变形不协调，进一步导致了管廊结构受拉破坏。

图 9  管廊结构损伤因子云图

Figure 9.  Damage factors of utility tunnel structure

下载: 全尺寸图片 幻灯片

3.   结语

以苏通GIL综合管廊为研究背景，建立精细化三维综合管廊及内部结构有限元模型，考虑上覆水的影响，分析不同地震动输入下综合管廊地震反应特性，得出以下主要结论：

（1）综合管廊外部衬砌在地震动作用下，拱腰与拱肩受力较大，且随着地震动强度的增大而增大。管片轴力及剪力受地震动输入方向的影响较小，轴力最大值出现在拱底和右拱脚位置，剪力最大值位于右拱腰附近。

（2）相同工况条件下，综合管廊拱顶张开量一般大于拱底，纵向地震动对管廊环间张开量的影响较大，当地震动达到峰值附近时，张开量增长明显，并在一定范围内波动。

（3）随着地震动强度的增加，综合管廊内部结构损伤加剧。当地震动峰值加速度为0.3 g时，管廊结构内部出现了严重损伤且变形较明显，总体上横向地震动对管廊结构破坏起主要作用。