引言

随着“21世纪海上丝绸之路”和“海洋强国”战略的提出与实施，我国加大了对海洋领域的开发和利用，大量海洋工程结构与基础工程正在快速建设发展中。我国的南海、渤海、台湾海峡等沿海地区位于环太平洋地震带西侧，海洋工程结构面临着严峻的地震灾害威胁（李小军，2006）。为保证海洋工程的地震安全性，首先需要对海域场地的地震动进行准确的预测。已有地震动预测手段主要有3种，即根据地震动衰减关系预测、理论解析预测以及数值模拟预测。然而，由于海域地震记录的缺失以及海域土层信息的不完善，导致海域地震动衰减关系确定缺乏相应的数据。相较于陆域场地，海域场地表层通常具有较厚的软土层（李小军等，2020；Hu等，2020），无法直接参照陆域地震动衰减关系。此外，众多研究人员将海域地形简化为规则模型，通过理论方法求解海域地震动响应（Brekhovskikh，1980；朱镜清，1988；Okamoto等，1999；柯小飞等，2019），但此类方法仅适用于规则地形（如水平成层场地），无法直接用于复杂地形的海域场地地震动模拟。

随着计算机技术的发展以及数值计算方法的成熟，加之近海地震台网建设的推进，使得利用数值模拟方法求解复杂海域地震动响应成为可能。Link等（2009）推导出流体-固体-声场相互作用的二维有限元格式。通过对Suruga海湾的地震动模拟，Nakamura等（2012）认为，海域不规则地形以及海水层对S波尾波的幅值和持时有显著影响。Li等（2017）考虑了海水层和饱和土层的影响，得到了成层海域场地的地震动传递函数，随后以南加州近海场地作为算例，并将数值解与实际地震动进行对比，其结果表明海水层会影响靠近海水P波共振频率附近的海域场地竖向运动。Dhakal等（2017）借助6个海底地震台站记录，对日本Sagam海湾的场地非线性效应进行了分析，结果表明，与陆地台站记录相比，海域场地更易进入非线性，且更易过滤掉地震波中的高频部分。王笃国等（2021）基于南海海域实测钻孔，分析不同地震动输入下海域场地地表峰值加速度和特征周期的变化规律，认为海域场地地震动长周期成分较为丰富。可见，相对于陆域场地地震动，海域场地的地震动模拟需将上覆饱和土层以及海水层纳入考虑。此外，相较于一维模型以及二维模型，三维海域模型能更准确的模拟海底地震动（Wang等，2020；Takemura等，2020；Bao等，2022）。相应的，三维模型海域地震动数值模拟的计算量也急剧增长，常规的单进程计算已无法满足计算需求。为了提高计算效率，充分利用计算机性能，Okamoto等（2010）、Liu等（2017）、Maeda等（2017）分别使用GPU并行、CPU主从通信并行、CPU缓存通信并行对三维场地进行地震动模拟。但即便采用并行计算，受计算机硬件以及计算时长的限制，现有三维大规模复杂海域场地地震动研究中，仍只能对频率低于0.5 Hz的地表运动进行较为准确的模拟（Okamoto等，2017；Takemura等，2019；Oba等，2020；Takemura等，2021）。若想将模拟频率向高频拓展，需要更为精细的海域模型，这也意味着更多的单元数量和更大的计算量。而以上求解复杂海域地震动的数值算法中，多将海水、饱和土、基岩的介质方程离散，随后进行整体求解或分区耦合求解，这使得针对不同介质需调用不同的求解器以获得该材料的响应，降低了计算效率。因此，除了使用并行技术提高计算效率以外，还需从算法的角度提高效率。

针对上述问题，为实现高效三维海域地震动模拟，本文采用流固耦合统一计算框架（陈少林等，2019a，2019b）模拟海水-海床相互作用，避免将海水、饱和多孔介质、基岩放在多个求解器求解，亦无需在界面交界处进行不同求解器之间的数据传递，提高计算效率；采用透射边界用于模拟无限域的影响（Liao等，1984），以自由场作为输入（柯小飞等，2019），实现地震波输入时海域场地地震响应的模拟。通过自编程实现建模、自由场分析和大规模海域场地地震响应并行分析的自动化计算，用户仅需输入相应的模型信息、入射波时程以及进程划分情况即可，使用方便。本文以日本东京湾为例，在超级计算机上模拟了脉冲SV波从模型底部垂直入射时的地震动响应。

