SV 波超临界角斜入射时层状地基地震动输入在ABAQUS 中的实现

张 季，谭灿星，叶国涛，巴振宁，梁建文

(1. 华东交通大学土木建筑学院，江西，南昌 330013；2. 天津大学建筑工程学院，天津 300354)

近年来，有限元数值方法广泛应用于场地效应和地下结构地震响应研究中，有限元模型的地震动输入是有限元数值方法求解相关问题的关键。受地震波传播路径、传播介质复杂性的影响，工程场地和工程结构在实际地震中往往受到非垂直入射的地震波作用，对于大型或长距离空间结构来讲，非垂直入射的地震波所产生的空间非一致地震作用将对结构地震响应产生显著的影响[1-2]，合理的地震动输入须充分考虑地震波入射角度的影响。

SV 波是地震动的一种主要波动形式，SV 波斜入射时存在临界角情况，超临界角入射时，土层将产生不均匀波，近地表质点将以Rayleigh 面波形式运动，并伴随更大的竖向地震作用，因此，准确模拟SV 波超临界入射下的地震动输入，对揭示任意角度斜入射SV 波对工程结构地震响应的影响规律至关重要。

根据地震作用施加方式的不同，地震动输入分为振动方法和波动方法两类，振动方法是“仅考虑SV 波垂直向上入射情况的一种特例，而波动方法可以进一步考虑弹性地基和地震波斜入射”[3]。波动方法往往通过求解自由场地震响应来实现，如文献[4 - 6]提出的通过求解自由场响应获得与地震动输入相对应的等效节点力，然后将等效节点力施加到人工边界上从而使人工边界处应力和位移与自由波场对应位置相一致。自由场的求解有解析法和数值法 2 种，解析法只能求解弹性均匀半空间问题，对于层状地基情况还需借助各种数值方法。

对于层状地基斜入射平面P 波、SV 波情况，不少学者已提出了相应的自由场数值算法。李山有等[7]利用显式有限差分法的内节点位移计算公式和基于水平成层介质波动传播特点建立的相邻节点间自由场运动的关系式，给出了计算水平成层半空间自由场的时域方程。刘晶波等[8]提出了一种弹性水平成层半空间中平面内波动斜入射时自由波场时域计算的一维化有限元方法。尤红兵等[9]利用土层和半空间的精确动力刚度矩阵，提出了地震波斜入射时水平层状场地的非线性地震反应计算方法。赵密等[10-11]提出一种精确模拟基岩半空间辐射阻尼的人工边界条件，并建立了计算平面波斜入射下成层半空间自由场地反应的一维化有限元时域解法。王笃国等[12-13]基于等效线性化频域数值求解方法和二维应变空间状态理论，提出了一种地震波斜入射时二维水平成层介质的非线性地震反应分析方法。卓卫东等[14]在文献[8]算法基础上，采用有限差分方法，提出了P-SV 波斜入射时有阻尼成层介质自由波场的一维化时域算法。最近，张奎等[15]提出了水下地基场地的斜入射地震动输入时域方法，赵密等[16]提出了阶梯地形成层场地的斜入射地震动输入方法。

可以看出，自由场数值计算方法可分为时域方法[7-8,10-11,14-15]、频域方法[9,12-13]2 类。从在有限元方法中实现层状地基斜入射地震动输入而言，时域计算方法应用更为广泛，但是正如文献[17 - 18]指出的一样，时域方法尚不能模拟SV波超临界角入射情况。而对于频域方法的地震动输入，据笔者所知，目前也尚未见到SV 超临界角地震动输入在有限元数值分析中的应用(如频域方法文献[12 - 13]明确表示：“本文仅对 SV 波入射角小于临界角的情况进行研究”)。因此，SV 波超临界角入射的地震动输入问题对于有限元数值模拟而言是一个亟待解决的问题。

