作为海洋大国，我国经济社会发展与海域资源开发密不可分，推动海洋经济高质量发展已成为建设海洋强国的战略核心。在众多工程结构中，隧道凭借其显著的社会效益与经济价值被广泛应用，其中沉管隧道作为水下重要工程设施，需应对复杂地质条件与高频地震活动的双重考验。

与陆域隧道抗震研究相比，目前针对海底隧道的抗震研究相对较少，受限于原型监测的困难，模型试验和数值计算已经成为海底隧道抗震研究的主要方法。禹海涛等^[1-4]基于自主研发的多台阵振动台联合试验平台，系统揭示了超长沉管隧道在非一致地震动激励下的动力响应规律和损伤机制。陆耀波等^[5]通过非均匀场地条件下沉管隧道振动台模型试验，研究了土体动力特性参数空间变异性对结构地震响应变化的影响。程新俊等^[6]对饱和砂土中的多段式沉管隧道进行了振动台试验研究，并总结了沉管隧道抗震的不利位置。KIM等^[7]通过一系列振动台试验，研究了不同回填材料性质和覆盖土层厚度的浅埋隧道的地震反应。

海底隧道地震响应的数值分析方法主要有梁-弹簧模型、质点-弹簧模型等基于位移的简化方法以及整体动力时程法。赵密等^[8]建立了基于规范位移模式的隧道纵向反应加速度方法，对不同场地条件下的隧道纵向地震响应规律进行了研究。CHEN等^[9]以汕头海底盾构隧道为工程背景，综合考虑了竖向、横向和纵向地震动的耦合作用，完成了长大海底隧道纵向抗震性能评估，并提出了盾构隧道的减隔震措施。刘晶波等^[10-13]基于波动理论与土-结构耦合效应，构建了适用于地下隧道纵向地震动力响应评估的整体式反应位移法理论，实现了纵向地震动激励下的隧道轴向变形响应分析。整体动力时程法能够综合考虑复杂的土层分布，土体和结构的非线性性能，并模拟纵向及横向地震动作用下隧道响应与土-结构相互作用。李心熙等^[14]对沉管隧道进行了精细化建模，重点探讨了地震动作用下隧道内力和变形响应对刚度变化的敏感性，揭示了刚度变化对隧道结构动力响应规律的影响。LI等^[15]将一维时域有限元法计算得到的自由场波动推广到三维，对隧道结构进行了非同步地震动作用下的精细非线性数值分析。CHENG等^[16-17]实现了应力场与渗流场的动态耦合分析，揭示了复杂环境荷载作用下海底隧道的动态响应机理。ZHAO等^[18-19]在Biot动力固结方程中引入新的残余孔隙水压力源项，揭示了沉管隧道周围海床的动态响应及液化机理。赵凯等^[20]对海底盾构隧道-竖井连接区进行了三维精细化建模，系统研究了海底盾构隧道的地震响应及端部效应，并提出了减隔震措施。

失稳上浮是沉管隧道的主要地震失效模式之一，已有研究较少涉及，为深入开展海水和代表性海床中隧道地震失稳机理的研究，本文基于渤海典型场地条件，建立了考虑海床-隧道动力相互作用的有限元数值分析模型。该模型能够考虑土骨架与孔隙水两相介质相互作用的动力有效应力和模拟宽应变范围海床土动力特性探讨了地震动作用下沉管隧道动力响应规律及其土体液化失稳上浮机理。

1   可液化海床有效应力分析方法

1.1   控制方程

在Biot动力固结方程框架内^[21]，考虑饱和土体中土颗粒与孔隙水的相互作用，建立动力有效应力分析方法：

式中：$\nabla$为拉普拉斯算子；上标“∙”为对时间求导；p为超孔压；n_s为孔隙率；k为渗透系数；u为土单元位移矢量；ε为土单元应变张量，角标e和p分别代表弹性应变张量和塑性应变张量；G为动剪切模量；ν为泊松比；γ_w为孔隙流体重度；K_f为流体压缩系数。

