地震P波作用下饱和土中考虑剪切与转动变形的衬砌结构动力响应分析
张鸿1, 2，高谦1，徐斌2，徐满清2，陈学嘉2
(1. 北京科技大学 土木与环境工程学院，北京，100083；
2. 南昌工程学院 土木与建筑工程学院，江西 南昌，330099)
摘要：对隧洞周边土体采用饱和土Biot理论，隧洞衬砌采用考虑剪切和转动变形的曲线梁振动理论，分析地震P波作用下饱和土体中圆形隧洞衬砌的动力响应问题。对于饱和土体中的散射波场采用波函数展开法求解，对曲线梁的振动控制微分方程采用一般化的微分求积法(GDQM)求解。由饱和土体与衬砌接触处的位移协调条件，采用最小二乘法确定波函数未知系数项。计算结果表明：当入射波频率较低时，衬砌结构的入射面与背对面的动力响应几乎是对称的；随着入射频率的增加，衬砌结构的入射面与背对面不再具有对称性，且衬砌结构入射面的动力响应要大于背对面的动力响应。
关键词：隧洞衬砌；饱和土；微分求积法(DQM)；曲线梁；波函数展开法
中图分类号：TU 435 文献标志码：A 文章编号：1672-7207(2014)06-1943-09
Dynamic response of tunnel lining considering both rotary and shear deformations embedded in saturated soil with seismic P wave
ZHANG Hong1, 2, GAO Qian1, XU Bin2, XU Manqing2, CHEN Xuejia2
(1. School of Civil and Environmental Engineering, University of Science and Technology Beijing, Beijing 100083, China;
2. Department of Civil Engineering, Nanchang Institute of Technology, Nanchang 330099, China)
Abstract: Dynamic response of a circular tunnel lining embedded in the saturated soil subjected to seismic P wave was investigated. The surrounding medium of the tunnel was considered as a saturated porous medium and described by Biot’s theory. The tunnel lining was treated as curved beams described by a curved beam theory, which can account for both rotary and shear deformation. The scattered wave field in the porous medium was obtained by the wave function expansion method. The differential equations governing the vibration of curved beam were discretized by the general differential quadrature (GDQ) method. The wave function coefficients in the truncated series can be obtained by using the displacement compatibility conditions between the porous medium and the tunnel lining. The displacement compatibility condition between the porous medium and the tunnel was implemented by boundary least square collocation method. The results show that for low-frequency incident waves, the responses of the liner segments at the incident side and the shelter side are almost symmetrical. With the increase of frequency, the responses of the liner pieces at the incident side and the shelter side are not symmetrical anymore. The responses of the incident side tunnel are larger than those of the shelter side tunnel.
Key words: tunnel lining; saturated soil; general differential quadrature (GDQ) method; curved beam; wave function expansion method
地下工程结构不可避免遭受地震波的作用，因此，评估地下结构在地震荷载作用下的动力响应是抗震设计所关心的问题之一[1-2]。当隧洞周边岩石、土的地质条件较差而难以承担上部各类荷载压力时，工程中一般采用衬砌结构承担隧洞围压的压应力，因此，在地下工程抗震设计中，必须分析隧洞衬砌在地震波作用下的动力响应，一般采用解析法与数值方法对隧洞在地震波作用下的动力响应进行分析。Chen[3]采用波函数展开法分析了圆形隧洞中的线弹性衬砌在SH波作用下的动力响应。Shi等[4]采用映射法结合波函数展开法，对非圆形隧洞中的线弹性衬砌在各向异性介质中受SH波的动力响应进行了分析。Moore等[5]采用系列映射法，对地震荷载作用下的双线隧洞在无限介质空间中的动力响应进行了分析。数值方法中边界元法及有限元法的运用较常见，如Stamos等[6-7]将线弹性衬砌及其周边的介质土体分为不同的计算域，利用衬砌与周边介质土体接触处的变形协调条件，建立统一的边界元积分方程，分析衬砌的动力响应。对于衬砌周边的地基土体，傅方等[8-9]采用弹性土体模型。人们认识到地震荷载作用下地基失效的一个重要原因是：土体中的液相会造成应力、应变和位移的累积变化及突变，采用弹性土体模型难以反映该工程现象，因此，采用水、土二相耦合的饱和土体模型更接近于实际情况。周香莲等[10-11]采用复变函数法，利用Biot波动理论[12-14]，对饱和土中的圆形衬砌结构的散射和动应力集中问题进行了求解。Liang等[15]采用波函数展开法，分析了圆弧形饱和凹陷地形对SV波的散射问题。对于隧洞衬砌结构，赵武胜等[16-17]采用弹性体Lame方程描述其动力响应，未考虑衬砌拉伸和转动影响。衬砌结构不仅在沿压缩P波入射方向会产生较大的压缩变形，而且在与入射P波垂直方向必然存在剪切变形。Eisenberger等[18]指出剪切和转动变形对曲线梁的高阶自振频率有较大的影响。目前还未见地震波作用下饱和土体中考虑拉伸与转动变形的衬砌动力响应的报道。对于环形结构，一般采用曲线梁振动方程，如Walz等[19]采用有限元法(FEM)进行求解，但计算收敛速度较慢[20]。为求解曲线梁振动方程，与采用低阶线性插值函数的有限元法不同的微分求积法(DQM)，在计算域内采用离散节点的高阶线性插值函数，因此，能够得到更高精度的微分方程解[21]，微分求积法(DQM)也因此广泛应用于高阶微分方程的求解中[21]。但当离散节点较多时，微分求积法(DQM)难以确定加权系数，Shu等[22]对此提出了离散节点位置及数量可以是任意的一般化微分求积法(GDQM)。基于此，本文作者对饱和土体中沿隧洞轴线方向垂直入射的地震P波动力响应问题，将隧洞周边土体、衬砌结构分别利用饱和土Biot动力理论与Friedman等[23]提出的曲线梁理论来描述，且饱和土体中波场散射采用波函数展开法求解，曲线梁的振动控制方程采用GDQM法求解，分析地震P波作用下衬砌结构的动力响应。
1 饱和土与衬砌的控制方程及其解
1.1 饱和土体控制方程与一般解
根据饱和土体Biot理论[12-14]，土骨架、水二相耦合的饱和土体本构方程为
(1)
式中：λ和μ为Lame常数；α和M为Biot参数；σij为土体总应力；εij为土骨架上的应变张量；ui，wi和p分别为饱和土体土骨架、水相位移及孔压；
，
。
饱和土体中孔隙流体沿i方向的流量qi为
(2)
由土骨架及孔隙流体位移所表示的饱和孔隙介质的运动方程为
(3)
式中：ρb和ρf分别为土体和孔隙流体密度；η和k分别为饱和土体黏滞系数与渗透系数；
；α∞为饱和土体弯道系数；
为孔隙率；“*”为卷积符号；K(t)为高频时与流体惯性有关的黏滞系数修正 项[14, 24]。
为得到Biot方程的一般解，定义时间t与频域ω之间的Fourier变换对为
(4)
式中：上标“^”表示频域内的量。利用直角坐标系下的Helmholtz分解法，考虑饱和土体中存在2种压缩P波(P1波、P2波)，土骨架位移可采用势函数表示为
(5)
式中：
和
分别与P1波和P2波相对应的频域内标量势函数；
(k=1，2，3)为矢量势函数；eijk为Levi-Civita符号，且矢量势函数(k=1，2，3)满足如下条件：
(6)
由文献[25]可知饱和土体的控制方程中不仅有土骨架、水相的位移矢量，而且包括孔隙水压力，且孔隙水压力满足如下关系：
(7)
式中：Af和As为常量，由Biot控制方程确定。
由式(1)，(3)及式(5)可得如下Helmholtz方程：
(8)
式中：kf，ks和kt分别为与饱和土体中P1波、P2波及S波相对应的复波数，具体表达式见文献[26]。
在极坐标系下，饱和土体位移、应力及孔压可由势函数，及
表示为
(9)
1.2 考虑剪切和转动变形的曲线梁理论
本文中隧洞中的衬砌结构为无限长，且沿轴向方向为均质各向同性弹性体，入射P波与隧洞轴向方向垂直，如图1所示。因此，隧洞衬砌结构及其周边土体相互作用可作为平面应变问题，衬砌结构的振动采用曲线梁振动方程描述，而隧洞周边土体采用饱和土体Biot理论来描述。
根据文献[23]，考虑拉伸、剪切及转动变形的曲线梁，在图2所示曲线坐标系(x，y，z)中有效横截面轴力N(x，t)、剪力Q(x，t)及弯矩M(x，t)为
(10)
式中：u，v和分别为轴向、径向位移及转角；E和G分别为梁的弹性模量及剪切模量；A和I分别为横截面面积及转动惯性矩；kc为曲线梁剪力修正系数。
对如图2所示环形衬砌，其外半径R为定值，则dx=Rdθ。由文献[23]可知，曲线梁的运动方程为
(11)
图1 入射波作用下饱和土体中隧洞衬砌结构计算模型
Fig. 1 Calculation model of tunnel liner structure embedded in porous medium subjected to seismic waves
图2 曲线梁及坐标系
Fig. 2 Curved beam with attached curvilinear coordinate system subjected to seismic waves
1.3 曲线梁振动方程的微分求积法(GDQ)
设函数f(x)在闭区间(
)可采用多项式近似表示为
(12)
式中：xj为闭合域[a，b]中的系列配置点，且
，
，
。由文献[22]可知，函数f(x)对变量x的n阶导数在配置点xi值可采用闭合域内所有配置点的函数值的加权线性之和来表示，即
(13)
式中：
为f(x)对变量x的n阶导数的加权系数，且
(14)
对于中心角为Θ的曲线梁，任意θ进行无量纲化为ζ=θ/(2π)，则ζ∈[0，1]。对θ进行无量纲化后，应用式(13)，则频域内的式(11)可表示为
(15)
N个配置点在区间[0，1]的值为
；i=1，2，…，N。
由GDQ法同样可得任意配置点K处曲线梁有效横截面轴力、剪力及弯矩分别为
(16)
采用N个配置点将环形曲线梁离散成N-1段，边界配置点记为第1和第N点，域内配置点为i=2，3，…，N-1点，由式(15)整理后可得
(17)
式中：Ub和Ud分别表示边界节点、域内节点的位移；Fb和Fd分别表示作用在边界节点、域内节点的外力。
考虑到边界配置点第1和第N点重合，由此可得
(18)
将式(18)代入式(17)，可得
(19)
对式(19)两端乘以矩阵
，由此可得域内离散节点的位移为
(20)
式中：
；
；
。
2 饱和土体与衬砌协调方程的最小二乘法解
饱和土体中存在衬砌结构，在地震P波的入射下，饱和土体中会产生2种波场，即入射场和散射场，因此，土体中的总波场为
(21)
式中：上标I和S分别表示入射场和散射场。
入射波场满足Helmholtz方程式(8)，因此，散射波场也同样满足Helmholtz方程，则在极坐标系下，散射场的势函数可采用无穷级数表示为
(22)
式中：an，bn和cn为未知的复系数；
表示第二类Hankel函数。
若势函数的无穷级数项采用2Ns+1项截断，则饱和土体中的总应力、位移分别为
(23)
(24)
根据衬砌外表面与饱和土体接触处的应力平衡及式(24)，再由式(17)可得衬砌的位移为
(25)
饱和土体与衬砌结构在接触点处的变形应协调一致，因此，配置点ζi处饱和土体位移应与式(25)中衬砌相关点的位移相等。由于对势函数的无穷级数项采用有限项数截断，因此，衬砌离散配置点处的变形协调条件比未知复系数an，bn和cn的多，拟采用最小二乘法确定未知复系数an，bn和cn。若衬砌结构的配置点数为N，则在配置点处土体位移与衬砌位移差的平方为
(26)
将式(24)和(25)代入式(26)，然后，将
对未知复系数an，bn和cn求导，且导数为0，则可近似保证边界接触点的变形协调，即
；-Ns≤k≤Ns (27)
由此可得与未知复系数an，bn和cn有关的3×(2Ns+1)个复线性方程，解此方程组可得到未知系数。
3 算例分析
3.1 衬砌结构的GDQ法验算
为验证文中GDQ法对衬砌结构离散的正确性，分析两端固定、横截面为圆形、半径 r环形曲线梁的自振频率。可令式(20)的右端项为0，即可得到曲线梁各阶自振频率。利用GDQ法，曲线梁离散的配置点数为N=60，则u1=0，v1=0，
；u60=0，v60=0，
=0。曲线梁的几何、材料参数为：A=4×10-4 m2，R=10r，E=2.1×1011 Pa，ρ=7.85×103 kg/m3，泊松比ν=0.3，kc=1/1.1。
将曲线梁中心角分别为π/3，2π/3和π时，环形曲线梁前4阶无量纲化自振模态频率
与文献[16，27]中的结果相比较(见表1)，从计算结果可知GDQ法计算的正确性。
3.2 计算收敛性分析
为了确定波场散射无穷级数的收敛性，考察饱和土体中的闭合环形衬砌结构在入射P波作用下的动力响应，如图1所示。饱和土体参数为：μ=3.0×109 Pa，λ=1.0×109 Pa，ρs=2.5×103 kg/m3，ρf=1.0×103 kg/m3，=0.3，α=0.9，M=5.0×109 Pa，η=1.0×10-3 Pa·s，k= 1.0×10-12 m2。衬砌结构的材料及几何参数为：R=6.0 m，厚度h=0.3 m，E=3.5×1010 Pa，ρ=2.0×103 kg/m3，泊松比ν=0.3，kc=0.909。入射P波沿轴方向，频率ω=2.0π×10 s-1，振幅为
，入射波可采用势函数表示为
，根据文献[10]，可采用极坐标系下的无穷级数表示为
(其中，Jn(*)为第一类Bessel函数)。采用GDQ法，衬砌配置点数为N=100。图3所示为势函数的无穷级数项采用Ns为3，8和9项截断时，衬砌的无量纲化轴力
、剪力
及弯矩
分布情况。从图3可知：随级数项的增大，衬砌的轴力、剪力及弯矩都收敛，而且轴力收敛较快。
表1 有固定边界圆形横截面的曲线梁的无量纲化自振频率与已知文献结果的比较
Table 1 Comparison of present results with known results for non-dimensional natural frequency of a circular beam with uniform cross-section and with fixed–fixed boundary condition
3.3 饱和土体中衬砌结构动力响应
考察不同频率的P波沿x轴方向入射时，饱和土体中闭合衬砌结构的动力响应，饱和土体参数为：
μ=1.0×107 Pa，λ=4.0×107 Pa，ρs=2.5×103 kg/m3，ρf=1.0×103 kg/m3，=0.3，α=0.95，M=5.0×108 Pa，η=1.0×10-3 Pa·s，k=1.0×10-13 m/s。衬砌结构的材料及几何参数为：R=6.0 m，厚度h=0.3 m，E=1.0×1010 Pa，ρ=2.0×103 kg/m3，泊松比ν=0.3，kc=0.909。图4(a)~(f) 所示分别为入射P波频率f=ω/2π=2，20，50 Hz时，衬砌的无量纲化位移
，
，转角
，轴力
，剪力
及弯矩
分布情况。从图4可见：
(1) 入射波的频率对衬砌的动力响应有较大的影响，并且衬砌不同部位的动力响应也不相同。当入射波频率较低时，衬砌结构的入射面(90°~180°~270°范围内)与背对面(90°~0°~270°范围内)的动力响应几乎是关于90°~270°半径线对称的。随入射频率的增加，
衬砌结构的入射面与背对面不再沿90°~270°半径线呈对称性，且衬砌结构的入射面的动力响应要大于背对面的动力响应，高频入射时，该现象更明显。
图3 衬砌无量纲化轴力、剪力及弯矩随不同级数项Ns的分布规律
Fig. 3 Distribution law of dimensionless axial force, shear force and bending moment of liner with different numbers of dimension series
图4 不同入射P波频率作用下隧洞衬砌动力响应
Fig. 4 Dynamic response of tunnel liner subject to incident P wave with different frequencies
(2) 在与入射P波方向垂直的90°~270°半径线方向，衬砌轴向位移最大，而在入射P波的0°~180°半径方向，有较大的径向位移。
(3) 当入射P波频率较低时，衬砌轴力最大值在90°与270°处，随入射P波频率增大，轴力最大值在入射波的迎面90°~270°范围内。衬砌剪力最大值在120°~150°及210°~240°的半径线之间，而衬砌的弯矩最大值出现在与入射波的0°~180°半径方向。
4 结论
(1) 采用GDQM法可以分析曲线梁自振特征频率，而且与散射波场的级数展开法结合分析衬砌动力响应时，计算具有收敛性。
(2) 饱和土体中隧洞衬砌结构在地震入射P波作用下，衬砌结构中不仅会产生轴力，而且还有剪力及弯矩。
(3) 入射波频率对衬砌的动力响应有较大的影响，并且衬砌不同部位的动力响应也不相同。
(4) 当入射波频率较低时，衬砌结构的入射面与背对面的动力响应几乎是对称的。随入射频率的增加，衬砌结构的入射面与背对面不再具有对称性，且衬砌结构的入射面的动力响应要大于背对面的动力响应，高频入射时，该现象更明显。
参考文献：
[1] 张赛民, 陈灵君, 刘海飞, 等. 基于模型分块交叉移动的学习型模拟退火算法地震反演[J]. 中南大学学报(自然科学版), 2012, 43(3): 1040-1046.
ZHANG Saimin, CHEN Linjun, LIU Haifei, et al. Seismic inversion based on simulated annealing of learning and cross-moving with divided block model[J]. Journal of Central South University(Science and Technology), 2012, 43(3): 1040-1046.
[2] 陈振宗, 陈福胜, 周永川. 隧道受震破坏调查分析与修复管理案例[J]. 岩石力学与工程学报, 2004, 23(S2): 4816-4822.
CHEN Zhenzong, CHEN Fusheng, ZHOU Yongchuan. Investigation and reparation of an earthquake-damaged tunnel[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(S2): 4816-4822.
[3] Chen Y L. The analysis of elastic liner in a cylindrical tunnel subjected to SH-waves[J]. Journal of the Chinese Institute of Engineers, 1980, 3(1): 21-29.
[4] Shi S X, Han F, Wang Z Q, et al. The interaction of plane SH waves and non-circular cavity surfaced with lining in anisotropic media[J]. Applied Mathematics and Mechanics, 1996, 17(9): 855-867.
[5] Moore I D, Guan F. Three-dimensional dynamic response of lined tunnels due to incident seismic waves[J]. Earthquake Engineering and Structural Dynamics, 1996, 25(4): 357–369.
[6] Stamos A A, Beskos D E. 3-D seismic response analysis of long lined tunnels in half-space[J]. Soil Dynamics and Earthquake Engineering, 1996, 15(2): 111-118.
[7] Kattis S E, Beskos D E, Cheng A H D. 2D dynamic response of unlined and lined tunnels in poroelastic soil to harmonic body waves[J]. Earthquake Engineering and Structural Dynamics, 2003, 32(1): 97-110.
[8] 傅方, 赵成刚, 李伟华, 等. SV波斜入射下局部地形对隧道地震响应的影响[J]. 北京交通大学学报, 2012, 36(6): 79-84.
FU Fang, ZHAO Chenggang, LI Weihua, et al. Influence of local topographic on seismic response of tunnels subjected to obliquely incident SV waves[J]. Journal of Beijing Jiaotong University, 2012, 36(6): 79-84.
[9] 李忠献, 王洪龙, 李宁. 考虑接头力学特性的盾构隧道地震响应分析[J]. 地震工程与工程振动, 2012, 32(6): 166-173.
LI Zhongxian, WANG Honglong, LI Ning. Seismic response analysis of shield tunnel considering mechanical characteristics of joints[J]. Journal of Earthquake Engineering and Engineering Vibration, 2012, 32(6): 166-173.
[10] 周香莲, 周光明, 王建华. 饱和土中圆形衬砌结构对弹性波的散射[J]. 岩石力学与工程学报, 2005, 24(9): 1572-15761.
ZHOU Xianglian, ZHOU Guangming, WANG Jianhua. Scattering of elastic wave by circular cavity with lining in saturated soil[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(9): 1572-15761.
[11] 姜领发, 王建华, 周香莲. 半空间饱和土中圆形衬砌对弹性压缩波的散射[J]. 岩土力学, 2008, 29(2): 315-320.
JIANG Lingfa, WANG Jianhua, ZHOU Xianglian. Scattering around a circular lining in saturated poroelastic half-space under dilatational waves[J]. Rock and Soil Mechanics, 2008, 29(2): 315-320.
[12] Biot M A. Theory of propagation of elastic waves in a fluid-saturated porous solid: I.Low frequency range[J]. The Journal of the Acoustical Society of America, 1956, 28(2): 168-178.
[13] Biot M A. Theory of propagation of elastic waves in a fluid-saturated porous solid:II.Higher frequency range[J]. The Journal of the Acoustical Society of America, 1956, 28(2): 179-191.
[14] Biot M A. Mechanics of deformation and acoustic propagation in porous media[J]. Journal of Applied Physics, 1962, 33(4): 1482-198.
[15] Liang J, You H. Dynamic stiffness matrix of a poroelastic multi-layered site and its Green`s functions[J]. Earthquake Engineering and Engineering Vibration, 2004, 3(2): 273-282.
[16] 赵武胜, 何先志, 陈卫忠. 盾构隧道地震响应分析方法及工程应用[J]. 岩土力学, 2012, 33(8): 2415-2421.
ZHAO Wusheng, HE Xianzhi, CHEN Weizhong. Method for analyzing seismic response of shield tunnel and its application[J]. Rock and Soil Mechanics, 2012, 33(8): 2415-2421.
[17] 晏成明, 胡建平, 杨勇. 衬砌厚度对地下交叉隧道地震响应的影响[J]. 水利与建筑工程学报, 2012, 10(3): 174-177.
YAN Chengming, HU Jianping, YANG Yong. Impact of lining thickness on the seismic response of underground cross tunnel[J]. Journal of Water Resources and Architectural Engineering, 2012, 10(3): 174-177.
[18] Eisenberger M, Efraim E.In-plane vibrations of shear deformable curved beam[J]. International Journal for Numerical Methods in Engineering, 2001, 52(11): 1221-1234.
[19] Walz J E, et al. Accuracy of finite element approximation to structural problems[R]. NASA Technical Note TN D-5728, 1970. Washington, D.C., United States, 1970: 5-8.
[20] Prathap G, Babu C R. An isoparametric quadratic thick curved beam element[J]. International Journal for Numerical Methods in Engineering, 1986, 23(9): 1583-600.
[21] Bellman R E, Casti J. Differential quadrature and long term integration[J]. Journal of Mathematical Analysis and Applications, 1971, 34(2): 235-238.
[22] Shu C, Richards B E. Application of generalized differential quadrature to solve two dimensional incompressible Navier Stokes equations[J]. International Journal for Numerical Methods in Fluids, 1992, 15(7): 791-798.
[23] Friedman Z, Kosmatka J B. An accurate two-node finite element for shear deformable curved beams[J]. International Journal for Numerical Methods in Engineering, 1998, 41(3): 473-498.
[24] Johnson D L, Koplik J, Dashen R. Theory of dynamic permeability and tortuosity in fluid-saturated porous-media[J]. Journal of Fluid Mechanics, 1987, 176: 379-402.
[25] Bonnet G. Basic singular solutions for poroelastic medium in the dynamic range[J]. The Journal of the Acoustical Society of America, 1987, 82(5): 1758-62.
[26] ZHANG H, GAO Q, XU B. Scattering wave field around a cavity with circular cross-section embedded in saturated soil using the boundary element method[J]. Journal of Central South University of Technology, 2013, 20(11): 3296-3304.
[27] Yildirim V. A computer program for the free vibration analysis of elastic arcs[J]. Computers and Structures, 1997, 62(3): 475-485.
(编辑 陈灿华)
收稿日期：2013-06-06；修回日期：2013-08-10
基金项目：国家自然科学基金资助项目(50969007，51269021)；江西省自然科学基金资助项目(20114BAB206012)；江西省教育厅科技项目(GJJ11252，GJJ12629)
通信作者：徐斌(1971-)，男，湖北孝感人，博士，副教授，从事饱和土体动力学问题研究；电话：13361622340；E-mail：xmq418@163.com
