P－SV波斜入射时成层半空间自由场的时域算法①

赵 密，杜修力，刘晶波，韩 强

（1．北京工业大学城市与工程安全减灾教育部重点实验室，北京 100124；2．清华大学土木工程系，北京 100084）

0 引言

地震作用下土与结构动力相互作用分析中，首先需要计算地震自由波场作为波动输入［1］。在自由场计算中，入射地震波可以简化为竖直入射或者倾斜入射的平面体波，场地可以简化为弹性均匀半空间或者弹性成层半空间。在一般工程结构的地震反应分析中，可以采用地震波竖直入射的均匀半空间或者成层半空间场地模型，该问题是空间一维的，能够方便地在时域内求解。对于高坝、核电站、大型地下结构等重大工程结构，则需要考虑地震波倾斜入射引起的场地运动非一致变化对结构地震反应的影响。地震波斜入射时，均匀半空间自由场可以根据波动传播规律在时域内计算［2－4］；而成层半空间自由场只能采用解析方法在频域内求解［5］，为了获得时域反应需要多次执行快速傅里叶变换，因而需要大量的计算时间和存储空间。

为了在时域内计算地震波斜入射时成层半空间自由波场，刘晶波和王艳［7－9］提出了一种简便易行的空间一维化计算方法。该方法采用矩形有限单元进行网格虚拟划分，竖直方向采用满足波动精度要求的网格尺寸；水平方向的网格尺寸等于时间步长与水平视波速的乘积。基底设置黏性人工边界条件［10］模拟基岩半空间的辐射阻尼，并将波动输入转化为等效荷载施加在人工边界结点上。将集中质量有限元法和时间中心差分相结合建立结点运动方程组。然后根据斯奈尔定律，将水平方向相邻结点的运动用该结点相邻时刻的运动表示，从而将求解结点运动的空间二维方程组化为空间一维方程组。求解此方程组，即得到自由波场中竖向一列结点的运动。最后根据行波传播的特点，可方便地确定整个空间的自由波场。这样，计算地震波斜入射时水平成层半空间自由场的空间二维问题转化为简单的空间一维问题。

对于出平面SH波问题［6］，修正的黏性边界能够完全吸收外行波，精确模拟基岩半空间的辐射阻尼；然而对于平面内P－SV波问题［7－8］，黏性边界及其修正形式均不能完全吸收外行波，形成的虚假反射波将污染自由场反应。地震波倾斜入射角越大，计算结果与真实波场之间的误差越大。实际上，黏性边界［10］、黏弹性边界［11－12］、多次透射公 式［13］、高精度人工边界条件［14－15］均不能精确和高效地模拟该问题中基岩半空间的辐射阻尼。本文利用平面内波动问题中外行波是传播方向已知的平面P波和SV波这一波动特性，提出一种适用于该问题的精确人工边界条件，应用该条件替代黏性边界，建立PSV波斜入射时成层半空间自由场的一维化时域算法。

1 人工边界和地震荷载

地震波（平面P－SV波）斜入射时成层半空间波动问题如图1所示。水平成层介质位于半空间上，成层介质表面满足应力为零的自由边界条件，各层介质交界面处满足位移和应力连续条件。第j（1…N）层弹性均匀介质的厚度为hj，密度为ρi，λi拉梅常数为 和Gj。下部弹性均匀半空间的密度为ρ，拉梅常数为λ和G。平面P波或者SV波从半空间无限远处向右上方倾斜入射，入射角为θ。在计算中，将下部半空间截去，此时第N层介质与半空间的交界面为人工边界，本节建立能够精确模拟下部半空间辐射阻尼的人工边界条件，给出人工边界处入射波的等效荷载。

图1 地震波斜入射时成层半空间波动问题Fig．1 Wave propagation problem in layered half space under seismic wave of oblique incidence．

1．1 精确的人工边界条件

半空间内向无穷远处传播的外行波由平面P波和平面SV波构成，两者的势函数可以分别表示为

当P波入射时，θP＝θ，可由式（3）确定θS；当SV 波入射时，θS＝θ，可由式（3）确定θP。位移的势函数表达式为

其中，上标撇号表示波形函数的导数。

对式（5）求时间偏导数，得到速度向量为

其中，变量上方的点号表示时间导数。弹性均匀介质的应力－位移本构关系为

对于外行波，将位移式（5）代入式（7）得到应力向量为

联立式（6）和（8），消去波形函数得

其中，阻抗矩阵为

在人工边界处，式（9）是精确的人工边界条件，它由边界速度计算半空间作用于成层介质的应力。有限元计算中，该边界条件可以采用集中离散方式实现。可以看到，对于不同的人工边界结点，本文边界条件是空间解耦的；但对于每个结点的两个自由度，本文边界条件是耦联的，这一点与黏性边界不同，正是耦联特性保证本文边界条件能够精确地同时吸收外行P波和SV波。

1．2 人工边界上的地震荷载

本节给出尚未施加人工边界条件时，入射波在人工边界处的等效荷载。

