P-SV波斜入射时有阻尼成层介质自由波场的一维化时域算法

卓卫东，高智能，谷 音(.福州大学 土木工程学院， 福建 福州 3506；2.福州大学 福建省土木工程多灾害防治重点实验室， 福建 福州 3506)

P-SV波斜入射时有阻尼成层介质自由波场的一维化时域算法

卓卫东1，2，高智能1，谷 音1，2
(1.福州大学 土木工程学院， 福建 福州 350116；2.福州大学 福建省土木工程多灾害防治重点实验室， 福建 福州 350116)

为了研究土层介质阻尼对自由波场计算的影响，在一维化时域算法基础上，采用有限差分方法，推导建立了考虑阻尼影响的水平成层弹性介质平面内自由波场求解的显式的数值逐步法公式，并采用MATLAB程序语言编制了相应的数值计算程序，进行平面内P-SV波斜入射情形下均匀弹性土层的算例分析。结果表明：阻尼对水平成层弹性介质平面内的自由波场有着重要的影响，与无阻尼均匀土层相比，P波斜入射下有阻尼均匀土层的水平和竖向位移幅值分别可减小53.1%和26.5%，SV波斜入射下水平和竖向位移幅值分别可减小48.2%和31.7%；阻尼使均匀土层自由表面的位移幅值明显减小，而中、下部节点的位移幅值减小相对较小。由此可见，土层介质阻尼对成层半空间平面内自由波场一维化计算的影响不可忽视。

成层介质；自由波场；一维化时域算法；土层介质阻尼；逐步法

在考虑地震作用下土-结构动力相互作用问题时，常需要计算自由波场以获得地震波动输入[1]。地震波动输入是否合理，直接影响计算结果的准确性。对一般工程结构，地震波可假设为竖直入射，此时自由波场的求解是一维问题，容易在时域内实现[2]；然而，对地铁车站、桥梁、大坝等大型大跨结构，地震波斜入射引起的地面运动非一致变化将对结构的地震响应产生重要影响[3-6]，此时需要考虑地震波的空间变异性。

在实际工程中，由于土层分布的复杂性，对自由波场的解析求解几乎是不可能的[7]。因此，采用有限元等数值方法研究复杂土层介质中地震波动问题，已经成为主要的研究方向。针对水平成层弹性半空间在地震波斜入射情形下的自由场计算问题，李山有等[8]根据波动传播水平视波速不变且已知的特点，建立了入射侧边界节点自由场波动计算的精确内插公式，并将其与计算内节点位移的显式差分公式相结合，进而提出了简化的时域数值模拟方法。刘晶波等[9-11]进一步提出一种一维化的时域算法，将水平成层弹性半空间在SH波和P-SV波斜入射下二维自由波场的计算问题简化为时域内的一维问题求解，且具有与二维有限元数值解同样的精度。赵密等[12]在刘晶波等[10]算法的基础上，提出一种模拟基岩半空间辐射阻尼的人工边界条件，采用该人工边界条件代替粘性边界条件，可获得更高的计算精度。

已有关于水平成层弹性半空间自由波场的简化时域算法，均未考虑土层介质阻尼的影响[8-12]。然而，土层介质阻尼可能对地震波斜入射下土层介质的自由波场有重要的影响。本文在刘晶波等[10]算法的基础上，考虑土层介质阻尼，推导建立P-SV波斜入射下水平成层弹性半空间自由波场求解的数值逐步法公式，并编制相应的数值计算程序；通过算例分析，讨论土层介质阻尼对P-SV波斜入射下水平成层弹性半空间自由波场的影响。

1 成层介质中平面波动问题的一维化

在土层介质自由场计算中，通常将土层介质假设为水平成层弹性半空间。按照刘晶波-王艳算法[10]，将水平成层弹性半空间划分为如图1所示的有限元离散化模型；图1中，竖向网格尺寸为Δy，可取为满足精度要求的任意值；水平方向网格尺寸为Δx，其取值需满足下式：

Δx=cx·Δt

(1)

其中：cx为P波或SV波的水平视波速，Δt为时间步长。

(2)

其中，M、C和K分别为离散化模型的质量矩阵、阻尼矩阵和刚度矩阵。

图1 水平成层半空间有限元模型

