单仓矩形综合管廊地震响应三维有限元分析及试验验证

吴 立，苏 勇，陈 军，韩佳欣，毛庆威，王丹生

(1. 湖北省城建设计院股份有限公司， 湖北 武汉 430051；2. 华中科技大学 土木与水利工程学院， 湖北 武汉 430074)

地下综合管廊是一种典型的浅埋型地下空间结构，它能够容纳排水管道、电信电缆、燃气管道等多种生命管线，可很大程度地提升城市管线管理效率。综合管廊设计寿命长，安全性要求高，一旦在强震作用下发生严重破坏，就会使城市的用水、用电、通信陷入瘫痪，造成极大的经济损失与社会影响。如1995年的日本阪神地震就造成大量的地面建筑及地下结构的破坏，经济总损失高达1000亿美元[1]。1999年的集集地震、2008年的汶川地震也使得当地许多地下结构遭受了不同程度的破坏。近年来，我国地下空间开发和地下结构建设方兴未艾，据不完全统计，我国地铁隧道运营里程达5000多km，综合管廊已建设里程约1700 km，今后几年地下综合管廊建设规模将有望达到1万km。因此，对我国地下结构开展深入的抗震研究显得更加重要。

同济大学李杰教授团队是国内较早对管廊抗震性能进行研究的团队，他们设计了一个整体管廊缩尺模型振动台试验，对单室矩形管廊在地震作用下的特征点加速度响应、钢筋应变以及结构内力等方面进行了数值与试验研究[2]。重庆大学的仉文岗等[3]利用振动台试验对双仓矩形地下综合管廊缩尺模型最不利时刻的弯矩特征进行了研究，并利用有限元软件ABAQUS建立了一个土体 - 管廊结构的二维模型进行数值模拟。重庆大学的冯立等[4]则通过振动台试验，在考虑接缝影响的基础上，研究了单室矩形综合管廊在不同种类、不同峰值加速度地震波激励下的地震响应。蒋录珍等[5]采用数值分析的方式，研究了不同的土体材料、结构特性以及土体与结构之间的接触摩擦等因素对结构地震响应的影响。从上述研究可知，目前国内学者已经开展了部分管廊结构地震响应分析的有限元数值模拟研究，但为了提高模型的收敛性和减少计算成本，大多局限在二维有限元建模与分析，且均未考虑管廊内部钢筋的影响。

本文采用有限元软件ABAQUS，对配置有钢筋的单室矩形综合管廊进行了三维有限元地震响应分析，研究了在不同强度地震波作用下管廊特征点的应力响应、加速度响应、管廊整体损伤情况以及内部钢筋的应力响应，并通过已有的振动台试验结果来验证本文管廊有限元模型以及数值分析的可靠性，最后依据计算结果为实际管廊工程的抗震设计提出了一些合理可行的建议。

1 单室矩形管廊三维有限元模型

1.1 土体 - 管廊三维有限元模型参数

本文依据史晓军等的管廊振动台试验[6]建立的土体 - 管廊三维有限元模型整体尺寸为4 m×1.8 m×2.37 m，其中综合管廊部分尺寸为0.6 m×0.6 m×1.3 m，壁厚为37.5 mm。本文管廊结构建模过程中考虑了钢筋影响，相比于素混凝土综合管廊模型更能反映实际情况中综合管廊的真实受力状态。土体、混凝土、钢筋的材料参数分别见表1，2，3。

表1 土体材料参数

表2 混凝土材料参数

表3 钢筋材料参数

为了让振动箱与土体一致受力，本文模型将振动箱与土体设置为了一个整体。起减小箱壁反射效应作用的聚氯乙烯泡沫板部分采用无限单元作为人工边界，这样可以在一定程度上减弱其边界效应，从而更好地模拟出实际工程中地下结构所处的半无限空间体。振动箱壁的质量施加到土体与无限元单元接触部位，模拟由振动箱产生的惯性力[7]。地震波采用El-Centro波，以加速度输入方式进行激励，其峰值加速度为3.26 m/s2。土体 - 管廊整体模型尺寸、管廊模型尺寸、钢筋笼模型、模型网格剖面图、El-Centro波加速度时程曲线分别如图1～5所示。