1．2．1 平面P波入射情况入射平面P波的势函数为

如果规定入射波合速度的正方向与入射波传播方向一致，则在人工边界和y轴交点处已知的入射P波合速度时程可以表示为

联立式（12）和（13），消去波形函数得

对于入射波，将式（11）代入式（4），然后代入式（7），得入射波作用下半空间作用于成层介质的应力向量

联立式（13）和（15），消去波形函数得

1．2．2 平面SV波入射情况

入射平面SV波的势函数为

其中，上标i表示入射波；是任意波形函数。

如果规定入射波合速度的正方向为入射波传播方向顺时针旋转90°角的方向，则在人工边界和y轴交点处已知的入射SV波合速度时程可以表示为

对于入射波，将式（17）代入式（4），然后代入式（7），得入射波作用下半空间作用于成层介质的应力向量为

联立式（19）和（21），消去波形函数得

2 有限元结合中心差分的一维化算法

对于半空间上部的N层介质，采用矩形四结点双线性平面应变有限单元进行虚拟网格划分，沿y轴方向的单元尺寸Δy为，沿x轴方向的单元尺寸为Δx，时间步长为 。x＝0列结点p时刻的集中质量有限元方程为

人工边界处波场满足叠加原理，可以分解为外行波场和入射波场，x＝0处人工边界结点的荷载和速度可以分别写为

时间导数的中心差分公式为

成层半空间中的平面波服从斯奈尔定律，即全部平面波的水平方向视波速相等。如果取水平方向的有限元尺寸Δx等于时间步长Δt与水平视波速cx的乘积，即Δx＝cxΔt。则根据斯奈尔定律，x＝0列结点p＋1时刻的位移可以用x＝－Δx列结点p时刻的位移表示，x＝0列结点p－1时刻的位移可以用x＝Δx列结点p时刻的位移表示，即

这样，空间二维问题能够被转化为简单的空间一维问题。将式（29）～（32）代入式（26），整理得到x＝0列结点从前两时刻位移递推当前时刻位移的线性代数方程组

式（33）的初始值为

逐个时刻的求解式（33），可以得到自由波场中竖向一列结点的运动。根据波动水平方向视速度相等的特点，可以方便地确定整个空间的自由波场。

竖直方向有限元尺寸Δy应满足波动精度要求。确定时间步长Δt的数值稳定性条件，以及方法的精度分析见文献［7－9］。

3 数值算例

为了验证本文改进方法的精度和有效性，分析半空间基岩表面覆盖单层土问题，该算例也见文献［9］的3．8节。土层厚度为50m，质量密度为1 000 kg／m3，P 波 波 速 为 346m／s，SV 波 波 速 为 200 m／s。基岩的质量密度为1 500kg／m3，P波波速为866m／s，SV波波速为500m／s，人工边界截去地表深度100m以下的基岩，即保留厚度50m的基岩作为有限元区域。在人工边界和y轴交点处入射波位移时程如图2所示。深度方向的空间步距为5 m，满足数值稳定性条件的时间步长为0．005s。

图2 入射波位移时程Fig．2 Time history of displacement of incident wave．

首先研究平面P波斜入射情况。文献［7，9］的研究表明，入射角越大，应用黏性边界的自由场算法的精度越低，限于篇幅，本文仅研究60°角大角度入射的情况。图3给出P波60°角入射时地表（坐标原点）的位移时程结果，其中理论解是采用频域传递矩阵法［5］结合快速傅里叶变换得到的。为了更清楚地观察本文方法相对于刘晶波－王艳（Liu－Wang）方法［7－9］的改进效果，图3也给出了两种方法计算结果相对于理论解的误差，即每种方法结果与理论解的差的绝对值。可以看到，新方法的计算结果与理论解吻合良好，而Liu－Wang方法的计算误差相对较大，说明本文人工边界条件能够进一步提高P波大角度入射时平面内自由场一维化时域算法的精度。

当平面SV波斜入射时，本文仅考虑入射角小于临界角的情况，即自由波场中仅含有均匀平面波的情况。与P波入射类似，本文研究SV波30°角入射的情况。图4给SV波30°角入射时地表（坐标原点）的位移时程结果，其中理论解同样是采用频域传递矩阵法［5］结合快速傅里叶变换得到。图4也给出了本文新方法计算结果和Liu－Wang方法［8－9］计算结果相对于理论解的误差。可以看到，新方法的计算结果与理论解吻合良好，而Liu－Wang方法的计算误差相对较大，说明本文人工边界条件能够进一步提高SV波大角度入射时平面内自由场一维化时域算法的精度。

图3 P波60°角入射时地表位移时程Fig．3 Time history of displacement on ground surface under P wave of 60°angle incidence．

4 结语

图4 SV波30°角入射时地表位移时程Fig．4 Time history of displacement on ground surface under SV wave of 30°angle incidence．