1.   基本理论

1.1   广义饱和多孔介质运动的数学模型

饱和多孔介质是一种孔隙率β介于0到1之间的材料，而理想流体和理想弹性固体可视为孔隙率分别为1和0的特殊饱和多孔介质。基于此，我们将固体、流体、饱和多孔介质统一视为广义饱和多孔介质（图1）。

图 1  介质关系示意图

Figure 1.  Schematic diagram of media relations

下载: 全尺寸图片 幻灯片

根据Biot多孔介质理论，饱和多孔弹性介质的基本微分方程表示如下（Biot，1956，1962）：

固相平衡方程：

式中，${{\boldsymbol{L}}_{\text{s}}}$、${{\boldsymbol{L}}_{\text{w}}}$为微分算子矩阵；${\boldsymbol{\sigma}} '$为固相有效应力矢量；$\tau$为平均孔隙水压力，以受拉为正；$P$为孔隙水压力，以受压为正；$ {\boldsymbol{U}} $、$ {\boldsymbol{u}} $分别代表液相和固相的位移，${\boldsymbol{\dot U}}$、${\boldsymbol{\dot u}}$分别代表液相和固相速度，${\boldsymbol{\ddot U}}$、${\boldsymbol{\ddot u}}$分别代表液相和固相加速度；$\rho_{\rm{s}}$、$\rho_{\rm{w}}$分别代表固相和液相的密度；$\beta$为孔隙率，$b=\beta^2\mu_0 /k_0$，$k_0$为流体渗透系数，$\mu_0$为流体粘度系数；$E_{\rm{w}}$为流体的体积模量；$e_{\rm{s}}$和$e_{\rm{w}}$分别代表固相和液相的体积应变。

假设$\beta {\text{ = }}1$、${\rho _{\rm{s}}}{\text{ = }}0$、${\mu _0}{\text{ = }}0$、$G{\text{ = }}0$，此时饱和多孔介质中无固相介质，即为纯流体状态。式（1）、（2）中的${\boldsymbol{\sigma}} '$和$b$均为0，式（1）等号两侧均为0，式（2）退化为理想流体方程$ - {\boldsymbol{L}}_{\rm{w}}^{\rm{T}}P = {\rho _{\rm{w}}}{\boldsymbol{\ddot U}} $。同理，假设$\beta {\text{ = 0}}$、${\rho _{\rm{w}}}{\text{ = }}0$、${\mu _0}{\text{ = }}0$、${E_{\rm{w}}}{\text{ = }}0$，此时饱和多孔介质中无液相介质，即为纯固体状态。式（1）、（2）中的$P$和$b$均为0，式（2）等号两侧均为0，式（1）退化为弹性波动方程$ {\boldsymbol{L}}_{\rm{s}}^{\rm{T}}{\boldsymbol{\sigma '}} = {\rho _{\rm{s}}}{\boldsymbol{\ddot u}} $。因此，通过对指定参数的修改，广义饱和多孔介质方程可同时描述弹性固体、饱和多孔介质、理想流体。

边界条件由Dirichlet边界条件、Neumann边界条件、人工边界条件共同组成。

Dirichlet边界条件：

Neumann边界条件：

式中，$\tilde {\boldsymbol{u}}$和$\tilde {\boldsymbol{U}}$分别代表另一广义饱和多孔介质边界处对应节点固相和液相的位移，$\tilde {\boldsymbol{\sigma}}$和$\tilde P$分别代表另一广义饱和多孔介质边界处对应节点固相应力和孔压，${\boldsymbol{n}}$是沿边界外法线方向的方向向量，${\hat {\boldsymbol{n}}}$是由方向导数组成的矩阵。为了模拟无限域的波辐射效应对计算区域的影响，在模型边界处施加了透射边界。该人工边界条件通过前一时刻内节点的位移来确定此时边界节点的位移，因此可以归类为Dirichlet边界条件。

