0.   引言

近年,中国海底隧道建设发展迅速,建成或拟建的达20多条。海底隧道通常需在海岸设置竖井,安装隧道内部排水、供电、通风换气、防灾救援、内外气压调节等重要设施。受强震袭击时,海底隧道与竖井连接部位易发生应力集中和变形过大问题,是抗震的薄弱部位。由于现场调查困难,强震引起的海底盾构隧道破坏的报道极少,现有研究结果也存在争议。1985年墨西哥Michoacán大地震,报道了墨西哥城盾构污水隧道遭受严重的结构破坏^[1]。赵武胜等^[2]报道了1985年墨西哥大地震和1995年阪神地震中盾构隧道与竖井连接部位的震害现象,如管片错动、环间螺栓剪断、拱顶轴向裂缝、漏水等。振动台模型试验^[3]、数值模拟^[4]均表明：盾构隧道–竖井连接部位管环张开量显著大于远离连接部位张开量。Duran等^[5]研究发现：地震作用下竖井的扭转对临近盾构隧道管环的变形和应变存在显著影响。震害经验与动力分析表明,隧道结构的地震破坏是由过大的地基变形而非惯性力引起的^[1,6]。因此,对海底盾构隧道与竖井连接部位进行专门的地震反应分析与抗震设计是非常必要的。

盾构隧道与竖井连接部位地震反应特性的现有研究,未考虑管环间的非连续变形特征,难以分析管环截面应力和张开量的分布情况；采用的材料本构模型过于简单,难以反映材料非线性特性对盾构隧道与竖井连接部位地震反应的真实影响。此外,地震动为单向、双向输入,未见三向地震动激励下盾构隧道与竖井连接处地震反应分析。因此,对海底隧道与竖井连接部位进行三向非线性地震反应分析是必要的。

地下结构的纵向地震反应分析,通常采用实体有限元模型的动力时程分析法和简化的土弹簧–空心梁单元模型的拟静力广义反应位移法。广义反应位移法具有简单且计算量小的特点,近年在盾构隧道、地下管廊纵向地震反应分析中广泛应用^[7-9],且纳入《地下结构抗震设计标准》^[10],但未见采用广义反应位移法分析盾构隧道与竖井连接处地震反应特性的研究。

综上,本文以汕头海湾海底盾构隧道–竖井连接部位为案例,采用动力时程分析法和广义反应位移法进行平行分析,考虑海床土和混凝土的材料非线性、三向（纵轴向+横向+竖向）和双向（纵轴向+竖向）地震动输入,探讨了海底盾构隧道–竖井连接部位的地震反应特性与损伤机理。该研究结果可为盾构隧道与竖井连接部位的抗震设计提供有益的参考依据。

1.   数值分析模型

1.1   三维精细化有限元模型

汕头海湾超大直径海底盾构隧道的管环外径14.5 m,内径13.3 m,每节管环长2 m；竖井长50 m,宽25 m,壁厚1.2 m,埋深30 m。竖井区域海床土沿隧道纵轴线的地质剖面如图1所示。三维动力时程分析模型示意图见图2：纵向长225 m,横向宽180 m；取剪切波速≥500 m/s且其下不存在更低波速的界面作为地震基岩面,深度取至100 m。盾构隧道每节管环长2 m,单孔取50节。土体和结构均采用拉格朗日实体单元模拟,单元类型为C3D8R。对场地网格尺寸做精细化处理,竖井、区间隧道四周土体的网格尺寸1 m；其余部分的土体,竖向网格尺寸随深度增大,取1～3 m,水平向网格尺寸随离竖井距离增大而增大,取1～5 m；竖井、区间隧道网格尺寸1 m。土单元数为496624,结构单元数为35036,总的自由度数为1789431。结构和土体之间定义接触面,取摩擦系数0.3。采用梁单元模拟钢筋,通过Embedded方法嵌入混凝土内部。模型侧边和底部设置黏弹性人工边界^[11]。

图  1  沿隧道纵向轴线断面地质示意图

Figure  1.  Geological sketch of section along longitudinal axis of tunnel

下载: 全尺寸图片 幻灯片

图  2  三维动力时程分析的精细化有限元模型

Figure  2.  Refined finite element model for three-dimensional dynamic time-history analysis

下载: 全尺寸图片 幻灯片

1.2   基于壳–梁单元的广义反应位移法模型