赵凯等^[22]引入等效循环应力比ESR^[23]量化复杂应力路径，综合考虑了动应力路径对塑性剪切体应变发展规律的影响，重新定义了复杂动应力路径下塑性剪切-体应变增量的规准化模型：

式中：$\overline q$为土体等效循环应力；$\sigma_{\rm{p}}'$为初始平均有效应力；σ_z竖向应力；σ_x为水平向应力；T为动荷载周期；|q(t)|为应力路径上任意一点与原点的距离。

将复杂动应力路径引起的塑性-剪切体积应变增量进行归准化，并作为超孔隙水压力累积源项引入Biot动力固结方程框架中，其表达式为：

式中：${k_1}$、${k_2}$均为与土性相关的拟合参数；${\rm ESR}_{{\text{t}}}$为体积门槛循环应力比。

1.2   本构模型

将Maing法则修正的Davidenkov骨架曲线和能够描述土体剪切破坏和张拉破坏的复合Mohr-Coulomb屈服准则结合，分别描述弹性范围内的应力-应变曲线和海床土体接近液化时及液化后的塑性大变形^[24-27]，土骨架的循环应力-应变关系如图1所示。

图  1  土骨架动力应力-应变曲线

Figure  1.  Schematic illustration of dynamic stress-strain curve of soil skeleton

下载: 全尺寸图片 幻灯片

土体在循环荷载作用下的软化现象通过模量退化模型来描述，并以残余孔隙水压力作为退化参数^[28]：

式中：$G_{\max }^t$为t时刻最大动剪切模量；$\overline p$为残余孔隙水压力；r_u为残余孔压比。

1.3   有效应力方法的实现

基于FLAC3D软件提供的二次开发平台，实现了该有效应力分析方法子程序的开发，子程序实现流程如图2所示。计算过程中，子程序首先根据传入的初始有效应力场对土体的初始剪切模量进行修正；随后子程序获取FLAC3D主程序中计算出的应变增量dε_ij、应变张量ε_ij、应力张量σ_ij和总孔隙水压力p，并判别土单元是否发生应力转向。若没有发生应力转向：计算当前的等效剪应变并更新当前时变剪切模量；若发生了应力转向：更新应力转向点$({\gamma _{{\text{rev}}}},{\tau _{{\text{rev}}}})$，根据式(4)和式(5)计算当前加卸载周期内的等效循环应力比ESR，并根据式(6)和式(7)计算剪切-体积耦合应变增量$\varepsilon _{jj}^{\rm p}$，并代入式(2)中计算残余孔隙水压力$\overline p$，随后根据式(8)计算退化后的最大动剪切模量。在获得当前时变剪切模量G^t后计算土单元试探应力dσ_ij，若应力F($\sigma _{ij}^{{t}} + {\rm{d}}{\sigma _{ij}}$)超过屈服面则根据Mohr-Coulomb准则对骨架曲线进行修正，之后根据当前应力状态通过Biot动力固结方程式(1)和式(2)计算振荡孔隙水压力$\tilde p$。最后子程序将计算得到的土单元应力张量和孔隙水压力返回主程序，并进行下一步计算。为防止数值计算失稳，设置了平均有效应力最小阈值P_t，即土体液化后仍具有一定的残余剪切强度。

图  2  动力有效应力方法计算流程

Figure  2.  Flow chart of dynamic effective stress method

下载: 全尺寸图片 幻灯片

1.4   方法验证

对不同应力路径下饱和南京细砂不排水循环剪切试验进行数值模拟，以验证本文采用的有效应力方法的准确性和合理性。试验中试样相对密度D_r=50%，饱和重度γ=19.2 kN/m²，初始有效围压σ=100 kPa，加载循环应力比CSR=0.15，单元测试模型计算参数见表1。

