局部褶皱对层状地基马蹄形孔洞散射的影响分析

李志远，钟 红，胡志强，林 皋

(1. 海岸与近海工程国家重点实验室，大连理工大学，大连 116024；2. 工程抗震研究所，建设工程学部水利工程学院，大连理工大学，大连 116024；3. 中国水利水电科学研究院工程抗震研究中心，北京 100048)

多次震害调查表明，地下结构对地震波有散射、折射、衍射等作用，显著改变了地表的地震响应。水平岩层在构造运动作用下会产生弯曲，形成局部褶皱。现有研究多假定场地为均质或水平成层场地，忽略了局部复杂地形对地下结构散射场的影响。因此研究局部褶皱对层状地基中地下结构散射作用的影响对地震安全性评价具有重要指导意义。

天然的工程场地常被假定为水平成层分布。Thomson 和Haskell[1 − 2]提出的传递矩阵算法，可以提供分层介质中弹性波动控制方程的系统解。为了解决求解时指数溢出的问题，Kausel 和Roesset[3]重新推导了力和位移关系，并建立了精确刚度法，Chen[4]定义了全局广义的反射-透射矩阵。但是，这类算法只适用于求解具有规则边界条件的水平成层地基，对包含不规则边界的场地不能直接使用。Chen[5]利用模态综合法研究了不规则层状地基中Love 波的传播特性，Lee 等[6 − 7]采用加权余量法求解了体波和面波在不规则层状地基的传播特性。值得指出的是以上研究假定不规则土层均限于有限区域，远场仍为水平成层场地，且不包含孔洞等局部构造。Luco 和Apsel[8]基于传递矩阵法提出了一种层状地基的格林函数，在此基础上得到了移动荷载在层状地基的动力响应[9]，利用层状地基的移动荷载的格林函数，构造了一种间接边界元法，可以将求解三维结构的动力响应。利用这种间接边界元法，Barro 和Luco[10]分析了SV 波、P 波入射时，层状地基中埋置孔洞的对地表幅值的响应。他们还分析了层状半空间存在衬砌结构的隧洞时，入射SV 波、P 波，孔洞衬砌的动应力集中系数[11]。Wuttke 等[12]结合离散波数法和边界元法给出了一种两步法用于分析层状地基中波动传播。国内学者尤红兵和梁建文[13]利用间接边界元法求解了层状地基中洞室对SV 波的散射。巴振宁等[14]利用间接边界元法分析了场地动力特性对衬砌隧道地震反应的影响。赵密等[15]提出一种土-结构相互作用的高效方法，提高了深厚土层地下结构地震分析的计算效率。

值得指出的是以上文献求解散射场是通过对散射波和自由场波动的叠加而得到，这种方法将地下结构与周围土体作为一个整体分析，散射波的求解通常需要对散射体施加边界条件，当结构复杂时处理起来会变的异常复杂，因此以上研究的散射体均假设为圆形、椭圆形等。考虑多种复杂条件，如包含局部褶皱层状地基中马蹄形孔洞散射波的高精度解几乎没有。

本文基于Lin 等[16]提出散射场分析模型，通过一条人工边界将层状介质中复杂边界条件的波动散射问题转化为子结构波动响应求解的问题。这种转化对于线弹性介质是准确的，没有引入任何假设，当辐射边界条件和波动输入为高精度解时，复杂边界的散射问题也是高精度的。基于本文方法，讨论局部褶皱对层状场地中马蹄形孔洞散射场的影响。

1 模型及计算方法

1.1 子结构法控制方程

图1 为典型的土-结构相互作用体系，虚线ABCD 包围的区域为近场复杂结构，由虚线ABCD一直延伸到无穷远处的部分为远场，虚线ABCD交界面。本文只列出关键公式，详细推导见文献[17]。频域中子结构法的控制方程为：

图1 土-结构相互作用体系示意图Fig.1 Soil-structure interaction system

式中： K 和 M分别为通过有限元计算得到刚度阵和质量阵； ω为圆频率。注意到当近场结构为规则开挖土体时，整个体系退化为具有规则边界的自由场，利用式(1)对自由场体系进行分析，可以得到如下关系：

将式(3)代入式(1)中可得：