本文在Wolf 给出的弹性层状地基和半空间精确的动力刚度矩阵[19]基础上，结合傅里叶时频变换推导SV 波任意角度斜入射下的地震动输入等效节点力计算公式，重点探讨SV 波超临界角入射下频域刚度矩阵法求解层状地基地震动输入在ABAQUS 有限元软件中应用的有效性和准确性。同时，本文还将频域刚度矩阵法与等效线性化方法结合，通过引入二维应变空间状态理论来改进文献[9]基于一维理论建立的等效线性化方法，以使频域刚度法进一步适用于求解二维非线性层状地基的斜入射地震动输入。

1 频域刚度矩阵法

自由场模型如图1 所示，假定P 波、SV 波从基岩面以θ 角度入射，θ 为波阵面与水平x 轴的夹角。自由波场的求解思路如下：首先，采用层状地基和其下半空间精确的动力刚度矩阵建立自由场整体动力刚度矩阵；然后，由基岩入射来波确定基岩露头运动并以之计算施加于自由场的外荷载向量；接着，采用高斯消元法求解该荷载向量作用下自由场运动方程，即求得自由场响应；最后，结合傅里叶时频变换和等效线性化方法，进一步求解自由场时域响应解。本文将上述自由场求解办法称之为频域刚度矩阵法。

图 1 自由场模型示意图Fig.1 Schematic plots of free field model

1.1 土层和半空间动力刚度矩阵

频域内土层的动力刚度矩阵KL可表示为：

式中：方程左边项P、R 分别表示作用在土层界面的水平向和竖向荷载，下标1、2 分别表示荷载作用于土层上表面和土层下表面；i 为单位复数；刚度矩阵为KL对称矩阵，其系数kij参见文献[19]；u、w 分别表示土层界面处水平向和竖向位移，其下标1、2 分别表示土层上表面、下表面位移。

半空间的刚度矩阵为KR可表示为：

式中：方程左边项P0、R0分别表示作用在半空间表面的水平向和竖向荷载，刚度矩阵KR为对称矩阵，其系数rij参见文献[19]；u0、w0分别表示半空间表面水平向和竖向位移，即基岩露头运动。

1.2 荷载向量与自由场响应的求解

以ASV0和AP0作为入射SV 波和P 波的来波幅值，对于均匀半空间，该来波幅值作用下的半空间地表运动由下式确定：

集整土层刚度矩阵KL和其下半空间刚度矩阵KR得自由场整体刚度矩阵K，结合外荷载向量QT，通过高斯消元法即可求得自由场各土层界面的位移响应UT：

式(5)得到的位移是土层界面上的位移，土层内坐标为(x, z)处的位移由下式确定：

式中：cP和cS分别为土的压缩波速和剪切波速；ζ 为阻尼比。

记G、λ 为拉梅常数，土层内坐标为(x, z)处的应力为：

上述所得位移、应力均为频域解，故还需通过傅里叶逆变换将频域解转换至时域解。

1.3 等效线性化方法

等效线性化方法的原理，就是采用一组等效的剪切模量和阻尼比来代替不同应变幅值下的剪切模量和阻尼比，从而将非线性问题转变成线性问题求解[20]。等效线性化算法一般通过迭代方法进行，每一次迭代运算均是对地基进行一次线弹性地震反应计算，本文以频域刚度矩阵法建立线弹性地震反应程序，然后将频域刚度矩阵法与等效线性化相结合建立自由场的非线性地震反应程序，程序流程如图2 所示。由于自由场受斜入射P 波或SV 波作用，场地的受力状态实际上是二维平面应变状态，此时，土层内任意一点的最大剪应变γmax由该点的应变状态εx、εz、εxz确定：

根据等效线性化方法原理，等效剪应变γeff取为最大剪应变γmax时程峰值的0.65 倍：

2 粘弹性人工边界与地震动输入

2.1 粘弹性人工边界设置

图 2 等效线性化方法流程图Fig.2 Flow chart of equivalent linear analysis method

