在近断层区域，桥梁、大坝、隧道等大型建/构筑物承受地震动惯性力和断层错动的耦合破坏作用下，大型建/构筑物与周围土体的动力相互作用对结构本身的地震反应有显著影响。尤其是对于结构与围岩相互作用的地下结构而言，地层的错动变形在地震破坏效应中往往起着关键性作用^{[1 − 3]}。因此，考虑土-结构相互作用的抗震分析方法研究具有重要意义^{[4 − 5]}。

数值模拟是分析土-结构动力相互作用问题的主要方法^{[6 − 10]}。基于数值模拟方法进行土-结构相互作用的地震反应分析时，关键在于对半无限介质中地震波传播的模拟。为了提升计算效率，需要从半无限土体中截取有限尺寸的计算模型，并利用人工边界技术模拟有限尺寸模型的地震波动。近场波动问题中常用的人工边界可分为全局人工边界和局部人工边界两类^{[11 − 12]}。全局人工边界可以精确模拟半无限地基的全部数理条件，但是由于全局人工边界的频域求解方法复杂，难以应用于大型结构动力分析。局部人工边界是一种近似的人工边界方法，因为数值解法简单且具有时空解耦性而应用广泛。粘弹性边界是一种应力型局部人工边界^{[13 − 14]}，因为具有理论明确、应用方便、没有边界失稳问题等优点，得到了广泛的应用和发展。目前，杜修力等^[12]与刘晶波等^[15]提出的两种粘弹性边界形式获得了较高的认可度。完美匹配层也是一种局部人工边界，因其具有较好的吸收任意入射角和频率的出射波并且在边界处没有任何的反射^{[16 − 19]}，在地震波传播的数值模拟方法上得到了广泛的应用。

结构动力响应的准确模拟取决于地震作用的选择和作用形式。在传统的土-结构动力相互作用的抗震分析方法中，通常将地震动简化为平面波并直接作用在岩土地基上，并通过计算入射波的时间延迟模拟近场波动^[20]。该类方法的局限性在于无法模拟震源效应和断层位错产生的永久地层变形，而且将地表记录到的地震动简化为地底平面波的方式依旧存在争议。赵伯明等^[21]根据位错理论^{[22 − 23]}和位移人工边界提出地下工程的抗位错分析方法，但该方法基于静态位错理论，无法模拟近场波动。因此，只有明确“震源-路径-结构”的地震波传播的完整动力过程，才能合理地确定地震的荷载形式。

基于震源的地震动模拟是地震学研究的基本工具，其研究思路是通过定义震源模型和介质模型来模拟从震源到路径的地震波传播过程。在模拟地震波传播的各种数值方法中，交错网格有限差分法因为数学概念直观和并行计算效率高而得到广泛的应用。交错网格可以在同一网格的不同节点上计算应力和速度分量，从而提升了计算的效率和精度^[24]。有限差分法(Finite Difference Method，FDM)的原理是利用差商代替偏导数近似地获得差分方程，然后通过求解相应的差分方程获得波动方程的近似解^[25]。MADARIAGA^[26]首先提出在时间和空间上具有二阶精度的交错网格FDM模拟地震地面运动时程，LEVANDER^[27]随后将其空间上的精度扩展到四阶。目前，已经有多种基于交错网格FDM模拟地震波传播的算法，例如GMS^[28]、FDMPI^[29]、FDSim^[30]和OpenSWPC^[31]等。

基于震源的地震动模拟已经广泛应用于震源反演、强地面运动预测和地震灾害预测等，但其应用于土-结构动力相互作用的抗震研究还比较匮乏。本文结合粘弹性人工边界和地震波传播的数值模拟算法，提出一种基于有限断层震源的土-结构动力相互作用抗震分析方法。在此基础上，本文介绍了该方法的实现步骤，并通过自由场模型算例和土-结构相互作用模型算例验证该方法的有效性。

1   粘弹性边界

粘弹性边界是一种局部应力型边界条件。根据有限尺寸模型外域的本构关系，该人工边界将计算域的单侧外源波表示为边界上的应力状态。粘弹性边界是时空解耦的，该人工边界可以认为是截断边界节点上的弹簧-阻尼器系统，如图1所示。粘弹性边界通过将自由场运动转化为边界节点上的等效节点力来实现地震波的输入，人工边界面上l节点i方向的等效节点力${f_{li}}$的计算公式为^{[32 − 33]}：

式中：$K$和$C$分别为对应l节点i方向的弹簧刚度和阻尼系数；${A_l}$为l节点的控制面积(见图1)；$u_{li}^{\text{f}}$、$v_{li}^{\text{f}}$和$\sigma _{li}^{\text{f}}$分别为截断边界处l节点i方向的自由场位移、速度和应力。