根据行波传播特点以及Snell定律，可以将y轴相邻列节点的运动用y轴上节点的运动表示，最终使P-SV波斜入射下水平成层弹性半空间自由波场的计算问题就转化为一维问题求解，具体公式推导详见文献[10]，这里不再赘述。

2 考虑介质阻尼的一维化时域算法

2.1 考虑介质阻尼的数值逐步法的计算公式

首先，采用刘晶波等[10]算法的思路，将P-SV波斜入射下水平成层弹性半空间二维自由波场的计算问题转化为一维问题。其次，采用数值方法，在时域内直接求解式(2)所列的运动方程。在数值计算中可取m=0，即先求得y轴上各节点的位移后，再根据行波传播特点，依次确定水平成层弹性半空间中的自由波场。以下，采用有限差分方法，推导建立考虑土层介质阻尼时求解待求列节点(y轴上节点)位移的计算列式。

根据式(2)，在考虑土层介质阻尼时，y轴上任一内节点(0,n)在pΔt时刻的运动方程为：

(3)

其中，M0,n和M0,0分别为节点(0,n)和节点(0,0)的集中质量；Ci,j、Ki,j分别为节点(0,n)和节点(i,j)之间的阻尼系数和刚度系数。

y轴与自由表面交界处的节点(0,0) 在pΔt时刻的运动方程为

(4)

在采用数值方法计算水平成层弹性半空间的自由波场时，显然必须从半无限介质中切取有限的计算区域，并在区域边界引入合适的人工边界条件。目前已发展了多种人工边界[13-22]，不失一般性，本文假定人工边界为黏性边界。对y轴与黏性边界交界处的节点(0,N)，其在pΔt时刻的运动方程为：

(5)

(6)

其中：ρ为介质质量密度；cs和cp分别为介质剪切波速和压缩波速；θ为地震波入射角；u0(xB，yB，t)和σ0(xB，yB，t)分别为入射波在黏性边界节点B上产生的位移矢量和应力矢量。

式(3)～式(5)所列的运动方程中的加速度项可利用中心差分法近似计算：

(7)

(8)

依据行波传播特点，同时将式(7)和式(8)代入式(3)～式(5)中，并简记u0,n=un，M0,n=Mn，经整理得到如下矩阵形式的方程：

(9)

式(9)中，

(10)

(11)

从式(10)和式(11)可见，式(9)所列的矩阵方程左边的系数矩阵是稀疏的三对角矩阵，右边的向量仅与边界节点输入的等效荷载以及待求列节点(y轴上节点)在pΔt时刻及其前一时刻(p-1)Δt的位移有关，故只要给定边界节点等效荷载以及初始时刻各节点的位移和速度值，即可通过逐步法求解上述方程组，得到y轴上各节点在时域内的位移解，进而确定全部自由波场。

2.2 计算步骤及程序实现

由式(9)所列的矩阵方程可见，本文基于有限差分方法建立的考虑介质阻尼的水平成层弹性半空间自由波场的数值逐步法的计算公式是显式的。具体计算步骤如下：

(1) 计算图1所示的有限元离散化模型的质量矩阵M和刚度矩阵K，其中质量矩阵M采用集中质量法计算，并假定阻尼矩阵C；

(2) 根据黏性人工边界条件，确定人工边界节点等效荷载时程；

(12)

(6) 根据行波传播特点确定计算区域内各节点的自由波场。

根据上述计算步骤，基于MATLAB语言编制了相应的计算程序。限于篇幅，这里没有给出计算程序的源代码。

2.3 算法的稳定性和计算精度

本文所建立的P-SV波斜入射下考虑介质阻尼的水平成层弹性半空间自由波场的一维化时域算法，其稳定性条件仍与刘晶波等[10]算法的稳定性条件相同，即要求：

(13)

为保证一维化时域算法的计算精度，刘晶波和王艳[10]建议，有限差分网格尺寸应满足以下离散化准则：

(14)

其中：Δ为网格尺寸；λ为波长。

3 算例分析

以下，选择物理性质均匀的土层进行算例分析。均匀土层的几何参数和物理参数取值参考文献[10]，具体如表1所列；其有限元离散化模型见图2，其中竖向网格尺寸Δy取为5 m，水平方向网格尺寸Δx根据式(1)确定。数值计算中，取P波波速为866 m/s，SV波波速为500 m/s；为满足稳定性条件，时间步长取为0.005 s。