近年来，粘弹性人工边界由于其吸收散射波的有效性和在商用软件中的易于实现性，已广泛应用于场地和地下结构的地震反应分析研究中。在ABAQUS 有限元软件中，粘弹性人工边界可通过在计算模型的人工边界上施加并联的接地弹簧和阻尼器元件来实现。本文将粘弹性人工边界的弹簧、阻尼元件相关参数取为[21]：

式中：K、C 分别为弹簧、阻尼器元件的刚度系数和阻尼系数，下标bn、bt 分别为该元件施加于边界的法向和切向；ρ 为介质密度；r 为散射源至人工边界的距离；参数A、B 通过数值实验获得，本文根据文献[21]建议，取A=0.8，B=1.1。在粘弹性人工边界的施加过程中，还需要将边界节点的影响区域面积考虑到刚度系数和阻尼系数参数中。

2.2 地震动输入方法

采用刘晶波教授[4]提出的地震动输入方法，通过在人工边界节点上施加等效节点力的方式实现地震动输入，等效节点力应使人工边界上的位移和应力与原自由场相同。施加在有限元模型左、右和底边界的等效节点力分别为：

3 方法验证

3.1 模型1：弹性均匀半空间模型

首先，以平面 SV 波任意角度斜入射二维弹性均匀半空间为算例，通过频域刚度矩阵法计算有限元模型的地震动输入，来验证频域刚度矩阵法对任意角度斜入射SV 波地震动求解的适用性。半空间弹性模量E 为2 GPa，泊松比ν 为0.3333，密度ρ 为2000 kg/m3，此时SV 波的临界角为30°。选择图3 所示的脉冲波作为输入波，该脉冲波的解析表达式为：

其中，A=1，T=0.5，z(a)=a3H(a)，H(a)为Heaviside函数(当a<0 时，H(a)=0；当a≥0 时，H(a)=1)。

图 3 入射波Fig.3 Incident wave

采用ABAQUS 有限元建立宽1600 m、高600 m的均匀半空间有限元模型，模型如图4 所示。按照网格尺寸不大于剪切波最小波长的1/10 划分单元，本例中，计算截止频率取10 Hz，故单元网格尺寸按照6 m 长度设置。取散射波源为自由表面中心点(即B 点)，在模型左、右和底边界分别设置粘弹性人工边界。采用基于本文频域刚度矩阵法的自编Fortran 程序计算人工边界节点上的位移、速度和应力响应，然后在ABAQUS有限元软件中通过Python 程序读取这些响应值、计算等效节点力并施加于各侧人工边界节点上，最后提交Job 计算有限元模型的动力响应。

图 4 有限元模型示意图Fig.4 Schematic plots of finite element model

分别计算SV 波入射角为1°～90°时地表B 点的位移响应，提取各入射角度下B 点的水平向位移幅值ut和竖向位移幅值wt，将其与文献[19]结果进行对比，如图5 所示，可以看出，ABAQUS结果与文献[19]结果整体上非常吻合。必须指出的是，当入射角大于临界角30°时，ABAQUS 所得水平向位移略低于文献结果，究其原因，可能是由于在粘弹性人工边界上仍然产生了轻微的波反射现象，因为粘弹性人工边界理论是基于柱面波波动理论建立的，当SV 波超临界角入射时，近地表面将产生非均匀波，这种非均匀波作用下质点以面波形式运动，粘弹性人工边界对这种非均匀波的吸收效果可能不如对平面波的吸收效果好。

图 5 地表位移幅值随入射角的变化曲线Fig.5 Curves of surface displacement amplitude with different incident angle

