基于FK-FE 混合方法的位错点源作用下全过程结构地震反应模拟

巴振宁，鲁世斌，付继赛，梁建文，芦 燕

(1.中国地震局地震工程综合模拟与城乡抗震韧性重点实验室，天津 300354；2.天津大学土木工程系，天津 300354)

准确可靠的地震动输入是计算分析工程结构地震反应的前提[1]。目前工程应用中大多采用地面运动预测方程(GMPEs)来获取地震动，但该方法基于遍历假设，使得震源参数、传播路径和场地条件对地震动时空分布的影响无法得到精确反应[2]。因此，建立一种基于物理的全过程分析方法(包含震源破裂、地震波场传播、土-结构相互作用和结构地震反应)，对精确结构地震反应分析具有重要意义。然而，由于全过程分析涉及地球物理、地震学、地震工程等多个学科，存在地壳层到工程结构的多尺度跨越问题，目前仍是一项难题。

近年来，诸多学者对全过程模拟进行了有益探索，并提出一些全过程模拟方法，按照建模方法主要可分为3 类：直接分析法、弱耦合两步法和强耦合两步法。其中，直接分析法是指将震源、传播路径和结构建立在同一模型中直接进行动力反应分析，例如：MAZZIERI 等[3]利用非连续伽辽金技术和谱元法，将区域尺度场地和工程结构建立在同一模型中进行分析，然而由于SEM中常用的显式时间推进的条件稳定性，无法广泛应用于梁、板和墙等构件的精细模拟；BRUN 等[4]利用混合异步时间积分器(HATI)，将SEM 和有限元(FEM)耦合在同一模型，并应用于大坝非线性动力分析，该方法可直接考虑从震源到结构的全过程反应，但由于数值方法的计算量随频率分辨率的16 倍增长，需消耗大量计算资源；弱耦合两步法是指先建立区域尺度模型求得地表地震动，然后直接在结构底部输入所得地震动进行分析，例如：KRISHNAN 等[5]利用谱元法(SEM)求得南加州地区在两次7.9 级设定地表作用下的地震波场，然后直接施加在一栋18 层钢框架结构底部进行动力分析，该方法操作流程较为简单，但未考虑土-结构相互作用的影响；强耦合两步法指将计算区域分为区域尺度模型(包含震源和传播路径)和局部区域模型(包含局部场地和结构)，首先对区域尺度模型进行求解，然后利用区域缩减法将所求得的地震波场施加在局部区域的边界，进而对局部区域进行动力分析，例如：ICHINMURA 等[6-7]、MCCALLEN 等[8]、ZHANG 等[9]和 张 磊 等[10]学者，分别采用摄动法、有限差分(FDM)、SEM、UCSB 等方法计算地震波场在区域尺度模型中的传播，采用FEM 计算局部区域动力反应，该方法基于区域缩减[11-13]方法，相对于整体模拟法更为便捷，且具有可考虑土-结构相互作用的优点，应用相对较广，但计算效率仍受限于区域尺度模型的最高识别频率和计算能力。

针对区域尺度模型中地震波场的传播问题，笔者前期建立了一种基于精确动力刚度矩阵[14-15]的频率波数域(FK)半解析方法[16-17]，该方法具有精度高、计算快和有效频带宽(可达0 Hz～20 Hz)等优点，可有效克服上述计算瓶颈。因此基于以上现状，本文以强耦合两步法为框架，将FK 方法和FE 方法相结合，建立了从震源到工程结构全过程地震反应分析的FK-FE 混合方法。同时，为便于使用，笔者进一步将FK 方法二次开发至大型商用有限元分析软件ABAQUS 中形成可视化插件，实现了仅在ABAQUS 平台下即可完成位错点源作用下全过程地震反应分析的功能，极大提高了建模和计算效率。

本文主要工作如下：首先，对FK 方法和FK-FE混合方法的理论进行阐述；然后，对FK-FE 混合方法在ABAQUS 软件中的主要操作步骤、二次开发流程和插件功能进行详细介绍；在验证该混合方法正确性的基础上，进一步以Benchmark 框架模型(美国SAC Steel Project 中的9 层抗弯钢框架结构)为例建立全过程分析模型，探讨了地壳层速度结构、和近地表场地条件对钢框架结构地震反应的影响。

