地下结构地震响应的计算模型1)

林 皋

(大连理工大学海岸与近海工程国家重点实验室，大连116024)(大连理工大学工程抗震研究所，大连116024)

地下结构地震响应的计算模型1)

林 皋2)

(大连理工大学海岸与近海工程国家重点实验室，大连116024)(大连理工大学工程抗震研究所，大连116024)

地震时地下结构在围岩的约束作用下发生变形，其动态特性与地面结构有很大不同.自二十世纪七八十年代以来，地下结构抗震设计与研究取得了很大进展.总的看来，工程设计中普遍采用的计算方法与设计导则大都建立在比较简单假定的基础上，实际的岩土介质条件都是十分复杂的.地下结构抗震研究的近期成果则表现在对地下结构动力分析中的波动散射问题提出了波函数展开法以及边界积分方程方法等多种计算方法.但计算相对复杂，在工程设计中的推广应用有一定困难.此文致力于地下结构计算模型的改进，使之具有良好的计算精度与效率，又便于工程应用.为此，提出了一种地下结构抗震响应分析的新的计算模型.模型具有较广泛的适应性，可以进行河谷、孔洞、地下铁道、隧洞等地下结构的散射与绕射分析.对于复杂层状的地质条件，提出了格林函数求解简便而有效的方法.数值算例论证了方法的精度和效率.

地下结构，地震响应，波的散射，格林函数，层状半空间地基

引言

地下工程结构在城市生产和生活中发挥着重要作用.地下结构承担着城市公共设施的许多重要功能，在通讯、能源供应(供电、供气)、供水排水、交通运输(地下铁道、海底隧道)、国防和民防(防空、防爆)等许多方面，发挥着现代工业生产和城镇生活中大动脉的作用.地下结构还常用作储油库、储气罐；地下结构也较常用来建设地下电站的厂房.随着国民经济和城市建设的发展，地下工程结构的高度和跨距也愈来愈大，这对地下结构的抗震设防提出了更高的要求.

历史上多次大地震中地下结构的破坏造成城市供电、供水、通讯的中断，引起国民经济的重大损失.著名事例[1]如：1906年美国旧金山大地震(震级 8.2)，由于供水中断造成火灾蔓延，破坏区域达12.2km2(4.7平方哩)，大火燃烧数日，所造成的损失占地震总损失的80%；1923年日本关东大地震(震级8.2)，在东京80%的城区中引起火灾，并有25条铁路隧洞受有破坏，横滨市由于排水总管破坏，除火灾外又引起次生洪水灾害.由于供电、供水等地下管线的破坏对民众的健康和安全关系重大，因此，这些地下工程结构又称为生命线工程.

20世纪80年代以前，国际上还缺乏地下结构或生命线工程抗击地震荷载的设计规范化条款.对一些重要的地下管道，日本采用震度法(地震系数法)进行设计，这并不能反映地下结构所遭受地震作用的实际情况.美国旧金山大地震(1906)、日本关东大地震 (1923)、新潟地震 (1964)、宫城县地震(1978)和日本海中部地震(1983)中管道等地下结构遭到很大程度的破坏，给城市的生产和生活造成了很大影响，促进了地下结构抗震研究的发展.这一时期(主要是二十世纪七八十年代)发表了大量研究论文[2]，涵盖地下结构地震震害和地震响应观测、影响因素、地下结构地震反应分析和计算模型的建立、设计要领和设计准则等各个方面.与此同时，日本进行了地下结构地震响应的观测和模型试验研究，其中比较有代表性的是田村和久保等的工作[37].

1977年 8月美国土木工程师学会生命线地震工程学科技委员会召开了地震工程学专业会议[2]，1981年第7届世界地震工程会议期间举行了生命线地震工程学最新进展委员会的会议.对地下结构或生命线工程抗震研究的情况进行了总结分析，深化了地下结构抗震特性的认识.

根据会议总结和会议相关论文的论述[27]，地下结构的地震响应具有以下特点：

(1)大多数现场观测资料表明，地震时地下管线和水下隧道无论轴向和侧向均随围岩发生着相同的运动.

(2)观测到的地下管线和水下隧道的轴向应变远较弯曲应变显著.管道转弯处的弯曲变形与直线段处的弯曲变形具有相同的量级.隧道的轴向刚度相对较大，对围岩变形起一定的限制作用.

(3)地下管线和隧道地震时产生的惯性力，对结构自身的地震响应只产生非常小的影响.

(4)地下管线的存在对围岩地震动的特性和扰动只有微弱的影响.

本文只讨论由地震波传播所引起的地下结构的地震响应分析.

1 地下结构地震响应计算模型的发展

由于地震时地下结构受围岩介质的约束作用而发生变形，引起震害，因此地下结构的地震响应分析应包括以下三方面的内容：