图6 给出了SV 波入射角为15°、35°和60°时ABAQUS 有限元方法和频域刚度矩阵方法得到的地表B 点位移时程曲线，图中刚度矩阵法曲线由前述自编Fortran 程序所得，可以看出，两种方法所得时程曲线是非常吻合的。继续以图6 分析三个入射角下B 点的运动状态，当入射角为15°时，水平向位移和竖向位移同时到达最大值，说明二者是同相位的运动；当入射角为35°时，水平向位移落后于竖向位移；当入射角为60°时，竖向位移落后于水平向位移，这说明，当SV 波超临界角入射时，地表点的水平向位移和竖向位移存在相位差。实际上，正是由于这种相位差的存在，才使得SV 波超临界角入射时地表质点的运动呈椭圆型面波的运动形态。图7 进一步给出了三种入射角度下基于ABAQUS 结果的地表B 点运动轨迹图，可以看出，当SV 波15°入射时，质点的运动轨迹为一直线；当SV 波35°入射时，质点的运动轨迹近似为一顺时针旋转的椭圆；当SV 波60°入射时，质点的运动轨迹近似为一逆时针旋转的椭圆。这些现象表明，对于均匀半空间模型来讲，采用本文刚度矩阵方法计算有限元模型的地震动输入并进行有限元数值计算，准确再现了平面SV 波超临界角入射时地表质点的面波型运动形态。

图 6 入射角为15°、35°和60°时地表位移时程(模型1)Fig.6 Time histories of surface displacement when incident angle is 15°, 35° and 60° (for model 1)

3.2 模型2：弹性均匀半空间上单一土层模型

下面，以平面 SV 波斜入射作用下弹性均匀半空间上单一土层为算例，来进一步讨论本文方法的计算精度。场地参数如表1 所示，根据表中参数易知，土层中，剪切波波速与压缩波波速之比为1∶2，基岩中二者之比也为1∶2，并且土层压缩波波速和基岩剪切波波速相等。为了便于分析，设置地震波在土层底面以下300 m 处输入，建立1800 m 宽、600 m 高的有限元模型，仍以图3所示的脉冲波作为输入波，但计算时间截取至5 s。分别采用频域刚度矩阵法和有限元方法，计算入射角为15°、30°和60°时场地的位移响应，其中60°为超临界角入射情况(SV 波临界角为30°)。

图 7 入射角为15°、35°和60°时地表位移轨迹线Fig.7 Trajectory lines of surface displacement when incident angle is 15°, 35° and 60°

表 1 基岩上卧单一土层场地参数Table 1 Site profile of a single soil layer lying on bedrock

图8 给出了两种方法所得位移时程曲线的比较，可以看出，无论是入射角小于、等于或大于临界角，ABAQUS 有限元法所得位移与刚度矩阵法所得位移均吻合很好。图9 进一步给出了t=2.56 s时场地的位移幅值云图，云图反映了此刻土层和基岩中的波阵面走势，不难看出，土层中的波阵面比其下基岩的波阵面更为复杂，土层中的波阵面存在不同方向波阵面叠加情况，而基岩中的波阵面，除入射波波阵面(图中0-1 区段)外，其余波阵面均方向一致、彼此平行。

图 8 入射角为15°、30°和60°时地表位移时程(模型2)Fig.8 Time histories of surface displacement when incident angle is 15°, 35° and 60° (for model 2)

图 9 入射角为60°时场地位移幅值云图(t=2.56 s)Fig.9 Cloud picture of the magnitude of displacement when incident angle is 60° (t=2.56 s)

由于波阵面的形成由波在土层界面的反射和折射定律(即Snell 定律)所决定，因此，根据Snell定律可以从理论上推出任意时刻的波阵面图。图10、图11 为通过Snell 定律推导出的入射角为60°时t=2.56 s 时刻SV 波和P 波的波阵面示意图。以波阵面与水平线的夹角定义入射波、反射波、折射波的角度(注意波阵面与波的传播方向是相互垂直的)，对图中界面1 来说，此时SV 波从基岩进入土层中，Snell 定律可表示为：

式中：θ 为SV 波的入射角；θR、αR分别为基岩中SV 波、P 波的反射角；θL、αL分别为通过界面1 进入土层中的SV 波、P 波折射角；式中分母项为土层和基岩的压缩波或剪切波波速，其中，以下标S 表示剪切波速、下标P 表示压缩波速，以上标R 表示基岩、上标L 表示土层。