1 FK-FE 混合方法

图1 给出了全过程地震反应求解分析的示意图，包括图1(a)所示的位错点源-传播路径-建筑结构整体模型和图1(b)所示的求解步骤。具体求解时将整体模型划分为STEP Ⅰ 所示的不含建筑结构的层状半无限空间模型和STEP Ⅲ 所示的有限元模型。

图1 全过程结构地震反应求解分析示意图Fig.1 Schematic of end-to-end structural seismic analysis

层状半无限空间模型自下而上包括半无限空间、地壳层及近地表覆盖土层，用于求解位错点源作用下层状半无限空间的地震波传播；有限元分析模型包括人工边界子结构、内部域土体和建筑结构，用于求解STEP Ⅰ 求得的地震波场输入下建筑结构的地震反应，该模型可考虑土-结构相互作用的影响。此外，在人工边界子结构外层节点设置黏弹性边界用于吸收有限元模型中建筑结构产生的外行散射波

本文建立FK-FE 混合方法的总体步骤如下：

STEP Ⅰ：采用FK 方法计算层状半无限空间(不含结构)人工边界子结构内外层节点位置位移时程。

STEP Ⅱ：将 STEP I 获得的位移时程输入人工边界子结构(杀死有限元模型内部单元)，采用有限元方法求得人工边界子结构外层节点等效节点力时程。

STEP Ⅲ：将 STEP II 获得的等效节点力时程施加在人工边界子结构外层节点，采用有限元方法对整个有限元模型(激活内部单元)精细分析求得建筑结构的地震反应。

以下分别对FK 方法模拟位错点源作用下地震波场和FK-FE 混合方法实现原理进行介绍。

1.1 FK 方法：位错点源作用下地震波场模拟

位错点源作用下场地某处位移响应表示为[18]：

式中：ui为x点沿i方向的位移分量；t为时间；∗为卷积运算符；g(x,t)为位错点源作用下动力格林函数；STF(t)为震源时间函数；Mpq为地震矩张量的六种位态分量(p,q=x,y,z)，在笛卡尔坐标系中为一对沿p方向作用相反、沿q方向分开的力偶。

地震矩张量位态分量Mpq可表示为与走向、倾角和滑动角相关的形式[19]：

式中：M0为位错点源地震矩； ϕ 、 δ 和 λ分别为位错点源的走向、倾角和滑动角。

位错点源格林函数的表达式为[17,20]：

式中： (Rkm,Skm,Tmk)为面谐基矢量；i 为虚数单位；k为水平波数；ω为圆频率；下标r、 θ和z为柱坐标系下x点位移分量方向；m为方位角模数。详细公式推导可参考文献[17，20]。

式(3)中位移幅值分量ui可表示为：

式中：Jm为第一类m阶贝塞尔函数；um和wm分别为平面内水平和竖向位移幅值；vm为平面外位移幅值；um、vm和wm具体表达式见文献[17]。

1.2 人工边界子结构法：FK-FE 混合方法实现

FK-FE 混合方法的建立主要包括两个关键技术：有限元模型地震波场输入和有限元模型边界处理。在地震波场输入方面，采用人工边界子结构法[21-22]输入FK 方法得到的地震波场；在边界处理方面，采用三维黏弹性动力人工边界[23-24]来吸收有限元模型中建筑结构产生的散射波。以下对人工边界子结构法和三维黏弹性动力人工边界进行介绍。

1.2.1 人工边界子结构法

参考文献[21]，输入地震波场定义为自由波场，且可以转化为在边界上施加的等效节点力，等效节点力的确定只和自由波场(图1 中STEP I)有关。因此，建立如图2 所示的有限元模型和相对应的自由场有限元模型。将模型按节点位置划分为人工边界子结构内、外层节点和剩余内部节点。

图2 人工边界子结构法求解示意Fig.2 Schematic for the substructure of artificial boundary method solution

按照节点位置，将两模型的运动方程以分块矩阵的形式表示：

式中：M为质量矩阵；C为阻尼矩阵；K为刚度矩阵；u、u˙ 和u¨分别为位移、速度和加速度；Peff为施加在人工边界子结构外层节点的等效节点力；下标i、b 和e 分别代表内部域节点、人工边界子结构内层节点和外层节点，上标“0”表示自由场。