图1 土体 - 管廊整体模型/m

图2 矩形管廊模型/m

图3 钢筋笼模型

图4 整体模型网格剖面

图5 El-Centro波加速度时程曲线

在本文建立的三维模型中，土体计算区域采用C3D8单元，边界区域采用CIN3D8单元，综合管廊混凝土部分采用C3D8单元，钢筋采用梁桁架T3D2单元；所涉及到的相互作用包括综合管廊与土体的摩擦接触，纵筋与箍筋为绑定约束，整个钢筋笼则是以内置区域的方式嵌入混凝土中。

1.2 材料本构模型1.2.1 混凝土

混凝土采用混凝土塑性损伤(Concrete Plastic Damage，CDP)模型。ABAQUS软件中的CDP模型是基于Lubine等[8,9]提出的材料损伤应变等效原理发展而来。这种模型假设材料在发生损伤后，其应变等效为无损伤应变，因此可以得出损伤后混凝土的应力应变关系式为：

(1)

σ=E0(1-D)ε

(2)

式中：σ为正应力；ε为应变；E0为无损混凝土的弹性模量；Ed为损伤混凝土的弹性模量；D为损伤参数。GB50010-2010《混凝土结构设计规范》[10]中单轴受力混凝土的损伤参数D是基于名义应变，而ABAQUS输入的损伤因子d是基于非弹性应变，因此需要进行转换。根据田连波等[10]提出的损伤因子取值理论，当混凝土出现损伤后再受力时，其产生的弹性势能与相同荷载下无损伤时的混凝土弹性势能在形式上相同，此时只需将应力转变为等效应力即可，或者直接采用出现损伤时的弹性模量即可[11]。则有：

(3)

(4)

结合GB50010-2010《混凝土结构设计规范》可得：

E0(1-D)=E0(1-d)2

(5)

从而得出转换公式：

(6)

1.2.2 土体

土体选用Mohr-coulomb模型，其遵循的准则为剪切破坏准则：

τf=c+σtanφ

(7)

式中：τf为土体的剪应力；c，φ分别为粘聚力与内摩擦角。Mohr-coulomb模型的剪切屈服面函数为：

F=Rmcq-ptanφ-c=0

(8)

式中：q为等效压应力；p为Mises等效应力；Rmc为控制屈服面在π平面的形状，其表达式为：

(9)

(10)

式中：Θ为极偏角；r为第三偏应力不变量。

1.2.3 钢筋

由于一般结构破坏时，钢筋的应变尚未进入强化阶段，使用理想弹塑性模型就已经足够精确，因此本文中钢筋的本构模型采用理想弹塑性模型，其应力 - 应变曲线如图6所示。

图6 钢筋应力 - 应变曲线

塑性材料本构方程为：

(11)

式中：σs为钢筋应力；εs为钢筋应变；Es为钢筋的弹性模量；fy为钢筋的屈服强度代表值；εy为与fy相对应的钢筋屈服应变，本文取0.15%。

2 单室矩形管廊地震响应分析

2.1 管廊特征点应力响应分析

本节选取峰值加速度为0.1g，0.2g，0.4g的El-Centro波三种工况，对管廊端部截面特征点应力响应进行研究。特征点分布如图7所示，在整个地震作用过程中各特征点第一主应力和第三主应力的最大值分别见表4，5。

图7 管廊端部截面特征点分布

表4 特征点第一主应力最大值 MPa

由以上两表可以看出：

(1)依照GB50010-2010《混凝土结构设计规范》，C30混凝土的抗压强度设计值为14.3 N/mm2，抗拉强度设计值为1.43 N/mm2，在三种工况下管廊特征点的压应力均远小于14.3 N/mm2，故管廊结构不会出现压裂的情况；而结构的拉应力除了0.1g地震波工况中各特征点的拉应力均小于1.43 N/mm2以外，另外两种工况下管廊结构中均有多个特征点的拉应力值超过1.43 N/mm2，这说明在0.2g，0.4g地震波作用下混凝土因拉应力过大而产生较大塑性变形甚至局部破坏。

