0.   引言

近年来地震频发，对地下结构和工程场地的安全性构成了严重威胁，地下结构抗震分析逐渐成为当前热点问题^[1]。数值模拟方法作为此问题的重要研究手段之一，因其适用性强、能够模拟复杂结构形式及场地条件等优势而被广泛采用。目前地下结构抗震分析常用的数值方法包括有限元、边界元及相应耦合方法等，但在实际应用中仍存在不足，如地震动输入及非自由场散射效应^[2]、半无限地基场地模拟及边界效应^[3]等因素均会影响数值模拟的计算精度和效率^[4]，并存在计算尺度、计算效率、计算精度之间的制约，因此，如何精确、高效地模拟地下结构与地层系统的地震响应特征仍是亟待解决的关键问题。

为了模拟半无限介质中地下结构的地震响应，有限元方法通常基于构建局部人工边界来处理无限域和地震动输入问题，其中基于等效节点力的黏（弹）性边界方法应用最为广泛。刘晶波等^{[3, 5]}率先提出基于黏（弹）性边界等效节点力的有限元模型地震动输入方法，实现了地下结构抗震分析的波动输入模拟；谭辉^[6]基于多尺度分析思路，采用人工边界子结构法对大尺度震源-复杂场地-结构模型进行了研究；赵密等^[7]基于黏弹性边界提出一种深埋地下结构高效时程分析方法；金丹丹等^[4]通过构建大尺度复杂地形真实场地研究了复合场地非线性地震响应特征。由于上述人工边界的精度与等效地震荷载计算、波源距离、波动传播方向等因素有关，故其边界位置选取问题仍待完善。针对该问题，Bielak等^[8]将外源激励转化为等效荷载施加在虚拟截断边界区域，提出区域缩减法(DRM)。随后，Wang等^[9]基于DRM研究了层状地基斜入射地震动输入问题；Zhang等^[10]基于DRM实现了局部复杂结构和地形的多尺度地震动模拟；Kontoe等^[11]和胡丹等^[12]将DRM与基于等效节点力的黏弹性边界方法进行对比，认为DRM方法能精确描述等效地震力，并且具有截断边界不受限制、所需计算模型规模更小等优势。

虽然有限元方法（FEM）在地下结构抗震分析中应用较广，但在处理大尺度工程场地模型动力问题时仍存在计算效率问题。相比而言，边界元方法（BEM）具有降维和同时满足无限远辐射条件的优势而被广泛应用于复杂地形条件的地震模拟^[13]。然而，由于在边界处将波动方程转化为边界积分方程而构造的满秩矩阵，因此对于复杂地下结构的建模计算成本较高^[14]。鉴于边界元在处理无限域问题方面的优势，有限元-边界元耦合逐渐成为解决地下结构地震响应问题的有效分析方法^[15-16]。但目前基于BEM或FEM的传统耦合方法尚需构建统一矩阵，破坏了有限元的带状矩阵特性，故导致计算效率较低，而且分区耦合的迭代收敛性问题会影响计算精度^[17]。为此，目前仍迫切需要建立能够合理模拟地下结构与工程场地耦合体系地震响应特征的精确、高效分析方法。

本文针对地下结构地震响应模拟在计算尺度、计算精度和计算效率方面的实际需求，旨在提出一种地下结构-地层耦合系统地震响应分析的高效混合数值模拟方法。相比于已有研究，该方法可以更全面考虑大尺度工程场地和复杂地下结构模型（如多尺度问题），基于区域缩减法实现自由场/非自由场运动的精确输入，且在保证计算精度的前提下对地下结构地震响应进行快速参数化分析。方法基于区域缩减法将边界元和有限元相融合的核心思想，将地下结构-地层整体模型进行划分，即近场地层及结构为内域，采用有限元模拟，远场地层条件为外域，采用边界元求解外域自由场或地形影响下的非自由场地震动，通过DRM构造矩阵方程将外域动力响应转化为等效地震荷载，并施加在内域有限元模型截断边界处，进而实现对内域中地下结构的地震响应分析。通过二组典型算例分析，验证了该方法的可靠性和高效性。

1.   模型和方法