(1)围岩介质，特别是复杂岩土介质中自由场的地震波动特性.由于地下结构的存在对总体地震波动场的扰动一般较小，因此，可首先研究地下结构建造以前自由场的波动特性.

(2)地震时地下结构周围的地震波场特性.包括自由场和散射场两部分.重点需要计算的是由于地下结构不同的几何形状和动力刚度特性引起的散射波场特性.

(3)地震时地下结构和周围岩土介质间的动力相互作用特性.由于地下结构的质量和刚度与周围岩土介质相比是一个微量，因此地下结构自身的惯性力对其地震变形只起很微小的作用.但是地下结构的刚度，尤其是轴向刚度仍然对围岩的地震变形发生足够的反作用，这使地震时地下结构和围岩介质间产生有差异的地震变形.

二十世纪七八十年代以来，地下结构的抗震分析取得了很大的进展，但是，对上述三方面的研究深度，仍然不能使人满意.下文将对计算模型、规范标准以及文献中较广泛采用的有限元数值计算方法等方面，根据笔者的体会，作简要阐述.

自1972年Trifunac[8]发表关于半圆形河谷在SH平面波作用下的散射效应的论文以来，文献中关于地下结构地震响应的计算模型方面发展了很多数值计算方法.主要针对不同形状的河谷、冲积层、地下孔洞、地下管道、地下隧道等在SH波、SV波、P波和瑞利波作用下的散射效应提出了许多计算模型与方法.具体可参见有关文献阐述[912].比较有代表性并获得较广泛应用的数值计算方法主要有波函数展开法和边界积分方程方法两种.

波函数展开(WFE)法主要适用于均质空间或半空间散射问题的求解[1316].大多数求解的为二维散射问题，三维问题的求解局限于轴对称散射[17]，个别的也求解了非轴对称散射问题[18].

边界积分方程(BIE)方法又称边界元方法，其特点是只需进行散射体表面的离散，使问题降阶一维；同时满足无穷远处的辐射条件.但需获得计算域的基本解，增加了问题求解的困难，并增加了计算工作量.边界积分方程方法又可区分为直接BIE方法和间接BIE方法两种[1925].

将边界积分方程方法与Haskell和Thomson[26-27]提出的传递矩阵方法或Kausel[28]提出的刚度矩阵方法相结合，可以求解层状地基的格林函数，从而求解层状半空间的散射问题.

文献中也有的将边界积分方程法与格林函数的离散波函数展开法相结合进行散射问题的求解[2930].一般情况下，边界积分方程不适于求解复杂不均匀介质的散射问题，而有限元法则具有较广泛的适应性，对计算域介质的不均质问题处理非常方便，将两者相结合可以发挥各自的优越性.近场采用有限元法分析，远场辐射边界则采用边界积分方程求解[31].同理，也可以将有限元法与波函数展开法结合求解散射问题[32].

近年来，中国学者将前苏联学者Muskhelishvili所提出的求解弹性力学问题的复变函数法[33]应用于求解散射问题也取得了一定效果[3435].此外，文献中还提出了一些其他方法，本文将不再赘述.

二十世纪七八十年代以后，一些国家的抗震设防标准，主要是土木工程的抗震规范和核电站的抗震规范,包含了地下结构的抗震设计准则.但相对于房屋和桥梁等地面结构来说，条文较简略，难以满足实际抗震设计的要求[3638].

现有规范中关于地下结构的抗震分析方法，主要针对地下埋设管道的计算.对于无限长管，假设管道的地震变形与均匀无限介质中波传播产生的变形完全相同[3940]，据此计算管的轴向应变与弯曲应变.这一做法在核电结构等规范[41]中应用较为普遍.

一般规范标准中，很少有考虑复杂地基介质中的波传播特性对地下结构地震响应影响的规定.只有日本沉埋隧道的抗震设计规程提出了考虑沿线地基地质条件变化的地震动输入计算模型[38].其中用多质点弹簧阻尼体系来模拟围岩介质中地震波的传播.

文献中应用有限元方法进行地下结构的抗震分析较普遍.将地下结构与围岩介质划定一范围作为计算域，用有限元离散进行分析.这种方法可以考虑地下结构的各种几何形状和围岩介质的不同特性.但是由于缺乏有关地下结构抗震设计的导则和标准，因此，在网格的划分和布置方面比较自由，不同研究者的计算结果离散性较大.此外，在计算域范围、边界条件的选择等方面都含有一定的任意性.例如，侧面边界有的采用能量传递边界，有的采用黏性边界，有的则采用水平向可自由滑行的边界；底面边界有的采用黏性边界，有的则采用一定深度处的固定边界.地震动输入的模型也有一定的任意性.总体上，对自由场地震动的分析，散射场的计算以及地下结构与围岩介质的动力相互作用分析方面是否与实际情况相符都缺乏必要的检验，计算精度及其可靠性缺乏论证.将这样的计算结果用来进行地下结构的抗震设计与抗震安全评价是难以令人满意的.