通过三维自由场非线性地震反应计算（图3（a））,提取结构–土层界面处指定节点的位移时程,并施加于模型结构,称为广义反应位移法。基于Timoshenko梁理论,岸侧明挖的箱型隧道简化为箱型梁、海底盾构隧道简化为环型梁,采用空心的B31梁单元模拟；采用S4R壳单元模拟竖井；竖井两侧与隧道接触的面上节点自由度耦合至面中心,中心点与明挖隧道箱型截面梁单元进行自由度绑定,与盾构隧道环型截面梁单元采用连接单元连接,如图3（b）所示。自由场地震反应计算的土单元数为241200；壳、梁单元数为1692,整体模型自由度数为763896。在梁单元的上、下、左和右四个方向及壳单元的外侧面设置土弹簧,共1348个弹簧支点,每个支点各设置沿隧道的纵向、横向（切向）及法向3个土弹簧,共4044个土弹簧。将从自由场目标点中提取的4044条位移时程,如图3中的点a～j所示,输入到土弹簧支点,对隧道–竖向系统进行拟静力计算,即可获得隧道–竖井连接部位的地震反应。

图  3  基于Timoshenko壳–梁单元的广义反应位移法模型

Figure  3.  Model for generalized response displacement method based on Timoshenko shell-beam elements

下载: 全尺寸图片 幻灯片

1.3   材料模型与参数

（1）土体

采用基于Davidenkov骨架曲线、不规则加卸载准则构造的Non-Masing本构模型^[9,12-13]描述地震作用土体应力–应变关系,以自定义子程序模块嵌入Abaqus软件。Davidenkov骨架曲线表达式为

其中,

式中 G_max=ρv_s²为最大剪切模量,ρ为土体密度,v_s为土体剪切波速；A,B和${\gamma }_0$为土的试验参数。钻孔波速剖面分布见图4；试验给出的各类土动剪切模量比G/G_max、阻尼比λ与剪应变$\gamma$的关系曲线^[9],如图5所示；各类土的本构模型参数列于表1。

图  4  典型钻孔S波和P波波速剖面图

Figure  4.  Typical velocity profiles of S- and P-waves

下载: 全尺寸图片 幻灯片

图  5  海底隧道场址典型土的G/G_max–γ和λ–γ试验曲线

Figure  5.  G/G_max–γ and λ–γ curves of typical soils for subsea tunnel site

下载: 全尺寸图片 幻灯片

表  1  海床土的本构模型参数

Table  1.  Parameters of constitutive model for subsea soils

土层	A	B	${\gamma }_{\text{0}}$/(10-4)
淤泥	1.03	0.40	5.1
淤泥质土	1.01	0.39	6.7
中粗砂	1.20	0.37	7.4
粉质黏土	1.12	0.41	9.0
花岗岩	1.30	0.40	10.0

下载: 导出CSV 

| 显示表格

（2）混凝土

采用混凝土损伤模型（CDP）描述循环荷载下混凝土的不可逆累积塑性变形和刚度损失^[14]。CDP模型需定义塑性参数与损伤参数,塑性参数描述材料屈服面变化形式,损伤参数利用拉、压应力–应变关系和损伤因子描述材料加卸载特性。隧道、竖井的混凝土强度等级C60,相应的损伤模型参数列于表2；根据比例应变法,求得拉、压损伤因子^[15]。

表  2  混凝土损伤模型参数

Table  2.  Parameters of concrete from plasticity tests

弹性模量Ec /GPa	泊松比μ		初始屈服应力σco /MPa	抗压强度σcu /MPa	抗拉强度σtu /MPa	膨胀角Ψ/(°)	偏心率δ	fbo/fco	Kc	黏性系数u
36	0.2	33.27		47.57	3.88	38.0	0.1	1.16	2/3	0.0005
注：fbo / fco为双轴与单轴极限抗压强度比；Kc为拉伸子午面和压缩子午面上的第二应力不变量之比。

下载: 导出CSV 

| 显示表格

（3）连接螺栓

沿海底盾构隧道纵向会发生压缩、拉伸和弯曲变形,如图6所示。受拉时,环间拉力由螺栓承担,接头处轴向抗拉刚度为所有螺栓刚度之和；受压时,混凝土管环承担压力,接头轴向刚度为管环受压刚度；受弯时,拉力由螺栓承担,压力由混凝土管环承担。相邻管环间设置42根M42型螺栓,管环间防水允许张开量限值为15 mm,螺栓的物理力学参数见文献[9]。

图  6  管环间荷载传递机理示意图

Figure  6.  Schematic diagram of load transfer mechanism of ring segment joints

下载: 全尺寸图片 幻灯片

