圆形隧道横断面地震响应简化解析方法的数值验证

赵 密，邵伟昂，黄景琦，杜修力，王 媛

(1. 北京工业大学城市与工程安全减灾教育部重点实验室，北京 100124；2. 北京科技大学土木与资源工程学院，北京 100083)

隧道作为我国山区的生命线工程，在国家基础设施建设中占有重要作用。由于隧道埋置于土体或岩体介质中，相较于地上结构，以往的研究者认为隧道结构具有较好的抗震性能。但在近几次大地震中，如Kobe(1995)、Chichi(1999)和Wenchuan(2008)等地震，很多的隧道结构遭受严重破坏，造成了大量的人员伤亡及财产损失[1−3]。因此，在隧道结构安全性设计中除需考虑周围岩土体的静力荷载外，还需考虑地震带来的附加地震荷载。对隧道设计者而言，简单有效的隧道内力及变形计算方法对隧道抗震设计尤为重要[4−5]。

目前隧道抗震设计主要有以下三种方法：数值分析方法、模型试验方法和简化解析方法。数值分析方法包括动力时程分析方法[6− 7]及拟静力简化分析方法(如整体式反应位移法[8]、Pushover分析方法[9]等)，其优点为可考虑众多复杂条件，但对设计者的数值计算理论和软件操作能力要求较高；模型试验[10]通过相似比设计可很好的反应真实情况，但周期长、费用高，较难在实际工程设计中大规模采用；相较于前两种方法，简化解析方法由于能够给出隧道地震内力及变形的计算公式，可快速对隧道地震反应进行计算，因此在隧道结构的初步设计中得到广泛应用。Hoeg[11]提出了地震荷载作用下圆形隧道轴力、弯矩和位移的解析解；Einstein 等[12]推导了无穷远处剪应力作用下圆形地下结构的应力和变形；Wang[13]考虑了土-结构相互作用，分别给出了圆形和矩形地下结构的内力求解公式；Penzien[14]基于梁单元假设，提出了一种矩形隧道和圆形隧道内力及变形的解析方法；Bobet[15− 16]考虑隧道处于饱和沉积物中并受到快速加载的地震荷载，分别针对排水与不排水条件推导了圆形隧道内力的计算公式。另外，Park 等[17]通过引入接触刚度系数，推导了能够考虑土体与地下结构不同接触条件下的静力解析解。

由于能够简便快捷的计算隧道衬砌的抗震设计内力，简化的解析方法常在初步设计中被用于估算隧道的内力与变形。然而上述各简化解析方法的求解公式并不完全相同，不同解析方法对隧道内力与变形的计算结果甚至出现相差较大的情况。因此，需对此类简化解析公式的准确性进行必要的对比验证。Hashash 等[18]指出在土体与结构接触完全滑动的情况下，Wang 和Penzien 的方法计算结果一致，而当土体与结构处于完全绑定的情况下，Wang 的计算结果明显高于Penzien 的计算结果；Sedarat 等[19]在无滑移条件下通过将Wang和Penzien 的解析解与拟静力有限元的数值结果进行对比，得出Penzien 的解析解严重低估了隧道轴力；Kontoe 等[20]同样通过拟静力有限元数值模拟，考虑不同的接触状态及不同的土-结构相对刚度，对Wang、Penzien、Bobet 和Park 等四种方法进行了较为系统的验证。在国内，晏启祥等[21]、宾佳[22]对简化解析方法也进行了相关的验证。上述验证工作通常采用拟静力有限元模型对简化解析方法进行数值验证，将地震荷载等效为均匀的剪应变或剪应力，并忽略结构惯性效应及隧道对地震波的散射作用。而实际地震作用下的隧道结构响应为一种动力响应过程，需考虑地震传播的时空效应、结构的惯性效应及散射场作用[23−24]。另外，简化解析方法假定结构埋置于无限岩土体中，未考虑实际工程中地表对浅埋结构的影响。因此，采用理论上相对较为完备的近场波动有限元方法，开展不同场地土类别及不同埋深下简化解析方法对隧道内力预测精度的验证，给出简化解析方法的适用条件及计算精度具有重要的工程意义。