当自由场有限元模型的人工边界子结构单元网格划分与实际有限元模型一致时：

结合式(5)～式(7)可得等效节点力：

根据式(8)即可得到最终施加于人工边界子结构外层节点的等效节点力的表达式：

具体推导过程可见文献[21 - 22]。由式(9)可知，等效节点力的求解只与人工边界子结构内外层节点的运动有关，而与有限元模型内部域节点无关。因此，采用1.1 节中描述的FK 方法得到层状半无限空间中人工边界子结构位置处的节点(内外层节点)动力响应后，即可通过人工边界子结构法获得相应等效节点力(外层节点)，进而完成内部域的FE 分析。

值得注意的是，运用式(9)计算等效节点力时需要获得人工边界子结构外层节点的位移、速度和加速度，但在实际应用中有限元软件可自动求解速度和加速度，因此只需给定节点位移即可求得等效节点力。

1.2.2 三维黏弹性动力人工边界

人工边界子结构外层节点法向与切向弹簧和阻尼器参数参考式(10)和式(11)：

式中：CN和CT分别为阻尼器法向和切向阻尼系数；KN和KT分别为弹簧法向和切向刚度；G为剪切模量；ρ 为质量密度；cP为压缩波速；cS为剪切波速；R为散射波源到人工边界的垂直距离；αN和αT分别表示法向和切向弹簧的修正系数，参考文献[25]分别取为1.34 和0.67。

2 FK-FE 混合方法在ABAQUS 中的二次开发

为实现FK-FE 混合方法在工程实际问题中的可视化应用，通过ABAQUS 提供的二次开发接口将上述工作开发为具有可视化图形界面的程序插件。以下对FK-FE 混合方法在软件中的主要步骤及插件二次开发流程简要说明。

2.1 FK-FE 混合方法在ABAQUS 软件实现主要步骤

FK-FE 混合方法在ABAQUS 软件中的实现主要步骤如图3 所示，具体操作流程如下：

图3 FK-FE 混合方法在软件中的主要实现步骤Fig.3 Implement steps of the FK-FE hybrid method in software

Step 1：建立土-结构相互作用有限元模型；

Step 2：采用生死单元法，杀死土-结构相互作用有限元模型内部单元，建立人工边界子结构；

Step 3：利用Python 脚本，提取人工边界子结构所有节点位置坐标；

Step 4：将位置坐标传递给FK 方法自编Fortran 程序，完成相应节点位移时程计算；

Step 5：利用Python 脚本计算三维黏弹性动力人工边界条件参数，并赋予人工边界子结构外层节点；

Step 6：修改.inp 文件，输入人工边界子结构两层节点位移时程，计算得到等效节点力Peff；

Step 7：修改.inp 文件，施加等效节点力Peff于人工边界子结构外层节点；

Step 8：激活内部单元，完成位错点源作用下土-结构相互作用有限元分析，求得建筑结构地震响应。

2.2 二次开发流程

本插件开发主要分为FK 方法Fortran 程序在ABAQUS 中的接入和ABAQUS-GUI 程序代码编写 两 部 分。 其 中， FK 方 法Fortran 程 序 在ABAQUS 中的接入实现了层状半无限空间地震波场在有限元模型中的输入，ABAQUS-GUI 程序代码实现了FK-FE 混合方法在ABAQUS 中的手动操作到内核自动执行程序和可视化界面的转化。

需要注意的是，由于不同版本ABAQUS 所对应的内置Python 解释器不同，进行编译和相应numpy 库下载时，需处理好版本兼容问题。本文所用ABAQUS 版本为6.14-4，Python 解释器版本为2.7，NumPy 库版本为NumPy 1.11.2。

二次开发设计思路如图4 所示，具体流程如下：

图4 二次开发设计思路Fig.4 Design ideas of secondary development

Step 1：将FK 方法Fortran 子程序写为.F90 格式，其中.F90 代码中Module 封装的变量用Common重新定义为全局变量，合并修改后的.F90 子程序，封装为子函数(Subroutine)；