式(4)即为子结构法的控制方程。由式(4)可以看出，欲求复杂场地的动力响应，需要计算三部分内容：自由场交界面处的动力刚度和波动位移；近场子结构和开挖土体的动力刚度；控制方程的求解。

1.2 层状地基辐射边界

地下隧洞等结构沿轴向长度远大于截面尺寸，因此可以考虑为二维平面应变问题。建立直角坐标系，沿水平向右为x 轴正向，垂直向下为z 轴正向。层状地基的材料参数为：拉梅常数 λ和G，密度为 ρ，泊松比ν。忽略体力，用位移表示的波动方程为：

引入应力作为对偶变量可将方程化简为一阶常微分方程：

式中：V=[ uxux−σxz−σzz]T； H为水平波数κ、频率 ω和地基材料参数的函数[18]。任选一土层分析，定义土层[za,zb]，上层面为 za，下层面为zb，土层的厚度为ς，无外力作用时，土层两端位移向量ub和应力向量 pa，满足如下方程：

相邻土层在层面处满足位移连续和应力连续条件。相邻土层1[za,zb]和土层2[zb,zc]都可表示成对偶方程式(7)的形式，利用层面间位移和应力的连续条件，消去层面b 处的位移和应力向量，可得到合并后的土层c[za,zc]的对偶方程，得到系数矩阵间的关系为：

采用精细积分计算，土层厚度可以取的非常小(如τ=ς/220)，这时初始区段的系数矩阵可通过Taylor 级数展开获得[19]。当土层厚度非常小时，F趋于单位阵，式(8)不再适用。因此针对微层间的合并时，式(8)做如下修改：

其中，L=(I+G1Q2)−1。

由初始土层进行合并，可得到任意层面间的对偶关系。下面采用子结构法求解频域-波数域中的整体刚度矩阵。首先对层状地基在竖直方向进行离散，在计算点处设置土层层面，其次对相邻计算点所在土层进行合并，以区段[za,zb]为例，可得到层面间的对偶方程(式(7))。为便于使用子结构法，将对偶方程转化成刚度形式：

式中：

或者写成紧凑形式：

式中： P 表示外荷载的向量； K为区段的刚度矩阵；U 为位移向量。

最后是整体刚度矩阵的组装。在这种情况下，全局载荷矢量对应于界面处的外部作用力，整体刚度的组装可应用结构动力学子结构法的定理和技术。整体的刚度矩阵通过在层状体系的每个“节点”(层面)处叠加相邻土层的贡献来获得：

式中： R∞为半空间层面处的动力刚度[19]；n 为土层个数。

1.3 自由场波动输入

如图2 所示，无限半空间上卧n−1 个均匀弹性层，有n 层层面，图中 ρj、 λj、Gj表示第j 层的密度和拉梅常数。每层单独建立坐标系，坐标原点设在该层的顶面上，水平向右为x 轴正向，竖直向下为z 轴正向。

图2 层状地基中波动传播Fig.2 Wave propagation in a layered half-space

平面内右行简谐波的一般解可写成：

式中： φ、 ψ分别为SV 波、P 波的位势函数； κ1为水平向波数； κ2S、 κ2P分别为SV 波和P 波的竖直向波数； ES、 EP分别为SV 波、P 波上行波幅值；FS、 FP分别为SV 波、P 波下行波幅值。位移满足如下关系：

将式(14)代入式(15)可得各层面处的位移表达式；由位移-应变关系式，以及应力-应变关系式可以得到各层面处的应力表达式。引入应力-位移状态向量，整理后可表示为：

式中：

顶层的波幅向量可表示为：

式中：RSP、 RSS为SV 波的反射系数；RPS、 RPP为P 波的反射系数，可通过地表自由面的边界条件求得[20]。

下面建立半空间和顶层波幅向量间的关系，并求解顶层和半空间的波幅向量。由第2 节可知，层状体系可以通过精细积分进行高效合并。利用精细积分合并层状体系，最终得到界面zn−1和界面 z0间的对偶方程：

式中，F1n、G1n、Q1n为合并层状体系后得到的系数矩阵，结合式(16)可得半空间和顶层的波幅向量关系：

其中：

式(30)写成紧凑形式：QEE=[amk]， (m,k=1,2,···,4)。将式(20)代入式(23)可得：