表1 均匀土模型参数

假定土层介质阻尼为Rayleigh阻尼，则可按下式计算其阻尼矩阵C：

C=αM+βK

(15)

其中，α、β为比例系数，它们由下式计算：

(16)

式(16)中，ξ为均匀土层的阻尼比，这里取为0.1；ωi和ωj分别为均匀土层两个特定的自振圆频率，这里分别取为第1阶和第2阶自振圆频率。

利用ABAQUS软件建立均匀土层的有限元模型，通过动力特性分析，得到该均匀土层前两阶的自振频率：ω1=8.8645 rad/s，ω2=10.0280 rad/s；将其代入式(16)，求得α=0.9412，β=0.01058；将α、β的数值代入式(15)，得到均匀土层的阻尼矩阵C。

图2 均匀土层的有限元离散模型

[10]，假定入射P波(SV波)为持时0.5 s的单位脉冲，其位移时程如图3所示。

3.1 P波斜入射下均匀土层的自由波场分析

利用所编制的计算程序，计算了P波以不同入射角入射情形下均匀土层的自由波场；本文分析所取的入射角变化范围为0°～90°，并以15°为间隔。选取y轴上自由表面节点A1、中部节点B1以及底部边界节点C1作为观测点(见图2)，对计算结果进行分析讨论。由于无阻尼情形下一维化算法的计算结果已得到验证[10]，本文主要讨论考虑土层介质阻尼情形下一维化算法的计算结果。

图3 入射波位移时程

图4绘出了P波以30°角斜入射情形下，分别采用一维化算法及有限元法计算得到的有阻尼均匀土在3个观测点处的水平位移和竖向位移响应(限于篇幅，这里没有给出P波以其它角度斜入射情形下的对比分析，下同)。从图4中可以发现，两种数值方法计算得到的3个观测点的位移时程曲线规律一致；与有限元法计算结果相比较，一维化算法计算得到的位移时程存在“振幅衰减”和“位移超前”现象，这一现象可能与一维化算法采用了差分列式有关。总的来看，一维化算法的计算结果可以满足工程需要。

图5绘出了P波以30°角斜入射情形下，采用一维化算法计算得到的无阻尼和有阻尼均匀土在3个观测点处的水平位移及竖向位移响应。从图5中可以发现，无阻尼与有阻尼均匀土层的位移时程曲线波形一致；与无阻尼均匀土层相比，有阻尼均匀土层在各个时刻的位移响应绝对值均减小了，尤其是在自由表面观测点A1处，减小幅度最大。根据计算结果，自由表面观测点A1处的水平和竖向位移在0.34 s附近达到最大值；与无阻尼均匀土层相比，有阻尼均匀土层在自由表面观测点A1处的最大水平位移从1.11 cm(无阻尼时)减小至0.73 cm，减小幅度达34.2%；最大竖向位移绝对值从1.67 cm(无阻尼时)减小至1.24 cm，减小幅度达25.7%。可见，P波斜入射下土层介质阻尼对自由波场有重要的影响。

图6绘出了P波以不同入射角入射情形下，无阻尼和有阻尼均匀土层在表面观测点A1处的位移幅值随入射角的变化情况。从图6中可以发现，在入射角小于60°时，无阻尼和有阻尼均匀土层的水平位移幅值均随着入射角的增大而增大；在入射角超过60°时，水平位移幅值则随着入射角的增大而减小；而竖向位移幅值均随着入射角的增大而减小。与无阻尼均匀土层相比，相同入射角下有阻尼均匀土层的水平和竖向位移幅值均减小了；在入射角为75°时，两者水平位移幅值相差最大，在入射角为0°时，两者竖向位移幅值相差最大。计算结果表明，在P波入射角为75°时，有阻尼均匀土层在自由表面观测点A1处的水平位移幅值从1.62 cm(无阻尼时)减小至0.76 cm，减小幅度达53.1%；在P波入射角为0°时，有阻尼均匀土层在自由表面观测点A1处的竖向位移幅值从2.00 cm(无阻尼时)减小至1.47 cm，减小幅度达26.5%。

图4 一维化算法与有限元法计算结果比较

图5 P波30°斜入射时无阻尼和有阻尼均匀土层的位移响应比较

图6 观测点A1位移幅值随P波入射角的变化曲线