地下空间的利用发展很快，适应各种需要的地下结构的高度、跨距等都逐渐增加，这对地下结构的抗震设计提出了新的要求.提高地下结构抗震研究的水平势在必行.现有的地下结构地震观测的资料还比较缺乏，并具有一定的局限性，多偏重于可能发生的地震破损形态，例如，沿隧道全线或断面内若干测点的变形和加速度等，对地震作用下不同地下结构的变形特点尚缺乏深入认识.用以进行新形势下地下结构的地震响应分析还有待进一步深入研究.

2 建议的地下结构地震响应分析的基本计算模型

本文的目的在于寻求一种理论上容易理解，计算不复杂，具有较好的计算精度，同时又便于在实际工程中应用的新方法，以便为地下结构抗震设计导则和规程的制订提供一定的参考.

下面从土与结构动力相互作用 (dynamic soilstructure interaction)的基本理论出发来建立地下结构地震响应分析的方程[42].应用子结构分析方法，土与结构动力相互作用的计算方程可表示如下(图1)

式中,ω为激励频率,ζ表示滞回阻尼系数.如需考虑高阶阻尼影响时，S的表达式可参见文献[43],进行地震作用分析时Ps=0，结构和无限地基的相互作用力Rb可表示为

图1 结构与无限地基系统Fig.1 Structure-unbounded soil system

将式(2)和式(3)代入式(1)可得

在图 1(b)中，将地基被挖去的部分表示为结构，则式(5)中矩阵Sbb-SbsSSsb表示地基被挖去部分的刚度(指被挖去部分凝聚于边界b上的刚度).这样式(5)左方括号内的矩阵+即表示地基未开挖前b边界的动力刚度；而即表示地基未开挖前b边界的自由场地震动[44].

图2 散射体的计算模型示意图Fig.2 Computational model for surface and subsurface irregularities and inhomogeneities

需要指出，式(8)的应用需要根据情况作出相应的变化.例如对于冲击层填充的河谷(图3)，计算公式相应地转换为

图3 填充河谷散射体Fig.3 Scattering by alluvial valley

图4 散射场计算边界的选择Fig.4 Selection of the outer boundaries for scattering analysis

可以指出，这里所建立的计算公式，既可适用于二维结构，也可适用于三维结构；既可适用于均质半

3 数值验算

本文将以均质半无限空间中的散射问题为主，通过算例检验计算模型的有效性.为计算均质半无限空间的动力刚度,比例边界有限元方法(SBFEM)[45]是有效的方法.SBFEM为半解析半数值性的求解方法，其径向可获得解析解，环向则达到有限元的计算精度.由于只需要进行边界面的离散，问题的维数降低一阶；并可自动满足无限远处的辐射条件.因此SBFEM兼有FEM和BEM的优点，也不需要基本解，计算简便.

3.1 半圆形河谷的散射

如图5所示,河谷边界ABC,计算边界DEFG，相应计算域大小为4a×2a(a为河谷半径).设围岩介质材料特性参数为：质量密度ρ0=1.5，剪切波速cs0=1，泊松比ν0=1/3；河谷中填充材料参数为：质量密度ρ1=1，剪切波速cs1=0.5，泊松比ν0=1/3.无量纲频率 η =(ωa)/(πcs)，波长 λ =2/η.计算采用式(9).网格离散如图5.各子域边界采用了3节点二次单元离散，共144单元，288节点.相应DEFG边界设48单元，97节点.入射波从底部垂直输入，计算结果精度高，与Luco和de Barros的解[46]相符性良好 (图 6).Luco和 de Barros，以及 Dravinski和Mossessian[47]采用的都是间接边界积分方程(IBIE)方法.

图5 半圆形河谷的计算图形Fig.5 Schematic view of semi-circular valley

图6 半圆形沉积河谷散射问题的检验Fig.6 Verificatio of the scattering by semi-circular valley

图6 半圆形沉积河谷散射问题的检验(续)Fig.6 Verificatio of the scattering by semi-circular valley(continued)

3.2 地下洞室衬砌的地震响应