其中：

当P 波入射时，EPn=1， ESn=0；当SV 波入射时，ESn=1， EPn=0。求解式(24)可以得到顶层的上行波的幅值和半空间中下行波幅值。将顶层上行波系数和反射系数带入式(20)可得到顶层的波幅向量HE1，将顶层波幅向量代入式(16)可得顶层的状态向量，其他各层的状态向量可由精细积分方法合并得到两端边值的对偶关系，最终可求解。本文提出的层状自由场的求解方法，可称为改进传递矩阵法，与传统传递矩阵法的主要区别在微层传递矩阵，即式(23)，其优点在于基于区段混合能的表示，使 Q和 G对于层厚很小时也不会发生数值病态[21]，保证了求解的顺利进行。同时，改进传递矩阵法的便于推广到各向异性层状地基波动问题的求解。

2 方法验证

计算均质半空间中垂直入射SV 波时地表位移幅值的响应，与文献[11]的结果对比验证本文方法的精度。如图3 所示，均质半空间存在浅埋圆柱形孔洞，均质半空间的密度ρ0=2200 kg/m3，泊松比ν=0.333，剪切波速cS0=600 m/s。圆柱形孔洞半径r=6.0 m，埋置深度为H=1.5r=9.0 m。虚线ABCD 包括的部分为近场区域4r×3.5r。定义位移响应幅值与入射波幅值的比值为无量纲位移幅值。垂直入射SV 波，入射波频率为η=ωr/πcS0=2.0，图4 为圆柱形孔洞场地地表的位移幅值，计算结果与Luco 和Barros[11]采用间接边界元法的计算结果对比，两者吻合性良好。

图3 均质半空间埋置圆形孔洞Fig.3 Schematic view of the shallow embedded circular cavity

图4 浅埋圆柱形孔洞场地垂直入射SV 波地表位移幅值Fig.4 Surface displacement generated from scattering by a shallow embedded circular cavity for vertically incident SVwave with particle motion along x-axis

3 算例与分析

地质运动作用下，褶皱是一种常见的工程地形，本节讨论层状半空间与上卧土层交界面为正弦波形状时的对马蹄形孔洞散射作用的影响。本文假定不规则土层交界面为有限区域，远场仍为水平成层地基。

计算模型参数为：厚度为H=4r的单层土体，置于半空间之上；马蹄形孔洞的穹顶半径为r=10 m，边墙高度为0.785r=7.85 m，埋置深度D=1.5r=15 m。半无限空间的材料参数：质量密度ρ0=2000 kg/m3，泊松比ν0=0.35，剪切波速cS0=700 m/s，阻尼比ζ0=0.02；上卧土层的密度ρ1=2000 kg/m3，泊松比ν1=0.35，阻尼比ζ1=0.05，剪切波速cS1=350 m/s2。定义无量纲频率η=ωr/(πcS1)=1.0。近场场地采用有限元计算近场模型，最大的单元尺寸为最小波长λmin=2r/ηmax=r的1/8。

如图5 所示，讨论两种不同的工况，工况A 是褶皱埋入半空间中，工况B 是褶皱埋入土层中。实际工程中褶皱可能是各种形状，本文为方便描述，讨论正弦形褶皱，褶皱形状分布满足如下方程：

图5 局部褶皱层状地基中马蹄形孔洞Fig.5 Horseshoe-shaped cavity embedded in a layered halfspace with local folds

褶皱的中心点与马蹄形孔洞的中心在同一竖直线上，褶皱的分布范围为x=−1.5r 和x=1.5r之间。上卧土层的厚度为H=4r，马蹄形孔洞的埋置深度D=2.0r。入射SV 波，入射波频率η=2，图6 为工况A 和工况B 地表的位移幅值，作为参考，水平成层场地(即不考虑褶皱地形)的地表位移幅值也绘制在图中。保持参数不变，入射频率η=2的P 波，计算结果绘制在图7 中。从图6 中可以看出，SV 波入射时局部褶皱对地表位移幅值影响明显；在地表−r

