A COMPUTATIONAL MODEL FOR SEISMIC RESPONSE ANALYSIS OF UNDERGROUND STAUCTURES
-
摘要: 地震时地下结构在围岩的约束作用下发生变形，其动态特性与地面结构有很大不同。自二十世纪七八十年代以来，地下结构抗震设计与研究取得了很大进展。总的看来，工程设计中普遍采用的计算方法与设计导则大都建立在比较简单假定的基础上，实际的岩土介质条件都是十分复杂的。地下结构抗震研究的近期成果则表现在对地下结构动力分析中的波动散射问题提出了波函数展开法以及边界积分方程方法等多种计算方法。但计算相对复杂，在工程设计中的推广应用有一定困难。此文致力于地下结构计算模型的改进，使之具有良好的计算精度与效率，又便于工程应用。为此，提出了一种地下结构抗震响应分析的新的计算模型。模型具有较广泛的适应性，可以进行河谷、孔洞、地下铁道、隧洞等地下结构的散射与绕射分析。对于复杂层状的地质条件，提出了格林函数求解简便而有效的方法。数值算例论证了方法的精度和效率。Abstract: Under earthquake excitation, the deformation of underground structures is restricted by the surrounding soil media. The dynamic behavior of it displays quite differently 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 structures is proposed. The model is versatile to deal with wave scattering and diffraction 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 effective technique is developed for the evaluation of Green's functions. Numerical examples are provided to validate the accuracy and efficiency of the proposed approach.
-
Keywords:
- underground structure /
- seismic response /
- wave scattering /
- Green's function /
- layered half-space
-
引言
地下工程结构在城市生产和生活中发挥着重要作用.地下结构承担着城市公共设施的许多重要功能，在通讯、能源供应(供电、供气)、供水排水、交通运输(地下铁道、海底隧道)、国防和民防(防空、防爆)等许多方面，发挥着现代工业生产和城镇生活中大动脉的作用.地下结构还常用作储油库、储气罐；地下结构也较常用来建设地下电站的厂房.随着国民经济和城市建设的发展，地下工程结构的高度和跨距也愈来愈大，这对地下结构的抗震设防提出了更高的要求.
历史上多次大地震中地下结构的破坏造成城市供电、供水、通讯的中断，引起国民经济的重大损失.著名事例[1]如：1906年美国旧金山大地震(震级8.2)，由于供水中断造成火灾蔓延，破坏区域达12.2 km2 (4.7平方哩)，大火燃烧数日，所造成的损失占地震总损失的80%；1923年日本关东大地震(震级8.2)，在东京80%的城区中引起火灾，并有25条铁路隧洞受有破坏，横滨市由于排水总管破坏，除火灾外又引起次生洪水灾害.由于供电、供水等地下管线的破坏对民众的健康和安全关系重大，因此，这些地下工程结构又称为生命线工程.
20世纪80年代以前，国际上还缺乏地下结构或生命线工程抗击地震荷载的设计规范化条款.对一些重要的地下管道，日本采用震度法(地震系数法)进行设计，这并不能反映地下结构所遭受地震作用的实际情况.美国旧金山大地震(1906)、日本关东大地震(1923)、新潟地震(1964)、宫城县地震(1978) 和日本海中部地震(1983) 中管道等地下结构遭到很大程度的破坏，给城市的生产和生活造成了很大影响，促进了地下结构抗震研究的发展.这一时期(主要是二十世纪七八十年代)发表了大量研究论文[2]，涵盖地下结构地震震害和地震响应观测、影响因素、地下结构地震反应分析和计算模型的建立、设计要领和设计准则等各个方面.与此同时，日本进行了地下结构地震响应的观测和模型试验研究，其中比较有代表性的是田村和久保等的工作[3-7].
1977年8月美国土木工程师学会生命线地震工程学科技委员会召开了地震工程学专业会议[2]，1981年第7届世界地震工程会议期间举行了生命线地震工程学最新进展委员会的会议.对地下结构或生命线工程抗震研究的情况进行了总结分析，深化了地下结构抗震特性的认识.
根据会议总结和会议相关论文的论述[2-7]，地下结构的地震响应具有以下特点：
(1) 大多数现场观测资料表明，地震时地下管线和水下隧道无论轴向和侧向均随围岩发生着相同的运动.
(2) 观测到的地下管线和水下隧道的轴向应变远较弯曲应变显著.管道转弯处的弯曲变形与直线段处的弯曲变形具有相同的量级.隧道的轴向刚度相对较大，对围岩变形起一定的限制作用.
(3) 地下管线和隧道地震时产生的惯性力，对结构自身的地震响应只产生非常小的影响.
(4) 地下管线的存在对围岩地震动的特性和扰动只有微弱的影响.
本文只讨论由地震波传播所引起的地下结构的地震响应分析.
1. 地下结构地震响应计算模型的发展
由于地震时地下结构受围岩介质的约束作用而发生变形，引起震害，因此地下结构的地震响应分析应包括以下三方面的内容：
(1) 围岩介质，特别是复杂岩土介质中自由场的地震波动特性.由于地下结构的存在对总体地震波动场的扰动一般较小，因此，可首先研究地下结构建造以前自由场的波动特性.
(2) 地震时地下结构周围的地震波场特性.包括自由场和散射场两部分.重点需要计算的是由于地下结构不同的几何形状和动力刚度特性引起的散射波场特性.
(3) 地震时地下结构和周围岩土介质间的动力相互作用特性.由于地下结构的质量和刚度与周围岩土介质相比是一个微量，因此地下结构自身的惯性力对其地震变形只起很微小的作用.但是地下结构的刚度，尤其是轴向刚度仍然对围岩的地震变形发生足够的反作用，这使地震时地下结构和围岩介质间产生有差异的地震变形.
二十世纪七八十年代以来，地下结构的抗震分析取得了很大的进展，但是，对上述三方面的研究深度，仍然不能使人满意.下文将对计算模型、规范标准以及文献中较广泛采用的有限元数值计算方法等方面，根据笔者的体会，作简要阐述.
自1972年Trifunac[8]发表关于半圆形河谷在SH平面波作用下的散射效应的论文以来，文献中关于地下结构地震响应的计算模型方面发展了很多数值计算方法.主要针对不同形状的河谷、冲积层、地下孔洞、地下管道、地下隧道等在SH波、SV波、P波和瑞利波作用下的散射效应提出了许多计算模型与方法.具体可参见有关文献阐述[9-12].比较有代表性并获得较广泛应用的数值计算方法主要有波函数展开法和边界积分方程方法两种.
波函数展开(WFE)法主要适用于均质空间或半空间散射问题的求解[13-16].大多数求解的为二维散射问题，三维问题的求解局限于轴对称散射[17]，个别的也求解了非轴对称散射问题[18].
边界积分方程(BIE)方法又称边界元方法，其特点是只需进行散射体表面的离散，使问题降阶一维；同时满足无穷远处的辐射条件.但需获得计算域的基本解，增加了问题求解的困难，并增加了计算工作量.边界积分方程方法又可区分为直接BIE方法和间接BIE方法两种[19-25].
将边界积分方程方法与Haskell和Thomson [26-27]提出的传递矩阵方法或Kausel[28]提出的刚度矩阵方法相结合，可以求解层状地基的格林函数，从而求解层状半空间的散射问题.
文献中也有的将边界积分方程法与格林函数的离散波函数展开法相结合进行散射问题的求解[29-30].一般情况下，边界积分方程不适于求解复杂不均匀介质的散射问题，而有限元法则具有较广泛的适应性，对计算域介质的不均质问题处理非常方便，将两者相结合可以发挥各自的优越性.近场采用有限元法分析，远场辐射边界则采用边界积分方程求解[31].同理，也可以将有限元法与波函数展开法结合求解散射问题[32].
近年来，中国学者将前苏联学者Muskhelishvili所提出的求解弹性力学问题的复变函数法[33]应用于求解散射问题也取得了一定效果[34-35].此外，文献中还提出了一些其他方法，本文将不再赘述.
二十世纪七八十年代以后，一些国家的抗震设防标准，主要是土木工程的抗震规范和核电站的抗震规范, 包含了地下结构的抗震设计准则.但相对于房屋和桥梁等地面结构来说，条文较简略，难以满足实际抗震设计的要求[36-38].
现有规范中关于地下结构的抗震分析方法，主要针对地下埋设管道的计算.对于无限长管，假设管道的地震变形与均匀无限介质中波传播产生的变形完全相同[39-40]，据此计算管的轴向应变与弯曲应变.这一做法在核电结构等规范[41]中应用较为普遍.
一般规范标准中，很少有考虑复杂地基介质中的波传播特性对地下结构地震响应影响的规定.只有日本沉埋隧道的抗震设计规程提出了考虑沿线地基地质条件变化的地震动输入计算模型[38].其中用多质点弹簧阻尼体系来模拟围岩介质中地震波的传播.
文献中应用有限元方法进行地下结构的抗震分析较普遍.将地下结构与围岩介质划定一范围作为计算域，用有限元离散进行分析.这种方法可以考虑地下结构的各种几何形状和围岩介质的不同特性.但是由于缺乏有关地下结构抗震设计的导则和标准，因此，在网格的划分和布置方面比较自由，不同研究者的计算结果离散性较大.此外，在计算域范围、边界条件的选择等方面都含有一定的任意性.例如，侧面边界有的采用能量传递边界，有的采用黏性边界，有的则采用水平向可自由滑行的边界；底面边界有的采用黏性边界，有的则采用一定深度处的固定边界.地震动输入的模型也有一定的任意性.总体上，对自由场地震动的分析，散射场的计算以及地下结构与围岩介质的动力相互作用分析方面是否与实际情况相符都缺乏必要的检验，计算精度及其可靠性缺乏论证.将这样的计算结果用来进行地下结构的抗震设计与抗震安全评价是难以令人满意的.
地下空间的利用发展很快，适应各种需要的地下结构的高度、跨距等都逐渐增加，这对地下结构的抗震设计提出了新的要求.提高地下结构抗震研究的水平势在必行.现有的地下结构地震观测的资料还比较缺乏，并具有一定的局限性，多偏重于可能发生的地震破损形态，例如，沿隧道全线或断面内若干测点的变形和加速度等，对地震作用下不同地下结构的变形特点尚缺乏深入认识.用以进行新形势下地下结构的地震响应分析还有待进一步深入研究.
2. 建议的地下结构地震响应分析的基本计算模型
本文的目的在于寻求一种理论上容易理解，计算不复杂，具有较好的计算精度，同时又便于在实际工程中应用的新方法，以便为地下结构抗震设计导则和规程的制订提供一定的参考.
下面从土与结构动力相互作用(dynamic soil-structure interaction)的基本理论出发来建立地下结构地震响应分析的方程[42].应用子结构分析方法，土与结构动力相互作用的计算方程可表示如下(图 1)
$$\begin{array}{l} \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{M}}_{ss}}} & {{\mathit{\boldsymbol{M}}_{sb}}}\\ {{\mathit{\boldsymbol{M}}_{bs}}} & {{\mathit{\boldsymbol{M}}_{bb}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\ddot u}}_s^t}\\ {\mathit{\boldsymbol{\ddot u}}_b^t} \end{array}} \right\} + \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{C}}_{ss}}} & {{\mathit{\boldsymbol{C}}_{sb}}}\\ {{\mathit{\boldsymbol{C}}_{bs}}} & {{\mathit{\boldsymbol{C}}_{bb}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\dot u}}_s^t}\\ {\mathit{\boldsymbol{\dot u}}_b^t} \end{array}} \right\} + \\ \quad \quad \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{K}}_{ss}}} & {{\mathit{\boldsymbol{K}}_{sb}}}\\ {{\mathit{\boldsymbol{K}}_{bs}}} & {{\mathit{\boldsymbol{K}}_{bb}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{u}}_s^t}\\ {\mathit{\boldsymbol{u}}_b^t} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{P}}_s}}\\ { - {\mathit{\boldsymbol{R}}_b}} \end{array}} \right\} \end{array}$$ (1) 式中， ${\pmb M}$ , ${\pmb C}$ , ${\pmb K }$ 分别表示结构的质量、阻尼和刚度矩阵； $\ddot {\pmb u}^t$ , $\dot{\pmb u}^t$ , ${\pmb u}^t$ 表示结构的绝对加速度、速度和位移；下标 $b$ 表示结构和无限地基交界面上的自由度.下标 $s$ 则表示结构其余部分的自由度； ${\pmb P}_s $ 表示结构上的外加荷载； ${\pmb R}_b $ 表示结构和无限地基的相互作用力.由于 ${\pmb R}_b $ 为激励频率的函数，方程(1) 一般在频率域进行求解.对于有限域的结构动力刚度 ${\pmb S}$ 一般可将其表示为如下的低阶频率展开式
$$\mathit{\boldsymbol{S}}\left( \omega \right) = \mathit{\boldsymbol{K}}\left( {1 + 2\zeta {\rm{i}}} \right) - {\omega ^2}\mathit{\boldsymbol{M}}$$ (2) 式中, $\omega $ 为激励频率, $\zeta $ 表示滞回阻尼系数.如需考虑高阶阻尼影响时， ${\pmb S}$ 的表达式可参见文献[43], 进行地震作用分析时 ${\pmb P}_s ={\bf 0}$ ，结构和无限地基的相互作用力 ${\pmb R}_b $ 可表示为
$${\mathit{\boldsymbol{R}}_b}\left( \omega \right) = \mathit{\boldsymbol{S}}_{bb}^g\left( \omega \right)\left( {\mathit{\boldsymbol{u}}_b^t\left( \omega \right) - \mathit{\boldsymbol{u}}_b^g\left( \omega \right)} \right)$$ (3) 式中， ${\pmb S}_{bb}^g \left(\omega \right)$ 表示界面 $b$ 处无限地基的动力刚度，为激励频率的函数； ${\pmb u}_b^g \left(\omega \right)$ 表示结构不存在时，界面 $b$ 处的散射场地震动.为简便起见，下文式中的 $\left(\omega \right)$ 将予以省略.
将式(2) 和式(3) 代入式(1) 可得
$$\left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{S}}_{ss}}} & {{\mathit{\boldsymbol{S}}_{sb}}}\\ {{\mathit{\boldsymbol{S}}_{bs}}} & {{\mathit{\boldsymbol{S}}_{bb}} + \mathit{\boldsymbol{S}}_{bb}^g} \end{array}} \right]\left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{u}}_s^t}\\ {\mathit{\boldsymbol{u}}_b^t} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{l}} {\;\;{\bf{0}}}\\ {\mathit{\boldsymbol{S}}_{bb}^g\mathit{\boldsymbol{u}}_b^g} \end{array}} \right\}$$ (4) 从式(4) 的上、下两式中消去 ${\pmb u}_s^t $ ，可有
$$\left( {{\mathit{\boldsymbol{S}}_{bb}} - {\mathit{\boldsymbol{S}}_{bs}}\mathit{\boldsymbol{S}}_{ss}^{ - 1}{\mathit{\boldsymbol{S}}_{sb}} + \mathit{\boldsymbol{S}}_{bb}^g} \right)\mathit{\boldsymbol{u}}_b^t = \mathit{\boldsymbol{S}}_{bb}^g\mathit{\boldsymbol{u}}_b^g$$ (5) 在图 1(b)中，将地基被挖去的部分表示为结构，则式(5) 中矩阵 ${\pmb S}_{bb} -{\pmb S}_{bs} {\pmb S}_{ss} ^{ -1}{\pmb S}_{sb} $ 表示地基被挖去部分的刚度 ${\pmb S}_{bb}^e $ (指被挖去部分凝聚于边界 $b$ 上的刚度).这样式(5) 左方括号内的矩阵 ${\pmb S }_{bb}^e + {\pmb S}_{bb}^g $ 即表示地基未开挖前 $b$ 边界的动力刚度；而 ${\pmb u}_b^t $ 即表示地基未开挖前 $b$ 边界的自由场地震动 ${\pmb u}_b^f$ [44].
亦即
$$\mathit{\boldsymbol{S}}_{bb}^f\mathit{\boldsymbol{u}}_b^f = \mathit{\boldsymbol{S}}_{bb}^g\mathit{\boldsymbol{u}}_b^g$$ (6) 式中
$$\left. {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{S}}_{bb}^f = \mathit{\boldsymbol{S}}_{bb}^e + \mathit{\boldsymbol{S}}_{bb}^g}\\ {\mathit{\boldsymbol{S}}_{bb}^e = {\mathit{\boldsymbol{S}}_{bb}} - {\mathit{\boldsymbol{S}}_{bs}}\mathit{\boldsymbol{S}}_{ss}^{ - 1}{\mathit{\boldsymbol{S}}_{sb}}} \end{array}} \right\}$$ (7) 方程(6) 表示无限地基的边界 $b$ 在未开挖前的自由场地震动 ${\pmb u}_b^f $ 与开挖后该处散射场地震动 ${\pmb u}_b^g $ 的关系. Wolf在1985年发表了这一表达式[45], 并被许多文献广泛应用于结构与地基的动力相互作用分析，但仅限于地上结构.这一思想可以有效地应用于地下结构的散射分析，如图 2(a)和图 2(b)所示.根据式(7)，散射场的计算公式为
$$\left( {\mathit{\boldsymbol{S}}_{bb}^f - \mathit{\boldsymbol{S}}_{bb}^e} \right)\mathit{\boldsymbol{u}}_b^g = \mathit{\boldsymbol{S}}_{bb}^f\mathit{\boldsymbol{u}}_b^f$$ (8) ![]()
图 2 散射体的计算模型示意图Figure 2. Computational model for surface and subsurface irregularities and inhomogeneities式中， ${\pmb u}_b^g $ 代表散射波场的位移.这样将波的散射问题的计算转化为无限地基自由场地震动 ${\pmb u}_b^f $ 与动力刚度矩阵 ${\pmb S}_{bb}^f $ 的计算.由于这时边界条件比较规则(不存在空洞和不规则嵌入体)， ${\pmb S}_{bb}^f$ 和 ${\pmb u}_b^f $ 的求解将比 ${\pmb S}_{bb}^g$ 和 ${\pmb u}_b^g $ 的求解更为简便易行.
需要指出，式(8) 的应用需要根据情况作出相应的变化.例如对于冲击层填充的河谷(图 3)，计算公式相应地转换为
$$\left( {\mathit{\boldsymbol{S}}_{bb}^f - \mathit{\boldsymbol{S}}_{bb}^e + \mathit{\boldsymbol{S}}_{bb}^ * } \right)\mathit{\boldsymbol{u}}_b^g = \mathit{\boldsymbol{S}}_{bb}^f\mathit{\boldsymbol{u}}_b^f$$ (9) 式中， ${\pmb S}_{bb}^\ast $ 代表填充河谷部分的动力刚度矩阵； ${\pmb S}_{bb}^e $ 代表河谷采用与周围岩土介质相同的材料构成时的动力刚度矩阵.
在推导式(8) 时，矩阵 ${\pmb S}_{bb}^g $ 和位移矢量 ${\pmb u}_b^g $ 或 ${\pmb S}_{bb}^f $ 和 ${\pmb u}_b^f $ 的建立均基于散射体边界 $abc$ 或 $abcd$ 的结点自由度(参见图 4(a)和图 4(b))而形成，这只是一种选择.实际上，散射场边界可以有多种选择.例如，可选为稍远离散射体的边界 $ABCD$ (参见图 4).后者的边界比较规则，计算既相对简便，也有利于提高计算结果的精度[42].这时再利用式(4) 和式(6) 进行子结构 $ABCD$ 的动力响应分析，就可以容易地求得散射体边界 $abc$ 或 $abcd$ 的散射位移 ${\pmb u}_b^g $ ，也可进行有衬砌的地下结构(例如隧洞)的地震响应计算.这样，图 4中所表示的两类不同问题的散射都可根据 $ABCD$ 边界的 ${\pmb S}_{bb}^g $ 和 ${\pmb u}_b^g $ 而求得.
可以指出，这里所建立的计算公式，既可适用于二维结构，也可适用于三维结构；既可适用于均质半无限空间中的地下结构的地震响应分析，也可适用于复杂地基中地下结构的地震响应分析.由于均质半无限空间中动力刚度矩阵的建立和自由场地震动的求解比较方便，因此下文中数值验算将主要针对二维均质半空间的算例进行.复杂地基中动力刚度矩阵的建立需要求解地基内部的格林函数.下文将进一步介绍求解的基本思想和主要的计算公式.
3. 数值验算
本文将以均质半无限空间中的散射问题为主，通过算例检验计算模型的有效性.为计算均质半无限空间的动力刚度, 比例边界有限元方法(SBFEM)[45]是有效的方法. SBFEM为半解析半数值性的求解方法，其径向可获得解析解，环向则达到有限元的计算精度.由于只需要进行边界面的离散，问题的维数降低一阶；并可自动满足无限远处的辐射条件.因此SBFEM兼有FEM和BEM的优点，也不需要基本解，计算简便.
3.1 半圆形河谷的散射
如图 5所示, 河谷边界 $ABC$ , 计算边界 $DEFG$ ，相应计算域大小为 $4a×2a$ ( $a$ 为河谷半径).设围岩介质材料特性参数为：质量密度 $\rho _0 =1.5$ ，剪切波速 $c_{\rm s0} =1$ ，泊松比 $\nu _0 =1/3$ ；河谷中填充材料参数为：质量密度 $\rho _1 =1$ ，剪切波速 $c_{\rm s1} =0.5$ ，泊松比 $\nu _0 =1/3$ .无量纲频率 $\eta = (\omega a)/({\rm{ \mathsf{ π} }}{c_{\rm{s}}})$ ，波长 $\lambda = 2 / \eta $ .计算采用式(9).网格离散如图 5.各子域边界采用了3节点二次单元离散，共144单元，288节点.相应 $DEFG$ 边界设48单元，97节点.入射波从底部垂直输入，计算结果精度高，与Luco和de Barros的解[46]相符性良好(图 6). Luco和de Barros，以及Dravinski和Mossessian[47]采用的都是间接边界积分方程(IBIE)方法.
3.2 地下洞室衬砌的地震响应
研究地下隧道衬砌的地震响应. de Barros和Luco[48]求解了这一问题，采用间接边界积分方程方法模拟外围半空间土介质的作用，再结合Donnell的壳体理论计算隧洞衬砌外半径的波动响应.子域括号内的数字代表该子域相似中心的坐标.隧道结构如图 7(a)所示.隧道埋深 $H = 8.33r_1 = 7.573a$ .衬砌材料质量密度 $\rho _0 = 2.24×10^3$ kg/m $^3$ ，弹性模量 $E_0 = 1.6×10^{10}$ N/m $^2$ ，泊松比 $\nu _0 =0.2$ .围岩地基介质材料质量密度 $\rho _1 = 7.665×10^3$ kg/m $^3$ ，弹性模量 $E_1 = 6.9×10^8$ N/m $^2$ ，泊松比 $\nu _0 =0.45$ .位移按 $r =a$ 进行规格化.假设SV波垂直入射，无量纲频率 $\eta = (\omega a)/({\rm{ \mathsf{ π} }}{c_{\rm{s}}})$ =0.132.计算的衬砌响应如图 8所示.可见与文献解的相符性较好.
![]()
图 8 SV波垂直入射作用下隧道衬砌外表面的动力响应Figure 8. Dynamic response at the outer radius of the tunnel liner subjected to vertically incident SV wave4. 复杂地基条件下地下结构的地震响应分析
复杂地基条件下地下结构抗震分析的关键在于地基动力刚度 ${\pmb S}_{bb}^f $ 的求解.这是一个相对复杂的技术课题.在大坝、核电厂、桥梁和房屋的结构与地基动力相互作用分析中，常将无限地基简化为均质地基进行处理.但实际的天然地基则常常十分复杂.为此，在大坝与核电厂的抗震分析中，对图 9所示的3种典型的不均质地基，我们分别提出了地基动力刚度的计算模型.
其中的水平分层地基比较常见，下文将重点研究.对于分块不均质地基，当各块之间可以用同一相似中心加以描绘时(图 9(b))，提出采用比例边界有限元的处理方法[49]，对近场含不规则不均质捕虏体等复杂地基(图 9(c))则提出采用分步阻尼影响抽取法[50]进行求解.所提出的方法虽然是结合结构与地基动力相互作用分析进行阐述，但可适用于地下结构的抗震分析.
对各向同性和非各向同性层状半空间不均质地基的动力刚度及格林函数的求解，作者研究组从2012年起提出了基于积分变换的方法，并逐步得到完善[43, 51-53]，并可适用于横观各向同性的层状介质地基的求解.
下面介绍其基本思想和主要方程.
4.1 层状介质的波动方程
在圆柱坐标系下，层状介质(图 10)中任一层以位移所表示的波动方程具有如下形式
$$\left. {\begin{array}{*{20}{l}} {\left( {\lambda + 2G} \right)\frac{{\partial {\theta _t}}}{{\partial r}} - \frac{{2G}}{r}\frac{{\partial {\omega _z}}}{{\partial \theta }} + 2G\frac{{\partial {\omega _\theta }}}{{\partial z}}}\\ {\quad \quad \rho \frac{{{\partial ^2}{{\tilde u}_r}}}{{\partial {t^2}}}}\\ {\left( {\lambda + 2G} \right)\frac{{\partial {\theta _t}}}{{r\partial \theta }} + 2G\frac{{\partial {\omega _z}}}{{\partial r}} - 2G\frac{{\partial {\omega _r}}}{{\partial z}} = }\\ {\quad \quad \rho \frac{{{\partial ^2}{{\tilde u}_\theta }}}{{\partial {t^2}}}}\\ {\left( {\lambda + 2G} \right)\frac{{\partial {\theta _t}}}{{\partial z}} - \frac{{2G}}{r}\frac{\partial }{{\partial r}}\left( {r{\omega _\theta }} \right) + \frac{{2G}}{r}\frac{{\partial {\omega _r}}}{{\partial \theta }} = }\\ {\quad \quad \rho \frac{{{\partial ^2}{{\tilde u}_z}}}{{\partial {t^2}}}} \end{array}} \right\}$$ (10) 式中
$$\left. {\begin{array}{*{20}{l}} {{\theta _t} = {{\tilde \varepsilon }_r} + {{\tilde \varepsilon }_\theta } + {{\tilde \varepsilon }_z} = }\\ {\quad \quad \frac{1}{r}\frac{{\partial \left( {r{{\tilde u}_r}} \right)}}{{\partial r}} + \frac{1}{r}\frac{{\partial {{\tilde u}_\theta }}}{{\partial \theta }} + \frac{{\partial {{\tilde u}_z}}}{{\partial z}}{\rm{ }}}\\ {{\omega _r} = \frac{1}{2}\left( {\frac{1}{r}\frac{{\partial {{\tilde u}_z}}}{{\partial \theta }} - \frac{{\partial {{\tilde u}_\theta }}}{{\partial z}}} \right)}\\ {{\omega _\theta } = \frac{1}{2}\left( {\frac{{\partial {{\tilde u}_r}}}{{\partial z}} - \frac{{\partial {{\tilde u}_z}}}{{\partial r}}} \right)}\\ {{\omega _z} = \frac{1}{2}\left( {\frac{1}{r}\frac{{\partial \left( {r{{\tilde u}_\theta }} \right)}}{{\partial r}} - \frac{1}{r}\frac{{\partial {{\tilde u}_r}}}{{\partial \theta }}} \right)} \end{array}} \right\}$$ (11) 其中， $\lambda$ , µ为拉梅常数； $\tilde {u}_r $ , $\tilde {u}_\theta $ , $\tilde {u}_z$ 为 $r$ , θ, $ z$ 方向的位移.
对式(10) 进行Fourier-Bessel变换，得出波数域的表达式.变换式如下
$$\begin{array}{l} {\mathit{\boldsymbol{u}}_n}\left( {k,z,\omega } \right) = \\ \;\quad \quad {\mathit{\boldsymbol{a}}_n}{\rm{ }}\int_0^{ + \infty } {r{\mathit{\boldsymbol{C}}_n}} {\rm{ }}\int_0^{2{\rm{ \mathsf{ π} }}} {{\mathit{\boldsymbol{T}}_n}} \int_{ - \infty }^{ + \infty } {\mathit{\boldsymbol{\widetilde u}}\left( {r,\theta ,z,t} \right){{\rm{e}}^{ - {\rm{i}}\omega t}}{\rm{d}}t{\rm{d}}\theta {\rm{d}}r} {\rm{ }} \end{array}$$ (12) 相应的反变换式为
$$\begin{array}{l} \mathit{\boldsymbol{\tilde u}}\left( {r,\theta ,z,t} \right) = \\ \quad \quad \;\frac{1}{{2{\rm{ \mathsf{ π} }}}}\int_{ - \infty }^{ + \infty } {{{\rm{e}}^{{\rm{i}}\omega t}}} \sum\limits_{n = 0}^\infty {{\mathit{\boldsymbol{T}}_n}} {\rm{ }}\int_0^{ + \infty } {k{\mathit{\boldsymbol{C}}_n}{\mathit{\boldsymbol{u}}_n}\left( {k,z,\omega } \right){\rm{d}}k{\rm{d}}\omega } {\rm{ }} \end{array}$$ (13) 式中 $k$ 代表波数，其余系数的表达式为
$$\left. {\begin{array}{*{20}{l}} {{a_n} = \frac{1}{2}\;\;(n = 0)}\\ {{a_n} = \frac{1}{{\rm{\pi }}}\;\;(n \ne 0)} \end{array}} \right\}$$ (14) $${\mathit{\boldsymbol{T}}_n} = {\rm{diag}}\left[ {\begin{array}{*{20}{l}} {\left( {\begin{array}{*{20}{l}} {\cos n\theta }\\ {\sin n\theta } \end{array}} \right)} & {\left( {\begin{array}{*{20}{l}} {\cos n\theta }\\ {\cos n\theta } \end{array}} \right)} & {\left( {\begin{array}{*{20}{l}} {\cos n\theta }\\ {\sin n\theta } \end{array}} \right)} \end{array}} \right]$$ (15) $${\mathit{\boldsymbol{C}}_n}(kr) = \left[ {\begin{array}{*{20}{l}} {\frac{1}{k}{{\rm{J}}_n}{{(kr)}_{,r}}} & {{\rm{ }}\frac{n}{{kr}}{{\rm{J}}_n}(kr){\rm{ }}} & {\quad 0}\\ {\frac{n}{{kr}}{{\rm{J}}_n}(kr)} & {\frac{1}{k}{{\rm{J}}_n}{{(kr)}_{,r}}} & {\quad 0}\\ {\quad 0} & {\quad 0} & { - {{\rm{J}}_n}(kr)} \end{array}} \right]$$ (16) 式(15) 中的两个子式，上面的子式对应于荷载与变形相对于 $X$ 轴对称的情况；下面的子式则为荷载与变形相对于 $X$ 轴成反对称的情况.其中 ${\rm J}_n(kr)$ 为柱贝塞尔函数. Green函数计算时，针对圆形单元面上作用的均布水平向荷载和竖向荷载，可选择 $n = 1$ 和 $n =0$ 两种情况(图 11)与之对应.当圆形单元半径趋于极限 $\Delta r \to 0$ 时，均布荷载转化为集中水平力和集中竖向力作用的情况.即，在实际应用时，主要可限于 $n =1$ 和 $n =0$ 的两种情况.
按式(12) 进行Fourier-Bessel变换后得出频率-波数域内解耦的SV-P波动方程(变量为 $u_r.u_z $ )和SH波动方程(变量为 $u_\theta $ ) [51-52]，也可参见文献[28].
$$\begin{array}{l} \left[ {\begin{array}{*{20}{l}} \mu & {\quad 0}\\ {\,0} & {\lambda + \mu } \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{{\mu ''}_r}}\\ {{{\mu ''}_z}} \end{array}} \right] + k\left[ {\begin{array}{*{20}{l}} {\quad 0} & { - \lambda + \mu }\\ {\lambda + \mu } & {\quad 0} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{{\mu '}_r}}\\ {{{\mu '}_2}} \end{array}} \right] - \\ \quad \left( {{k^2}\left[ {\begin{array}{*{20}{l}} {\lambda + 2\mu } & {\,0}\\ {\quad 0} & \mu \end{array}} \right] - \rho {\omega ^2}} \right)\left[ {\begin{array}{*{20}{l}} {{\mu _r}}\\ {{\mu _2}} \end{array}} \right] = {\bf{0}} \end{array}$$ (17) $$\mu {{u''}_\theta } - \left( {{k^2}\mu - \rho {\omega ^2}} \right){u_\theta } = 0$$ (18) 对横观各向同性介质中的波动，解耦方程同样成立，但相应系数改变如下
$$\begin{array}{l} \left[ {\begin{array}{*{20}{l}} {{d_{44}}} & {\;0}\\ {\;0} & {{d_{33}}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{{\mu ''}_r}}\\ {{{\mu ''}_z}} \end{array}} \right] + k\left[ {\begin{array}{*{20}{l}} {\quad 0} & { - \left( {{d_{13}} + {d_{44}}} \right)}\\ {{d_{13}} + {d_{44}}} & {\quad \quad 0} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{{u'}_r}}\\ {{{u'}_z}} \end{array}} \right]\\ \quad \left( {{k^2}\left[ {\begin{array}{*{20}{l}} {{d_{11}}} & 0\\ 0 & {{d_{44}}} \end{array}} \right] - \rho {\omega ^2}} \right)\left[ {\begin{array}{*{20}{l}} {{u_r}}\\ {{u_z}} \end{array}} \right] = {\bf{0}}\\ \end{array}$$ (19) $${d_{44}}{{u''}_\theta } - \left( {{k^2}{d_{66}} - \rho {\omega ^2}} \right){u_\theta } = 0$$ (20) 式中， ${u}'_r, {u}”_r $ 分别表示 ${{\rm{d}} u_r }/ {{\rm{d}} z}$ ， ${{\rm{d}}^2u_r }/ {{\rm{d}} z^2}$ .方程(17)、(18) 或方程(19)、(20) 均可表示为以下统一的形式
$$\begin{array}{l} \mathit{\boldsymbol{K}}_{22}^m{\left( {{u^m}} \right)^{\prime \prime }} + \left( {\mathit{\boldsymbol{K}}_{21}^m - \mathit{\boldsymbol{K}}_{12}^m} \right){\left( {{\mathit{\boldsymbol{u}}^m}} \right)^\prime } - {\rm{ }}\\ \quad \left( {\mathit{\boldsymbol{K}}_{11}^m - \rho {\omega ^2}{\mathit{\boldsymbol{I}}^m}} \right){\mathit{\boldsymbol{u}}^m} = {\bf{0}}\qquad (m = 1,2){\rm{ }} \end{array}$$ (21) 式中, 上标 $m =1$ 代表SV-P波动， $m =2$ 代表SH波动. ${\pmb I }^m$ 为单位阵， $m =1$ 为2阶单位阵， $m =2$ 为1阶单位阵.
$${\mathit{\boldsymbol{u}}^1} = {\left[ {{u_r}\;\;{u_z}} \right]^{\rm{T}}},\;\;\;{\mathit{\boldsymbol{u}}^2} = {u_\theta }$$ (22) 对各向同性介质(方程(17) 与方程(18))
$$\left. {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{K}}_{11}^1 = {k^2}\left[ {\begin{array}{*{20}{l}} {\lambda + 2\mu } & 0\\ {\quad 0} & \mu \end{array}} \right],\mathit{\boldsymbol{K}}_{22}^1 = \left[ {\begin{array}{*{20}{l}} \mu & {\quad 0}\\ 0 & {\lambda + 2\mu } \end{array}} \right]}\\ {\mathit{\boldsymbol{K}}_{21}^1 = k\left[ {\begin{array}{*{20}{l}} {\quad 0} & 0\\ {\lambda + \mu } & 0 \end{array}} \right],\mathit{\boldsymbol{K}}_{12}^1 = \left[ {\begin{array}{*{20}{l}} 0 & {\lambda + \mu }\\ 0 & {\quad 0} \end{array}} \right]}\\ {\mathit{\boldsymbol{K}}_{11}^2 = {k^2}\mu ,\;\;\mathit{\boldsymbol{K}}_{22}^2 = \mu ,\;\;\mathit{\boldsymbol{K}}_{21}^2 = \mathit{\boldsymbol{K}}_{12}^2 = {\bf{0}}} \end{array}} \right\}$$ (23) 对于横观各向同性介质
$$\left. {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{K}}_{11}^1 = {k^2}\left[ {\begin{array}{*{20}{l}} {{d_{11}}} & {\;0}\\ {\;0} & {{d_{44}}} \end{array}} \right],\mathit{\boldsymbol{K}}_{22}^1 = \left[ {\begin{array}{*{20}{l}} {{d_{44}}} & {\;0}\\ {\;0} & {{d_{33}}} \end{array}} \right]}\\ {\mathit{\boldsymbol{K}}_{12}^1 = k\left[ {\begin{array}{*{20}{l}} {\;0} & {\;{d_{13}} + {d_{44}}}\\ {\;0} & {\;\quad 0} \end{array}} \right],\mathit{\boldsymbol{K}}_{21}^1 = k\left[ {\begin{array}{*{20}{l}} {\;0} & {\;0}\\ {\;{d_{13}} + {d_{44}}} & {\;0} \end{array}} \right]}\\ {\mathit{\boldsymbol{K}}_{11}^2 = {k^2}{d_{66}},\;\;\mathit{\boldsymbol{K}}_{22}^2 = {d_{44}},\;\;\mathit{\boldsymbol{K}}_{21}^2 = \mathit{\boldsymbol{K}}_{12}^2 = {\bf{0}}{\rm{ }}} \end{array}} \right\}$$ (24) 为了进行微分方程(21) 的求解，引入如下定义的应力矢量 ${\pmb p}^m$ 作为位移 ${\pmb u}^m$ 的对偶变量
$${\mathit{\boldsymbol{p}}^1} = - {\left[ {{\tau _{rz}}\;\;{\sigma _z}} \right]^{\rm{T}}},\;{\mathit{\boldsymbol{p}}^2} = - {\tau _{\theta z}}$$ (25) 可以证明 ${\pmb p}^m$ 满足关系式[54]
$${\mathit{\boldsymbol{p}}^m} = - ({\mathit{\boldsymbol{K}}_{22}}({\mathit{\boldsymbol{u}}^m})' + \mathit{\boldsymbol{K}}_{21}^m{\mathit{\boldsymbol{u}}^m})\;\;\;(m = 1,2)$$ (26) 以下的求解过程对 $m=1$ 和 $m= 2$ 都适用，所以在下文中可将 $m$ 省略.
利用式(26)，微分方程(21) 可化为如下对耦形式
$${\mathit{\boldsymbol{u}}^\prime } = \mathit{\boldsymbol{Au}} + \mathit{\boldsymbol{Dp}},\quad {\mathit{\boldsymbol{p}}^\prime } = - \mathit{\boldsymbol{Bu}} - \mathit{\boldsymbol{Cp}}$$ (27) 矩阵 ${\pmb A}, {\pmb B}, {\pmb C}, {\pmb D}$ 和矩阵 ${\pmb K}_{ij} \ (i, j = 1, 2)$ 的变换关系如下
$$\left. {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{A}} = - \mathit{\boldsymbol{K}}_{22}^{ - 1}{\mathit{\boldsymbol{K}}_{21}},\;\;\mathit{\boldsymbol{B}} = {\mathit{\boldsymbol{K}}_{11}} - {\mathit{\boldsymbol{K}}_{12}}\mathit{\boldsymbol{K}}_{22}^{ - 1}{\mathit{\boldsymbol{K}}_{21}} - \rho {\omega ^2}\mathit{\boldsymbol{I}}}\\ {\mathit{\boldsymbol{C}} = {\mathit{\boldsymbol{K}}_{12}}\mathit{\boldsymbol{K}}_{22}^{ - 1} = - {\mathit{\boldsymbol{A}}^{\rm{T}}},\;\;\mathit{\boldsymbol{D}} = - \mathit{\boldsymbol{K}}_{22}^{ - 1}}\\ {\mathit{\boldsymbol{B}} = {\mathit{\boldsymbol{B}}^{\rm{T}}},\;\;\mathit{\boldsymbol{D}} = {\mathit{\boldsymbol{D}}^{\rm{T}}}{\rm{ }}} \end{array}} \right\}$$ (28) 在状态空间，对耦微分方程(27) 可以联立求解
$$\mathit{\boldsymbol{V'}}{\rm{ = }}\mathit{\boldsymbol{HV}}$$ (29) 式中
$$\mathit{\boldsymbol{V}} = \left\{ {\begin{array}{*{20}{l}} \mathit{\boldsymbol{u}}\\ \mathit{\boldsymbol{p}} \end{array}} \right\},\mathit{\boldsymbol{H}} = \left[ {\begin{array}{*{20}{l}} \mathit{\boldsymbol{A}} & {\;\mathit{\boldsymbol{D}}}\\ { - \mathit{\boldsymbol{B}}} & { - {\mathit{\boldsymbol{A}}^{\rm{T}}}} \end{array}} \right]$$ (30) 齐次一阶线性常微分方程(29) 的解为指数函数.
$$\mathit{\boldsymbol{V}} = \exp \left( {\mathit{\boldsymbol{Hz}}} \right)\mathit{\boldsymbol{c}}$$ (31) 采用钟万勰提出的精细积分方法[55]可以将解进行精确计算.
对层状半空间中的任意一层，或其中的一层，设层厚为 $\eta = z_b -z_a $ ，则层上、下两端的位移和应力将成立如下关系
$${\mathit{\boldsymbol{V}}_b} = \exp \left( {\mathit{\boldsymbol{H}}\eta } \right){\mathit{\boldsymbol{V}}_a}$$ (32) 或写为
$${\mathit{\boldsymbol{V}}_b} = \mathit{\boldsymbol{T}}{\mathit{\boldsymbol{V}}_a}$$ $$\left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}_b}}\\ {{\mathit{\boldsymbol{p}}_b}} \end{array}} \right\} = \left[ {\begin{array}{*{20}{l}} {{T_{uu}}} & {{T_{up}}}\\ {{T_{pu}}} & {{T_{pp}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}_a}}\\ {{\mathit{\boldsymbol{p}}_a}} \end{array}} \right\}$$ (33) 式中
$$\begin{array}{l} \mathit{\boldsymbol{T}} = \exp \left( {\mathit{\boldsymbol{H}}\eta } \right) \approx \mathit{\boldsymbol{I}} + \mathit{\boldsymbol{H}}\eta + \frac{1}{{2!}}{\left( {\mathit{\boldsymbol{H}}\eta } \right)^2} + \\ \quad \quad \frac{1}{{3!}}{\left( {\mathit{\boldsymbol{H}}\eta } \right)^3} + \frac{1}{{4!}}{\left( {\mathit{\boldsymbol{H}}\eta } \right)^4} + \cdots \end{array}$$ (34) 为达到高精度的计算效果， ${\pmb T}$ 按如下方式进行计算[55]
$$\mathit{\boldsymbol{T}} = \exp \left( {\mathit{\boldsymbol{H}}\eta } \right) = {\left[ {\exp \left( {\mathit{\boldsymbol{H}}\eta /b} \right)} \right]^b} = {\left[ {\exp \left( {\mathit{\boldsymbol{H}}\tau } \right)} \right]^b}$$ (35) 式中， $\tau = \eta/ b$ ， $b$ 可取为任意整数.当取 $b = 2^N, N = 20$ 时，这相当于将 $\eta $ 分成 $2^N = 1 048 557$ 个微细薄层用式(34) 计算，从而可以达到任意希望的精度，而计算工作量不大.文献[55]中给出了计算 ${\pmb T}$ 的递推公式，应用方便.
为了进行层的结合和计算Green函数，解也宜表示成如下对耦变量的形式
$$\left. {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}_b} = \mathit{\boldsymbol{F}}{\mathit{\boldsymbol{u}}_a} + \mathit{\boldsymbol{G}}{\mathit{\boldsymbol{p}}_b}}\\ {{\mathit{\boldsymbol{p}}_a} = - \mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{u}}_a} + {\mathit{\boldsymbol{F}}^{\rm{T}}}{\mathit{\boldsymbol{p}}_b}} \end{array}} \right\}$$ (36) 式中矩阵 ${\pmb F}, {\pmb G}, {\pmb Q}$ 也由 ${\pmb T}$ 求出
$$\left. {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{Q}} = \mathit{\boldsymbol{T}}_{pp}^{ - 1}{\mathit{\boldsymbol{T}}_{pu}},\;\;\mathit{\boldsymbol{G}} = {\mathit{\boldsymbol{T}}_{up}}\mathit{\boldsymbol{T}}_{pp}^{ - 1}}\\ {{\mathit{\boldsymbol{F}}^{\rm{T}}} = \mathit{\boldsymbol{T}}_{pp}^{ - 1}} \end{array}} \right\}$$ (37) 对于水平分层介质，采用积分变换方法，可以求得波动方程在波数域的简便形式(式(17) 和式(18) 或式(19) 和式(20)).进一步引入对偶变量，又可将其化为一阶齐次线性常微分方程(29)，其解可采用精细积分方法达到任意希望的精度.最后将解表示成对耦形式(36)，可便于分层之间的结合和Green函数的求解，说明如下.
4.2 相邻层的结合
当将相邻的两层：层1 $\left[{z_a, z_b } \right]$ 和层2 $\left[{ z_b, z_c } \right]$ 进行结合时，应用式(36)，可有
$$\left. {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}_b} = {\mathit{\boldsymbol{F}}^1}{\mathit{\boldsymbol{u}}_a} + {\mathit{\boldsymbol{G}}^1}{\mathit{\boldsymbol{p}}_b}}\\ {{\mathit{\boldsymbol{p}}_a} = - {\mathit{\boldsymbol{Q}}^1}{\mathit{\boldsymbol{u}}_a} + {{\left( {{\mathit{\boldsymbol{F}}^1}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{p}}_b}}\\ {{\mathit{\boldsymbol{u}}_c} = {\mathit{\boldsymbol{F}}^2}{\mathit{\boldsymbol{u}}_b} + {\mathit{\boldsymbol{G}}^2}{\mathit{\boldsymbol{p}}_c}}\\ {{\mathit{\boldsymbol{p}}_b} = - {\mathit{\boldsymbol{Q}}^2}{\mathit{\boldsymbol{u}}_b} + {{\left( {{\mathit{\boldsymbol{F}}^2}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{p}}_c}} \end{array}} \right\}$$ (38) 从式(38) 中消去 ${\pmb u}_b $ 和 ${\pmb p}_b $ ，即可求得合并后新层 $c \left[{ z_a, z_c } \right]$ 的力与位移关系
$$\left. {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}_c} = {\mathit{\boldsymbol{F}}^c}{\mathit{\boldsymbol{u}}_a} + {\mathit{\boldsymbol{G}}^c}{\mathit{\boldsymbol{p}}_c}}\\ {{\mathit{\boldsymbol{p}}_a} = - {\mathit{\boldsymbol{Q}}^c}{\mathit{\boldsymbol{u}}_a} + {{\left( {{\mathit{\boldsymbol{F}}_c}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{p}}_c}} \end{array}} \right\}$$ (39) 式中
$$\left. {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{F}}^c} = {\mathit{\boldsymbol{F}}^2}{{({\rm{I}} + {\mathit{\boldsymbol{G}}^1}{\mathit{\boldsymbol{Q}}^2})}^{ - 1}}{\mathit{\boldsymbol{F}}^1}}\\ {{\mathit{\boldsymbol{G}}^c} = {\mathit{\boldsymbol{G}}^2} + {\mathit{\boldsymbol{F}}^2}{{({{({\mathit{\boldsymbol{G}}^1})}^{ - 1}} + {\mathit{\boldsymbol{Q}}^2})}^{ - 1}}{{({\mathit{\boldsymbol{F}}^2})}^H}}\\ {{\mathit{\boldsymbol{Q}}^c} = {\mathit{\boldsymbol{Q}}^1} + {{({\mathit{\boldsymbol{F}}^1})}^H}{{({{({\mathit{\boldsymbol{Q}}^2})}^{ - 1}} + {\mathit{\boldsymbol{G}}^1})}^{ - 1}}{\mathit{\boldsymbol{F}}^1}} \end{array}} \right\}$$ (40) 4.3 Green函数的计算
层状半空间内部点Green函数的计算是地下结构抗震分析的核心部分.文献中广泛采用的刚度矩阵方法(stiffness matrix method)需要进行大型矩阵求解[28].这里提出对偶变量的转换方法，计算简捷方便，同时又准确高效.计算图形如图 12和图 13所示.图中源点(source)表示力的作用点，接受点(receiver)表示位移的计算点.已知源点和接受点的位置，可将层状半空间划分为3个子集，应用4.2节中阐述的方法可以分别建立各个子集形如式(36) 所示的对耦变量方程，并相应求出其 ${\pmb F }_i, {\pmb Q }_i, {\pmb G }_i $ 矩阵( $i=u, m, l$ ，其中 $u$ 表示上部， $m$ 表示中部， $l$ 表示下部)方程. Green函数的计算，基本上可划分为图 13所示的6种工况.应用式(36)，即可求得各种工况的Green函数计算的矩阵表达式[53].
![]()
图 13 层状半空间不同工况的Green函数计算Figure 13. Computation of Green's function for various cases of boundary conditions of multilayered half-space下面给出不同工况的Green函数的矩阵表达式.下标 $a, b$ 表示源点和接受点所在位置.
(1) 工况(a)和(b)层状地基表面点的Green函数
工况(a)下部为刚性地基
$${\mathit{\boldsymbol{u}}_a} = - {({\mathit{\boldsymbol{Q}}_l} + \mathit{\boldsymbol{F}}_l^{\rm{T}}\mathit{\boldsymbol{G}}_l^{ - 1}{\mathit{\boldsymbol{F}}_l})^{ - 1}}{\mathit{\boldsymbol{p}}_a}$$ (41) 工况(b)下部为弹性半空间
$${\mathit{\boldsymbol{u}}_a} = {[ - {\mathit{\boldsymbol{Q}}_l} + \mathit{\boldsymbol{F}}_l^{\rm{T}}{({\mathit{\boldsymbol{R}}_\infty } - {\mathit{\boldsymbol{G}}_l})^{ - 1}}{\mathit{\boldsymbol{F}}_l}]^{ - 1}}{\mathit{\boldsymbol{p}}_a}$$ (42) (2) 工况(c)和(d)下部固定的层状地基内部点的格林函数
工况(c)
$${\mathit{\boldsymbol{u}}_b} = \mathit{\boldsymbol{B}}_1^{ - 1}{\mathit{\boldsymbol{A}}_1}{\mathit{\boldsymbol{p}}_a}$$ (43) 式中
$$\left. {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{A}}_1} = {\mathit{\boldsymbol{F}}_m}({\mathit{\boldsymbol{G}}_u} + {\mathit{\boldsymbol{F}}_u}\mathit{\boldsymbol{Q}}_u^{ - 1}\mathit{\boldsymbol{F}}_u^{\rm{T}})}\\ {{\mathit{\boldsymbol{B}}_1} = \mathit{\boldsymbol{I}} + {\mathit{\boldsymbol{G}}_m}({\mathit{\boldsymbol{Q}}_l} + \mathit{\boldsymbol{F}}_l^{\rm{T}}\mathit{\boldsymbol{G}}_l^{ - 1}{\mathit{\boldsymbol{F}}_l})} \end{array}} \right\}$$ (44) 工况(d)
$${\mathit{\boldsymbol{u}}_a} = \mathit{\boldsymbol{B}}_2^{ - 1}{\mathit{\boldsymbol{A}}_2}{\mathit{\boldsymbol{p}}_b}$$ (45) 式中
$$\left. {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{B}}_2^{ - 1} = ({\mathit{\boldsymbol{G}}_u} + {\mathit{\boldsymbol{F}}_u}\mathit{\boldsymbol{Q}}_u^{ - 1}\mathit{\boldsymbol{F}}_u^{\rm{T}})}\\ {{\mathit{\boldsymbol{A}}_2} = {\mathit{\boldsymbol{Q}}_m}\mathit{\boldsymbol{F}}_m^{ - 1}{{({\mathit{\boldsymbol{Q}}_l} + \mathit{\boldsymbol{F}}_l^{\rm{T}}\mathit{\boldsymbol{G}}_m^{ - 1}{\mathit{\boldsymbol{F}}_l})}^{ - 1}} + }\\ {\quad {\mathit{\boldsymbol{Q}}_m}\mathit{\boldsymbol{F}}_m^{ - 1}{\mathit{\boldsymbol{G}}_m} + \mathit{\boldsymbol{F}}_m^{\rm{T}}} \end{array}} \right\}$$ (46) (3) 工况(e)和(f)下部为弹性半空间的层状地基内部点的格林函数
工况(e)
$${\mathit{\boldsymbol{u}}_b} = \mathit{\boldsymbol{B}}_3^{ - 1}{\mathit{\boldsymbol{A}}_3}{\mathit{\boldsymbol{p}}_a}$$ (47) 式中
$$\left. {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{A}}_3} = {\mathit{\boldsymbol{A}}_1}}\\ {{\mathit{\boldsymbol{B}}_3} = \mathit{\boldsymbol{I}} + {\mathit{\boldsymbol{G}}_m}\left[ {{\mathit{\boldsymbol{Q}}_l} - \mathit{\boldsymbol{F}}_l^T{{({\mathit{\boldsymbol{R}}_\infty } - {\mathit{\boldsymbol{G}}_l})}^{ - 1}}{\mathit{\boldsymbol{F}}_l}} \right]} \end{array}} \right\}$$ (48) 工况(f)
$${\mathit{\boldsymbol{u}}_a} = \mathit{\boldsymbol{B}}_4^{ - 1}{\mathit{\boldsymbol{A}}_4}{\mathit{\boldsymbol{p}}_b}$$ (49) 式中
$$\left. {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{B}}_4} = {\mathit{\boldsymbol{B}}_2}}\\ {{\mathit{\boldsymbol{A}}_4} = {\mathit{\boldsymbol{Q}}_m}\mathit{\boldsymbol{F}}_m^{ - 1}{{\left[ {{\mathit{\boldsymbol{Q}}_l} - \mathit{\boldsymbol{F}}_l^{\rm{T}}({\mathit{\boldsymbol{R}}_\infty } - {\mathit{\boldsymbol{G}}_l}){\mathit{\boldsymbol{F}}_l}} \right]}^{ - 1}} + }\\ {\quad {\mathit{\boldsymbol{Q}}_m}\mathit{\boldsymbol{F}}_m^{ - 1}{\mathit{\boldsymbol{G}}_m} + \mathit{\boldsymbol{F}}_m^{\rm{T}}} \end{array}} \right\}$$ (50) 所有矩阵的阶数对SV-P波为2 $\times$ 2，对SH波为1 $\times$ 1，计算很简便.最后可将波数域各种工况格林函数的计算，写成如下形式的统一的表达式(参见式(17) 和式(18) 或式(19) 和式(20), $u_r$ , $u_z$ 和 $u_B$ 是解耦的，可分别求解
$$\left. {\begin{array}{*{20}{l}} {\left\{ {\begin{array}{*{20}{l}} {{u_r}(k,\omega )}\\ {{u_z}(k,\omega )} \end{array}} \right\} = \left[ {\begin{array}{*{20}{l}} {{F_{rr}}} & {{F_{rz}}}\\ {{F_{zr}}} & {{F_{zz}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{l}} {{p_r}(k,\omega )}\\ {{p_z}(k,\omega )} \end{array}} \right\}}\\ {m = 1} \end{array}} \right\}$$ (51) $${u_\theta }(k,\omega ) = {F_{\theta \theta }}P\theta (k,\omega )$$ (52) 对图 11表示的格林函数计算的荷载形式，例如，半径为 $\Delta r$ 的圆形单元上分布的 $Z$ 向和 $X$ 向均布荷载 $\tilde {p}_z $ 和 $\tilde {p}_x $ ，其相应的波数域的表达式为
$$\begin{array}{l} {p_z}(k,\omega ) = \frac{1}{{2{\rm{ \mathsf{ π} }}}}\int_{r = 0}^{\Delta r} {r( - {{\rm{J}}_0}(kr))} {\rm{ }}\int_0^{2{\rm{ \mathsf{ π} }}} {{{\tilde p}_z}{\rm{d}}\theta {\rm{d}}r = } \\ \quad \quad - \frac{{{{\tilde p}_z}\Delta r}}{k}{{\rm{J}}_1}(k\Delta r) \end{array}$$ (53) $$\left\{ {\begin{array}{*{20}{l}} {{p_r}(k,\omega )}\\ {{p_\theta }(k,\omega )} \end{array}} \right\} = \frac{1}{\mathit{\boldsymbol{\pi }}}\int_{r = 0}^{\Delta r} {\mathit{\boldsymbol{F}}(kr)dr = \frac{{{{\tilde p}_x}\Delta r}}{k}\left\{ {\begin{array}{*{20}{l}} 1\\ 1 \end{array}} \right\}{{\rm{J}}_1}(k\Delta r)} $$ (54) 式中
$$\begin{array}{l} F(kr) = \frac{1}{k}\left[ {\begin{array}{*{20}{l}} {r{{\rm{J}}_1}{{(kr)}_{,r}}} & {{{\rm{J}}_1}(kr)}\\ {{{\rm{J}}_1}(kr)} & {r{{\rm{J}}_1}{{(kr)}_{,r}}} \end{array}} \right]\int_0^{2\mathit{\boldsymbol{\pi }}} {\left[ {\begin{array}{*{20}{l}} {\cos \theta \quad \quad }\\ {\quad \quad - \sin \theta } \end{array}} \right]} .\\ \quad \quad \left\{ {\begin{array}{*{20}{l}} {{{\tilde p}_x}\cos \theta }\\ { - {{\tilde p}_x}\sin \theta } \end{array}} \right\}{\rm{d}}\theta {\rm{ }} \end{array}$$ 在上述式中，令 $\Delta r \to 0$ ， ${\rm{\pi }} \Delta r^2 \tilde {p}_z \to \tilde {P}_z $ 或 ${\rm{\pi }} \Delta r^2 \tilde {p}_x \to \tilde {P}_x $ 就可得出集中力 $\tilde {P}_z$ 或 $\tilde {P}_x $ 作用的表达式.再将式(53)、式(54) 代入式(51)、式(52) 并利用(13) 进行反变换，就可得出圆柱坐标系物理域中格林函数的表达式，可统一写成
$$\left\{ {\begin{array}{*{20}{l}} {{u_r}(r,\theta ,\omega )}\\ {{u_\theta }(r,\theta ,\omega )}\\ {{u_z}(r,\theta ,\omega )} \end{array}} \right\} = \left[ {\begin{array}{*{20}{l}} {{G_{rx}}} & {{G_{ry}}} & {{G_{rz}}}\\ {{G_{\theta x}}} & {{G_{ry}}} & 0\\ {{G_{zx}}} & {{G_{ry}}} & {{G_{zz}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{l}} {{{\tilde p}_x}}\\ {{{\tilde p}_y}}\\ {{{\tilde p}_z}} \end{array}} \right\}$$ (55) 图 11中表示了 $\tilde p_z$ 和 $\tilde p_x$ 作用的求解， $\tilde p_y$ 的求解与 $\tilde p_x$ 类似.下面通过算例来检验其计算效果.
4.4 数值检验
研究均质半无限地基上表面弹性层中半圆形河谷的散射(图 14).散射计算需要求解层状半空间中的格林函数.计算中介质材料特性参数假设为：均质半空间 $\rho _0 = 4 / 3$ ， $r_0 = 0.333$ ， $c_{\rm s0} = 2$ ；表面弹性层 $\rho _1 = 1$ ， $r_1 = 0.333$ ， $c_{\rm s1} = 1$ ；河谷中冲积层 $\rho _2 = 0.667$ ， $r_2 = 0.333$ ， $c_{\rm s2} = 0.5$ .几何尺寸：河谷半径 $r$ ，表面弹性层厚 $h = 1.5r$ .无量纲频率 $\eta = (\omega a)/({\rm{ \mathsf{ π} }}{c_{\rm{s}}})$ .
![]()
图 14 半无限地基上表面弹性层中的半圆形河谷的散射Figure 14. Scattering of a semi-circular valley embedded in the surface layer overlying elastic half-space在SV波垂直入射情况下半圆形充填河谷地表位移幅值计算结果如图 15所示，可见与文献[56]的相符性很好.
![]()
图 15 SV波垂直入射情况下层状半空间中半圆形河谷地表位移幅值Figure 15. Amplitude of displacement response of the semi-circular valley embedded in layered half-space under vertically incident SV wave5. 地下结构地震响应分析的其他有关问题
5.1 地下结构地震响应的时域分析
地下结构地震响应的时域计算应用运动方程(1)，关键是相互作用力 ${\pmb R} = {\pmb S}_{bb}^g {\pmb u}_b^g $ 或 ${\pmb R } = {\pmb S }_{bb}^f {\pmb u}_b^f $ 的计算.由于无限域中 ${\pmb R }$ 是激励频率的函数，需将其转换成时域的表达式，如式(51) 所示(假设初值条件 ${\pmb u} (t = 0) = {\bf 0 }$ , $\dot{\pmb u} (t = 0) = {\bf 0 }$ )
$$\mathit{\boldsymbol{R}}(\omega ) = {\mathit{\boldsymbol{S}}^\infty }(\omega )\mathit{\boldsymbol{u}}(\omega )$$ (56) $$\mathit{\boldsymbol{R}}(t) = {\rm{ }}\int_0^t {{\mathit{\boldsymbol{S}}^\infty }(t - \tau )\mathit{\boldsymbol{u}}(\tau )dt} {\rm{ }}$$ (57) 式中， ${\pmb S}^\infty (\omega)$ 代表无限域的动刚度 ${\pmb S}_{bb}^g (\omega)$ 或 ${\pmb S}_{bb}^f (\omega)$ ； ${\pmb S}^\infty (t)$ 代表单位位移脉冲响应函数
$${\mathit{\boldsymbol{S}}^\infty }(t) = \frac{1}{{2{\rm{ \mathsf{ π} }}}}\int_{ - \infty }^\infty {{\mathit{\boldsymbol{S}}^\infty }(\omega ){{\rm{e}}^{{\rm{i}}\omega t}}{\rm{d}}\omega } $$ (58) 式(58) 只具有形式上的意义.由于 ${\pmb S}^\infty (\omega)$ 非平方可积，含有奇异性，式(58) 不能直接用于脉冲响应函数的计算[57]，需进行必要的处理，将另文阐述.
5.2 地下结构地震响应的简化计算
地下结构也可采用简化的数值方法进行初步的抗震分析.以图 16所示的地下结构为例来加以说明.图中虚线范围以内为计算域，域内可采用有限元等数值方法进行分析.根据第2节计算模型的推导，虚线所表示的边界应满足相互作用力的条件 ${\pmb R} = {\pmb S}_{bb}^g {\pmb u}_b^g $ 或 ${\pmb R} = {\pmb S }_{bb}^f {\pmb u}_b^f $ .也可选择较大的计算域，并设立能量传递边界，通过边界进行自由场的地震动输入.文献中一般的能量传递边界主要建立在均质无限地基假定的基础上，不适于复杂地基条件下地下结构的地震响应分析.如采用时间、空间非耦合的局部人工透射边界，计算准确度并不理想，需要采用高阶局部人工边界，但高阶局部边界计算精度的提高是以计算自由度的增加为代价的.我们提出的层状地基耦合性(或全域性)时域传递边界[58]，不仅计算精度高，而且计算效率也高，可供参考.
![]()
图 16 地下结构抗震分析的简化模型Figure 16. A simplified model for seismic analysis of under-ground structures6. 结束语
提出一种基于结构与地基动力相互作用方程的地下结构抗震分析模型，可以很方便地进行地震波散射场与地下结构和围岩动力相互作用的计算.模型具有广泛的适用性.可以考虑层状半空间等复杂地基的影响.计算准确、高效.本文以二维均匀半空间中地下结构的地震响应为例，阐述地下结构抗震计算模型的设想.但是对复杂层状半空间地下结构的地震响应只提出了计算方法和相关方程还来不及展开，以后将陆续进行阐述.
致谢: 博士生李志远、韩泽军和课题组李建波、胡志强等在本文完成中发挥了重要作用, 在此表示感谢. -
![]()
图 2 散射体的计算模型示意图
Figure 2. Computational model for surface and subsurface irregularities and inhomogeneities
![]()
图 8 SV波垂直入射作用下隧道衬砌外表面的动力响应
Figure 8. Dynamic response at the outer radius of the tunnel liner subjected to vertically incident SV wave
![]()
图 13 层状半空间不同工况的Green函数计算
Figure 13. Computation of Green's function for various cases of boundary conditions of multilayered half-space
![]()
图 14 半无限地基上表面弹性层中的半圆形河谷的散射
Figure 14. Scattering of a semi-circular valley embedded in the surface layer overlying elastic half-space
![]()
图 15 SV波垂直入射情况下层状半空间中半圆形河谷地表位移幅值
Figure 15. Amplitude of displacement response of the semi-circular valley embedded in layered half-space under vertically incident SV wave
-
[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 https://www.researchgate.net/publication/229781958_A_review_of_the_response_of_buried_pipelines_under_seismic_excitations
[2] Wang LRL, O'Rourke MJ. Overview of buried pipelines under seismic loading. Journal of Technical Councils of ASCE, 1978, 104(TC1):121-130 http://cedb.asce.org/CEDBsearch/record.jsp?dockey=8414
[3] Okamoto S, Tamura C. Behavior of subaqueous tunnel during earthquakes. International Journal of Earthquake Engineering & Structural Dynamics, 1973, 1(3):253-266 doi: 10.1002/eqe.4290010306/references
[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 doi: 10.1007/978-4-431-54892-8_5.pdf
[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 doi: 10.1002/eqe.4290010307/epdf
[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 doi: 10.1002/eqe.4290230507/full
[10] Wong HL. Effects 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 effects on strong ground motion. Soil Dynamics and Earthquake Engineering, 1987, 6(2):124-132 doi: 10.1016/0267-7261(87)90022-4
[12] Von Thun JL. Earthquake Engineering and Soil Dynamics Ⅱ-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 https://www.researchgate.net/publication/229600850_Response_of_a_circular_cylindrical_shell_to_disturbances_in_a_half-space
[14] Wong KC, Shah AH, Datta SK. Dynamic stresses and displacements in a buried tunnel. Journal of Engineering Mechanics, 1985, 111(2):218-234 doi: 10.1061/(ASCE)0733-9399(1985)111:2(218)
[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 doi: 10.1002/eqe.4290180309/citedby
[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 https://www.researchgate.net/publication/229758333_Transient_Scattering_of_Elastic_Waves_by_Dipping_Layers_of_Arbitrary_Shape_Part_1_Antiplane_Strain_Model
[17] Lee VW. Three-dimensional diffraction of plane P, SV & SH waves by a hemispherical alluvial valley. International Journal of Soil Dynamics and Earthquake Engineering, 1984, 3(3):133-144 doi: 10.1016/0261-7277(84)90043-3
[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 https://www.researchgate.net/publication/229876733_Scattering_of_plane_harmonic_SH_SV_P_and_rayleigh_waves_by_non-axisymmetric_three-dimensional_canyons_A_wave_function_expansion_approach
[19] Wong HL, Jennings PC. Effects of canyon topography on strong ground motion. Bulletin of the Seismological Society of America, 1975, 65(5):1239-1257 http://authors.library.caltech.edu/48928/1/1239.full.pdf
[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 doi: 10.1002/eqe.4290201003/abstract
[21] Wong HL. Effect of surface topography on the diffraction of P, SV, and Rayleigh waves. Bulletin of the Seismological Society of America, 1982, 72(4):1167-1183 https://www.researchgate.net/publication/252396578_Effect_of_surface_topography_on_the_diffraction_of_P_SV_and_Rayleigh_waves
[22] Sánchez-Sesma FJ, Campillo M. Diffraction 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 http://www.bssaonline.org/content/81/6/2234.short
[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 http://bssa.geoscienceworld.org/content/84/4/1169
[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 doi: 10.1002/eqe.4290240109/references
[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 doi: 10.1002/eqe.4290230507/full
[26] Haskell NA. The dispersion of surface waves on multilayered media. Bulletin of the Seismological Society of America, 1953, 43(1):86-103 http://oai.dtic.mil/oai/oai?verb=getRecord&metadataPrefix=html&identifier=AD0638148
[27] Thomson WT. Transmission of elastic waves through a stratified solid medium. Journal of Applied Physics, 1950, 21(2):89-93 doi: 10.1063/1.1699629
[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 diffraction of SH waves by an irregular surface. The Journal of the Acoustical Society of America, 1985, 77(1):1-5 doi: 10.1121/1.392258
[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 https://www.researchgate.net/publication/277838914_Time-domain_response_of_a_semi-circular_canyon_for_incident_SV_P_and_Rayleigh_waves_calculated_by_the_discrete_wavenumber_boundary_element_method
[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 http://bssa.geoscienceworld.org/content/ssabull/77/5/1784.full.pdf
[32] Shah AH, Wong KC, Datta SK. Diffraction 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, Noorhoff, 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 doi: 10.1016/0165-2125(82)90025-7
[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 doi: 10.1016/j.soildyn.2013.02.007
[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 http://bssa.geoscienceworld.org/content/64/5/1551.full.pdf
[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 finite element method in structural dynamics. International Journal for Numerical Methods in Engineering, 2009, 77(8):1139-1171 doi: 10.1002/nme.v77:8
[44] Lin G, Han Z, Li J. Soil-structure interaction analysis on anisotropic stratified medium. Géotechnique, 2014, 64(7):570-580 doi: 10.1680/geot.14.P.043
[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 https://www.researchgate.net/publication/227813746_Three-dimensional_response_of_a_layered_cylindrical_valley_embedded_in_a_layered_half-space
[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 https://www.researchgate.net/publication/277796192_Scattering_of_plane_harmonic_P_SV_and_Rayleigh_waves_by_dipping_layers_of_arbitrary_shape
[48] de Barros FCP, Luco JE. Seismic response of a cylindrical shell embedded in a layered viscoelastic half-space. Ⅱ:Validation and numerical results. Earthquake Engineering & Structural Dynamics, 1994, 23(5):569-580 https://www.researchgate.net/publication/230372924_Seismic_response_of_a_cylindrical_shell_embedded_in_a_layered_viscoelastic_half-space_II_Validation_and_numerical_results
[49] Lin G, Du J, Hu Z. Earthquake analysis of arch and gravity dams including the effects of foundation inhomogeneity. Frontiers of Archite https://www.researchgate.net/publication/227275536_Earthquake_analysis_of_arch_and_gravity_dams_including_the_effects_of_foundation_inhomogeneity
b50 Yin X, Li J, Wu C, et al. ANSYS implementation of damping solvent stepwise extraction method for nonlinear seismic analysis of large 3-D structures. Soil Dynamics & Earthquake Engineering, 2013, 44(1): 139-152 http://industry.wanfangdata.com.cn/dl/Detail/NSTLQK?id=NSTLQK_NSTL_QKJJ0230330084
[51] 林皋, 韩泽军, 李建波.层状地基任意形状刚性基础动力响应求解.力学学报, 2012, 44(6):1016-1027 doi: 10.6052/0459-1879-12-101 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(inChinese)) doi: 10.6052/0459-1879-12-101
[52] 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(inChinese))
[53] 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
[54] 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)
[55] Zhong W, Williams FW, Bennett PN. Extension of the wittrick-williams algorithm to mixed variable systems. Journal of Vibration & Acoustics, 1997, 119(3): 334-340
[56] Zhong WX. On precise integration method. Journal of Computational & Applied Mathematics, 2004, 163(1): 59-78
[57] de Barros FCP, Luco JE. Amplification of obliquely incident waves by a cylindrical valley embedded in a layered half-space. Soil Dynamics & Earthquake Engineering, 1995, 14(3): 163-175
[58] Wolf JP, Song C. Finite-element Modelling of Unbounded Media. Chichester: Wiley, 1996
下载:
计量
- 文章访问数: 1191
- HTML全文浏览量: 219
- PDF下载量: 823