3.2 SV波斜入射下均匀土层的自由波场分析

利用所编制的计算程序，计算了SV波以不同入射角入射情形下均匀土层的自由波场。考虑到数值计算的稳定性要求，这里所取的入射角变化范围为0°～30°，并以5°为间隔。同样选取y轴上自由表面节点A1、中部节点B1以及底部边界节点C1作为观测点，对计算结果进行分析讨论。

图7绘出了SV波以30°角斜入射情形下，采用一维化算法计算得到的无阻尼和有阻尼均匀土在3个观测点处的水平位移及竖向位移响应。从图7中可以发现，SV波斜入射下，无阻尼与有阻尼均匀土层的位移时程曲线同样波形一致；与无阻尼均匀土层相比，有阻尼均匀土层在各个时刻的位移响应绝对值均减小了，尤其是在自由表面观测点A1处，减小幅度最大。根据计算结果，自由表面观测点A1处的水平和竖向位移在0.42 s附近达到最大值；与无阻尼均匀土层相比，有阻尼均匀土层在自由表面观测点A1处的最大水平位移从1.66 cm(无阻尼时)减小至0.86 cm，减小幅度达48.2%；最大竖向位移绝对值从1.04 cm(无阻尼时)减小至0.71 cm，减小幅度达31.7%。可见，SV波斜入射下土层介质阻尼对自由波场也有重要的影响。

图8绘出了SV波以不同入射角入射情形下，无阻尼和有阻尼均匀土层在表面观测点A1处的位移幅值随入射角的变化情况。从图8中可以发现，无阻尼和有阻尼均匀土层的水平位移幅值均随着入射角的增大而减小，而竖向位移幅值均随着入射角的增大而增大。与无阻尼均匀土层相比，相同入射角下有阻尼均匀土层的水平和竖向位移幅值均减小了；在入射角为30°时，两者水平和竖向位移幅值相差最大。计算结果表明，在SV波入射角为30°时，有阻尼均匀土层在自由表面观测点A1处的水平位移幅值从1.66 cm(无阻尼时)减小至0.86 cm，减小幅度达48.2%；竖向位移幅值从1.04 cm(无阻尼时)减小至0.71 cm，减小幅度达31.7%。

4 结 语

本文基于刘晶波等算法，提出了考虑土层介质阻尼影响的水平成层弹性半空间平面内自由波场的一维化时域算法，并编制了相应的数值计算程序。通过算例分析，讨论了土层介质阻尼对P-SV波斜入射下水平成层弹性半空间自由波场的影响。分析结果表明：

图7 SV波30°斜入射时无阻尼和有阻尼均匀土层

图8 观测点A1位移幅值随SV波入射角的变化曲线

(1) 介质阻尼对水平成层弹性半空间平面内的自由波场有重要的影响。与无阻尼均匀土层相比，P波和SV波斜入射下有阻尼均匀土层的水平位移幅值分别可减小53.1%和48.2%，竖向位移幅值分别可减小26.5%和31.7%。

(2) 与无阻尼均匀土层相比，有阻尼均匀土层自由表面的位移幅值减小较大，而中部和底部节点的位移幅值减小相对较小。

参考文献：

[1] Wolf J P. Dynamic Soil-structure Interaction[M]. Englewood Cliff: Prentice-Hall, 1985.

[2] 廖振鹏.工程波动理论导论[M].2版.北京:科学出版社,2002.

[3] Wolf J P, Obernhuber P. Effects of horizontally traveling waves in soil-structure interaction[J]. Nuclear Engineering and Design, 1979,57(2):221-244.[4] 李山有，廖振鹏,周正华.大型结构地震反应数值模拟中的波动输入[J].地震工程与工程振动,2001,6(2):1-5.[5] 潘旦光,楼梦麟,范立础.多点输入下大跨度结构地震反应分析研究现状[J].同济大学学报,2001,29(10):1213-1219.

[6] 杜修力,陈 维,李 亮,等.斜入射条件下地下结构时域地震反应分析初探[J].震灾防御技术,2007,2(3):290-296.

[7] 傅淑芳,刘宝诚.地震学教程[M].北京:地震出版社,1991.

[8] 李山有,王学良,周正华.地震波斜入射情形下水平成层半空间自由场的时域计算[J].吉林大学学报(地球科学版),2003,33(3):372-376.