根据“等效连续化模型”的概念,假设：管环为不考虑管片拼接影响的弹性体；管环边侧各点为从控制点,边侧中心点为主控制点,采用多点约束（multi- point constraints）将所有从控制点自由度耦合至主控制点,使用连接单元（connector element）作为管环间的连接弹簧,如图7所示。通过定义连接单元的轴向刚度K_u和转动刚度$K_{\theta }$来模拟管环间的荷载传递。连接单元轴向、转动刚度取值见文献[9]。

图  7  管环多点约束及连接单元示意图

Figure  7.  Schematic diagram of multi-point constraints and connection elements of ring segments

下载: 全尺寸图片 幻灯片

（4）土弹簧

通过土弹簧方式将土体的变形作用于隧道、竖井,以模拟土–结构相互作用。土弹簧系数为^[10]

式中 k_l,k_t和k_n分别为沿隧道（或竖井）的纵向（竖向）、横向（切向）和法向的土弹簧系数；G为应变水平相容的土体动剪切模量,通过自由场地震反应计算确定。

1.4   输入地震动

根据汕头海湾隧道场址地震安全评价结果,50 a超越概率63%（小震）、10%（中震）和2%（大震）的基岩峰值加速度（peak bedrock acceleration, PBA）,近似可取0.1g,0.2g,0.4g,其中,g为重力加速度。由于缺乏场址内的强震记录,根据历史地震资料,选取国外地震台网中3次强震的基岩加速度记录作为基岩输入地震动,详细信息列于表3。图8为地震记录的加速度时程及其对应傅氏谱图。Darfield记录主频段集中在低频,Iwate记录的高频成分发育,而Kumano记录的频谱较宽。这3次地震的加速度记录的频谱特性存在显著差异。三向输入方向设定为：隧道纵向为NS向记录、横向为EW向记录,竖向为UD向记录。双向输入方向设定为：隧道纵向为NS向记录,竖向为UD向记录。为考虑输入地震动影响,将原始记录的NS向峰值加速度PGA调幅为0.1g（小震）,0.2g（中震）和0.4g（大震）,EW向、UD向的记录等比调幅。

表  3  基岩输入地震动的原始地震记录信息

Table  3.  Information of original earthquake recordings for input bedrock motions

地震名称	台站	震级Ms	震中距/km	分量	PGA/g	持时D5-95 /s	卓越频率fp /Hz
Darfield	Page Road Pumping(PRP) Station	7.1	53.53	NS	0.230	22.29	0.34
				EW	0.198	21.98	0.38
				UD	0.318	15.08	11.22
Iwate	IWTH27	7.2	56.85	NS	0.347	20.44	8.54
				EW	0.237	20.32	4.38
				UD	0.128	22.4	14.21
Kumano	KUMANO	6.5	68	NS	0.014	19.08	3.52
				EW	0.013	20.14	3.93
				UD	0.012	20.26	0.46

下载: 导出CSV 

| 显示表格

图  8  基岩输入地震动加速度时程及其傅氏谱（规准到PGA = 0.1g）

Figure  8.  Time histories of acceleration and Fourier spectra of input bedrock motions (0.1g level produced by scaling)

下载: 全尺寸图片 幻灯片

2.   环间张开量

环间张开量是评价盾构隧道抗震性能的重要指标,过大的张开量会导致发生漏水事故。基于对称性,取B通道为研究对象,对竖井和B通道连接接头8个测点位置进行分析,如图9所示。本文将盾构隧道第一环与竖井墙面相连环缝处简称为隧道–竖井接头。定义峰值张开量为左管环右端点位移时程与右管环左端对应点位移时程之差的绝对值最大值；环间张开量极值为相邻两环各测点的峰值张开量的最大值。

图  9  隧道–竖井接头处的监测位置

Figure  9.  Monitoring positions at junction of tunnel and shaft

下载: 全尺寸图片 幻灯片

2.1   接头横截面的峰值张开量

图10给出了三维动力时程分析法计算的隧道–竖井接头各测点峰值张开量。不同频谱的小震作用时,峰值张开量在接头的环形截面上具有对称性；随着PBA增大,各测点的峰值张开量增大,呈现出环形非对称的特征,外拱肩部位的最大。这说明采用连接单元模拟螺栓,可以较合理地反映峰值张开量的环形变形特征。究其原因：地震时隧道–竖井接头受其四周土体变形的影响,且海洋土的强非线性特性随地震动强度的增加也愈加强烈。大震、Darfield记录激励时,外拱肩部位的峰值张开量为17.32 mm,大于防水限值,存在安全隐患；其余工况均未超出防水限值。总体上,峰值张开量较大的位置多为拱顶和外拱肩部位。这是由于横波垂直或近似垂直于隧道轴线时,圆形隧道会发生椭圆变形^[16]；另外,隧道的上覆土质较软,隧道与土层的刚度差异大,地震作用下拱顶管片与土层之间的相对滑移较大。建议在设计海底隧道时,对拱顶、外拱肩采取特殊的抗震措施,以减小张开量。