刘晶波和王艳［6－9］建立了一种在时域内计算地震波斜入射时成层半空间瞬态自由场的数值方法。本文对该方法进行了改进，提出一种新的精确人工边界条件替代黏性边界，解决了平面P－SV波大角度入射时由于底部吸收边界条件精度低而导致整套自由场时域算法精度降低的问题。本文考虑自由波场中仅含有均匀平面波的情况，自由波场中含有非均匀平面波和面波的波动问题有待深入研究。

［1］ Wolf J P．Dynamic Soil－Structure Interaction［M］．Englewood Cliffs：Prentice－Hall，1985．

［2］ 徐海滨，杜修力，赵密，等．地震波斜入射对高拱坝地震反应的影响［J］．水力发电学报，2011，30（6）：159－165．

XU Hai－bin，DU Xiu－li，ZHAO Mi，et al．Effect of Oblique Incidence of Seismic Waves on Seismic Responses of High Arch Dam［J］．Journal of Hydroelectric Engineering，2011，30（6）：159－165．

［3］ 苑举卫，杜成斌，刘志明．地震波斜入射条件下重力坝动力响应分析［J］．振动与冲击，2011，30（7）：120－126．

YUAN Wei－ju，DU Chen－bin，LIU Zhi－ming．Time－domain Seismic Response for Gravity Dam to Obliquely Incident and Seismic Waves［J］．Journal of Vibration and Shock，2011，30（7）：120－126．

［4］ 郜新军，赵成刚，刘秦．地震波斜入射下考虑局部地形影响和土结动力相互作用的多跨桥动力响应分析［J］．工程力学，2011，28（11）：237－243．

GAO Xin－jun，ZHAO Cheng－gang，LIU Qin．Seismic Response Analysis of Multi－span Viaduct Considering Topographic Effect and Soil－structure Dynamic Interaction Based on Inclined Wave［J］．Engineering Mechanics，2011，28（11）：237－243．

［5］ 傅淑芳，刘宝诚．地震学教程［M］．北京：地震出版社，1991．

FU Shu－fang，LIU Bao－cheng．Seismology Tutoria［M］l．Beijing：Seismic Press，1991．

［6］ 刘晶波，王艳．成层半空间出平面自由波场的一维化时域算法［J］．力学学报，2006，38（2）：219－225．

LIU Jing－bo，WANG Yan．A 1－D Time－domain Method for 2－D Wave Motion in Elastic Layered Half－space by Antiplane Wave Oblique Incidence［J］．Chinese Journal of Theoretical and Applied Mechanics，2006，38（2）：219－225．

［7］ 刘晶波，王艳．成层介质中平面内自由波场的一维化时域算法［J］．工程力学，2007，24（7）：16－22．

LIU Jing－bo，WANG Yan．A 1DTime－domain Method for Inplane Wave Motion of Free Field in Layered Media［J］．Engineering Mechanics，2007，24（7）：16－22．

［8］ LIU Jing－bo，WANG Yan．A 1DTime－domain Method for inplane Wave Motions in A Layered Half－space［J］．Acta Mechanica Sinica，2007，23：673－680．

［9］ 王艳．非一致地震动场数值方法研究及在结构动力分析中的应用［D］．北京：清华大学，2007．

WANG Yan．Research on the Numerical Method for Asynchronous Seismic Wave Motions and Its Application in Dynamic Analysis of Structures［D］．Beijing：Tsinghua University，2007．

［10］ Lysmer J，Kulemeyer R ．Finite Dynamic Model for Infinite Media［J］．Journal of Engineering Mechanics，ASCE，1969，95：759－877．

［11］ LIU Jing－bo，LÜYan－dong．A Direct Method for Analysis of Dynamic Soil－structure Interaction Based on Interface Idea［A］／／Proceedings of the Chinese－Swiss Workshop on Dynamic Soil－Structure Interaction［C］．Beijing：International Academic Publishers，1997．

［12］ 杜修力，赵密，王进廷．近场波动模拟的人工应力边界条件［J］．力学学报，2006，38（1）：49－56．

DU Xiu－li，ZHAO Mi，WANG Jin－ting．A Stress Artificial Boundary in FEA for Near－field Wave Problem［J］．Chinese Journal of Theoretical and Applied Mechanics，2006，38（1）：49－56．

［13］ 廖振鹏．工程波动理论导论（第二版）［M］．北京：科学出版社，2002．

LIAO Zhen－peng．Introduction to Wave Motion Theories for Engineering（Second Edition）［M］．Beijing：Science Press，2002．

［14］ DU Xiu－li，ZHAO Mi．A local Time－domain Transmitting Boundary for Simulating Cylindrical Elastic Wave Propagation in Infinite Media［J］．Soil Dynamics and Earthquake Engineering，2010，30（10）：937－946．

［15］ ZHAO Mi，DU Xiu－li，LIU Jing－bo，et al．Explicit Finite Element Artificial Boundary Scheme for Transient Scalar Waves in Two－dimensional Unbounded Waveguide［J］．International Journal for Numerical Methods in Engineering，2011，87（11）：1074－1104．