图 10 入射角为60°时SV 波波阵面示意图(t=2.56 s)Fig.10 Schematic diagram of wavefront of SV wave when incident angle is 60° (t=2.56 s)

图 11 入射角为60°时P 波波阵面示意图(t=2.56 s)Fig.11 Schematic diagram of wavefront of P wave when incident angle is 60° (t=2.56 s)

在界面1 上，此处θ=60°，根据式(13)知，入射SV 波将在界面1 上产生反射SV 波、折射P 波和 折 射SV 波，且θR=60°、θL=25.66°、αL=60°；由于折射P 波角度和入射SV 波角度相等，因此，折射P 波和入射SV 波在同一波阵面上(如图11所示)。波阵面的这一特征在图9 中的位移云图得到了准确的反映，图9 中入射SV 波波阵面0-1 与折射P 波波阵面1-2 均在一条直线上。

在界面2 上，此处入射波是P 波(见图11 中波阵面1-2)，P 波入射角度为αL=60°，根据式(13)知，在界面2 上将产生反射P 波(波阵面2-3)和反射SV 波(波阵面2-8)，P 波反射角为60°，SV 波反射角为25.66°。很明显，经过界面2 产生的这两个波阵面在图9 中也得到了准确的反映。

进一步分析界面3，此处入射波是由折射P 波1-2 在自由表面反射回来的P 波2-3，对应入射角度为60°，根据式(13)知，在界面3 上将产生反射P 波(波阵面3-4)、反射SV 波(波阵面3-10)和折射SV 波(波阵面3-16)，相应反射角、折射角分别为60°、25.66°和60°。不难发现，这三个波的波阵面在图9 位移云图中同样有相应的反映，但是由于波在多次反射折射中能量的分解，这三个波波阵面上的位移幅值均较小，因此位移云图中这三个波的波阵面显示比较微弱。

依次分析后续各界面，不难得出土层和基岩中各P 波和SV 波的波阵面走势(如图10、图11所示)。需要指出的是，土层中所有的P 波、SV波由土层进入基岩后，将只产生折射SV 波而无折射P 波，且折射SV 波的角度均为60°。特别地，对界面8 来说，此处既有入射P 波(波阵面6-8)，也有入射SV 波(波阵面2-8)，这两个入射波经过界面8 以后均产生角度为60°的折射SV 波，也就是说，波阵面8-17 是由2 个折射SV 波叠加而成。

由上述分析知，图9 位移云图所反映的波阵面图与Snell 波反射折射定理是完全相吻合的，因此，采用频域刚度矩阵法求解地震动输入能够准确模拟SV 波超临界角入射条件下层状地基的地震波场。

3.3 模型3：非线性层状地基模型

下面以非线性层状地基模型来验证本文基于刚度矩阵法和等效线性法的非线性地震动输入的计算精度。土层、基岩的弹性参数如表2 所示，对应的土层动模量曲线和阻尼比曲线如图12 所示，以加速度峰值调幅为0.1 g 的El Centro 波作为输入地震动，其加速度时程和傅里叶谱如图13所示。

首先，采用频域刚度矩阵方法和EERA 方法[22]计算SV 波垂直入射时场地的地表加速度响应，通过比较二者差别来验证频域刚度矩阵法等效线性化程序的正确性；然后将频域刚度矩阵法算得的地震动输入作用于ABAQUS 有限元模型，在ABAQUS中再计算有限元模型的地表加速度，将其与前两者结果进行对比，来验证有限元地震动输入的准确性；最后，再以SV 波超临界角入射为例，分别采用频域刚度矩阵法和ABAQUS 有限元方法分别计算地表加速度响应，对比二者响应结果，来验证超临界角入射条件下非线性层状地基的地震动输入的准确性。

表 2 土层参数(模型3)Table 2 Soil profiles for model 3