为了进一步研究局部褶皱对散射波的影响随频率变化的特性，图8 中绘制了SV 波入射时观测点P1 (x/r=1.5) 的位移频谱。从图8 可以看出，在低频段(η<0.75)局部褶皱对地表水平向和竖直向位移幅值的影响均较弱；随着频率的增加，局部褶皱对地表位移幅值的影响增强，且不同的局部褶皱几何构造对地表位移幅值影响的程度不同。从整体上观察位移频谱的分布情况，可以看出：在低频段(η<0.75)位移幅值有明显的共振现象，可以认为，此时局部褶皱和马蹄形孔洞对散射场影响均较弱，这时的波动响应接近自由场的波动响应；在0.75<η<1.5区段，位移幅值较小，也没有明显的共振现象，这时马蹄形孔洞是影响散射场分布的主要因素；在1.5<η时，三种工况的位移幅值分布显著不同，局部褶皱的影响明显，此时局部褶皱已经显著的改变了位移幅值峰值和共振频率。图9 为P 波入射时观测点P1 (x/r=1.5) 的位移频谱。与SV 波入射时相比，具有如下特点：局部褶皱对地表位移幅值影响变显著的起始频率更低，SV 波出现在η=0.75附近，P 波出现在η=0.5附近；局部褶皱对地表位移幅值的影响不是简单的随着频率的增加而增强，变化更加复杂，比如在1.5<η<1.75区间，三种工况的位移幅值比较接近，此频率区段局部褶皱对位移幅值的影响要低于频段1.0<η<1.5。

图6 SV 波入射时地表位移幅值Fig.6 Surface displacements caused by SV wave

图7 P 波入射时地表位移幅值Fig.7 Surface displacements caused by P wave

图8 SV 波入射时观测点P1 的位移频谱 (P1 (x/r=1.5))Fig.8 Displacement spectrums of observation point P1 caused by SV wave (P1 (x/r=1.5))

图9 P 波入射时观测点P1 的位移频谱 (P1 (x/r =1.5))Fig.9 Displacement spectrums of observation point P1 caused by P wave (P1 (x/r=1.5))

场地地震安全性评价中，位移峰值是一个重要参数。因此，图10 绘制了SV 波入射时地表各处位移峰值分布图，纵坐标表示频率在η ∈[0,2]区间中位移幅值的峰值。从图10 可以看出，工况A 的褶皱几何构造增大了地表水平向位移峰值，工况B 的褶皱几何构造降低了地表水平向位移峰值；对于竖直向位移幅值，在近马蹄形孔洞区域(−r

图10 SV 波入射时地表位移峰值Fig.10 Peak surface displacements caused by SV wave

图11 P 波入射时地表位移峰值Fig.11 Peak surface displacements caused by P wave

4 结论

本文提出了一种求解近区域包含复杂构造的层状场地波动响应的算法。基于子结构法建立了复杂场地散射问题的控制方程，将复杂场地散射问题转化具有规则边界条件层状地基的动力刚度求解和波动响应求解。这种转化对于线弹性介质是完全准确的，因此当具有规则边界条件层状地基的辐射问题和波动输入问题得到高精度解答时，复杂场地散射问题的解答也是高精度的。通过与现有文献结果对比，验证了本文方法的精度，并利用该方法，分析了SV 波和P 波入射时，局部褶皱对层状马蹄形孔洞场地的散射场影响。

研究结果表明，局部褶皱对马蹄形孔洞的散射场具有显著影响。在低频段局部褶皱对SV 波入射时马蹄形孔洞散射场的影响不明显，随着频率的增加，局部褶皱的影响增大；局部褶皱对P 波入射时马蹄形孔洞散射场的影响更加复杂，在某些高频段，也出现了局部褶皱对散射场影响不显著的情况。在SV 波或P 波入射时，在近马蹄形孔洞区域，局部褶皱对竖直向位移峰值影响明显，在距离马蹄形孔洞稍远的位置，局部褶皱对竖直向位移峰值几乎没有影响；整体分析可得，地表竖直向位移峰值受局部褶皱和马蹄形孔洞的共同影响；地表水平向位移峰值受褶皱地形的影响更显著。本文仅分析了两种简单的局部褶皱构造对散射场的影响，但本文方法可用于任意近区域复杂构造对地下结构散射场影响的分析。需要指出，本文假定孔洞和褶皱均满足平面应变理论，若考虑场地三维散射效应，需将控制方程扩展到2.5 维。