图  1  粘弹性边界的示意图

Figure  1.  Schematic diagram of viscoelastic boundary

下载: 全尺寸图片 幻灯片

本文选用杜修力等^[10]提出的粘弹性边界参数取值方法进行计算，三维粘弹性边界的弹簧-阻尼参数为：

式中：$K_{\mathrm{N}}$和$C_{\mathrm{N}}$分别为法向弹簧刚度系数和法向阻尼系数；$K\mathrm{_T}$和$C\mathrm{_T}$分别为切向弹簧刚度系数和切向阻尼系数；$\lambda$为拉梅常数；$G$为剪切模量；$R$为有限尺寸模型的几何中心到该节点所在边界面的距离；$\rho$为介质的密度；${c_{\text{p}}}$和${c_{\text{s}}}$分别为纵波波速和横波波速。

因此，只需要获得截断边界处所有节点的位移、速度和应力，便可通过粘弹性边界模拟有限尺寸模型的地震波动。

2   基于震源的地震动模拟

本文采用了空间上具有四阶精度、时间上具有二阶精度的交错网格有限差分法模拟断层破裂引起的地震动和同震变形。该算法基于广义Zener体(Generalized Zener Body, GZB)^[34]，并采用理想匹配层(Perfectly Matching Layer, PML)^[19]作为模型的吸收边界。本节仅介绍算法中的关键理论和相关参数，例如速度模型参数和震源参数。关于交错网格有限差分法的原理和粘弹性介质理论的细节可以参考MOCZO等^[30]、IGEL^[35]和ROBERTSSON等^[36]的研究。

在三维笛卡尔坐标系下，由加速度和应力描述的运动方程如下^[24]：

式中：$\rho$为密度；${\ddot u_p}$、${\sigma _{pq}}$和${f_p}$分别为质点加速度、应力张量和体力($p,{\text{ }}q = x,{\text{ }}y,{\text{ }}{\textit{z}}$表示对应方向分量)。各向同性的弹性介质中的应力分量由式(4)计算^[24]：

式中：$\lambda$和$\mu$为拉梅常数，在GZB粘弹性体中可以由多个不同的松弛时间和蠕变时间表示^[34]，并通过假设所有Zener体的松弛时间和蠕变时间之比为常数，实现P波品质因子(${Q_{\text{p}}}$)和S波品质因子(${Q_{\text{s}}}$)在地震频率范围(小于1 Hz)内近似为常数；${\delta _{pq}}$为Kronecker函数；$e$为应变，并根据式(5)计算^[24]：

式中，$u$为质点位移。

式(3)~式(5)可以使用交错网格有限差分法进行显式求解^[24]。在地表采用零应力边界条件，其余边界采用理想匹配层作为模型的吸收边界。

在基于震源的地震动模拟算法中，地壳速度模型可以是均匀、分层均匀和非均匀等形式。而地壳速度模型需要给定的关键参数如下：密度($\rho$)、纵波波速(${c_{\text{p}}}$)、横波波速(${c_{\text{s}}}$)、P波品质因子(${Q_{\text{p}}}$)和S波品质因子(${Q_{\text{s}}}$)。

有限断层震源可以详细描述断层尺度、产状、破裂过程和断层面位错的不均匀性。有限断层震源由断层面上的一系列等效体力$f$模拟，当断层破裂经过某个点源时，该点源都会辐射出假定的震源时间函数($\dot M(t)$)的地震波。断层的震源机制可以由双力偶源的三个参数(走向角${\phi _{\text{s}}}$、倾角${\delta _{\text{d}}}$和滑移角${\lambda _{\text{r}}}$)定义，如图2所示；震源的位置根据其相对坐标定义；地震的震级大小通过断层滑动量(D)和断层面积(S)来描述。断层的破裂过程可以由震源时间函数($\dot M(t)$)、起始时间(${T_0}$)和上升时间(${T_{\text{r}}}$)描述。

图  2  有限断层模型的几何示意图

Figure  2.  Geometric diagram of finite-fault model

下载: 全尺寸图片 幻灯片

3   分析流程

图3展示了基于有限断层震源的土-结构相互作用抗震分析方法的计算模型。首先截取包含结构和地基在内的有限元模型，并通过粘弹性边界模拟有限尺寸模型的地震波动。然后。根据给定的震源参数和速度模型参数模拟地震波在地壳中传播，并得到粘弹性边界上所有节点的应力、位移和速度。