Step 2：利用Python 语言开源库中的f2py.exe工具(f2py.exe 可通过下载与Python 对应版本的numpy 库得到，是一个Fortran 程序到Python 程序的接口转换工具)以及GFortran 和mingw32 编译器对.F90 程序进行编译，形成.pyd 文件(编译指令为“f2py -m Dislocation_Source -c Dislocation_Source _Abaqus.f90 –fcompiler = gFortran –compiler=mingw32”，封装的.pyd 模块名称为Dislocation_Source)；

Step 3：在内核执行程序中写入.pyd 模块的调用命令，关联FK 方法Fortran 程序与二次开发插件(调用命令为模块导入指令“import Dislocation_Source”和子函数调用指令“Dislocation _Source.Sub”，模块内子函数名称为Sub)；

Step 4：利用ABAQUS 软件自带的RSG 构造器进行插件程序设计，主要内容包括设计交互界面的“GUI”(图形界面文件)和负责关联脚本函数的“Kernel”(注册文件)两部分。

2.3 插件简介

开发完成的插件由图形界面文件、注册文件、内核执行程序、.pyd 文件以及.txt 文件组成。各组件功能说明见表1，插件界面展示如图5～图7 所示。

表1 二次开发过程各组件功能Table 1 The function of each component in the secondary development process

图5 建立人工边界子结构参数设置界面Fig.5 Parameter setting interface for building substructure of artificial boundary

图6 FK 法输入波场参数设置界面Fig.6 Parameter setting interface of calculating the input waves field with the FK method

图7 施加等效节点力参数设置界面Fig.7 Parameter setting interface of applying equivalent nodal force

开发完成后，将插件放入ABAQUS 的起始工作路径，即可在 “Plugin-ins”菜单中找到相应的插件启动按钮。插件控件功能说明详见表2。

表2 FK-FE 混合方法控件功能说明Table 2 Control function description of FK-FE hybrid method

3 方法验证

本节分别采用FK 方法和FK-FE 混合方法计算同一层状半空间模型在位错点源作用下的地震反应来验证FK-FE 混合方法的精度，采用FK-FE混合方法模拟2021 年漾濞6.4 级地震动，并与强震记录对比来验证方法的有效性。

3.1 FK-FE 混合方法精度验证

验证模型如图8 所示，场地的速度结构参数见表3[25]。内部域有限元模型尺寸为400 m×400 m×400 m，网格大小10 m。设置3 个观测台站位置坐标分别为A(0, 0, 0)、B(0, 100, 200)和C(100, 100, 300)。位错点源埋置在地表下4 km 处，地震矩M0=1022dyne.cm，走向、倾角、滑动角分别为135°、70°和30°。震源时间函数为Ricker 子波，设定中心频率f0为5 Hz，波峰位置为0.3 s，计算时间间隔为0.025 s，Ricker 子波时域和频域波形如图9 所示。

表3 验证模型场地速度结构[26]Table 3 Site velocity structure of validate model

图8 FK-FE 混合方法验证模型Fig.8 Verification analysis model of the FK-FE hybrid method

图9 震源时间函数Fig.9 Source time function

图10～图11 分别给出了采用两种方法得到的观测台站处的位移、加速度时程的相对误差分布，两种结果十分吻合，位移时程误差均在±0.5%以内，加速度时程误差均在±1%以内，检验了本文所建立的FK-FE 混合方法精度。

图10 FK 方法和FK-FE 混合方法位移时程对比Fig.10 Comparison of displacement time history for the FK method and the FK-FE hybrid method

图11 FK 方法和FK-FE 混合方法加速度时程对比Fig.11 Comparison of acceleration time history for the FK method and the FK-FE hybrid method

3.2 FK-FE 混合方法有效性验证

有效性验证模型与精度验证模型类似，以观测台站为原点建立内部域有限元模型，模型尺寸为100 m×100 m×100 m，网格大小2 m。参考漾濞6.4 级地震相关信息[27]设定震源参数，位错点源埋深7.5 km，地震矩M0=4.47×1025dyne.cm，走向、倾角、滑动角分别为135°、82°和-165°。震源时间函数取钟形函数(式12)，τ表示上升时间，取τ=0.5 s，计算时间间隔为0.025 s。漾濞地区一维地壳波速结构可参考文献[28]。