本文采用近场波动有限元方法对隧道地震响应进行模拟求解。其中，以黏弹性人工边界条件模拟无限域的辐射阻尼效应，以等效节点力方法模拟地震动的波动输入。最终，以数值模拟结果为标准解，对Wang、Bobet 和Park 等3 种经典简化解析方法在不同的场地土类别及不同埋深下的精度及适用条件进行研究。

1 圆形隧道横断面地震响应简化解析方法

当圆形隧道横断面受竖向入射剪切波作用时，简化解析方法通常假设隧道-围岩土系统的周边远场地层处于纯剪状态，如图1(a)所示。圆形隧道此时的变形为一种椭圆形的剪切变形，如图1(b)所示。此时隧道内力取决于远场剪应力的大小及围岩土-结构的相对刚度。为描述围岩土与结构之间的相对刚度，Wang[13]通过采用压缩比C和柔度比F两个参数对隧道结构与周围围岩土介质的相对刚度进行量化。隧道衬砌的内力与变形是压缩比C和柔度比F的函数。压缩比C和柔度比F的表达式为：

式中：下标m 和l 分别代表围岩土和隧道衬砌；E和ν分别表示介质的弹性模量与泊松比；r和t分别表示隧道半径及衬砌厚度；I为隧道衬砌的截面惯性矩。

图1 简化解析方法的分析模型Fig.1 Analytical model of simplified analytical method

对于围岩土与结构的接触条件，简化解析方法中往往假定为完全无滑移和全滑移两种接触条件，且两种接触条件预测内力的结果往往差异较大[20]。实际工程中围岩土介质与隧道衬砌接触面较为粗糙，不可能处于完全滑移状态，因此全滑移接触的假设与实际情况差距较大。杜修力等[25]研究结果表明在小震情况下接触面粗糙程度对地下结构横断面地震响应影响并不显著，因此本文仅讨论无滑移条件下简化解析方法的精度。

1.1 Wang 的解析解

Wang[13]给出了在完全滑移状态及无滑移状态

在无滑移条件下，Wang 未给出弯矩的公式，但建议可采用全滑移下的弯矩公式进行求解。

1.2 Bobet 的解析解

1.3 Park 的解析解

2 近场时域波动有限元方法

在近场时域波动有限元方法中，通常在模型截断边界处设置人工边界条件来模拟无限域的辐射阻尼效应。相对于黏性边界[26]与透射边界[27−28]，黏弹性边界[29− 30]能够较好地模拟无限域地基的弹性恢复性能和辐射阻尼效应，且在有限元软件中方便实现。另外，地震动作为外源荷载通常采用等效节点力的方法进行施加[31−32]。因此，考虑人工边界条件的近场时域波动有限元的控制方程可写为[24]：

式中：下标I和B分别表示内部节点和人工边界节点；U、U˙ 和U¨分别为位移、速度和加速度；M、C和K分别为质量矩阵、阻尼矩阵和刚度矩阵；FB为截断边界节点上的地震荷载向量；FI为内部节点上的荷载向量。子矩阵CB和KB为黏弹性人工边界的对角矩阵。对于人工边界上的特定节点l，对角矩阵[CB]l和[KB]l的对角项为黏弹性人工边界的弹簧-阻尼元件参数，其表达式为[29]：

3 简化解析方法的数值验证

3.1 时域波动有限元分析模型建立

采用通用有限元程序ABAQUS[33]，建立圆形隧道横断面的二维有限元分析模型，如图2 所示，模型宽度为120 m。为讨论隧道埋深的影响，隧道拱顶至地表的距离H分别设定为1 m、3 m、5 m、7 m、10 m、17 m、27 m、37 m、47 m、57 m、67 m 和77 m。隧道衬砌外半径为3 m，衬砌厚度为0.3 m。在模型底部和两侧施加黏弹性人工边界条件来模拟无限域的辐射阻尼效应。具体方法为在人工边界面所有单元节点上，分别在法向和切向上施加并联的弹簧单元与阻尼单元。弹簧单元与阻尼单元的参数通过式(17)进行求解。围岩土体介质及隧道衬砌均采用平面应变单元模拟。围岩土体介质与隧道衬砌的接触面采用完全绑定约束以模拟无滑移接触条件。