图3对比了本文方法和室内试验获得的超孔隙水压力时程曲线。由图3可见，本文方法较好地预测了超孔隙水压力的发展规律，与试验结果仍存在一定的偏差可能是因为不同工况在试验过程(如制样、饱和、加载过程等)中存在的不确定性造成的^[29-30]。

表  1  单元测试参数

Table  1.  Mechanical parameters of the soil element

初始剪切模量Gmax	泊松比ν	A	B	γ0/(×10−4)	k1	k2	ESRt	黏聚力	内摩擦角/(o)	抗拉强度
50 MPa	0.3	1.02	0.35	4	0.76	0.58	0.02	0	30	0
注：A、B和γ0为骨架曲线参数。

下载: 导出CSV 

| 显示表格

图  3  不同应力路径下超孔压时程预测值与实测值的对比

Figure  3.  Comparison between the predicted and measured excess pore pressure time-histories under different stress paths

下载: 全尺寸图片 幻灯片

2   计算模型建立

以渤海海域中的典型钻孔作为场地条件^[31-32]，建立沉管隧道-可液化海床相互作用模型，该场地钻孔示意及相对应的剪切波速如图4所示。

图  4  渤海典型钻孔示意

Figure  4.  A typical borehole at the Bohai Strait

下载: 全尺寸图片 幻灯片

计算模型以剪切波速≥500 m/s作为地震动输入基岩面，本文中的模型场地总厚度为100 m。模型中沉管隧道为总宽度L = 20 m的三跨结构，2个边跨宽8 m，中间跨宽4 m，沉管隧道高8 m，顶板、底板、立柱及侧墙的厚度均为1 m，结构顶板覆土深度为1 m。模型中海水的深度为10 m，通过在海床表面施加等效静水压力来模拟海水的作用。假定整个地震动作用过程中，海床均处于不排水状态。设定海床与沉管隧道为法向硬接触，切向摩擦接触，且接触面不透水。采用四节点等应变有限差分单元建立均匀网格布局，计算域网格在波传播方向的尺寸限定为土层地震波传播的截止频率所对应的最短波长的1/8~1/10^[9]，在本文模型中，海床土体网格尺寸为1 m×1 m，沉管隧道的网格尺寸为0.25 m×0.25 m，如图5所示。沉管隧道的计算参数见表2，海床土体的计算参数见表3。

图  5  模型网格示意图

Figure  5.  Grid layout of the tunnel-seabed system in detail

下载: 全尺寸图片 幻灯片

表  2  沉管隧道计算参数

Table  2.  Parameters of immersed tunnel

弹性模量E/(kN·m−2)	泊松比ν	重度ρ/(kg·m−3)	埋深h/m
3.5×107	0.18	2500	1

下载: 导出CSV 

| 显示表格

将模型通过静力固结分析得到的有效应力和静孔隙水压力设定为动力计算的初始条件，静力固结的具体实施步骤为：1)分级加载与边界约束：对模型施加标准重力荷载和孔隙水压力，并通过固定约束边界条件仅允许沿重力方向的变形；2) 应力场迭代修正：将每次求解后的有效应力和孔隙水压力重新加载至初始模型，再次进行计算；3)应力集中消除：重复步骤2直到结构和土体中的应力集中现象基本消失。进行上述步骤可以显著降低初始不平衡力对动力分析的影响。

在进行动力计算时，为有效限制计算域尺寸并抑制边界效应，场地两侧边界采用竖向约束与水平施加弹簧阻尼器的黏弹性边界组合，模型横向宽度取400 m，为沉管隧道宽度的20倍；对于岩土材料，临界阻尼比取值一般取在2%~5%，本文中临界阻尼比取值为5%。

表  3  海床计算参数

Table  3.  Parameters of seabed soil

类别	剪切模量/(kN·m-2)	泊松比	Davidenkov参数				孔压参数				Mohr-Coulomb模型参数
			A	B	γ0/(×10−4)		k1	k2	ESRt		黏聚力/kPa	内摩擦角/(°)	抗拉强度