图12 给出了FK-FE 混合方法模拟的2021 年6.4 级漾濞地震中2 个远场台站处地震动加速度时程与实测记录的对比结果。53SDX 为施甸台站，震中距127.88 km，53YRH 为仁和镇台站，震中距142.33 km。计算结果均为50 s 的加速度时程，曲线末端所标数值为地震时地面最大峰值加速度PGA，单位为cm/s2，黑色与红色时程曲线分别为强震记录和模拟结果。比较发现，本文模拟的地震动加速度时程在幅值、持时、波形上与强震记录对比良好，论证了FK-FE 混合方法的有效性。

图12 模拟地震动时程与实测记录对比Fig.12 The time-history comparison between the simulations and the measured records

4 算例分析——以9 层Benchmark钢框架模型为分析对象

本节将建立的FK-FE 混合方法应用于美国SAC Steel Project 的 9 层Benchmark 抗弯钢框架的全过程地震反应分析，探讨了地壳层速度结构与场地类别两个因素对该结构地震反应的影响。该Benchmark 钢框架已被国内外多位学者使用[29-30]，具有一定代表性与普适性，模型具体资料可参考文献[29]。

4.1 层状半无限空间模型

层状半无限空间模型场地速度结构参数见表4[31]。位错点源埋置在地表下18 km 处，全局坐标为(20 km, 20 km, 18 km)。震源参数设置如下：走向、倾角和滑动角分别为135°、70°和30°，地震矩M0=3.55×1026dyn.cm。震源时间函数同验证模型，中心频率f0为5 Hz，截止频率fmax约为15 Hz，波形如图9 所示。

表4 层状半无限空间模型场地速度结构[31]Table 4 Site velocity structure of layered infinite half-space model

4.2 有限元模型

土结构相互作用会显著影响框架结构的地震反应[32]。因此，本文运用ABAQUS 建立内部域土-结构相互作用有限元模型，如图13 所示。

图13 土-结构相互作用有限元模型 /mFig.13 Finite element model of soil-structure interaction

模型布置方面，内部域土体模型尺寸为200 m×200 m×23.65 m，网格大小为3 m～5 m，满足15 Hz波场计算精度；土体模型顶面为自由面，其余边界面设置黏弹性边界；上部Benchmark 钢框架外层灰色一榀框架采用连接器以MPC-pin 的形式铰接于主体结构；框架柱底和周围场地土进行绑定约束。

单元类型选取方面，场地土体采用三维实体单元(C3D8)，Benchmark 钢框结构采用梁单元(B31)。

材料选取方面，场地土为弹性，其各项参数根据实际工况中场地速度结构进行设置。钢材采用双线性等向强化本构模型，塑性阶段(屈服后)模量取为1/100 倍的弹性阶段弹性模量。有限元模型采用Raiyleigh 阻尼，本模型取场地和结构的第一阶和第二阶振型自振频率计算Raiyleigh 阻尼系数α 和β，其中土层阻尼比见表4，钢材阻尼比取0.03。

4.3 结构模型合理性分析

对9 层Benchmark 钢框架有限元模型(不含场地)进行模态分析，如表5 所示，将所得前四阶振型自振频率与文献[29]结果对比，自振频率误差均控制在±5%以内，验证了结构模型的合理性。

表5 钢框架模型自振频率与文献[29]结果对比Table 5 Results comparison of this paper and literature [29]for Natural frequency of steel frame model

4.4 地壳层速度结构对建筑结构地震反应的影响

地震波地壳层传播过程中受到行波效应以及非均匀地形的影响，会发生反射和折射，导致幅值和方向发生变化[33]。为探讨地壳速度结构对建筑结构地震反应的影响，以表4 速度结构为基准(表6 中速度结构1)，以莫霍面压缩波速(8.10 km/s)为上限，控制地壳层介质泊松比与密度不变，递增地壳各层剪切波速，得到表6 所示三种地壳层介质。本文重点关注了结构一侧一榀框架的动力响应，各响应输出点位置示意图14 如所示。

表6 三种地壳层速度结构Table 6 Three crust velocity structures

图14 结果展示位置示意图Fig.14 Schematic plots of position for result display