(2)当地震波强度较小时，管廊底板仍处于弹性变形阶段，拉应力比其他部位大；当地震波强度较大时，底板已经发生塑性变形，此时增大地震强度后拉应力值变化很小。

(3)随着地震波强度的增大，侧板跨中点a4的拉应力值增长幅度明显小于顶板和底板特征点，且始终小于C30混凝土的抗拉强度设计值，因此在实际工程中可适当减少侧板跨中处的配筋以降低费用。

2.2 地震结束时管廊受拉损伤分析

由上述分析我们可以看到，在0.1g地震波工况下管廊内部最大的拉应力和压应力较小，结构并未产生破坏，而在0.2g，0.4g地震波工况下管廊内部产生的拉应力较大，而压应力仍然远小于其抗压强度设计值，可见结构会发生受拉破坏。为了研究在0.2g，0.4g地震波作用后管廊的具体受拉损伤状态，本节主要对这两种工况中地震作用结束时的受拉损伤云图进行分析，受拉损伤云图如图8所示。

图8 两种工况地震波结束时刻受拉损伤云图

由图8可以看出：

(1)在两种工况中，地震作用结束时管廊底板的受拉损伤情况最为严重，尤其是在底板端部靠近腋角的局部红色区域，0.2g地震波工况受拉损伤系数最大为0.2188，0.4g地震波工况受拉损伤系数达到了0.4938，这些区域在地震作用下容易产生裂缝甚至破坏，因此在实际工程中可适当加密配筋；

(2)随着地震波强度的增大，管廊的各个角点和顶板的受拉损伤系数也出现了明显增长，尽管这些部位的损伤情况没有底板严重，但在设计抗震设防等级较大地区的管廊时，顶板和角点区域应适当进行加固。

2.3 地震结束时管廊内钢筋Mises应力分析

在分析了管廊受拉损伤情况之后，本节将进一步对管廊内部钢筋在地震结束时的应力进行分析，0.2g，0.4g地震波这两种工况中地震作用结束时的钢筋Mises应力云图如图9所示。

图9 两种工况地震波结束时刻受拉损伤云图

由图9可以看出：

(1)在0.2g地震波工况下，地震作用结束时钢筋最大应力为85.03 MPa；在0.4g地震波工况下，地震波作用结束时钢筋最大应力为299.5 MPa。而模型中钢筋的抗拉强度为716 MPa，说明这两种工况下的钢筋均未达到屈服，且仍有较大的承载潜力。

(2)在这两种工况下，底板部位的钢筋应力均明显大于其他部位，这是由于底板混凝土内拉应力较大甚至出现局部拉裂的情况，因此其内部钢筋的拉应力也会明显高于其他部位，这也与上述分析结果一致。

2.4 管廊特征点加速度响应分析

本节选取El-Centro波峰值加速度分别为0.1g，0.2g，0.4g的三种加载工况，对整个加载过程中管廊特征点的最大水平加速度进行分析。特征点的选取与2.1节相同。各特征点的最大水平加速度响应见表6，侧板特征点的最大水平加速度响应如图10所示。

表6 特征点最大水平加速度 ×g

图10 三种工况下侧板特征点最大水平加速度

由表6和图10可以看出：

(1)在以上三种地震作用工况中，同一水平高度的特征点最大水平加速度值基本相同，这与施加的地震波为水平方向有关。

(2)从左侧板和右侧板特征点最大水平加速度图可以看出，越接近底板的特征点最大水平加速度越大，说明越靠近底板的区域受到地震的影响越大，这与我们之前进行应力分析以及损伤分析的结果是一致的。

3 试验验证

为了验证上述地下综合管廊三维有限元建模及地震响应分析的可靠性，本文选取了同济大学史晓军等[6,13]的管廊振动台试验结果来进行检验。他们试验采用的试件、土体的尺寸及材料属性均与本文相同，试验中各部分传感器的布置详见文献[6]。