研究地下隧道衬砌的地震响应.de Barros和Luco[48]求解了这一问题，采用间接边界积分方程方法模拟外围半空间土介质的作用，再结合Donnell的壳体理论计算隧洞衬砌外半径的波动响应.子域括号内的数字代表该子域相似中心的坐标.隧道结构如图7(a)所示.隧道埋深H=8.33r1=7.573a.衬砌材料质量密度 ρ0=2.24×103kg/m3，弹性模量E0=1.6×1010N/m2，泊松比 ν0=0.2.围岩地基介质材料质量密度ρ1=7.665×103kg/m3，弹性模量 E1=6.9×108N/m2，泊松比 ν0=0.45.位移按r=a进行规格化.假设SV波垂直入射，无量纲频率η=ωa/(πcs)=0.132.计算的衬砌响应如图8所示.可见与文献解的相符性较好.

图7 均质半无限地基中的隧道衬砌Fig.7 Tunnel and its liner embedded in uniform half-space

图8 SV波垂直入射作用下隧道衬砌外表面的动力响应Fig.8 Dynamic response at the outer radius of the tunnel liner subjected to vertically incident SV wave

4 复杂地基条件下地下结构的地震响应分析

其中的水平分层地基比较常见，下文将重点研究.对于分块不均质地基，当各块之间可以用同一相似中心加以描绘时(图9(b))，提出采用比例边界有限元的处理方法[49]，对近场含不规则不均质捕虏体等复杂地基(图9(c))则提出采用分步阻尼影响抽取法[50]进行求解.所提出的方法虽然是结合结构与地基动力相互作用分析进行阐述，但可适用于地下结构的抗震分析.

图9 不均质地基的计算模型Fig.9 Computational model for inhomogeneous unbounded soil

对各向同性和非各向同性层状半空间不均质地基的动力刚度及格林函数的求解，作者研究组从2012年起提出了基于积分变换的方法，并逐步得到完善[43,5153]，并可适用于横观各向同性的层状介质地基的求解.

下面介绍其基本思想和主要方程.

4.1 层状介质的波动方程

在圆柱坐标系下，层状介质(图10)中任一层以位移所表示的波动方程具有如下形式

图10 水平层状半空间Fig.10 Horizontally layered strata overlying half-space

对式(10)进行Fourier-Bessel变换，得出波数域的表达式.变换式如下

相应的反变换式为

式中k代表波数，其余系数的表达式为

式(15)中的两个子式，上面的子式对应于荷载与变形相对于X轴对称的情况；下面的子式则为荷载与变形相对于X轴成反对称的情况.其中Jn(kr)为柱贝塞尔函数.Green函数计算时，针对圆形单元面上作用的均布水平向荷载和竖向荷载，可选择n=1和n=0两种情况(图11)与之对应.当圆形单元半径趋于极限Δr→0时，均布荷载转化为集中水平力和集中竖向力作用的情况.即，在实际应用时，主要可限于n=1和n=0的两种情况.

图11 Green函数计算荷载单元Fig.11 Loaded element for evaluation of Green’s function

按式(12)进行Fourier-Bessel变换后得出频率--波数域内解耦的SV-P波动方程(变量为ur.uz)和SH波动方程(变量为uθ)[5152]，也可参见文献[28].

对横观各向同性介质中的波动，解耦方程同样成立，但相应系数改变如下

式中,上标m=1代表SV-P波动，m=2代表SH波动.Im为单位阵，m=1为2阶单位阵，m=2为1阶单位阵.

对各向同性介质(方程(17)与方程(18))

对于横观各向同性介质

为了进行微分方程(21)的求解，引入如下定义的应力矢量pm作为位移um的对偶变量

可以证明pm满足关系式[54]

以下的求解过程对m=1和m=2都适用，所以在下文中可将m省略.

利用式(26)，微分方程(21)可化为如下对耦形式

矩阵A,B,C,D和矩阵Kij(i,j=1,2)的变换关系如下

在状态空间，对耦微分方程(27)可以联立求解

齐次一阶线性常微分方程(29)的解为指数函数.

采用钟万勰提出的精细积分方法[55]可以将解进行精确计算.

对层状半空间中的任意一层，或其中的一层，设层厚为η=zb-za，则层上、下两端的位移和应力将成立如下关系

或写为

式中

为达到高精度的计算效果，T按如下方式进行计算[55]

式中，τ=η/b，b可取为任意整数.当取b=2N,N=20时，这相当于将η分成2N=1048557个微细薄层用式(34)计算，从而可以达到任意希望的精度，而计算工作量不大.文献[55]中给出了计算T的递推公式，应用方便.

为了进行层的结合和计算Green函数，解也宜表示成如下对耦变量的形式

式中矩阵F,G,Q也由T求出

对于水平分层介质，采用积分变换方法，可以求得波动方程在波数域的简便形式(式(17)和式(18)或式(19)和式(20)).进一步引入对偶变量，又可将其化为一阶齐次线性常微分方程(29)，其解可采用精细积分方法达到任意希望的精度.最后将解表示成对耦形式(36)，可便于分层之间的结合和Green函数的求解，说明如下.