1.2   节点运动方程及求解

使用伽辽金法将式（1）、（2）离散，同时将边界条件以及式（3）中的相容方程代入，可以得到任意节点 i 的解耦运动平衡方程（陈少林，2019b）：

式中，$ {{\boldsymbol{M}}_{{\rm{s}}i}} $、$ {{\boldsymbol{M}}_{{\rm{w}}i}} $分别为集中到节点的固相和液相质量，$ {\boldsymbol{F}}_i^{\rm{s}} $和$ {\boldsymbol{F}}_i^{\rm{w}} $分别代表固相和液相的节点力，$ {\boldsymbol{T}}_i^{\rm{s}} $和$ {\boldsymbol{T}}_i^{\rm{w}} $分别代表固相和液相的节点粘滞阻力，$ {\boldsymbol{S}}_i^{\rm{s}} $和$ {\boldsymbol{S}}_i^{\rm{w}} $分别代表固相和液相的界面节点力。由于本文中地震响应是通过边界点的自由场位移输入，这里的$ {\boldsymbol{f}}_i^{\rm{s}} $和$ {\boldsymbol{f}}_i^{\rm{w}} $仅代表除地震荷载以外的固相和液相外力。

对于内部节点，固、液相界面力$ {\boldsymbol{S}}_i^{\rm{s}} $和$ {\boldsymbol{S}}_i^{\rm{w}} $均为0，对式（8）、（9）进行时步积分可得节点固、液相位移：

式中，$ \Delta t $为计算步距，上标带有p−1、p和p+1的物理量分别代表对应参数在$ t = (p - 1)\Delta t $、$t=p \Delta t$和$ t = (p + 1)\Delta t $时刻的值，$ m_i^{\rm{s}} $和$ m_i^{\rm{w}} $分别为集中到节点i的固相和液相的质量。

当节点i处于2种材料交界面时（图2），基于隔离体的概念，对式（8）、（9）进行时步积分可得（陈少林，2019b）：

图 2  界面力示意图

Figure 2.  Schematic diagram of interfacial force

下载: 全尺寸图片 幻灯片

式中，$\hat {\boldsymbol{u}}_i^{(p + 1)}$、$\hat {\boldsymbol{U}}_i^{(p + 1)}$分别代表未施加界面力时的固相和液相位移，$\Delta {\boldsymbol{u}}_{Ni}^{(p + 1)}$是由法向力${\boldsymbol{S}}_{Ni}^{{\rm{s}}p}$产生的固相位移，$\Delta {\boldsymbol{u}}_{Ti}^{(p + 1)}$是由切向力${\boldsymbol{S}}_{Ti}^{{\rm{s}}p}$产生的固相位移，$\Delta {\boldsymbol{U}}_{Ni}^{(p + 1)}$是由法向力${\boldsymbol{S}}_{Ni}^{{\rm{w}}p}$产生的液相位移，$\Delta {\boldsymbol{U}}_{Ti}^{(p + 1)}$是由切向力${\boldsymbol{S}}_{Ti}^{{\rm{w}}p}$产生的液相位移。

2种广义饱和多孔介质交界面处应满足如下界面连续条件（Deresiewicz等，1964）：

式中，$ {\tilde {\boldsymbol{\sigma}} _{N}} 、{ \tilde{\boldsymbol{\tau}}}、 {\tilde{\boldsymbol{ \sigma}} _{T}}、 {\tilde P}、 {\tilde {\boldsymbol{u}}_{N}}、 {\tilde {\boldsymbol{u}}_{T}}、{ \tilde \beta}、 { \tilde{\boldsymbol{U}}_{N}} $表示另一广义饱和多孔介质的对应参数，而含下标N的字母代表法向对应的物理量，含下标T的字母代表切向对应的物理量，$ \tau $表示平均孔隙水压力。由式（12）、（13），结合边界条件可得：