基于有限断层震源的土-结构相互作用抗震分析流程如图4所示。主要分为以下几个步骤：第一步，建立土-结构相互作用的有限元数值模型，并提取模型截断边界上所有节点的坐标；第二步，根据给定的震源参数和速度模型参数，基于地震动模拟算法获得截断边界上所有节点的应力、位移和速度；第三步，在有限元数值模型的边界上设置粘弹性人工边界，根据节点的应力、位移和速度计算等效节点力，从而模拟有限元模型的地震波动。

图  3  抗震分析模型

Figure  3.  Seismic analysis model

下载: 全尺寸图片 幻灯片

图  4  分析流程图

Figure  4.  Analysis procedure diagram

下载: 全尺寸图片 幻灯片

4   自由场模型算例

为验证基于有限断层震源的土-结构相互作用抗震分析方法的计算精度，本节选取均匀弹性自由场的有限元数值模型算例进行计算分析。图5展示了在地壳速度模型中建立的有限断层震源模型，均质速度模型的参数为：尺寸(长×宽×高)为100 km×100 km×60 km；密度$\rho = 2700{\text{ kg/}}{{\text{m}}^{\text{3}}}$；纵波波速${c_{\text{p}}} = 6000{\text{ m/s}}$；横波波速${c_{\text{s}}} = 3500{\text{ m/s}}$；P波品质因子${Q_{\text{p}}} = 200$；S波品质因子${Q_{\text{s}}} = 200$。

图  5  速度模型中的有限断层震源和有限元模型

Figure  5.  Finite-fault source and finite element model in velocity model

下载: 全尺寸图片 幻灯片

有限元模型的参数为：尺寸(长×宽×高)为20 km×20 km×3 km，空间网格尺寸为0.25 km；有限元模型几何中心在地壳速度模型中的坐标为(50 km, 50 km, 1.5 km)；模型顶部为自由地表，并在其余五个面(共10 401个节点)上设置粘弹性人工边界。如图5所示，地壳速度模型中z方向向下为正，为了方便计算，我们对粘弹性边界上节点的应力、位移和速度进行坐标变换，使得有限元模型中z方向向上为正。

基于地震矩理论^[37]和断层破裂参数公式^[38]，在逆断层作用下，对应的矩震级为${M_{\text{w}}} = 7.0$的典型震源参数为^[31]：尺寸(长×宽)为50 km×15 km；断层埋深5 km，断层面几何中心的水平投影坐标为(50 km, 50 km)；平均滑移量为2.0 m(单凹凸体位错分布模式)；走向角为${0^\circ }$；倾角为${45^\circ }$；滑移角为${90^\circ }$；有限断层震源沿走向和倾向用50×15个子断层表示(子断层尺寸1.0 km×1.0 km)。

断层的位错分布模式会显著影响有限元模型的位移分布，而破裂传播方式将影响同震变形场的发展。因此，在土-结构相互作用的地震响应研究中，断层的位错分布模式和破裂传播方式是不可忽视的因素。断层面上的位错分布可以根据地表的同震位移进行反演^[39]，也可以基于凹凸体震源模型^[40]预测断层的位错分布。本文的位错分布按照SOMERVILLE等^[41]的单凹凸体位错分布模式(见图6)：凹凸体居中分布，凹凸体面积占断层总面积的22%(15 km×11 km)，凹凸体上位错为$2.01\overline u$；背景区域面积占断层总面积的78%，背景区域的位错为$0.71\overline u$。子断层的震源时间函数采用Kupper单周期小波模式^[42]；所有的子断层的持续时间${T_{\text{r}}}$都设置为4 s。破裂从断层面的一个角落开始传播(见图6，破裂开始的子断层的起始时间${T_0}$为0.1 s)，并以2.5 km/s的恒定破裂速度在断层面上传播。通过调整一系列子断层的破裂起始时间${T_0}$来模拟断层上破裂的传播，而起始时间可以根据子断层的相对位置和破裂速度计算。

图  6  凹凸体位错分布模式

Figure  6.  Asperity distribution mode

下载: 全尺寸图片 幻灯片

根据速度模型参数和震源参数，基于地震动模拟算法计算截断边界上节点的应力、位移和速度，并通过粘弹性边界模拟有限元模型的地震波动。在均质速度模型和位错不均匀分布(单凹凸体位错分布模式)的条件下，有限元模型不同时刻的同震位移场如图7所示。地层内部的变形规律与人工边界处一致，位移的变化均匀合理，验证了本文方法使用粘弹性边界模拟地震波动的合理性。

图  7  不同时刻的同震位移

Figure  7.  Coseismic displacement at different times

下载: 全尺寸图片 幻灯片