图15～图16 分别给出了FK-FE 自由场模型中地表P点三分量加速度时程和反应谱结果。为便于观察，将反应谱中所有y向分量缩小5 倍、z向分量缩小25 倍。可以发现，地壳层速度结构软硬程度的变化对地震动加速度响应的影响较为明显，其中尤以水平向分量最为显著。随地壳层介质变硬，地震动响应出现时刻提前，幅值减小，持时更短，衰减更加迅速。且反应谱以短周期成分为主，幅值随介质变硬而变小，但包含的频率成分基本不变。例如：相较于速度结构1，速度结构3 工况下地表P点三分量加速度响应分别减小约72.6%、86.3%和43.8%；衰减为峰值幅值1/10所用时间分别减小2.30 s、1.95 s 和1.15 s；反应谱幅值减小约52.1%、45.0%和36.6%。

图15 自由场模型中地表P点加速度时程Fig.15 Acceleration time history of pointPon the surface for the model of free field

图16 自由场模型中地表P点加速度三分量加速度反应谱Fig.16 Three-component acceleration response-spectrum of pointPon the surface for the model of free field

图17 给出了FK-FE 土-结构相互作用模型中地表P点三分量加速度时程。比较土-结构相互作用模型与自由场模型地表加速度时程发现，考虑土-结构相互作用的影响后，地表加速度响应在三种地壳层速度结构下均有不同程度的放大。例如：速度结构1 中，地表三分量加速度分别放大1.16 倍、1.65 倍和2.19 倍，速度结构2 为1.18 倍、2.23 倍和2.54 倍，速度结构3 则是1.4 倍、2.18 倍和2.61 倍。由此可见，土-结构相互作用对较硬地壳介质下结构地震反应的影响更为显著。

图17 土-结构相互作用模型中地表P点加速度时程Fig.17 Acceleration time history of pointPon the surface for the model of soil-structure interaction

图18 分别给出了不含结构的FK-FE 自由场模型和含有结构的FK-FE 土-结构相互作用模型在不同地壳层速度结构下地表P点的三分量加速度反应谱，对应阻尼比为5%。

图18 两种模型地表P点三分量加速度反应谱Fig.18 Three-component acceleration response-spectrum of pointPon the surface for the two model

在土-结构相互作用的模型中，反应谱幅值同样会随地壳介质变硬而变小，这一现象在短周期成分较为明显，而对长周期成分影响较小。值得注意的是，考虑土-结构相互作用的影响后，加速度反应谱峰值对应周期增大，反应谱低频成分增多，且在土-结构相互作用模型中，峰值对应周期会随地壳介质的变硬而增大。相较速度结构1，速度结构3 工况下三分量加速度反应谱峰值对应周期分别提高了0.1 s、0.075 s 和0.05 s。

图19 给出了框架角柱柱顶处(图14)在三种地壳层速度结构条件下x与y方向的位移时程；图20给出了结构角柱柱顶产生最大位移时三种地壳层速度结构下钢框架变形云图，其中图20(a)～图20(c)为一榀框架(图14 所示位置框架)x向变形云图，图20(d)～图20(f)为框架整体变形云图。

图19 不同地壳层速度结构下顶层角柱位移时程Fig.19 Displacement time history of top corner pillar in different crust velocity structures

结果显示，随着地壳层剪切波速和压缩波速的增大，位移峰值点出现时刻逐渐前移，且峰值大小逐渐减小。例如：三种速度结构x向位移峰值分别出现在11.875 s、10.375 s 和9.25 s，最大相对位移分别为21.6 cm、19.8 cm 和16.7 cm，和地壳层速度结构1 相比，速度结构3 工况下建筑结构x向最大位移减小了22.7%；y向位移峰值分别出现在11.875 s、10.375 s 和9.275 s，最大相对位移分别为21.1 cm、19.4 cm 和18.1 cm，和地壳层速度结构1 相比，速度结构3 建筑结构y向最大位移减小了14.3%。

上述现象是由于在较硬的地壳层介质中地震波传播速度较快，P 波初至时刻提前，且不同地壳层速度结构对地震波放大作用也不同导致的。

此外，不同地壳层速度结构下，顶层角柱峰值位移时刻结构最大形变分布沿层高也有所不同。对于本算例所示的一榀框架，速度结构1 工况下结构最大形变出现在第3 层～4 层的位置；速度结构2 中出现在第4 层～5 层及第8 层～9 层位置；速度结构3 则在第2 层～3 层及第9 层位置。