图 12 剪切模量比G/Gmax、阻尼比ζ 与剪应变γ 关系曲线Fig.12 Shear modulus ratio G/Gmax and damping ratio ζ curves versus shear strain γ

图 13 El Centro 波加速度时程及其傅里叶谱Fig.13 Accelerogram and Fourier spectrum of El Centro wave

具体地，通过以下方式，将基于刚度矩阵法和等效线性法的非线性地震动输入作用于ABAQUS有限元模型：

1) 采用自编的频域刚度矩阵法Fortran 程序建立一维土层等效线性化地震反应计算模型，进行等效线性化迭代运算，迭代运算完成后，基于有限元模型边界节点坐标信息，输出边界节点处土层等效剪切模量、等效阻尼比以及非线性地震响应信息(包括位移、速度和应力响应)；

2) 根据所得土层等效剪切模量计算自由场一阶自振频率，并根据式(12)和边界节点非线性地震响应信息，计算粘弹性边界上的等效节点力；

3) 以步骤1)中所得土层的等效剪切模量和等效阻尼比作为ABAQUS 有限元模型中的土层参数，并通过Python 代码将等效节点力批量化施加到模型左、右、底边上，然后在ABAQUS 软件中进行动力时程计算。

需要指出的是，由于ABAQUS 有限元阻尼理论与频域方法复阻尼理论不同，本文在ABAQUS中采用瑞雷阻尼考虑土层的阻尼比，瑞雷阻尼C=αM+βK，M 为质量矩阵，K 为有限元刚度矩阵，参数α、β 与阻尼比ζ 的关系如下式：

本文参照文献[23]的做法，将式中ω1取为步骤2) 中所得自由场一阶自振频率，ω2取为输入地震动的卓越频率。

图14 给出了SV 波垂直入射时(θ=0°)刚度矩阵方法、EERA 方法和ABAQUS 方法所得地表加速度响应的比较，可以看出，三种方法结果均较为吻合(其中，刚度矩阵法与EERA 的吻合度更高，是因为二者均是频域方法)。因此，本文根据土层动力刚度矩阵建立的等效线性化程序是准确的，且由此程序输出的非线性地震动输入也是准确的。

图15 给出了SV 波大于临界角入射时(θ=38.6°)刚度矩阵方法和ABAQUS 方法所得地表水平向、竖向加速度响应的比较，可以看出，两种方法结果也均吻合很好。这说明，刚度矩阵法所得的非线性地震动输入对于SV 波大于临界角情况来说，也是适用的。

图 14 刚度矩阵方法、 EERA 方法与ABAQUS 方法所得地表加速度响应比较(θ =0°)Fig.14 Comparisons of surface acceleration computed by stiff matrix method, EERA method and ABAQUS (θ=0°)

图 15 刚度矩阵方法、ABAQUS 方法所得地表加速度响应比较(θ=38.6°)Fig.15 Comparisons of surface acceleration computed by stiff matrix method and ABAQUS (θ=38.6°)

最后还需要指出的是，本文频域刚度矩阵方法之所以适用于SV 波超临界角入射情况，是因为该方法所采用的刚度矩阵是弹性层状土的精确动力刚度矩阵；此外，本文仅讨论了SV 波入射情况，对任意角度入射P 波情况，本文方法也同样适用，限于篇幅，不再给出对比结果。

4 结论

本文采用频域刚度矩阵法，推导了SV 波任意角度斜入射下的地震动输入等效节点力计算公式，通过ABAQUS 有限元软件分别模拟了SV 波斜入射下弹性均匀半空间、弹性基岩上卧单一土层和非线性层状地基的地震波场，探讨了基于频域刚度矩阵法的任意角度斜入射地震动输入应用于ABAQUS 有限元软件的有效性和准确性。结果表明，采用频域刚度矩阵法可在ABAQUS 中实现层状地基任意角度斜入射SV 波的地震动输入，且方法具有很高的计算精度。

本文方法可进一步应用于SV 波超临界角入射的场地效应研究和地下空间结构地震响应分析研究中。