3.1 管廊与土体特征点加速度响应对比

选取文献[6]中1.0g地震荷载施加到本文建立的管廊有限元模型上，分别将有限元模型顶部输出点、底部输出点的加速度响应与文献[6]中的试验结果进行对比。输出点的位置及该处加速度响应时程对比分别如图11，12所示。

图11 土体与管廊输出点位置/mm

图12 不同输出点加速度响应结果对比

对比管廊有限元与试验加速度响应曲线可知：

(1)数值模拟的特征点加速度时程曲线与试验结果吻合较好，其相位同步，幅值相差不大。

(2)无论是数值模拟结果，还是试验结果，都反映出管廊特征点的加速度幅值与土体基本相同，这是由于管廊结构受周围土体的约束，其加速度响应应该服从于周围土体的加速度响应。

(3)数值模拟及试验结果还反映出，管廊底部的加速度峰值比管廊顶部的加速度峰值大，这说明在地震作用下管廊底部比顶部更容易受到破坏，这也验证了上述加速度响应分析结果的可靠性。

3.2 钢筋笼特征点应变响应对比

选取钢筋笼中部位置，将文献[6]中1.0g地震波工况下钢筋笼上24个输出点的应变数据通过Digitizer进行识别(识别误差不超过4%)，将识别的试验结果与数值模拟结果进行对比，钢筋笼各部分的应变结果对比如图13所示。

将数值模拟结果与试验结果进行对比可以看出：

(1)侧板模拟结果比顶板和底板模拟结果好，除了右侧内层钢筋的应变数据外，其余侧板上的输出点钢筋应变数据与试验结果走势一致，但应力幅值仍存在一定差异。

(2)钢筋笼外侧部分应变模拟结果好于内侧部分，其中左侧外层钢筋、底板外层钢筋的应变趋势与幅值均与试验结果吻合度较高。

(3)数值模拟的钢筋应变幅值普遍大于试验结果。

(4)试验结果和部分模拟效果较好的数值模拟结果均反映出，角点处的钢筋应变幅值大于跨中处的钢筋应变幅值，说明在角点处产生的应力更大，因此在实际工程设计管廊时可适当在角点处加密或加粗钢筋，也可以考虑将管廊角点处设计成弧形。

图13 应变模拟结果对比

3.3 误差分析

通过分析振动台试验的试验过程与数值模拟过程的差异，可以得到以下可能导致误差的原因：

(1)数值模拟中选取的土体模型参数无法完全表现出真实土体材料的复杂性，因此与试验中的实际土体性质存在一定差异。

(2)有限元模拟对振动台中层状剪切箱进行了一定的近似，与具体情况仍有一定的差别。

(3)原试验中部分传感器破坏无法正常工作，而且在试验中难免存在一定的试验误差，使得测试结果与理想结果有一定的偏差。

4 结 语

本文建立了考虑钢筋影响的单室矩形综合管廊三维有限元模型，对不同强度地震波作用下综合管廊结构的动力响应进行了研究，并利用一个已有的振动台试验结果验证了本文有限元数值分析的可靠性。综合已完成的工作，得出以下结论：

(1)管廊在三种不同强度地震波工况下的最大压应力均远小于混凝土抗压强度设计值，因此不会出现压裂现象；而在0.2g，0.4g地震波作用下均会因拉应力过大产生塑性变形甚至可能拉裂。

(2)管廊底板在地震作用下的损伤程度大于其他部位，其次是顶板及角点处，这些部位在实际工程中可适当加粗或加密配置钢筋，角点处可考虑进行平滑处理，如设计成弧形，以减小应力集中。

(3)在不同强度地震波工况下，管廊内部钢筋在地震作用结束时同样表现为底板部位拉应力较大，但仍处于弹性范围内。

(4)数值分析结果和试验结果进行对比后吻合程度较高，验证了数值分析方法的可靠性。