图  10  隧道–竖井接头各测点峰值张开量

Figure  10.  Peak opening widths of monitoring elements at junction of tunnel and shaft

下载: 全尺寸图片 幻灯片

低频成分丰富的Darfield记录激励时接头处的峰值张开量最大,高频成分丰富的Iwate记录激励时最小,频谱较宽的Kumano记录激励时居中。地震波传播通过松软海床土层时,会发生高频滤波、低频放大现象,而海底隧道对周围土层的位移具有依随性,隧道自身的动力特性对其影响较小。因此低频发育的地震记录作用时海底隧道管环的张开量尤为显著。

由图10可以看出：输入Darfield记录、Iwate记录和Kumano记录时,三向作用下计算的最大峰值张开量,中震时较双向震动计算的最大峰值张开量依次增大44.2%,76.7%和54.9 %；大震时较双向震动计算的最大峰值张开量依次增大52.9%,79.4%和59.1%。此外,双向震动时最大峰值张开量均发生在拱顶,三向震动时最大峰值张开量发生在外拱肩。因此,横向地震动对海底隧道安全具有严重不利影响,不容忽视。

2.2   环间张开量极值沿隧道纵轴线的分布

图11为广义反应位移法与动力时程分析法计算的环间张开量极值沿隧道纵轴线的分布。可以发现：基岩地震动特性对环间张开量极值的影响显著,但对沿隧道纵轴线分布的特征影响不大,不同地震记录激励时,其变化趋势相近：近接头处的环间张开量极值明显大于区间段,且接头处的环间张开量极值最大,在大震Darfield记录激励下,两种方法计算得的隧道–竖井接头处环间张开量极值>15 mm,在大震Kumano记录激励下,仅反应位移法计算值超限。

图  11  不同方法计算的环间张开量极值

Figure  11.  Maximum values of opening widths at ring intersegment calculated by different methods

下载: 全尺寸图片 幻灯片

不论PBA大小与基岩地震动的频谱特性差异,动力时程分析法和广义反应位移法计算的环间张开量极值沿隧道纵轴线分布的趋势有较好的一致性；另外,除极个别位置外,广义反应位移法计算的环间张开量极值大于动力时程法的计算值。两者的差异,可能存在以下原因：①广义反应位移法仅考虑了自由场土体的非线性效应,未考虑结构与四周土体动力相互作用引起的非线性效应；②广义反应位移法采用Timoshenko梁理论建模,该模型比三维动力时程法采用的实体模型具有更大刚度；③广义反应位移法是一种拟静力法,忽略了隧道周边向四周传播的波对周围土体的辐射阻尼,而三维动力时程法中由于土体阻尼的存在,其环间张开量极值曲线较壳–梁模型的平缓许多,削平了沿隧道纵轴线分布的部分峰值点。

另外,三向地震动作用时的环间张开量极值通常都大于双向地震动作用时的环间张开量极值,尤其在距离竖井约20 m以内,说明横向地震动对环间张开量存在显著的不利影响。究其原因：横向地震动激励使竖井产生较大扭转,加剧空间效应,弯扭作用下管环会产生更大的张开量。

3.   环间应力与地震损伤

3.1   接头环截面的峰值应力

图12为三向震动下动力时程分析法计算得到的隧道–竖井接头处各测点峰值应力（Mises应力时程的最大值）。在地震作用下管环间反复出现张开与闭合状态,峰值应力在接头的环形截面上呈X型分布,与竖轴共轭45°方向上的峰值应力较大。峰值应力最大值基本处于环形截面的外拱脚处,与峰值张开量最大值的位置呈现出对称性。因此,在三向地震动激励下,隧道–竖井接头部位会发生环形的扭转,当拱肩部位产生较大的张开量时,管环的变形会传递至另一端,致使拱脚处产生较大的挤压。由于管环混凝土的抗压刚度远大于连接螺栓的拉伸刚度,反而使拱脚部位的峰值应力大于拱肩部位的峰值应力。

图  12  三向震动下隧道–竖井接头各测点峰值Mises应力