4.2 相邻层的结合

当将相邻的两层：层1[za,zb]和层2[zb,zc]进行结合时，应用式(36)，可有

从式(38)中消去ub和pb，即可求得合并后新层c[za,zc]的力与位移关系

式中

4.3 Green函数的计算

层状半空间内部点Green函数的计算是地下结构抗震分析的核心部分.文献中广泛采用的刚度矩阵方法(sti ff ness matrix method)需要进行大型矩阵求解[28].这里提出对偶变量的转换方法，计算简捷方便，同时又准确高效.计算图形如图12和图13所示.图中源点(source)表示力的作用点，接受点(receiver)表示位移的计算点.已知源点和接受点的位置，可将层状半空间划分为3个子集，应用4.2节中阐述的方法可以分别建立各个子集形如式(36)所示的对耦变量方程，并相应求出其Fi,Qi,Gi矩阵(i=u,m,l，其中u表示上部，m表示中部，l表示下部)方程.Green函数的计算，基本上可划分为图13所示的6种工况.应用式(36)，即可求得各种工况的Green函数计算的矩阵表达式[53].

图12 层状半空间的Green函数Fig.12 Green’s function for multilayered half-space

图13 层状半空间不同工况的Green函数计算Fig.13 Computation of Green’s function for various cases of boundary conditions of multilayered half-space

下面给出不同工况的Green函数的矩阵表达式.下标a,b表示源点和接受点所在位置.

(1)工况(a)和(b)层状地基表面点的Green函数

工况(a)下部为刚性地基

工况(b)下部为弹性半空间

(2)工况(c)和(d)下部固定的层状地基内部点的格林函数

工况(c)

(3)工况(e)和(f)下部为弹性半空间的层状地基内部点的格林函数

工况(e)