式（19）、（20）具体展开式以及各参数含义详见（陈少林等，2019b）。由式（19）、（20）得到${\boldsymbol{S}}_{Ni}^{{\rm{s}}p}$和${\boldsymbol{S}}_{Ni}^{{\rm{w}}p}$，根据界面连续条件可得$\tilde {\boldsymbol{S}}_{Nk}^{{\rm{s}}p}$和$\tilde {\boldsymbol{S}}_{Nk}^{{\rm{w}}p}$，随后利用界面处固相位移连续条件，可得：

根据界面连续条件，可求得${\boldsymbol{S}}_{Tk}^{{\rm{s}}p}$。而界面处的液相切向力始终为0。求得界面力后，可由式（12）、（13）求得交界面上的节点位移。

对于模型边界处的节点，使用多次透射公式模拟外行波穿过人工边界的运动（Liao等，1984）：

式中，$u_{o{\rm{c}}}^{(p + 1)}$、$U_{o{\rm{c}}}^{(p + 1)}$是固/液相外行波（散射波）在$t = \left( {P + 1} \right)\Delta t$时刻的位移；o表示人工边界点，j表示人工边界相邻节点序号；N是透射阶数，$C_j^N$可表达为：

这种局部人工边界条件具有较高的通用性，与具体的波动方程无关，可直接用于饱和多孔介质中的波动问题。首先，通过传递矩阵方法（柯小飞等，2019）获得自由场响应，然后将式（22）分别应用于固相和液相的散射场，得到边界点在p+1时刻的散射位移。最后，将边界点的自由场位移与散射位移相加，可以得到边界点的总固/液相位移。陈少林等（2019a）将该计算方法应用于水平成层场地，并将相应计算结果与解析解进行对比，验证了该计算方法的精度。

1.3   逐元并行策略

大规模计算所需内存和计算量通常较大，如何提高计算效率、节省内存是大规模地震动模拟的关键。常规有限元计算中，通常先形成单元刚度阵、质量阵和阻尼阵，然后组装并存储模型整体刚度阵、质量阵和阻尼阵，对于大规模计算，所需内存较大，若采用隐式积分格式，需求解大型方程组，计算量大，效率较低。为此，本文采用逐元法计算，该方法根据单元尺寸以及材料属性将单元分类，仅需形成和存储各单元类的刚度阵、质量阵、阻尼阵，无需组装整体刚度阵、质量阵和阻尼阵，可大大节省内存；计算每个节点的反应时，只需通过该节点周围单元刚度阵和阻尼阵得到该节点的本构力，并通过集中质量及显式积分格式得到该节点的反应，无需求解大型方程组，效率较高，且便于并行计算。

在模拟海域场地地震动响应的过程中，需经过建模、自由场计算、海域场地地震波场计算，整体流程较为复杂，具体如下：（1）建模。当模型单元数量较大（千万级）时，现有商业有限元软件或建模软件的效率低下甚至无法输出。此外，若在建模过程中生成个别畸形单元便可能导致最终波场计算出现错误，甚至发散，需识别畸形单元和调整网格剖分，十分麻烦。因此，本文采用规则的六面体单元，通过自编程实现了建模-自由场计算-场地地震动求解的全流程自动化运算。通过读取用户输入的场地信息（如土层材料参数、高程信息等）和单元尺寸，用立方体将整体场地模型离散为正六面体单元，并赋予相应材料参数，计算和提取不同介质（流体、饱和土、基岩）界面的方向导数。规则六面体单元对于复杂地形的模拟不如四面体单元，但一般可通过加密网格来满足精度要求，且可保证网格质量，配合逐元法，仅需存储少量类别的单元阵即可完成波场模拟，极大节省了计算机内存，适用于大规模场地计算（Ichimura等，2007）。（2）自由场。根据模型边界区的土层信息，自动计算相应自由场，并输出为指定的二进制格式，供地震波模拟时读取使用。（3）地震波场模拟。按海域场地统一计算框架，采用并行计算，输出所需点的响应，如图3所示。本文通过自编程，实现了建模-自由场分析-三维复杂海域地震动模拟的全流程自动化，可保证建模网格质量，提高整体计算的稳定性和鲁棒性，方便使用。

图 3  逐元并行FEM地震波场模拟流程图

