引言

地铁、隧道等地下结构破坏会引起重大的安全事故及经济损失,因此人们对隧道等地下结构的安全性十分重视,针对地下结构动力分析的研究已有不少成果(韩大建等,1999Uenishi等,2000Choi等,2002Huo等,2003袁勇等,2016)。目前,研究沉管隧道纵向地震响应的手段有:①大型振动台试验和离心机试验(袁勇等,2016)。利用振动台等试验方法分析沉管隧道的地震响应具有较高的可信度,但一般用来模拟水平和竖向振动。②简化力学模型方法:梁-弹簧模型、质点-弹簧模型(或等效质点-弹簧模型),前者把隧道作为支撑在弹性地基上的连续梁,后者把隧道作为多质点系统,接头统一用弹簧模拟,用弹簧和阻尼来模拟土-结构之间的相互作用(韩大建等,1999)。这些简化方法基于众多假设条件,所以计算精确度有所欠缺。③采用ADINA、ABUQUS、ANSYS等软件进行仿真模拟。利用大型软件进行水域沉管隧道地震响应分析已得到广泛的应用,陈向红(2013)利用有限元软件分析竖向地震作用下,动水压力对浅埋海底隧道内力的影响;陈贵红(2002)分析了在竖向地震作用下,隧道的埋深、海床的地质条件等对沉管隧道的影响;李鹏(2013)利用ANSYS软件进行沉管隧道纵向地震响应分析,但并未考虑海水的作用;周晓洁等(2017)利用ABAQUS软件建立二维模型分析SV波斜入射沉管隧道地震响应。

本文利用ADINA软件,考虑地基的人工边界和海水-海床-隧道间的相互作用,建立三维沉管隧道模型,分析P波的入射角度对沉管隧道动力响应的影响。

1 流固耦合分析原理

本文所用流体为理想流体,其连续方程和Naiver-Stokes控制方程可表达为:

$ \frac{{\partial {u_i}}}{{\partial {x_i}}} = 0 $ (1)
$ \frac{{\partial {u_i}}}{{\partial t}} + \sum\limits_{j=1}^3 {{u_j}} \frac{{\partial {u_i}}}{{\partial {x_j}}} = {f_i} - \frac{{\partial p}}{{\rho \partial {x_i}}}\;\;\;\;\;\;i = 1, \;2, \;3 $ (2)

式中:u为水流速度(m/s);${f_i}$为质点体力(m/s2);p为压强(Pa);$\rho $为水的密度(kg/m3);ij为坐标方向。

根据流体的边界条件和流体与海床的耦合边界条件,利用泛函对流体的连续方程进行有限元离散,可以得到流体的振动控制方程(陈贵红,2002):

$ {{\bf{K}}_{\rm{L}}}{{\bf{ \pmb{\mathsf{ δ}} }}_{\rm{P}}} + {{\bf{M}}_{\rm{L}}}{{\bf{ \pmb{\mathsf{\ddot δ}} }}_{\rm{P}}} + {\rho _{\rm{L}}}{{\bf{S}}^{\rm{T}}}{{\bf{ \pmb{\mathsf{\ddot δ}} }}_{\rm{u}}} = 0 $ (3)

式中:KLML分别为流体的刚度矩阵(N/m)和质量矩阵(kg);ST为流固耦合影响矩阵的转置(m3);${{\bf{ \pmb{\mathsf{ δ}} }}_{\rm{P}}}$${{\bf{ \pmb{\mathsf{\ddot δ}} }}_{\rm{P}}}$为流体质点压力自由度位移(m)、加速度矩阵(m/s2);${{\bf{ \pmb{\mathsf{\ddot δ}} }}_{\rm{u}}}$为固体质点加速度矩阵(m/s2);${\rho _{\rm{L}}}$为流体密度(kg/m3)。

考虑海床与流体的耦合边界条件,同样对海床土与隧道进行有限元离散,可得其振动控制方程(陈贵红,2002):