图2 计算模型Fig. 2 Computation model

岩土体和隧道衬砌均采用线弹性本构模型，隧道衬砌密度为2500 kg/m3，泊松比为0.2，弹性模量为24.8 GPa。假定结构周围场地为均匀、各向同性介质。为研究场地中不同岩土体类别对简化解析方法精度的影响，参考《城市轨道交通结构抗震设计规范》(GB50909−2014)[34]，选取软弱土、中软土、中硬土、坚硬土或软岩和岩石等五种类别的岩土体，场地密度为2500 kg/m3，泊松比为0.25，每种岩土体类型选取2 组参数。岩土体介质的弹性模量、剪切波速及相应的柔度比见表1所示。

表1 场地岩土体参数Table 1 Parameters of surrounding ground

输入地震动为1995 年日本阪神地震中的Kobe波，其加速度时程曲线见图3 所示。截取能反映Kobe 波地震特性的前20 s 进行输入，加速度峰值6.114 m/s2。地震波作为剪切波从模型底部开始输入，传播方向沿竖向垂直传播。

图3 输入的Kobe 波加速度时程曲线Fig. 3 Acceleration time-history of the incident Kobe waves

3.2 简化解析方法的精度验证

3.2.1 隧道轴力Tmax

图4 给出了不同埋深及不同柔度比下隧道峰值轴力的数值模拟结果与各简化解析方法的计算结果。从图4 中的数值解可以看出，随着柔度比增加，即随着场地岩土体逐渐变硬，隧道轴力具有先增加后减小的规律。

三种简化解析方法具有相同的轴力预测值，且预测结果与数值模拟结果具有相同的变化趋势。对比数值结果与三种解析结果可以看出，在埋深较浅情况下(即H≤10 m)，小柔度比下的解析解与数值解差距较大，解析解严重高估了隧道轴力，但随着柔度比的增加，解析解与数值解的差距逐渐减小。对于埋深较深的情况(即H>10 m)，解析解的轴力小于数值结果，说明此时解析解低估了隧道轴力。另外，解析解在较小柔度比和较大柔度比情况下与数值解均吻合较为良好，但在中间柔度比下两者吻合程度稍差，即大埋深情况下，解析解与数值解吻合程度随柔度比增加具有先减小后增加的规律。

为定量描述简化解析方法预测轴力与数值解轴力之间的相对误差，现定义误差R如下：

式中：TmNax为数值解；TmAax为简化解析解。

由于Wang、Bobet 和Park 三种方法预测轴力相同，其误差值R也相同。表2 中给出了不同柔度比(即不同的场地土类别)及不同埋深下的轴力误差R。从表2 中可以看到，当隧道埋深较浅时，解析解预测的轴力误差往往较大，且误差大小随柔度比的增加而减小，即场地土越硬则预测误差越小。当隧道埋深较深时，解析解误差明显小于埋深较浅的情况。表2 中灰色值均为误差R超过15%的情况。若以误差小于15%为可接受的标准，简化解析方法在软弱土、中软土、中硬土、坚硬土或软岩和岩石五类场地土中用于预测隧道轴力的隧道埋深需分别大于5 m、3 m、3 m、1 m和1 m。

图4 不同埋深H 及不同柔度比F 下数值解与各简化解析解的轴力Fig. 4 Tunnel thrusts by numerical modelling and analytical solutions under different depths H and flexibility ratios F

表2 不同埋深H 及不同柔度比F 下简化解析解的轴力误差R /(%)Table 2 Errors R of these three analytical solutions for tunnel thrusts under different depths H and flexibility ratios F

3.2.2 隧道弯矩Mmax

不同埋深及不同柔度比下隧道的峰值弯矩如图5 所示。从图5 中的数值解可以看出，隧道处于软弱土和中软土中的弯矩明显大于隧道处于岩石场地中的情况。当埋深H≤17 m 和H≥67 m时，隧道轴力随柔度比的增加而减小。而当埋深处于17 m

图5 不同埋深H 及不同柔度比F 下数值解与各简化解析解的弯矩Fig. 5 Tunnel bending moments by numerical modelling and analytical solutions under different depths H and flexibility ratios F