本文基于区域缩减法将边界元和有限元相融合的核心思想，首先在地下结构-地层系统整体计算模型中构建截断边界${\mathit{\Gamma}}$，将计算模型划分为外部区域和内部区域，如图 1所示，其中，外部区域包括震源、自由场或地形条件（如河谷、山体、透镜体等局部不均匀地质构造）影响下的非自由场；内部区域包括复杂地下结构（群）及相邻地层。图中内域为${{\mathit{\Omega}} ^ - }$，外域为${{\mathit{\Omega}} ^ + }$，内域、外域及相邻边界面上的结点位移分别表示为${u}_{\text{i}}，{u}_{\text{e}}，{u}_{\text{b}}$。高效混合数值模拟方法主要包含以下2个步骤：

（1）采用边界元方法构建外域计算模型，首先给出自由场响应，其次考虑地形条件影响，在地形边界处将波动方程转化为边界积分方程构造散射场，求解地形效应引起的散射场响应，并将自由场和散射场相叠加，从而得到虚拟边界面${\mathit{\Gamma}}$处的自由场或地形影响下的非自由场动力响应。

（2）采用有限元方法构建内域计算模型，采用DRM构造内域和外域有限元矩阵方程并进行整合，给出等效荷载表达形式，即将外域动力响应转化为等效地震荷载，并施加在内域有限元模型截断边界处，最终实现内域中地下结构地震响应的参数化分析。

在该高效数值模拟方法中，第1步的目的是通过边界元法对外域计算模型进行动力分析，确定截断边界处的动力响应；第2步的目的是通过DRM将外域动力响应转化为内域的地震动输入，通过有限元法进行地下结构地震响应分析。

图  1  计算模型

Figure  1.  Computational model

下载: 全尺寸图片 幻灯片

本方法可以计算局部地形影响下的非自由场地震动，利用边界元计算大尺度场地条件下的动力响应，从而更大范围地缩减外域模型，提高地下结构动力分析的计算效率。以下主要针对所建立的高效混合数值模拟方法及计算流程进行详细阐述。

1.1   外域边界元求解

首先考虑自由场响应，即不存在地形条件影响。假设一个圆频率为$\omega$的平面SV波以$\theta$角从半无限空间入射，在直角坐标系中其波势函数可以表示为

式中：$k$为SV波波数；i为虚数单位。入射SV波将在地表产生反射P波和反射SV波，两者的波势函数分别表示为

式中：$a，b$分别为两种反射波系数；自由波场引起的应力为$\sigma _{ij}^{\text{f}}(i, j = x, y)$；位移为$u_i^{\text{f}}$，上标f表示自由场，具体表达式详见文献[18]。

当存在地形条件影响时，以河谷为例（其余情况类似），即考虑二维半无限空间域中的河谷地形，如图 2所示。采用边界积分方程法^[18]，基于单层位势理论在边界（场点）附近设定虚拟波源面（源点），散射波由虚拟波源的作用叠加而得，进而由边界条件建立方程求得虚拟波源密度。另外，结合本文模型及参考文献[19]中的点位配置特征，虚拟波源面半径取值0.4R~0.6R，波源点数取为隧道边界离散点数的一半，边界配点数满足每个波长7~10个点，即可保证较高的计算精度。针对本文地形影响下的非自由场研究对象，在河谷边界内引入一个虚拟波源面S₁，河道边界为S，来构造地层中地形条件引起的散射波。

图  2  外域BEM模型

Figure  2.  BEM model for outer domain

下载: 全尺寸图片 幻灯片

由于河谷地形存在，将弹性半无限空间波场分为自由场和散射场叠加。根据单层位势理论，通过在散射体内部虚拟面上的波源积分，可以模拟扩展到无穷远外部区域中的位移$u_i^{\text{s}}$、应力$\sigma _{ij}^{\text{s}}$，上标s表示地形引起的散射场，积分表达式为