所有矩阵的阶数对SV-P波为2×2，对SH波为1×1，计算很简便.最后可将波数域各种工况格林函数的计算，写成如下形式的统一的表达式(参见式(17)和式 (18)或式(19)和式(20),ur,uz和uB是解耦的，可分别求解

对图 11表示的格林函数计算的荷载形式，例如，半径为Δr的圆形单元上分布的Z向和X向均布荷载和，其相应的波数域的表达式为

在上述式中，令 Δr→0，πΔr2z→z或πΔr2x→x就可得出集中力z或x作用的表达式.再将式(53)、式(54)代入式(51)、式(52)并利用(13)进行反变换，就可得出圆柱坐标系物理域中格林函数的表达式，可统一写成

4.4 数值检验

研究均质半无限地基上表面弹性层中半圆形河谷的散射(图14).散射计算需要求解层状半空间中的格林函数.计算中介质材料特性参数假设为：均质半空间 ρ0=4/3，r0=0.333，cs0=2；表面弹性层 ρ1=1，r1=0.333，cs1=1；河谷中冲积层ρ2=0.667，r2=0.333，cs2=0.5.几何尺寸：河谷半径r，表面弹性层厚h=1.5r.无量纲频率η=ωa/(πcs).

在SV波垂直入射情况下半圆形充填河谷地表位移幅值计算结果如图15所示，可见与文献[56]的相符性很好.

图14 半无限地基上表面弹性层中的半圆形河谷的散射Fig.14 Scattering of a semi-circular valley embedded in the surface layer overlying elastic half-space

图15 SV波垂直入射情况下层状半空间中半圆形河谷地表位移幅值Fig.15 Amplitude of displacement response of the semi-circular valley embedded in layered half-space under vertically incident SV wave

5 地下结构地震响应分析的其他有关问题

5.1地下结构地震响应的时域分析

地下结构地震响应的时域计算应用运动方程(1)，关键是相互作用力 R=或R=的计算.由于无限域中R是激励频率的函数，需将其转换成时域的表达式，如式(51)所示(假设初值条件

式中，S∞(ω)代表无限域的动刚度(ω)或ω)；S∞(t)代表单位位移脉冲响应函数

式(58)只具有形式上的意义.由于S∞(ω)非平方可积，含有奇异性，式(58)不能直接用于脉冲响应函数的计算[57]，需进行必要的处理，将另文阐述.

5.2 地下结构地震响应的简化计算

地下结构也可采用简化的数值方法进行初步的抗震分析.以图16所示的地下结构为例来加以说明.图中虚线范围以内为计算域，域内可采用有限元等数值方法进行分析.根据第2节计算模型的推导，虚线所表示的边界应满足相互作用力的条件R=或R=.也可选择较大的计算域，并设立能量传递边界，通过边界进行自由场的地震动输入.文献中一般的能量传递边界主要建立在均质无限地基假定的基础上，不适于复杂地基条件下地下结构的地震响应分析.如采用时间、空间非耦合的局部人工透射边界，计算准确度并不理想，需要采用高阶局部人工边界，但高阶局部边界计算精度的提高是以计算自由度的增加为代价的.我们提出的层状地基耦合性(或全域性)时域传递边界[58]，不仅计算精度高，而且计算效率也高，可供参考.

图16 地下结构抗震分析的简化模型Fig.16 A simplifie model for seismic analysis of under-ground structures

6 结束语

提出一种基于结构与地基动力相互作用方程的地下结构抗震分析模型，可以很方便地进行地震波散射场与地下结构和围岩动力相互作用的计算.模型具有广泛的适用性.可以考虑层状半空间等复杂地基的影响.计算准确、高效.本文以二维均匀半空间中地下结构的地震响应为例，阐述地下结构抗震计算模型的设想.但是对复杂层状半空间地下结构的地震响应只提出了计算方法和相关方程还来不及展开，以后将陆续进行阐述.

致谢 博士生李志远、韩泽军和课题组李建波、胡志强等在本文完成中发挥了重要作用,在此表示感谢.

1 Ariman T,Gregory E Muleski.A review of the response of buried pipelines under seismic excitations.Earthquake Engineering&Structural Dynamics,1981,9(2)：133-152

2 Wang LRL,O’Rourke MJ.Overview of buried pipelines under seismic loading.Journal of Technical Councils of ASCE,1978,104(TC1)：121-130

3 Okamoto S,Tamura C.Behavior of subaqueous tunnel during earthquakes.International Journal of Earthquake Engineering&Structural Dynamics,1973,1(3)：253-266

4 Tamura C,Okamoto S,Kato K.Observations on dynamic strains of submerged tunnel during earthquakes.Bulletin of Earthquake Resistant Structure Research Center,University of Tokyo,1972,6：1-21

5 Tamura C.Dynamic behavior of submerged tunnel during earthquake.Transactions of the Japan Society of Civil Engineers,1975,7：197-200

6 Kubo K.Behavior of underground water pipes during an earthquake//Proceedings of the Fifth World Conference on Earthquake Engineering,Rome,Italy,1974,569-578

7 Kubo K,Katayama T,Ohashi A.Present state of lifeline earthquake engineering in Japan//The Current State of Knowledge of Lifeline Earthquake Engineering,Rome,Italy,1974,569-578

8 Trifunac MD.Scattering of plane SH waves by a semi-cylindrical canyon.Earthquake Engineering&Structural Dynamics,1972,1(3)：267-281

9 Luco JE,De Barros FCP.Seismic response of a cylindrical shell embedded in a layered viscoelastic half-space.I：Formulation.Earthquake Engineering&Structural Dynamics,1994,23(5)：553-567

10 Wong HL.E ff ects of surface topography and site conditions//Proc.of the EPRI Workshop on Strong Ground Motion Simulation and Earthquake Engineering Applications,Los Altos,California,April 30-May 3,1984

11 Sanchez-Sesma FJ.Site e ff ects on strong ground motion.Soil Dynamics and Earthquake Engineering,1987,6(2)：124-132

12 Von Thun JL.Earthquake Engineering and Soil Dynamics II—Recent Advances in Ground-Motion Evaluation.[s.l.]：ASCE,1988

13 El-Akily N,Datta SK.Response of a circular cylindrical shell to disturbances in a half-space.Earthquake Engineering&Structural Dynamics,1980,8(5)：469-477

14 Wong KC,Shah AH,Datta SK.Dynamic stresses and displacements in a buried tunnel.Journal of Engineering Mechanics,1985,111(2)：218-234

15 Eshraghi H,Dravinski M.Transient scattering of elastic waves by dipping layers of arbitrary shape part 1：Antiplane strain model.Earthquake Engineering&Structural Dynamics,1989,18(3)：397-415

16 Eshraghi H,Dravinski M.Transient scattering of elastic waves by dipping layers of arbitrary shape.Part 2：Plane strain model.Earthquake Engineering&Structural Dynamics,1989,18(3)：417-434

17 Lee VW.Three-dimensional di ff raction of plane P,SV&SH waves by a hemispherical alluvial valley.International Journal of Soil Dynamics and Earthquake Engineering,1984,3(3)：133-144

18 Eshraghi H,Dravinski M.Scattering of plane harmonic SH,SV,P and Rayleigh waves by non-axisymmetric three-dimensional canyons：A wave function expansion approach.Earthquake Engineering&Structural Dynamics,1989,18(7)：983-998

19 Wong HL,Jennings PC.E ff ects of canyon topography on strong ground motion.Bulletin of the Seismological Society of America,1975,65(5)：1239-1257

20 Zhang L,Chopra AK.Three-dimensional analysis of spatially varying ground motions around a uniform canyon in a homogeneous half-space.Earthquake Engineering&Structural Dynamics,1991,20(10)：911-926

21 Wong HL.E ff ect of surface topography on the di ff raction of P,SV,and Rayleigh waves.Bulletin of the Seismological Society of America,1982,72(4)：1167-1183

22 S´anchez-Sesma FJ,Campillo M.Di ff raction of P,SV,and Rayleigh waves by topographic features：A boundary integral formulation.Bulletin of the Seismological Society of America,1991,81(6)：2234-2253

23 Pedersen HA,Sanchez-Sesma FJ,Campillo M.Three-dimensional scattering by two-dimensional topographies.Bulletin of the Seismological Society of America,1994,84(4)：1169-1183

24 Luco JE,Wong HL,De Barros FCP.Three-dimensional response of a cylindrical canyon in a layered half-space.Earthquake Engineering&Structural Dynamics,1990,19(6)：799-817

25 Luco JE,De Barros FCP.Seismic response of a cylindrical shell embedded in a layered viscoelastic half-space.I：Formulation.Earthquake Engineering&Structural Dynamics,1994,23(5)：553-567

26 Haskell NA.The dispersion of surface waves on multilayered media.Bulletin of the Seismological Society of America,1953,43(1)：86-103

27 Thomson WT.Transmission of elastic waves through a stratifie solid medium.Journal of Applied Physics,1950,21(2)：89-93

28 Kausel E.Fundamental Solutions in Elastodynamics：A Compendium.New York,USA：Cambridge University Press,2006

29 Bouchon M.A simple,complete numerical solution to the problem of di ff raction of SH waves by an irregular surface.The Journal of the Acoustical Society of America,1985,77(1)：1-5

30 Kawase H.Time-domain response of a semi-circular canyon for incident SV,P,and Rayleigh waves calculated by the discrete wavenumber boundary element method.Bulletin of the Seismological Society of America,1988,78(4)：1415-1437

31 Mossessian TK,Dravinski M.Application of a hybrid method for scattering of P,SV,and Rayleigh waves by near-surface irregularities.Bulletin of the Seismological Society of America,1987,77(5)：1784-1803

32 Shah AH,Wong KC,Datta SK.Di ff raction of plane SH waves in a Half-space.Earthquake Engineering&Structural Dynamics,1982,10(4)：519-528

33 Muskhelishvili NI.Some Basic Problems of the Mathematical Theory of Elasticity.4th Edn.Netherlands,Noorho ff,New York,University of Groningen,1963(in Russian)

34 Liu D,Gai B,Tao G.Applications of the method of complex functions to dynamic stress concentrations.Wave Motion,1982,4(3)：293-304

35 Liu Q,Zhao M,Wang L.Scattering of plane P,SV or Rayleigh waves by a shallow lined tunnel in an elastic half space.Soil Dynamics and Earthquake Engineering,2013,49：52-63

36王汝梁.地下管线的震害、抗震验算、设计与措施.王珞珈等译.北京：地震出版社，2007(Wang Ruliang.Buried Lifeline Earthquake Engineering and Related Research.Wang LJ,et al.Trans.Beijing：Seismicity Press,2007(in Chinese))

37侯忠良.地下管线抗震.北京：学术书刊出版社，1990(Hou Zhongliang.Earthquake Resistance of Buried Pipelines.Beijing：Academic Press,1990(in Chinese))

38 Earthquake Engineering Committee.The Japan Society of Engineers.Earthquake Resistant Design for Civil Engineering Structures in Japan.1988

39 Sakurai A,Takahashi T.Dynamic stresses of underground pipelines during Earthquake//Proceedings of 4th World Conference on Earthquake Engineering,Chile,1969,811-895

40 Newmark NM.Problems in wave propagation in soil and rock//Proceedings of International Symposium in Wave Propagation and Dynamic Properties of Earth Materials,Albuquerque,New Mexico,1967,7-26

41 Yeh GCK.Seismic analysis of slender buried beams.Bulletin of the Seismological Society of America,1974,64(15)：1551-1562

42 Lin G,Li ZY,Li JB.A substructure replacement technique for the numerical solution of wave scattering problem(to be published)

43 Song C.The scaled boundary finit element method in structural dynamics.International Journal for Numerical Methods in Engineering,2009,77(8)：1139-1171

44 Lin G,Han Z,Li J.Soil–structure interaction analysis on anisotropic stratifie medium.G´eotechnique,2014,64(7)：570-580

45 Wolf JP.Dynamic Soil-Structure Interaction.[s.l.]：Prentice Hall Int.,1985

46 Luco JE,de Barros FCP.Three-dimensional response of a layered cylindrical valley embedded in a layered half-space.Earthquake Engineering&Structural Dynamics,1995,24(1)：109-125

47 Dravinski M,Mossessian TK.Scattering of plane harmonic P,SV,and Rayleigh waves by dipping layers of arbitrary shape.Bulletin of the Seismological Society of America,1987,77(1)：212-235

48 de Barros FCP,Luco JE.Seismic response of a cylindrical shell embedded in a layered viscoelastic half-space.II：Validation and numerical results.Earthquake Engineering&Structural Dynamics,1994,23(5)：569-580

49 Lin G,Du J,Hu Z.Earthquake analysis of arch and gravity dams including the e ff ects of foundation inhomogeneity.Frontiers of Architecture&Civil Engineering in China,2007,1(1)：41-50

50 YinX,LiJ,WuC,etal.ANSYSimplementationofdampingsolvent stepwise extraction method for nonlinear seismic analysis of large 3-D structures.Soil Dynamics&Earthquake Engineering,2013,44(1)：139-152

51林皋，韩泽军，李建波.层状地基任意形状刚性基础动力响应求解.力学学报，2012，44(6)：1016-1027(Lin Gao,Han Zejun,Li Jianbo.Solution of the dynamic response of rigid foundation of arbitrary shape on multi-layered soil.Chinese Journal of Theoretical and Applied Mechanics,2012,44(6)：1016-1027(in Chinese))

52 Lin G,Han Z,Li J.General formulation and solution procedure for harmonic response of rigid foundation on isotropic as well as anisotropic multilayered half-space.Soil Dynamics&Earthquake Engineering,2015,70：48-59

53 Lin G,Han Z,Lu S,et al.Wave motion equation and the Green’s function for transverse isotropic multi-layered half-space.(to be published)

54 Zhong W,Williams FW,Bennett PN.Extension of the wittrickwilliams algorithm to mixed variable systems.Journal of Vibration&Acoustics,1997,119(3)：334-340

55 Zhong WX.On precise integration method.Journal of Computational&Applied Mathematics,2004,163(1)：59-78

56 de Barros FCP,Luco JE.Amplificatio of obliquely incident waves by a cylindrical valley embedded in a layered half-space.Soil Dynamics&Earthquake Engineering,1995,14(3)：163-175

57 Wolf JP,Song C.Finite-element Modelling of Unbounded Media.Chichester：Wiley,1996

58 Lin G,Lu S,Liu J.Transmitting boundary for transient analysis of wavepropagationin layered media formulated based on acceleration unit-impulse response.Soil Dynamics&Earthquake Engineering,2016,90：494-509

A COMPUTATIONAL MODEL FOR SEISMIC RESPONSE ANALYSIS OF UNDERGROUND STAUCTURES1)