自由场模型的位移场在5 s~17 s逐渐增大；在17 s~25 s又逐渐减小；在30 s~40 s趋于稳定；50 s时的位移场可以认为是断层位错产生的地层永久变形。位移场在17 s附近达到变形峰值，该时刻取决于断层的破裂过程和地震波的传播。图8展示了有限元模型中沿断层倾向(y向)中线上地表各点(速度模型中的坐标：x=50 km；y=0 km)的位移时程。受到断层产状的影响，地表的竖向位移明显大于水平位移，且在断层上断线的地表投影附近(速度模型中的坐标y坐标为50 km~54 km)的位移变化最为剧烈。显然，地表不同站点的位移时程曲线不同，地震波的到达时间存在差异，位移达到峰值的时刻也不一致。因此，本文方法可以模拟大跨度结构抗震分析时的行波效应。此外，地震发生过程中的峰值位移与震后的地层永久位移存在显著差异，本文方法可以基于动力时程法分析断层位错产生的同震变形对结构的影响。

图  8  地表各点的竖向位移时程

Figure  8.  Vertical displacement time history at different sites on the ground

下载: 全尺寸图片 幻灯片

OKADA^{[22 − 23]}提出了半无限均匀弹性介质中矩形位错引起的位移场的解析表达式。为了验证本文方法的准确性，将有限元模型的地表永久位移的模拟结果与OKADA的理论解进行对比，如图9所示。水平变形的误差非常小，且主要分布在粘弹性人工边界处(见图9(a)和图9(b))，这可以归因于粘弹性边界的影响。竖向变形的误差主要分布在有限元模型中部，但误差值相对于变形值较小(见图9(c))。根据比较结果，模型解与解析解吻合较好，模型解的同震变形分布与解析解一致，相对误差较小。因此，本文基于有限断层震源的抗震分析方法具有良好的精度和准确性。

图  9  凹凸体位错分布模式下地表的同震位移 　/m

Figure  9.  Surface coseismic displacement under asperity dislocation distribution mode

下载: 全尺寸图片 幻灯片

5   土-结构相互作用模型算例

为进一步验证本文方法在土-结构相互作用系统的适用性，以成兰铁路穿越北川-映秀断裂的某隧道^[43]为例进行基于有限断层震源的抗震分析。建立如图10所示的围岩-隧道有限元模型。隧道顶板覆土厚度45 m，隧道断面简化为抗震效果最优的圆形，外径10 m，衬砌厚度500 mm，隧道的材料参数为：弹性模量为34.5 GPa，泊松比0.2，密度$\rho = 2500{\text{ kg/}}{{\text{m}}^{\text{3}}}$。当断层埋深较浅时，隧道衬砌的变形主要集中断层附近500 m范围内^[43]。为了分析断层位错产生的永久地层变形对隧道结构的影响，有限元模型长度(y向，即断层倾向)为2000 m，模型位于断层上、下盘范围内的长度均为1000 m(见图10)。考虑围岩对隧道结构的影响范围为3倍~5倍洞径，模型宽和高分别取为100 m和100 m。

图  10  围岩-隧道有限元模型 　/m

Figure  10.  Finite element model of surrounding rock and tunnel

下载: 全尺寸图片 幻灯片

龙门山断裂带宽30 km ~50 km，长度超过500 km，由近乎平行的3条逆冲兼走滑特征的后山断裂、北川-映秀断裂(中央断裂)和前山断裂组成。北川-映秀断裂是龙门山中央断裂带的主断裂，也是2008年汶川地震的发震构造之一。该断裂的倾角在近地表处为50°~80°，随着埋深增加逐渐变缓。本文对北川-映秀断裂的产状进行简化，其相关参数为：尺寸(长$\times$宽)为50 km×15 km；断层上断线埋深0.5 km，断层面几何中心的水平投影坐标为(40 km, 25 km)；滑移量为2.0 m，均匀位错分布模式；走向角为${0^\circ }$；倾角为${50^\circ }$；滑移角为${90^\circ }$(逆断层)；有限断层震源沿走向和倾向划分为50×15个点源(子断层尺寸为1.0 km×1.0 km)；震源时间函数采用Kupper单周期小波模式；所有子断层的起始时间${T_0}$为0.1 s；所有子断层的上升时间${T_{\text{r}}}$为4 s。