三种简化解析方法预测的隧道轴力随柔度比和埋深的变化与数值解具有相同的变化趋势。Bobet 方法与Park 方法预测的弯矩值相同，但小于Wang 方法的预测值。三种方法的预测弯矩均大于数值解，说明三种解析方法都高估了隧道弯矩，其中Wang 方法对弯矩预测的偏差最为严重。另外，对于同一埋深情况，随着柔度比增加三种简化解析方法的弯矩与数值模拟结果逐渐接近。在柔度比较大的情况下，如坚硬土场地和岩石场地，简化解析解与数值解基本重合。

为定量描述简化解析方法在预测弯矩时的误差，采用式(19)中定义的误差R进行描述。表3中给出了不同柔度比及不同埋深下的Wang 方法的弯矩误差R。表4 中给出了相应的Park 方法和Bobet 方法的弯矩误差R。从表3 中可以看出，Wang 方法的误差在大多数情况下均较大，仅在岩石场地中误差相对较小。因此，在实际设计中采用Wang 方法对弯矩进行计算并不合适。

表3 不同埋深H 及不同柔度比F 下各Wang 解析解的弯矩误差R /(%)Table 3 Errors R of Wang’s analytical solution for tunnel bending moments under different depths H and flexibility ratios F

相对于Wang 的方法，表4 中Bobet 方法和Park 方法的误差明显较小。与表2 中轴力的规律相似，浅埋情况下的弯矩误差大于深埋的情况。随着柔度比增加，即场地土逐渐变硬时，Bobet 与Park 的误差整体具有减小的趋势。在岩石类场地中，即使隧道拱顶距地表仅为1 m 时，弯矩误差仍然处于10%以内。表4 中灰色值均为误差超过15%的情况。同样以误差小于15%为可接受误差，简化解析方法在软弱土、中软土、中硬土、坚硬土和岩石等五类场土中采用Bobet 方法或Park 方法预测隧道弯矩，隧道埋深需要分别大于10 m、5 m、5 m、1 m 和1 m。

表4 不同埋深H 及不同柔度比F 下Bobet 与Park 解析解的弯矩误差R /(%)Table 4 Errors R of Bobet’s and Park’s analytical solution for tunnel bending moments under different depths H and flexibility ratios F

4 结论

本文采用时域波动有限元方法，研究了3 种经典简化解析方法对圆形隧道横断面地震内力计算的准确性。无滑移条件下，通过对不同场地土类别及不同埋深情况下的对比验证，得到如下结论：

(1) Wang、Bobet 与Park 三种简化解析方法对隧道轴力的预测值相同。在埋深较浅情况下(即H≤10 m)，小柔度比下的解析解与数值解差距较大，解析解严重高估了隧道轴力，但随着柔度比增加，解析解与数值解差距逐渐减小。在埋深较深情况下，解析方法预测轴力小于数值计算轴力，随柔度比的增加，解析解与数值解的吻合程度具有先减小后增加的规律；

(2) Bobet 方法与Park 方法预测的弯矩值相同，且小于Wang 方法预测的弯矩。三种方法预测弯矩值均大于数值解，说明三种解析方法都高估了隧道弯矩，其中Wang 方法的预测值偏差最为严重。另外，对于同一埋深情况，随着柔度比增加三种简化解析方法预测弯矩的精度逐渐增加；

(3) 综合考虑对隧道轴力与弯矩的预测精度，相对Wang 的方法，Bobet 与Park 方法更适合用于隧道内力计算。对于软弱土、中软土、中硬土、坚硬土和岩石等五种类别的场地土，隧道埋深分别大于10 m、5 m、5 m、1 m 和1 m 时，Bobet与Park 方法的轴力与弯矩误差均能小于15%。

本文数值模拟中围岩土及隧道结构均采用线弹性模型，而在强震情况下场地土及隧道结构往往会进入一定程度的非线性，对于非线性情况下的验证工作将更加具有工程价值。另外，本文针对不同场地土条件均采用同一条地震记录输入，而在实际地震中，不同场地类别的地震记录往往存在较大差异。因此，针对不同场地土类别，选用相应的地震记录作为输入地震动也将在后续工作中持续完善。