图21 给出了三种地壳层速度结构下结构层间位移角的分布情况。由图可见，随着地壳层剪切波速和压缩波速的增大，结构各层层间位移角逐渐减小，相较速度结构1 的位移响应，结构在地壳层速度结构3 场地条件下x向层间位移角减小为约30.5%～61.3%，y向层间位移角减小为约26.9%～62.0%。可见，地壳层速度结构的变化会显著影响结构层间位移角的大小。

图21 不同地壳层速度结构下最大层间位移角Fig.21 The maximum inter-story drift angle under different cases of crust velocity structures

以上数据表明，在抗震设计中，应考虑所在地区地壳层速度结构对建筑地震反应的影响，并对建筑结构进行合理规划。

4.5 不同场地类别对结构地震反应的影响

依据我国《建筑抗震设计规范》，进行建筑抗震设计时应考虑建筑场地对建筑抗震利害的影响，建筑场地类别的划分是较为重要的参考因素[33]。因此，本文以土层等效剪切波速为设计依据，设计了Ⅱ、Ⅲ、Ⅳ三种场地类别探究其对结构地震反应的影响。以表4 场地速度结构为基础，调整顶部覆盖层厚度为100 m，覆盖层参数见表7。震源参数设置同4.1 节。

表7 场地覆盖层介质物理参数Table 7 Physical parameter of site overburden layers

图22～图23 给出了钢框架顶层角柱上端在不同场地类别下x与y方向的位移时程结果以及层间位移角沿高度的分布情况。可以发现，当建筑结构位于较硬的场地覆盖层介质上时，产生的动力响应较小，随着场地覆盖层介质剪切波速和压缩波速的减小，结构的整体位移响应及结构各层层间位移角均逐渐增大。例如：Ⅱ类场地相比于Ⅳ类场地，x向最大位移由24.9 cm 增加至30.1 cm，幅值变化约17.2%；y向最大位移由25.6 cm 增加至32.9 cm，幅值变化约22.3%；x向最大层间位移角约增大16.9%～44.7%，y向最大层间位移角约增大20.8%～30.0%。

图22 不同场地类别下顶层角柱位移时程Fig.22 Displacement time history of top corner pillar under different cases of different site classes

图23 不同场地类别下最大层间位移角Fig.23 The maximum inter-story drift angle under different cases of different site classes

由此可见，结构在位错点源地震作用下，场地覆盖层介质的软硬程度会显著影响结构的地震反应。此外，观察位移时程可发现，场地类别对位移峰值出现时刻基本无影响。这是由于，虽然近地表覆盖层介质不同，但覆盖层厚度较薄(100 m)，因此对P 波初至时刻影响不大。

5 结论

本文基于区域缩减理论提出了一种FK-FE 混合方法。该方法采用FK 法计算位错点源作用下整体模型地震波场，运用FE 法模拟有限元分析域中土-结构相互作用，实现了从位错点源到结构的全过程物理模型分析。以9 层Benchmark 抗弯钢框架模型为应用实例，研究了地壳层速度结构与场地类别对结构地震反应的影响。得到主要结论如下：

(1) 本文建立的FK-FE 混合方法可模拟震源到建筑结构的全过程地震反应，能够有效解决宽频地震动作用下地壳层到工程结构的多尺度跨越问题，具有较高的求解精度。

(2) 所形成的FK-FE 混合方法ABAQUS 插件具有操作简便、应用性强的特点，可综合考虑震源和场地参数信息。该插件不仅可以对本文Benchmark抗弯钢框架进行全过程地震反应分析求解，对其他各种结构同样适用，例如：高层结构、大跨空间结构、地下结构和桥梁结构等。

(3) 位错点源作用下，地壳层速度结构和场地类别对上部结构的变形和层间位移角有显著影响。本文算例中，地壳层速度结构3(较硬)相比速度结构1(较软)，框架最大位移和层间位移角分别减小了22.7%和62.0%；Ⅱ类场地相较于Ⅳ类场地，框架最大位移和层间位移角分别减小了22.3%和44.7%。