[9] 刘晶波,王 艳.成层半空间出平面自由波场的一维化时域算法[J].力学学报,2006,38(2):219-225.

[10] 刘晶波,王 艳.成层介质中平面内自由波场的一维化时域算法[J].工程力学,2007,24(7):16-22.

[11] Liu Jing-bo, Wang Yan. A 1D time-domain method for inplane wave motions in a layered half-space[J]. Acta Mechanica Sinica, 2007,23(6):673-680.

[12] 赵 密,杜修力,刘晶波,等.P-SV波斜入射时成层半空间自由场的时域算法[J].地震工程学报,2013,35(1):84-90.

[13] Lysmer J, Kulemeyer R L. Finite dynamic model for infinite media[J]. Journal of Engineering Mechanics, ASCE, 1969,95(4):859-877.

[14] Wolf P J. A comparison of time-domain transmitting boundaries[J]. Earthquake Engineering and Structure Dynamics, 1986,14(14):655-673.

[15] Deeks A J, Randolph M F. Axisymmetric time-domain transmitting boundaries[J]. Journal of Engineering Mechanics, ASCE, 1994,120(1):25-42.

[16] 王振宇,刘晶波.成层地基非线性波动问题人工边界与波动输入研究[J].岩石力学与工程学报,2004,23(7):1169-1173.

[17] 刘晶波,王振宇,杜修力,等.波动问题中的三维时域粘弹性人工边界[J].工程力学,2005,22(6):46-51.

[18] 潘旦光,楼梦麟,董 聪.土层地震行波反应分析中侧向人工边界的影响[J].岩土工程学报,2005,27(3):308-312.

[19] 刘晶波,谷 音,杜义欣.一致粘弹性人工边界及粘弹性边界单元[J].岩土工程学报,2006,28(9):1070-1075.

[20] 杜修力,赵 密,王进廷.近场波动模拟的人工应力边界条件[J].力学学报,2006,38(1):49-56.

[21] 卢华喜,梁平英,尚守平.地基非线性波动问题中黏-弹性人工边界研究[J].岩土力学,2008,29(7):1911-1916.

[22] 张 波,李术才,杨学英,等.三维黏弹性介质人工边界研究[J].岩土力学,2009,30(11):3469-3475.

DOI:10.3969/j.issn.1672-1144.2016.06.007

A 1D Time-Domain Method for in-plane Wave Motion of Free Field in Layered Media with Damping under Obliquely Incident P-SV Waves

ZHUO Weidong1,2, GAO Zhineng1, GU Yin1,2

(1.CollegeofCivilEngineering,FuzhouUniversity,Fuzhou,Fujian350116,China；2.KeyLaboratoryforMultiDisasterPreventionandGovernanceofCivilEngineeringofFujianProvince,Fuzhou,Fujian350116,China)

In order to study the influence of layered-soil damping on the free field algorithm, based on a 1D time-domain method for in-plane wave motion of free field in layered media proposed by Liu and Wang, formulas of explicit time-stepping method to solve the in-plane wave motion of free field in layered media with damping are established by using finite difference method, and its numerical program has been developed with Matlab programming language. A case study of in-plane wave motion of free field in layered soil with damping under obliquely incident P wave and SV wave was carried out respectively. The results show that: the peak horizontal and vertical displacements in layered soil with damping can be reduced by 53.1% and 26.5%, respectively, under incident P wave, and the peak horizontal and vertical displacements can be reduced by 48.2% and 31.7%, respectively, under incident SV wave, comparing with layered soil without damping. The case study indicates that soil damping has significant influence on in-plane wave motion of free field in layered soil, and it can reduce greatly with the displacement amplitude of the free surface of the layered soil, and reduce slightly with the displacement amplitudes of the middle and lower nodes. Thus it can be seen that the influence of layered-soil damping on 1D time-domain method for the free field algorithm in layered half space by incident in-plane wave can not be ignored.

layered media; free field; 1D time-domain method; soil damping; time-stepping method

10.3969/j.issn.1672-1144.2016.06.004

2016-07-27

国家自然科学基金资助项目(51108088)

卓卫东(1966—)，男，福建莆田人，博士，教授，博导，主要从事桥梁抗震方面的研究工作。 E-mail: zhuowd@fzu.edu.cn

P315，TU311

A

1672—1144(2016)06—0018—07