Figure  12.  Maximum Mises stresses of monitoring elements at junction of tunnel and shaft under tri-directional ground motions

下载: 全尺寸图片 幻灯片

3.2   环间应力极值沿隧道纵轴线的变化

图13为三向激励下广义反应位移法与动力时程分析法计算的环间应力极值（每节管环的峰值Mises应力最大值）沿纵轴线的分布。类似于接头部位的峰时接头处的环间应力极值最大、高频成分丰富的Iwate记录激励时最小、频谱较宽的Kumano记录激励时居中,且环间应力极值随PBA的增大而增大。离竖井墙面约20 m范围内,环间应力极值随距离增大而显著减小,且广义反应位移法计算值大于动力时程分析法计算值。基岩地震动特性对离竖井约20～60 m范围的环间应力极值随距离的变化趋势有显著影响。总体而言,离竖井约60 m以外,两种方法计算的应力极值,收敛于一个稳定的较小值,两种方法得出的应力极值相近,沿隧道纵轴线的变化趋于一致。因此,对离竖井20 m范围内的管环,应施加必要的抗震增强措施,以降低应力集中引起的隧道管片破坏的可能性。

图  13  三向震动下不同方法计算的管环间Mises应力极值

Figure  13.  Peak Mises stresses at ring intersegment calculated by different methods under tri-directional ground motions

下载: 全尺寸图片 幻灯片

3.3   连接部位的地震损伤

图14给出了三向地震动激励下隧道–竖井接头处的地震损伤分布云图。刚度损伤因子SDEG取值介于0（蓝色：无损伤）和1（红色：完全损伤）之间,是应力状态与拉、压损伤因子的函数^[14],可综合评价混凝土拉压损伤状态。接头部位损伤程度随PBA增大而增大。由于竖井与盾构隧道的刚度差异显著,竖井墙体与盾构隧道间存在较大的纵向拉压变形。因此,在不同的PBA水平、频谱特性的地震动作用下,竖井墙的应力集中较之隧道管片更为严重,损伤更显著。PBA = 0.1g时,仅在竖井墙产生部分损伤；PBA = 0.2g时,低频丰富的Darfield记录及宽频带的Kumano记录激励下,竖井及距离竖井的第2～4环隧道管片的损伤程度及范围显著增大；PBA = 0.4g时,在Darfield记录和Kumano记录激励下,竖井及距离竖井的第4～6环隧道管片的混凝土损伤严重,竖井墙已大面积整体失效。与隧道–竖井接头相比,远离竖井的隧道管片损伤较轻。观察接头部位的损伤分布,可以发现：竖井墙与竖轴共轭45°方向上先出现损伤,损伤分布呈现X形,与图12中峰值Mises应力分布特征一致。

图  14  三向震动下隧道–竖井连接部位损伤分布云图

Figure  14.  Cloud diagram of damage indexes at junction of tunnel and shaft under tri-directional ground motions

下载: 全尺寸图片 幻灯片

4.   结论与讨论

通过海底盾构隧道–竖井连接部位非线性地震反应特性的三维动力时程分析和广义反应位移法的计算与分析,得到以下4点结论。

（1）多点约束和连接单元能有效地模拟盾构隧道管环间的变形特性；隧道–竖井接头拱顶与外拱肩处的峰值张开量最大；竖井墙与竖向共轭45°方向处的地震应力集中显著,且拱腰处的应力比拱肩处的稍大,竖井墙在该部位的地震损伤严重。

（2）由于海床软土场地高频滤波、低频放大的特性,低频丰富的地震动作用时,海底盾构隧道–竖井连接部位的地震反应更强烈。

（3）盾构隧道–竖井接头部位存在刚度突变,接头处的地震变形与应力反应远大于距离接头较远处隧道管环的地震变形与应力反应；横向地震动对接头处的应力与变形以及环间张开量的不利影响不容忽视。

（4）两种方法计算的接头处及管环间的地震变形与应力的基本特征一致,但广义反应位移法的计算值大于三维动力时程分析法的计算值。

对于分析模型,基于16CPUS的并行算法进行三维非线性动力时程分析,单次计算耗时约20 h；采用广义反应位移法进行三维自由场非线性地震反应分析+简化的壳–梁单元–土弹簧模型的拟静力计算,两次计算耗时约4 h。综上,采用广义反应位移法进行盾构隧道–竖井连接部位的非线性地震反应拟静力分析,可有效提高计算效率,但与三维非线性动力时程分析的计算结果相比,存在偏于安全的较大的误差。