地壳速度模型的长×宽×高为80 km×50 km×50 km，速度模型的密度$\rho$、纵波波速${c_{\text{p}}}$和横波波速${c_{\text{s}}}$参考中国川藏地区典型地壳分层模型参数^[43]，见表1。地壳的P波品质因子${Q_{\text{p}}}$和S波品质因子${Q_{\text{s}}}$可以基于BROCHER^[44]总结的横波速度与品质因子的经验关系计算。有限元模型尺寸为100 m×2000 m×100 m，其几何中心在速度模型中的坐标为(40 km, 29.82 km, 0.05 km)。有限元模型中的围岩参数应该与地壳第一层岩土的材料参数一致，经过换算，围岩的材料参数如下：弹性模量为6.36 GPa，泊松比为0.39，密度为$2000{\text{ kg/}}{{\text{m}}^{\text{3}}}$。

表  1  川藏地区典型地壳速度模型参数

Table  1.  Key elements of regional seismic velocity model in Sichuan-Tibet area

深度/km	纵波波速/ (km·s−1)	横波波速/ (km·s−1)	密度/ (kg·m−3)	P波品质 因子$Q_{\text{p}}$	S波品质 因子${Q_{\text{S}}}$
0~0.1	2.50	1.07	2000	153.21	76.61
0.1~2	3.80	2.20	2200	357.01	178.51
2~6	4.66	2.70	2500	486.05	243.02
6~16	6.20	3.60	3000	830.78	415.39
>16	6.56	3.80	3100	932.81	466.41

下载: 导出CSV 

| 显示表格

图11展示了断层破裂后25 s时隧道的整体位移，显然，隧道结构沿纵向(y向)发生了十分显著的不均匀变形。受到断层产状的影响，隧道的竖向(z向)位移较大。隧道在断层上盘(主动盘)处的整体变形较大，但沿纵向的相对错动量较小。在断层上断线附近，隧道发生了非常剧烈的相对变形，这或许将导致更加严重的震害。因此，针对近活动断层区域的隧道抗震分析中，隧道结构的永久变形和沿纵向的不均匀错动是不可忽视的因素。

图  11  隧道的整体变形图

Figure  11.  The deformation diagram of tunnel

下载: 全尺寸图片 幻灯片

图12展示了沿断层倾向(y向)不同位置的隧道拱顶的z向位移时程和z向加速度时程。隧道在上盘范围内的z向最大永久位移(断层破裂后25 s时)为1.77 m，z向加速度峰值达到了$0.93{\text{ m}} \cdot {{\text{s}}^{ - 2}}$。但是，在上盘范围内，隧道的相对变形较小，其z向相对错动量小于0.40 m。因此，受到主动盘抬升的影响，隧道在上盘范围内发生了显著的平移，并承受了较大的惯性力荷载。在断层上断线附近，隧道在纵向长度800 m(y=−400 m~400 m)范围内的z向相对错动量为1.17 m(见图12)。

图  12  隧道拱顶的竖向位移时程和竖向加速度时程

Figure  12.  Vertical displacement time history and vertical acceleration time history of the tunnel vault

下载: 全尺寸图片 幻灯片

因此，对于结构与围岩相互作用的地下结构而言，地层的错动变形在地震破坏效应中起着关键性作用。而本文的抗震分析方法可以对结构施加由地震动引起的惯性力和由断层位错产生的地层变形荷载。此外，该方法同样适用于跨断层桥梁、大坝等大型建筑物的抗震分析。

6   结论

本文将基于震源的地震动模拟方法引入土-结构动力相互作用分析的研究领域，提出一种考虑有限断层震源参数的抗震分析方法。新方法明确“震源-路径-结构”的地震波传播的完整动力过程，可以考虑断层产状、位错的分布形式和破裂的传播方式等因素对地震作用的影响。

基于有限断层震源的土-结构相互作用抗震分析方法的基本原理为：基于地震动模拟方法计算截断人工边界上的边界条件，并通过粘弹性人工边界模拟有限元模型的地震波动。实施流程为：首先，建立土-结构相互作用的有限元模型；然后，根据震源参数和速度模型参数，基于地震动模拟方法计算模型边界上所有节点的应力、位移和速度；最后，利用粘弹性人工边界模拟有限元模型的地震波动。

为验证本文提出的抗震分析方法的有效性，对自由场模型算例进行计算分析。考虑凹凸体位错分布的有限断层震源模型对地震波传播和同震变形的影响。自由场算例的模拟结果与OKADA位错理论的解析解吻合较好，证明本文提出的抗震分析方法具有较高的精度。

本文进一步验证了基于有限断层震源的抗震分析方法在土-结构相互作用模型(围岩-隧道有限元模型)的适用性。该方法可以对结构由地震动引起的惯性力和由断层位错产生的地层变形荷载。结果表明：对于结构与围岩相互作用的地下结构而言，地层的错动变形在地震破坏效应中起着关键性作用。此外，本文的方法同样适用于跨断层桥梁、大坝等大型建筑物的抗震分析。