Lin Gao2)
(State Key Laboratory of Coastal and O ff shore Engineering，Dalian University of Technology，Dalian 116024，China)(Institute of Earthquake Engineering，Dalian University of Technology，Dalian 116024，China)

Under earthquake excitation,the deformation of underground structures is restricted by the surrounding soil media.The dynamic behavior of it displays quite di ff erently from aboveground structures.Considerable progress has been made on the design and research of the underground structures since seventies-eighties of twentieth century.However,in the main,the widely used computational methods and guidelines in engineering design practice are based on rather simple assumptions,in reality,the actual soil conditions might be much more complex than ideally boundary conditions.Recent achievements of earthquake research on underground structures lie in the development of various computational methods for wave scattering problems of underground structures,such as the wave function expansion method,the boundary integral equation method etc.As the computation is somewhat complex,which impedes its application and dissemination in the engineering design practice.The author devotes himself to the improvement of the computational model for seismic analysis of underground structures,such that it achieves higher accuracy and efficiency,meanwhile it proves to be convenient for engineering design.To this end,a new model for seismic response analysis of underground structuresis proposed.The model is versatile to deal with wave scattering and di ff raction by canyons,subsurface cavities,subways and tunnels etc.In case of the presence of complex soil conditions like the layered half-space,a simple and e ff ective technique is developed for the evaluation of Green’s functions.Numerical examples are provided to validate the accuracy and efficiency of the proposed approach.

underground structure,seismic response,wave scattering,Green’s function,layered half-space.

TU93

：A

10.6052/0459-1879-16-301

2016–10–28 收稿，2017–01–25 录用,2017–03–02 网络版发表.

1)国家重点研发计划资助项目(2016YFB0201001).

2)林皋，中国科学院院士，教授，主要研究方向：水利工程与核电工程抗震.E-mail：gaolin@dlut.edu.cn

林皋.地下结构地震响应的计算模型.力学学报,2017,49(3)：528-542

Lin Gao.A computational model for seismic response analysis of underground stauctures.Chinese Journal of Theoretical and Applied Mechanics,2017,49(3)：528-542