粉土	7×104	0.30	1.03	0.50	1.8		0.43	0.93	0.02		0	30	0
粉质黏土	3×104	0.25	0.96	0.47	7.3		−	−	−		10	25	0
细砂	8×104	0.20	1.02	0.48	17.9		0.19	2.1	0.02		0	35	0

下载: 导出CSV 

| 显示表格

采用零厚度接触面单元模拟沉管隧道与周围海床土体之间的相对滑动及脱开等相互作用^[33-34]。海床土体与沉管隧道之间的相对滑动发生在接触面处剪应力超过摩擦阻力后，相对滑动条件可描述为：

式中：${\sigma '_{\text{n}}}$为接触面单元处的有效法向应力；$\mu$为沉管隧道与周围海床之间的摩擦系数，本文中取值为0.7。

由于缺乏渤海海峡的实际地震动记录，选取1995年日本神户地震中Port Island场地钻孔记录的Kobe波作为基岩输入地震动，将原始记录中的峰值加速度分别调幅为0.2 g和0.4 g，其调幅后的水平向加速度时程如图6所示。Kobe波的水平向加速度时程是不对称的，在负方向上的峰值加速度为正方向上的1.6倍。地震动由模型底部边界输入，为基岩面垂直向上传播的剪切波。

图  6  输入地震动加速度时程

Figure  6.  Input acceleration time history

下载: 全尺寸图片 幻灯片

3   结果分析

3.1   沉管隧道上浮

图7给出了PGA分别为0.2 g和0.4 g时，沉管隧道顶板两侧的上浮时程。沉管隧道两侧的上浮曲线发生了明显的分离，隧道右侧最终上浮量明显大于左侧，当PGA为0.2 g时沉管结构左右侧墙的上浮位移分别为0.08 m和 0.11 m，当PGA达到0.4 g时两侧的竖向位移分别为0.15 m和0.31 m，可认为沉管隧道发生了不均匀上浮。

图  7  沉管隧道不同位置竖向位移时程

Figure  7.  Time history of floating displacement of immersed tunnel

下载: 全尺寸图片 幻灯片

图8给出了沉管隧道周围海床表层震后残余竖向位移分布，以沉管隧道对称轴为中线，定义隆起为正，沉降为负，图8显示沉管隧道结构两侧海床表层呈现出显著差异沉降现象，当PGA=0.4 g时，沉管两侧海床土体最大沉降差达到了0.07 m。在海床-沉管隧道动力相互作用过程中，沉管侧边土体与沉管的相对滑动，产生向下的位移迫使沉管底部土体发生塑性变形隆起，成为沉管隧道失稳上浮的关键因素。当PGA=0.2 g时，沉管隧道左侧墙上方土体的位移为0.09 m，右侧墙上方土体的位移为0.13 m，位移差为0.04 m；当PGA增大至0.4 g时，左右两侧土体的竖向位移分别为0.14 m和0.3 m，位移差值达到了0.16 m。由此说明，随着输入地震动强度的增加，隧道左右侧的不均匀沉降现象愈加明显。

图  8  沉管隧道周围海床表面震后竖向残余位移

Figure  8.  Post-earthquake vertical residual displacement of seabed soil around immersed tunnel

下载: 全尺寸图片 幻灯片

3.2   沉管隧道-可液化海床相互作用

沉管隧道侧墙受到的摩擦阻力为沉管隧道抗浮力的重要组成部分，当沉管隧道所受的上浮力大于抗浮力时，沉管隧道将发生上浮失稳破坏。沉管两侧的侧摩阻力可由接触面上的法向有效应力计算得到：

式中：μ为海床与沉管隧道之间的摩擦系数；σ_iz为与沉管隧道侧墙接触土单元的法向有效应力；A_i为沉管隧道侧墙单元外表面面积；θ_i为沉管隧道侧墙外表面与竖直向上方向的夹角。

沉管隧道侧墙不同位置所受摩擦阻力存在一定的差异，选取侧墙中部及脚部作为分析点。图9给出了PGA为0.2 g和0.4 g情况下沉管隧道侧墙不同测点所受摩擦阻力的时程，地震动的作用下沉管隧道受到的摩擦阻力逐渐减小并最终稳定。沉管隧道两侧受到的摩擦阻力出现了明显的分离，沉管右侧所受摩擦阻力衰减快于左侧。当PGA=0.4 g时，与测点S1、S2和S4相邻海床均发生了液化，所受摩擦阻力基本消失。

由超孔隙水压力上升引起的浮力作用在沉管底板上，可由下式进行计算：

式中：P_i为与沉管隧道底板相接触的土单元的超孔隙水压力；B_i为沉管隧道底板单元外表面面积；η_i为沉管隧道底板外表面与竖直向下方向的夹角。

图  9  沉管隧道侧壁不同位置摩擦阻力时程

Figure  9.  Time histories of friction resistance at different positions on side wall of immersed tunnel

下载: 全尺寸图片 幻灯片

为分析沉管底板不同位置受到浮力的差异，图10给出了PGA为0.2 g时两侧墙底部测点位置处的浮力时程，沉管隧道底部由超孔隙水压力引起的浮力在地震动作用下逐渐增大，且底板右侧受到的超浮力的增长速度明显快于左侧。

图  10  沉管隧道结构底部浮力时程

Figure  10.  Time history of buoyancy at the bottom of immersed tunnel

下载: 全尺寸图片 幻灯片

图11给出了地震动作用30 s后沉管隧道周围海床的平均有效应力分布，平均有效应力的分布存在着明显的不对称性。可见沉管隧道中轴线右侧海床平均有效应力明显小于左侧，结合图9、图10可知，沉管隧道两侧的受力存在着明显的非对称性。

图  11  沉管隧道周围海床最终平均有效应力分布

Figure  11.  Final mean effective stress distribution of seabed around immersed tunnel

下载: 全尺寸图片 幻灯片

图12给出了沉管隧道两侧土单元P1和P2的剪应力-剪应变曲线，隧道两侧海床土体的剪应变水平产生了明显的差异，当PGA = 0.2 g和0.4 g时，隧道左侧海床土体的剪应变水平分别为0.15%和0.4%，远小于隧道右侧海床土体的剪应变水平的0.2%和1%。

图  12  沉管隧道结构侧面海床应力-应变滞回曲线

Figure  12.  Stress-strain hysteresis curve of seabed on the side of immersed tunnel

下载: 全尺寸图片 幻灯片

综合上述结果可知，地震作用下沉管隧道两侧的受力产生了明显的非对称性，沉管右侧所受的抗浮力衰减较左侧更快，且沉管隧道右侧的应变水平也较沉管隧道左侧更大。结构受力的非对称性与周围海床土体变形的非对称性进而引起了沉管隧道的不均匀上浮。

4   结论

本文基于“扩展Masing”法则和复合Mohr-Coulomb屈服准则构造了土体粘弹塑性应力-应变滞回曲线，以等效循环应力比ESR量化复杂应力路径描述剪切-体积应变增长规律，实现了饱和砂土液化过程的有效应力分析，研究了砂质海床沉管隧道的地震响应，得到以下结论：

(1) 本文方法能够合理地表征复杂应力路径下饱和砂土的孔隙水压力动态累积规律和液化过程。

(2) 输入地震动加速度时程的非对称性导致了沉管隧道两侧受力及海床土体变形的非对称性，进而引发了沉管隧道的不均匀上浮。

(3) 沉管隧道侧摩擦阻力的减小、底部浮力的增大和周围海床土体的大变形共同导致了沉管隧道的上浮，输入地震动峰值加速度增大了沉管隧道的上浮量，同时增大了其不均匀上浮程度。