$ {{\bf{K}}_{\rm{S}}}{{\bf{ \pmb{\mathsf{ δ}} }}_{\rm{u}}} + {{\bf{M}}_{\rm{S}}}{{\bf{ \pmb{\mathsf{\ddot δ}} }}_{\rm{u}}} - {\bf{S}}{{\bf{ \pmb{\mathsf{ δ}} }}_{\rm{P}}} = 0 $ (4)

式中:KSMS分别为固体的刚度矩阵和质量矩阵;S为流固耦合影响矩阵;${{\bf{ \pmb{\mathsf{ δ}} }}_{\rm{u}}}$${{\bf{ \pmb{\mathsf{\ddot δ}} }}_{\rm{u}}}$为固体质点位移、加速度矩阵。

因此,存在外部荷载时,流-固耦合的振动方程可以改写为:

$ \left[ {\begin{array}{*{20}{c}} {{{\bf{M}}_{\rm{S}}}}&0\\ {{\rho _{\rm{L}}}}&{{{\bf{M}}_{\rm{L}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{{\bf{ \pmb{\mathsf{\ddot δ}} }}}_{\rm{u}}}}\\ {{{{\bf{ \pmb{\mathsf{\ddot δ}} }}}_{\rm{P}}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{{\bf{K}}_{\rm{S}}}}&{ - {\bf{S}}}\\ 0&{{{\bf{K}}_{\rm{L}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\bf{ \pmb{\mathsf{ δ}} }}_{\rm{u}}}}\\ {{{\bf{ \pmb{\mathsf{ δ}} }}_{\rm{P}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\bf{R}}\\ 0 \end{array}} \right] $ (5)

式中R为固态受到除流体压力之外的荷载(N)。

目前,求解公式(5)的方法有迭代法(分离法)和直接耦合法(同步求解法),具体的求解步骤可参考文献(岳戈,2010)。

2 粘弹性人工边界

粘弹性人工边界由弹簧和阻尼器组成,由于忽略质量m,并把阻尼器的一端固定,原来的阻尼和刚度系数的精度有所降低,刘晶波等(2006)对粘弹性人工边界的参数进行了修正并给出了其取值范围。

切向边界的弹簧系数K和阻尼系数C分别为:

$ \left\{ {\begin{array}{*{20}{l}} {{C_{\rm{T}}} = \rho {c_{\rm{S}}}}\\ {{K_{\rm{T}}} = {\alpha _{\rm{T}}}\frac{G}{r}} \end{array}} \right. $ (6)

法向边界为:

$ \left\{ {\begin{array}{*{20}{l}} {{C_{\rm{N}}} = \rho {c_{\rm{P}}}}\\ {{K_{\rm{N}}} = {\alpha _{\rm{N}}}\frac{G}{r}} \end{array}} \right. $ (7)

式中:CNCT分别为法向阻尼与切向阻尼;KNKT分别为弹簧法向与切向刚度(N/m);r为人工边界到波源的距离(m);cPcS分别为P波和SV波波速(m/s);G为介质的剪切模量(Pa);$\rho $为介质质量密度(kg/m3);${\alpha _{\rm{N}}}$${\alpha _{\rm{T}}}$分别为法向与切向粘弹性人工边界参数,即为粘弹性人工边界的修改参数。

3 等效荷载输入

在近场波动分析中,由于使用人工边界来模拟无限域地基,当外部荷载从远场无限域进入近场有限域时,在人工边界处,需要考虑荷载的输入问题,即外源输入问题。目前处理外源输入问题主要有波场分解法和等效边界力法。人工边界处的总波场可分为内行场(下标r表示)和外行场(下标s表示),人工边界处位移和作用力可以分解为(赵密,2009):

${{\bf{u}}_{\rm{B}}} = {{\bf{u}}_{\rm{Br}}} + {{\bf{u}}_{\rm{Bs}}} $ (8)
$ {{\bf{F}}_{\rm{B}}} = {{\bf{F}}_{\rm{Br}}} + {{\bf{F}}_{\rm{Bs}}} $ (9)

式中:uBruBs分别表示内行场位移和外行场位移(m);FBrFBs分别是保证边界位移等于内行场所需要提供的抵抗人工边界条件的荷载和保证边界位移等于内行场所需要提供的抵抗有限域的荷载(N);其中FBs可以表示为:

$ {{\bf{F}}_{\rm{Bs}}} = - {\bf{K}}_{\rm{B}}^\infty {{\bf{u}}_{\rm{Bs}}} - {\bf{C}}_{\rm{B}}^\infty {{\bf{\dot u}}_{\rm{Bs}}} $ (10)

式中:${\bf{K}}_{\rm{B}}^\infty $${\bf{C}}_{\rm{B}}^\infty $为粘弹性人工边界的刚度矩阵和阻尼矩阵。

下面推导斜入射的FBr公式。本文只考虑xy平面内的斜入射,推导也只适用于平面斜入射。假设P波以α角从左侧进入有限域。


图 1 P波斜入射示意图 Fig. 1 Diagram of P-wave oblique incidence

利用文献(赵密,2009)中的方法推导出在边界处所施加的等效荷载。

左边界:

$ \begin{array}{*{20}{l}} {F_{lx}^{ - x}(t)A\frac{{\lambda + 2G{{\sin }^2}\alpha }}{{{c_{\rm{P}}}}}({{\dot u}_{\rm{P}}}(t - \Delta {t_1}) - {B_1}{{\dot u}_{\rm{P}}}(t - \Delta {t_2})) + }\\ {\;\;\;\;\;\;\;\;\;\;\;\;A\frac{{2G\sin \beta \cos \beta }}{{{c_{\rm{S}}}}}{B_2}{{\dot u}_{\rm{p}}}(t - \Delta {t_3}) + {K_{\rm{N}}}({u_{\rm{P}}}(t - \Delta {t_1})\sin \alpha - }\\ {\;\;\;\;\;\;\;\;\;\;\;\;{B_1}{u_{\rm{P}}}(t - \Delta {t_2})\sin \alpha + {B_2}{u_{\rm{P}}}(t - \Delta {t_3})\cos \beta ) + }\\ {\;\;\;\;\;\;\;\;\;\;\;\;{C_{\rm{N}}}({{\dot u}_{\rm{P}}}(t - \Delta {t_1})\sin \alpha - {B_1}{{\dot u}_{\rm{P}}}(t - \Delta {t_2})\sin \alpha + {B_2}{{\dot u}_{\rm{P}}}(t - \Delta {t_3})\cos \beta )} \end{array} $ (11)
$ \begin{array}{*{20}{l}} {F_{ly}^{ - x}(t)A\frac{{2G\sin \alpha \cos \alpha }}{{{c_{\rm{P}}}}}({{\dot u}_{\rm{P}}}(t - \Delta {t_1}) + {B_1}{{\dot u}_{\rm{P}}}(t - \Delta {t_2})) + }\\ {\;\;\;\;\;\;\;\;\;\;\;A\frac{{G({{\sin }^2}\beta - {{\cos }^2}\beta )}}{{{c_{\rm{S}}}}}{B_2}{{\dot u}_{\rm{p}}}(t - \Delta {t_3}) + {K_{\rm{T}}}({u_{\rm{P}}}(t - \Delta {t_1})\cos \alpha + }\\ {\;\;\;\;\;\;\;\;\;\;\;{B_1}{u_{\rm{P}}}(t - \Delta {t_2})\cos \alpha + {B_2}{u_{\rm{P}}}(t - \Delta {t_3})\sin \beta ) + }\\ {\;\;\;\;\;\;\;\;\;\;\;{C_{\rm{T}}}({{\dot u}_{\rm{P}}}(t - \Delta {t_1})\cos \alpha + {B_1}{{\dot u}_{\rm{P}}}(t - \Delta {t_2})\cos \alpha + {B_2}{{\dot u}_{\rm{P}}}(t - \Delta {t_3})\sin \beta )} \end{array} $ (12)

式中:$\Delta {t_1}$$\Delta {t_2}$$\Delta {t_3}$分别是左边界的直接入射的P波、地面反射P波和SV波相对于原点到达左边界每1点的延迟时间(假定左人工边界上的任意1点的坐标为(0,y));B1B2分别为反射P波和反射S波与入射P幅值的比值;α为P波入射和反射角、β为S波反射角。

前边界上的等效荷载表达为:

$ \begin{array}{*{20}{l}} {F_{lx}^{ - z}(t)={K_{\rm{N}}}({u_{\rm{P}}}(t - \Delta {t_4}){\rm{sin}}\alpha - {B_1}{u_{\rm{P}}}(t - \Delta {t_5})\sin \alpha + {B_2}{u_{\rm{P}}}(t - \Delta {t_6})\cos \beta ) + }\\ {\;\;\;\;\;\;\;\;\;\;\;\;{C_{\rm{T}}}({{\dot u}_{\rm{P}}}(t - \Delta {t_4}){\rm{sin}}\alpha - {B_1}{{\dot u}_{\rm{P}}}(t - \Delta {t_5})\sin \alpha + {B_2}{{\dot u}_{\rm{P}}}(t - \Delta {t_6})\cos \beta )} \end{array} $ (13)
$ \begin{array}{*{20}{l}} {F_{ly}^{ - z}(t)={K_{\rm{N}}}({u_{\rm{P}}}(t - \Delta {t_4})\cos \alpha + {B_1}{u_{\rm{P}}}(t - \Delta {t_5})\cos \alpha + {B_2}{u_{\rm{P}}}(t - \Delta {t_6})\sin \beta ) + }\\ {\;\;\;\;\;\;\;\;\;\;\;\;{C_{\rm{T}}}\left( {{{\dot u}_{\rm{P}}}(t - \Delta {t_4})\cos \alpha + {B_1}{{\dot u}_{\rm{P}}}(t - \Delta {t_5})\cos \alpha + {B_2}{{\dot u}_{\rm{P}}}(t - \Delta {t_6})\sin \beta } \right)} \end{array} $ (14)
$ F_{lz}^{ - z}(t)=A\frac{\lambda }{{{c_{\rm{P}}}}}({\dot u_{\rm{P}}}(t - \Delta {t_4}) - {B_1}{\dot u_{\rm{P}}}(t - \Delta {t_5})) $ (15)

式中:$\Delta {t_4}$$\Delta {t_5}$$\Delta {t_6}$分别为前后边界的直接入射的P波、地面反射P波和SV波的延迟时间。

后边界上的等效荷载可表达为:

$ \left\{ \begin{array}{l} F_{lx}^z(t) = F_{lx}^{ - z}(t)\\ F_{ly}^z(t) = F_{ly}^{ - z}(t)\\ F_{lz}^z(t) = - F_{lz}^{ - z}(t) \end{array} \right. $ (16)

底边界上的等效荷载可表达为:

$ F_{lx}^{ - y}(t) = A\frac{{2G\sin \alpha \cos \alpha }}{{{c_{\rm{P}}}}}{\dot u_{\rm{P}}}(t - \Delta {t_7}) + {K_{\rm{T}}}({u_{\rm{P}}}(t - \Delta {t_7})\cos \alpha ) + {C_{\rm{T}}}({\dot u_{\rm{P}}}(t - \Delta {t_7})\cos \alpha ) $ (17)
$ F_{ly}^{ - y}(t) = A\frac{{\lambda + 2G{{\cos }^2}\alpha }}{{{c_{\rm{P}}}}}{\dot u_{\rm{P}}}(t - \Delta {t_7}) + {K_{\rm{N}}}({u_{\rm{P}}}(t - \Delta {t_7})\cos \alpha ) + {C_{\rm{N}}}({\dot u_{\rm{P}}}(t - \Delta {t_7})\cos \alpha ) $ (18)

式中:$\Delta {t_7}$为底边界入射P波的延迟时间。

为了验证P波斜入射方法的模拟精度,假设输入P波的角度α,输入峰值为0.1g的El Centro波;土体尺寸$80{\rm{m}} \times 80{\rm{m}} \times 80{\rm{m}}$,网格尺寸为1m,土密度2000kg/m3,弹性模量为100MPa,泊松比0.33。计算所得的地面位移是输入位移的2倍(廖振鹏,2002)(图 2),说明地震波输入方法精确。


图 2 输入位移和地面位移的对比 Fig. 2 The comparison between input displacement motion and the ground displacement response
4 沉管隧道三维计算模型
4.1 工程背景

本文建立三维分析模型,分析P波与隧道横截面形成的斜入射角度对沉管隧道的影响参考港珠澳大桥沉管隧道的尺寸,高12m,宽38m,采用3个节段进行分析,每个节段长22m。计算条件为:横向长度118m,高度78m,隧道上覆土高度6m,水深10m。隧道和土体的参数见表 1。接头的力学曲线(刘鹏等,2014)见图 34,输入时程为0.1g的El Centro波。

表 1 土和隧道参数 Table 1 Parameters of soil and immersed tunnel

图 3 接头处的轴向位移与轴力的关系 Fig. 3 Axial displacement-force curve of the joint

图 4 接头处的切向位移与剪力的关系 Fig. 4 Tangential displacement-force curve of the joint

沉管隧道三维计算模型以三维八节点六面体划分,共74532个单元,59256个节点。土体和隧道采用实体单元,海水采用势流体单元。土体底部和侧面均为粘弹性边界,海水与土体接触边界为流固耦合边界(Fluid-Structure),海水侧面采用无限域边界(Fluid infinite region),海水表面采用自由表面(Free surface)。

4.2 结果及分析

从结果(图 814)可以得知:入射角为40°时,隧道应力峰值(包括剪应力峰值、正应力峰值)最大;入射角为0°—40°时,隧道的应力峰值逐渐增大,入射角为40°—60°时,隧道的应力峰值逐渐减小。


图 5 沉管隧道有限元模型 Fig. 5 Finite element model of immersed tunnel

图 6 中间管节截面示意图 Fig. 6 The section of middle tunnel element

图 7 El Centro波加速度曲线 Fig. 7 Acceleration time history of EL Centro wave

图 8 中间管节截面剪应力峰值图 Fig. 8 The peak shearing stress curves of middle tunnel element section

图 9 中间管节截面正应力峰值图 Fig. 9 The peak normal stress curves of middle tunnel element section

图 10 中间管节截面竖向位移峰值图 Fig. 10 The vertical peak displacement curves of middle tunnel element section

图 11 中间管节隔墙截面剪应力峰值图 Fig. 11 The peak shearing stress curves of middle tunnel partition wall

图 12 中间管节隔墙截面正应力峰值图 Fig. 12 The peak normal stress curves of tunnel partition wall

图 13 中间管节隔墙截面水平位移峰值图 Fig. 13 The horizontal peak displacement curves of middle tunnel partition wall

图 14 中间管节隔墙截面竖向位移峰值图 Fig. 14 The vertical peak displacement curves of middle tunnel partition wall

下面分别讨论顶板、底板和侧向板以及隔墙的应力峰值。

根据剪应力峰值图(图 8)发现,隧道截面的左侧(入射侧)剪应力峰值远大于隧道截面的右侧,隧道顶板(底板)和隔墙交接处以及截面4个角处及附近剪应力峰值相对较大。正应力峰值图(图 9)显示隧道顶板正应力峰值最大,其它地方正应力峰值相对较小,顶板的正应力峰值大约为底板的2倍。根据位移峰值图(图 10)可知,左侧板的位移峰值比底板和顶板大。

根据图 1114可知,左侧隔墙比右侧隔墙的应力峰值大,其中左侧为地震波斜入射一侧;隔墙中间剪应力峰值比两端大;隔墙上端正应力峰值比下端大。

5 结论

本文采用粘弹性边界和地震波以等效荷载输入的方法,利用ADINA软件建立三维模型,考虑海水和隧道、海床的流固耦合,分析P波斜入射对隧道截面内力的影响,得到的结论可为沉管隧道设计提供参考。

(1)P波斜入射角为0°—40°时,隧道的应力峰值逐渐增大,入射角为40°—60°时,隧道的应力峰值逐渐减小。

(2)隧道截面4个转角处及隔墙与顶板、底板的连接处的剪应力最大,所以隧道4角附近及隔墙与顶板、底板的连接处为薄弱位置。

(3)隧道顶板正应力最大,约为底板的2倍;隧道截面左侧位移远大于隧道截面右侧。

参考文献
陈贵红, 2002.沉管隧道抗震数值分析.成都: 西南交通大学.
陈向红, 2013.大型水下隧道与附属竖井的地震响应研究.北京: 北京交通大学.
韩大建, 周阿兴, 1999. 沉管隧道地震响应分析的等效质点系模型探讨[J]. 华南理工大学学报(自然科学版), 27(11): 108-114. DOI:10.3321/j.issn:1000-565X.1999.11.017
李鹏, 2013.饱和地基中隧道纵向地震反应的数值分析.北京: 清华大学.
廖振鹏, 2002. 工程波动理论导论[M]. 2版. 北京: 科学出版社.
刘晶波, 谷音, 杜义欣, 2006. 一致粘弹性人工边界及粘弹性边界单元[J]. 岩土工程学报, 28(9): 1070-1075. DOI:10.3321/j.issn:1000-4548.2006.09.004
刘鹏, 丁文其, 金跃郎, 等, 2014. 沉管隧道接头三维非线性刚度力学模型[J]. 同济大学学报(自然科学版), 42(2): 232-237.
袁勇, 禹海涛, 燕晓, 等, 2016. 超长沉管隧道多点振动台试验模拟与分析[J]. 中国公路学报, 29(12): 157-165. DOI:10.3969/j.issn.1001-7372.2016.12.020
岳戈, 2010. ADINA流体与流固耦合功能的高级应用[M]. 北京: 人民交通出版社.
赵密, 2009.近场波动有限元模拟的应力型时域人工边界条件及其应用.北京: 北京工业大学.
周晓洁, 张越宇, 何颖, 等, 2017. 地震SV波斜入射下沉管隧道的地震响应分析[J]. 地震工程学报, 39(4): 600-608. DOI:10.3969/j.issn.1000-0844.2017.04.0600
Choi J. S., Lee J. S., Kim J. M., 2002. Nonlinear earthquake response analysis of 2-D underground structures with soil-structure interaction including separation and sliding at interface. In: Proceedings of the 15th ASCE Engineering Mechanics Conference. New York: Columbia, 1-8.
Huo H. B., Bobet A., 2003. Seismic design of cut and cover rectangular tunnels-evaluation of observed behavior of Dakai station during Kobe earthquake. In: Proceedings of 1st World Forum of Chinese Scholars in Geotechnical Engineering. Shanghai: Tongji University Press, 456-466.
Uenishi K., Sakurai S., 2000. Characteristic of the vertical seismic waves associated with the 1995 Hyogo-ken Nanbu (Kobe), Japan earthquake estimated from the failure of the Daikai Underground Station[J]. Earthquake Engineering and Structural Dynamics, 29(6): 813-821. DOI:10.1002/(ISSN)1096-9845