式中：$\mathit{\boldsymbol{x}} \in S, \mathit{\boldsymbol{x'}} \in S'$；$b({x'_1}), c({x'_1})$分别为虚拟波源面$S'$上P波和SV波的波源密度，$G, T$分别为弹性半空间中的位移格林函数和应力格林函数。

由于边界元方法自动满足无限远处辐射条件的优势，为了获得虚拟波源的解，只需利用河谷表面的边界条件构造方程组即可。

此外，河谷边界S处应满足应力为零的条件：

式中：${t_i} = {\sigma _{ij}}{n_j}$表示边界力；${n_j}$为外法向量分量。

为了便于求解，还需对散射波场的积分方程进行离散化处理，结合边界条件，给出方程式矩阵：

采用最小二乘法求解方程组（7）从而得到虚拟波源密度$b, c$；再与格林函数矩阵T和G相乘，即可得到散射场中的应力和位移；然后和自由场结果叠加得到弹性波入射下半空间中任意位置处的动力响应解答；最后，借助快速傅立叶逆变换将频域解叠加，求解动力时程响应^[19]，即可得到内域与外域截断边界${\mathit{\Gamma}}$处的自由场响应或地形影响下的非自由场动力响应。

1.2   内域有限元求解

基于区域缩减法，将外域自由场或地形影响下的非自由场计算模型进行区域缩减（图 3），内域模型仅包括地下结构及相邻地层。

图  3  内域FEM模型

Figure  3.  FEM model for inner domain

下载: 全尺寸图片 幻灯片

不考虑阻尼效应，内域和外域的有限元格式动力控制方程分别为

式中，M为质量矩阵；K为刚度矩阵；上标$^\prime$表示外域；${F_{\text{b}}}$为内域和外域之间的相互作用力，${F_{\text{e}}}$为外域地震荷载。

将整体计算模型分解为自由场或考虑地形影响的非自由场和内部结构引起散射场的叠加，可表示为

其中，${\text{f'}}$表示自由场或地形影响下的非自由场，$s'$表示内域工程结构引起的散射场。

将内外域进行整合，得到模型整体平衡方程：

可以看出，式（11）左侧与内域和边界相关的项和式（8）是完全一致的，说明求解方程可以获得模型内域响应的精确解答。方程右侧的表达式含义是将外域远场作用转化为边界局部区域上的等效荷载。该等效荷载与荷载施加区域内的质量矩阵、刚度矩阵及外域动力响应相关，即

从式（12）可以看出，DRM方法通过构造节点荷载施加区域来重现地震波场的作用，即将1.1节通过边界元方法求得的自由场或地形影响下的非自由场地震响应代入式（12），即可获得内外域边界${\mathit{\Gamma}}$处的动力响应，从而实现外域向内域传递的等效荷载输入。内域地下结构引起的散射波，可由边界${\mathit{\Gamma}}$之外的模型吸收，如将外层模型设置为由多层具有高黏性阻尼的有限元单元组成的吸收层，从而避免反射波对内域动力响应产生影响。

1.3   混合方法实施

本文所提出的地下结构-地层系统地震响应分析的高效混合模拟方法计算流程如图 4所示。

图  4  数值方法实施

Figure  4.  Implementation of numerical method

下载: 全尺寸图片 幻灯片

具体实施步骤描述如下：

（1）采用边界元方法建立外域计算模型，按照地形条件设置波源点和场源点，在半无限域中构造自由场和散射场，得到平面波入射下弹性半空间中截断边界处的动力响应解答。

（2）依据工程场地与地下结构分析对象的实际需求，建立地层-结构体系的动力分析内域模型，采用DRM构造矩阵方程计算截断边界区域的等效输入荷载。

（3）修改内域模型的input计算文件，即将步骤（2）中的等效地震荷载施加在地层-结构模型对应的截断边界节点上实现外域地震荷载输入。

（4）选择合适的材料本构模型，进行地层-结构体系的地震动力响应及参数化分析。

需要注意的是，本方法在划分有限元网格和设置边界点时，应按照以下形式：有限元域沿边界${\mathit{\Gamma}}$共划分N个有限单元，边界元域沿边界${\mathit{\Gamma}}$的每个边划分N+1个边界单元，如图 5所示。令有限单元边界节点和边界元节点重合，即形成有限元域和边界元域的重合节点，实现外域地震动响应向内域的精确荷载传递，从而进行内域模型地震动力响应的参数分析。

图  5  混合方法单元节点分布

Figure  5.  Element node distribution of hybrid method

下载: 全尺寸图片 幻灯片

2.   数值算例

为了说明本文所建立的高效数值模拟方法用于地层-结构体系地震响应分析的可靠性和高效性，设计了两组典型算例，分别为无地形影响（简称标准工程）、考虑地形影响条件下的双线隧道地震响应分析，并分别采用本文方法与目前常用求解方法的计算结果进行对比分析，包括远置边界方法和传统黏弹性边界方法。

为便于后文分析，选取Ricker波和Chuetsu波作为地震动输入，其中Chuetsu波加速度和位移时程曲线如图 6所示，Ricker子波^[20]沿x-z平面内以SV波形式垂直入射。

式中，Ricker子波主频${f_{\text{p}}} = 3{\text{ Hz}}$，位移时程振幅为0.02 m，峰值时刻为1 s。

图  6  Chuetsu波加速度和位移时程曲线

Figure  6.  Time-history curves of accelerations and displacement of Chuetsu waves

下载: 全尺寸图片 幻灯片

2.1   标准工况双线隧道地震响应

本节标准工况假定为均匀半空间中双线圆形隧道地震响应，隧道衬砌外径为5.5 m，厚度为0.5 m，埋深为15 m，隧道间距为20 m，结构与地层接触边界为共节点约束，地层和隧道衬砌均假定为线弹性，其材料参数见表 1。

表  1  地层和衬砌的材料参数

Table  1.  Material parameters of strata and linings

材料参数	弹性模量/MPa	泊松比	密度/(kg·m-3)
衬砌	32500	0.20	2500
地层	100	0.25	2000

下载: 导出CSV 

| 显示表格

由于本标准算例不考虑地形影响，在外域计算时仅考虑自由场地震动输入即可。本文方法计算模型的宽度和高度分别取为60 m×30 m，单元数5199，见图 7。对于远置边界方法，由于受人工边界的影响，图 8给出了远置边界最小尺寸的验证，按照模型长宽比为3/2等比例缩小模型，验证结构响应的加速度峰值精度。图中横坐标为计算模型的宽度，左侧纵坐标表示结构和地层观测点的峰值加速度，右侧纵坐标表示结构和地层观测点的百分比误差。经过试算分析，对应的最小模型宽度（L）和高度（H）约为150 m×100 m，单元数18658，该尺寸下结构响应不受边界效应影响，并将其计算结果作为算例验证的参考解。

图  7  标准工况双线隧道

Figure  7.  Two tunnels under standard conditions

下载: 全尺寸图片 幻灯片

图  8  远置边界最小尺寸验证

Figure  8.  Validation of minimum size of far boundary

下载: 全尺寸图片 幻灯片

图 9，10分别给出了两种方法计算模型中地表观测点A和隧道结构上的观测点B处（位置见图 7）的位移动力响应时程。可以看出，本文方法得到的位移时程曲线与参考解吻合良好，表明本文方法的正确性。此外，图中还分别给出了两种地震动作用下位移结果与远置边界参考解的峰值误差，首先给出误差值定义$R = {{\left| {r{{(t)}_{\max }} - {r_0}{{(t)}_{\max }}} \right|} \mathord{\left/ {\vphantom {{\left| {r{{(t)}_{\max }} - {r_0}{{(t)}_{\max }}} \right|} {\left| {{r_0}{{(t)}_{\max }}} \right|}}} \right. } {\left| {{r_0}{{(t)}_{\max }}} \right|}}$，其中${r_0}{(t)_{\max }}$为峰值位移参考解，$r{(t)_{\max }}$为本文方法计算结果，图中${R_{\text{d}}}$表示本文方法的峰值误差。与远置边界参考解相比，不同地震动输入下的峰值误差均在8%以内，且单元数量减少72%。可见，本文方法在具有较高的计算精度前提下可以大幅提高计算效率。

图  9  Ricker地震动下位移时程

Figure  9.  Time histories of displacement from Ricker waves

下载: 全尺寸图片 幻灯片

图  10  Chuetsu地震动下位移时程曲线

Figure  10.  Time histories of displacement from Chuetsu waves

下载: 全尺寸图片 幻灯片

2.2   河谷地形双线隧道地震响应

本节研究均匀半空间中河谷地形影响下的双线圆形隧道地震响应，由于考虑均匀半空间地层中存在河谷地形影响，首先在外域计算中采用边界元给出河谷影响下的非自由场地震动进行输入。构建外域计算模型，其中河谷尺寸半径为30 m，外域河谷场地缩减后的有限元计算模型宽度和高度分别为60 m×30 m，见图 11，隧道和地层参数与2.1节一致。经验算，选取远置边界模型宽度和高度分别为300 m×200 m，其计算结果作为参考解，同时建立黏弹性边界模型，通过试算在保证计算精度的前提下模型最小尺寸确定为180 m×60 m，并与计算结果进行对比。

图  11  河谷地形双线隧道

Figure  11.  Two tunnels near river valley

下载: 全尺寸图片 幻灯片

图 12，13分别给出不同地震动条件下3种计算方法中观测点C处的位移时程反应谱。可以看出，黏弹性模型和本文方法在观测点C处的误差值在10%以内，结果均与参考解吻合良好，说明地形影响下的非自由场地震动可以实现精确输入。但相比参考解和黏弹性边界方法，计算模型分别缩减了97%和83%，可见，本文方法在保证精度要求的基础上，可以大范围缩减包括地形条件在内的外域计算模型的尺寸。此外，对于边界处地震波反射问题，本文方法相比黏弹性边界模型可以更好地吸收外行波，这是由于本文方法是在截断边界处进行地震输入，利用截断边界外的单元吸收外形波，而黏弹性边界需同时承担地震波输入和吸收外行波的作用，除了受到内部结构散射场的影响外，还与黏弹性边界模型尺寸大小有关。

图  12  Ricker波入射下观测点C处的位移时程

Figure  12.  Time histories of displacement at observation point C from Ricker waves

下载: 全尺寸图片 幻灯片

图  13  Chuetsu波入射下观测点C处的位移时程

Figure  13.  Time histories of displacement at observation point C from Chuetsu waves

下载: 全尺寸图片 幻灯片

表 2分别给出3种计算模型和方法下的计算耗时对比。计算机的基本配置如下：CPU：AMD Ryzen 7@3.2 GHz；内存：16.0 GB；CPU并行数量为4。选取地震波Chuetsu，采用隐式算法，时间步长为0.0051 s。可以看出：本文方法相比较远置边界参考解计算时间减少72%，相比基于等效节点力的黏弹性边界方法计算时间减少58%，说明本方法在保证较高计算精度的前提下极大提高了计算效率。

表  2  计算模型和耗时

Table  2.  Computational models and computational time

方法	单元数量	计算时间/s
远置边界	64089	3460
黏弹性边界	13218	1928
本文方法	5199	826

下载: 导出CSV 

| 显示表格

此外，在地震波斜入射问题的处理上，远置边界模型目前难以实现，黏弹性边界模型需要考虑不同侧边的时间延迟影响计算精度，而本文方法仅需要改变外域计算中的波势函数入射角度即可。以Chuetsu波为例，图 14，15分别给出SV波不同入射角情况下隧道结构和地表的动力响应，其中入射角度分别取为45°，60°和90°。可以看出：地表和衬砌结构动力响应均受入射角的影响，两观测点的水平和竖向位移放大效应发生改变。其中SV波垂直入射下，水平位移显著放大，随着入射角减小，地震动放大效应逐渐减弱。这是由于峡谷的存在，阻挡了地震波的传播，从而改变了结构周围的散射场分布，同时也影响了地震波的到达时间。

图  14  Chuetsu不同角度入射下观测点C处的位移时程

Figure  14.  Time histories of displacement at observation point C under different incident angles from Chuetsu waves

下载: 全尺寸图片 幻灯片

图  15  Chuetsu不同角度入射下观测点D处的位移时程

Figure  15.  Time histories of displacement at observation point D under different incident angles from Chuetsu waves

下载: 全尺寸图片 幻灯片

3.   结语

本文针对地下结构地震响应模拟的实际需求，提出了一种地下结构-地层系统地震响应分析的高效混合数值模拟方法，基于区域缩减法将边界元和有限元相融合的核心思想，将近场地层及结构作为内域采用有限元模拟，而采用边界元求解外域自由场或地形影响下的非自由场地震动，通过DRM构造矩阵方程将外域动力响应转化为等效地震荷载，并施加在内域有限元模型截断边界处，从而实现对内域中地下结构地震响应的高效模拟分析，并结合典型算例验证该方法的可靠性和高效性。结果表明，对于标准工况双线隧道地震响应，不同地震动输入下的峰值误差与远置边界参考解相比均在8%以内，验证了本文方法的可靠性，而对于河谷地形双线隧道地震响应，相比于参考解和黏弹性模型，计算模型分别缩减97%和83%，计算时间减少72%和58%，说明本文方法在保证计算精度的基础上，可以大范围缩减包含地形在内的外域模型尺寸，在具有较高的计算精度前提下提高计算效率。

本文所提出的方法还可以用于解决地震动不同角度入射问题，如斜入射地震动作用下地下结构的抗震设计与分析。需要说明的是，对于非均匀地形尺寸较小，且场地与结构较近的计算模型，为避免二次散射效应引起的误差，应将场地和结构一起建模，仅采用边界元计算自由场即可。而对于大尺度场地模型，尤其是大于结构尺寸2，3个数量级的情况，需采用边界元计算场地效应，以完成非自由场运动的精确输入和地下结构地震分析。此外，该方法还可以扩展到近场非线性或饱和场地地下结构地震响应分析，即考虑地层非线性和饱和地层中的流-固耦合作用对地下结构地震响应的影响，将另文分析。