Figure 3.  Flow chart showing the processes of seismic wave modeling using element-by-element parallel FEM

下载: 全尺寸图片 幻灯片

2.   数值算例

为验证上述算法及自编程序在大规模海域场地地震动模拟中的有效性，以东京湾为例模拟SV波由模型底部垂直入射时的情形，计算区域如图4黑色框线所示。模型水平尺寸为50 km×50 km，深度为0.8 km，由图可知东京湾东南方向存在部分山区，其余部分相对较为平整。地表高程及海底地形由Google Earth获得。Koketsu等（2009）根据波速和密度将东京湾附近区域的土层大致划分为4类，在120 km×120 km范围内给出各类土层相应参数以及各分界面深度的等高线图。本文计算区域涉及到其中2类土层、1个分界面。由于具体分界面数据难以获取，通过图5仅可知在计算区域范围内第一界面深度约为0.2～0.4 km，且由南向北逐渐变深。因此本文将该界面深度简化为由南向北从260 m至400 m线性递增。

图 4  东京湾及周边地形示意图

Figure 4.  Topography of the area around Tokyo bay

下载: 全尺寸图片 幻灯片

图 5  Koketsu等（2009）确定的第一界面深度

Figure 5.  The first interface depth determined by Koketsu et al.（2009）

下载: 全尺寸图片 幻灯片

假定模型仅由海水以及线弹性基岩构成，具体材料参数由表1给出。由模型底部垂直输入幅值为0.1 cm，脉冲宽度为1 s的位移脉冲波，时间步距0.0025 s，其位移时程、加速度时程以及相应傅里叶谱如图6所示。由图可知，输入位移波截止频率约为3 Hz，而对应的加速度频谱在高频区域仍有部分响应。为了满足计算精度的要求，有限元网格最大尺寸需满足（Moczo等，2000；杜修力，2009）：

表 1  东京湾材料参数（Koketsu 等, 2009）

Table 1.  Parameters of materials used in Tokyo bay （Koketsu et al., 2009）

材料	孔隙率/β	μ0	ρs /kg·m−3	Ρw/kg·m−3	ν	G /GPa	Ew /GPa	M /GPa	α	k0/μm2
海水	1	0	0	1000	0.020	0	2.25	2.25	1	1
基岩1	0	0	1850	0	0.437	0.666	0	—	0	0
基岩2	0	0	2080	0	0.395	2.080	0	—	0	0

下载: 导出CSV 

| 显示表格

图 6  输入脉冲

Figure 6.  Input pulse

下载: 全尺寸图片 幻灯片

式中，$ {\lambda _{\min }} $为离散网格模型中波动传播的最短波长，c_min为介质中的最小波速，f_max为波动问题数值模拟的截止频率。由于输入波截止频率为3 Hz以及该三维模型材料中最小剪切波为600 m/s，使用边长为20 m的正六面体单元对模型进行离散。这样不仅满足了计算精度要求，还极大减少了模型中的单元种类，节省内存的同时还可提高计算效率。

在模型四周和底部截断处施加人工边界用于吸收外行波。整体模型单元数量约为2.5亿，将模型划分为144个进程，采用64个核，每个进程约170万个单元，最终耗时38小时完成8000步计算。选取A、B、C作为观测点（图4），分别位于海底、平原、山顶。同时选取图4所示3个截面，给出基岩顶部的位移波场图。

图7为各观测点的位移时程、加速度时程以及加速度傅里叶谱。由图7（a）～（c）可知，受SV波由底部入射的影响，3个观测点x方向上的位移远大于另外2个方向，且能清晰观察到剪切波在不同介质间反射，第1个波峰到时也随着高程增加而延后。相较于处于平原的B点，处于山顶的C点受周边复杂地形的影响，各方向的位移响应均较为剧烈，这个特点在B、C点的加速度时程和傅里叶谱中均得到验证（图7（e）、图7（f）、图7（h）、图7（i））。值得注意的是，由图7（d）可知，A点的位移时程在10 s之后便较为稳定，但A点z方向的加速度却出现了明显抖动。这是因为A点处于海底，除了受底部入射的SV波影响，还受到海水的作用。SV波在前5 s 自模型底部向上入射，A点x方向上的位移和加速度均有响应，但由于海水不传递剪切作用，此时z方向几乎没有响应。与此同时，在近海岸处，海水受到岸边基岩影响产生位移。随着时间推移，由A点底部输入的SV波已经消散殆尽，而海水位移场传递到A点，由于海水海床法向连续，海水z向应力施加在基岩上，导致该点基岩出现z向的加速度。在整个计算过程中，我们假定材料为线弹性，因此无论地形如何，最终输出的加速度傅里叶谱的频谱特征与输入保持一致。

图 7  观测点的位移、加速度时程曲线和加速度傅里叶谱

Figure 7.  Displacement, acceleration time history curve and acceleration Fourier spectrum of observation points

下载: 全尺寸图片 幻灯片

图8为不同截面的位移波场图。截面1-1由海域和山区共同构成，因此整体波场较为复杂。对于海底区域，地势相对较为平整，且土层厚度最浅，因此该区域x方向位移在主响应结束后迅速归零，第一个波峰到时也最早。而在高程变化较为剧烈的山区，各点第一个波峰的到时存在差异，且当SV波到达自由表面后，出现了明显的散射现象。散射波又在不规则地形中多次反射，使得该区域位移波场出现较长时间的震荡（图8（a））。截面2-2是一个由平原-海域-平原组成的近似对称地形，波场几乎对称。与截面1-1相比，该截面较为规则，波场仅在地形变化相对剧烈的海岸区域发生局部轻微震荡。同样的，由于海域基岩层厚度最薄，入射波以及反射波的波峰到时均快于平原区域的（图8（b））。截面3-3是纯陆地区域，且地形变化缓慢，整体起伏较小，接近水平成层场地。因此，除了入射波与反射波以外，几乎看不出散射，整体波场也很快稳定（图8（c））。

图 8  截面位移波场图

Figure 8.  Displacement waveforms of cross sections

下载: 全尺寸图片 幻灯片

图9为不同时刻上层基岩顶部x方向位移的波场快照。基岩厚度最薄的海底区域最先响应，随后逐渐扩大至陆域，最后到达高程最高的山区。此外，由于底部基岩（基岩2）的剪切波速大于顶部基岩（基岩1），在该模型中，基岩2的高程由南向北逐渐递减，即在基岩厚度相同的区域，靠近南面的平均剪切波速更快，自模型底部入射的SV波也会更快到达基岩层表面。这导致在相同地形或相近高程区域（如海域）越靠近南面，SV波到时越早。

图 9  x方向位移的波场快照

Figure 9.  Snapshot of the wavefield with displacement in the x-direction

下载: 全尺寸图片 幻灯片

3.   结论

为实现三维大规模海域场地地震动模拟，本文从计算方法和策略2个方面对常规海域地震波场模拟进行改进。在计算方法方面，本文基于广义饱和多孔介质理论，将海水、海床放在同一框架中计算，避免不同求解器之间的数据传递；同时，采用集中质量显式有限元实现了方程解耦，从而提高计算效率。在计算策略方面，采用规则六面体单元建模，以及逐元并行策略，保证计算的鲁棒性，并提升计算效率，节省计算内存，使得大型场地计算成为可能。通过自编程，实现了建模-自由场计算-三维复杂海域地震动模拟的全自动化流程，实施方便。最后，以东京湾为例，模拟了SV波垂直入射的地震动场，验证了本文方法及自编程序的可行性和高效性。

本文方法仍存在以下不足：

（1）建模方面。本文采用规则单元，可减少单元类型和存储，但在模拟复杂地形方面效率不高；后续拟对局部地表采用不规则单元，其余部分采用规则单元，达到建模精度和存储之间的平衡。另外，本文算例中采用统一的单元尺寸，对于波速差异较大的介质，牺牲了计算量，在后续研究中，拟针对不同介质采用满足各自精度的单元尺寸，且不需过渡单元。

（2）材料本构方面。在本文算例中，将海床视作弹性介质，采用瑞利阻尼，未考虑强震作用下的材料非线性，后续拟采用更合理的本构模型。