引言

国内外大量震害观测资料及理论研究表明，沉积盆地对于进入其内部的地震波有显著的放大效应，进而导致盆地区域内震害加重（李雪强，2011）。如1985年墨西哥地震中，距震中超过400 km的墨西哥城由于其位于沉积盆地之上，发生的震害最为严重（Anderson等，1986）；1988年澜沧-耿马地震中距离震中200 km的施甸盆地出现Ⅶ度区高烈度异常（刘启方等，2013）；2008年汶川地震中位于渭河盆地中的西安市出现了Ⅴ度区中的Ⅵ烈度异常区（张宇翔等，2010）。鉴于此，沉积盆地地震动效应一直是地震学和地震工程学等领域研究人员关注的热点问题。

目前，众多学者对沉积盆地区域内的地震动进行了非线性研究。刘庆忠等（2015）采用Davidenkov动本构描述河漫滩-阶地土的非线性，将SH波作为地震动输入，采用ABAQUS软件计算分析了复杂场地地震动响应；董菲蕃等（2013，2014）分别采用等效线性和非线性（Martin-Seed-Davidenkov黏弹性动本构）方法对泉州盆地进行了地震动分析，通过改变基岩处的地震动强度，探究场地地形变化和土层不均匀性的地震动响应差异；金丹丹等（2012，2014）建立了福州盆地基岩场地二维模型，利用ABAQUS集群平台进行非线性模拟，土体动本构采用Martin-Seed-Davidenkov子程序计算，研究了盆地地形微变化和土层结构不均匀变化产生的地震动效应。

然而上述研究未考虑盆地内部饱和土体特性和近断层效应。滨海、河谷地带天然土体一般是饱和的，饱和场地地震波的传播与衰减特性同干土情况具有显著差别，同时面临更突出的液化、侧移、震陷等风险。近断层地震动由于断层距较小，受震源破裂过程的影响强烈，因此其地震动特点不同于远场地震动，如由断层滑动产生的永久地面位移、由破裂方向效应和“Fling step”效应引起的脉冲型地面运动等（王海云等，2006；刘启方等，2006）。考虑到近场地震动巨大的破坏力及与远场地震动的不同特点，刘中宪等（2017）、Liu等（2023）、Lu等（2024）对近断层沉积盆地地震动响应进行了模拟，给出了多种场地因素的影响，然而未考虑场地饱和特性，通过对不同震害情况的研究发现，饱和土层在地震中会发生液化，导致不可预知的次生震害。

本文利用ABAQUS软件对近断层附近的饱和沉积盆地进行了非线性模拟，建立了不同基岩面坡度的二维T型沉积盆地模型。利用二次开发平台，编写了基于Davidenkov骨架曲线的VUMAT子程序，深入讨论了不同盆地坡度、饱和特性、地震波强度作用下场地非线性反应。研究结果对于沉积盆地区域地震动估计及震害防御等均有重要意义。

1.   沉积盆地有限元模型

1.1   模型概况

沉积盆地模型如图1所示，外轮廓尺寸为500 m×200 m。由于沉积年代不同，盆地内部有3个沉积层，厚度分别为30、40、30 m，场地模型参数如表1所示。本文仅在盆地内部近地表土层内（第1层）考虑非线性本构。在场地地震动分析中，场地网格尺寸d与剪切波速$ {V}_{{\mathrm{s}}} $、频率$ f $有关，范围内应符合下式要求：

图 1  二维沉积盆地模型

Figure 1.  Two-dimensional sedimentary basin model

下载: 全尺寸图片 幻灯片

表 1  场地模型参数

Table 1.  Site model parameters

位置	层数	厚度/m	S波波速/（m·s−1）	网格类型	密度/（kg·m−3）	弹性模量E/MPa	静泊松比$ \mathrm{\vartheta } $	动泊松比$ {\mathrm{\vartheta }}_{\mathrm{d}} $
盆地内	第1层	30	205.8	CPE4 R	2 120	231	0.444	0.49
	第2层	40	300		1 800	453	0.4	—
	第3层	30	400		1 900	851	0.4	—
盆地外		200	600		2 000	1 920	0.334	—

下载: 导出CSV 

| 显示表格

盆地模型网格尺寸为5 m，在盆地内30 m厚度将网格尺寸加密，网格尺寸为2～3 m，可模拟的地震波频率范围为0～10 Hz。采用4种不同频谱特性的近断层地震波沿水平向输入场地。如图1所示，在盆地地表处设置8个测点，观察盆地内地表不同位置的地震动特性。测点①～④分别位于边坡1/2点及盆地底部脚点、1/4点、中点在地表的投影位置，测点⑤～⑦与测点①～③对称分布，测点⑧位于地表最左侧边缘处。

1.2   近断层地震波输入

为分析不同频谱特征的近断层地震波对场地非线性地震动的影响，选择汶川8.0级地震中江油含增51 JYH、绵竹清平51 MZK近断层速度脉冲型地震波及集集地震WGK-N-AT2、阪神地震KAK000-AT2加速度记录，以剪切波形式施加到场地模型。为对比分析地震波强度改变对场地非线性动力响应的影响，将地震波幅值曲线分别调整为0.1、0.2、0.3 g，调幅后的加速度时程曲线如图2所示，地震波反应谱归一化曲线如图3所示。

图 2  输入地震动调幅记录

Figure 2.  Input seismic amplitude modulation record

下载: 全尺寸图片 幻灯片

图 3  地震波加速度反应谱

Figure 3.  Seismic wave acceleration response spectrum

下载: 全尺寸图片 幻灯片

1.3   沉积盆地非线性模拟Davidenkov本构模型

ABAQUS软件中有丰富的内置材料模型，但这些模型并不能较好地处理岩土工程中的非线性动力问题，对常用的岩土非线性本构模型，需借助ABAQUS软件提供的用户材料子程序编译接口VUMAT或UMAT子程序，均支持FORTRAN语言编写，在模拟时主程序会自动调用编译的用户材料模型。本文开发的VUMAT子程序基本任务如下：①根据主程序传入的应变增量∆ε计算得到应力增量∆σ；②对计算过程中涉及到的状态变量SDV进行求解及更新。计算流程如图4所示。

图 4  VUMAT计算流程

Figure 4.  VUMAT calculation flowchart

下载: 全尺寸图片 幻灯片

Martin等（1982）提出了动剪应力-剪应变骨架曲线表达式：

式中，$ \tau $、$ \gamma $分别为剪应力和剪应变；$ {G}_{\max} $为初始动剪切模量；A、B和$ {\gamma }_{0} $为土的试验参数，此处的$ {\gamma }_{0} $并非参考剪应变幅值，仅做拟合参数使用。

将上述Davidenkov骨架曲线与Masing法则结合，可得Davidenkov应力-应变滞回曲线表达式（孙锐等，2004）：

式中，$ {\gamma }_{{\mathrm{c}}} $、${\tau }_{{\mathrm{c}}} $分别为加卸载转折点处的应变、应力。

应力-应变关系可表示为：

骨架曲线时变（$ \left(t+\Delta t\right) $时刻）切线剪切模量可通过式（2）求导得到：

滞回曲线时变切线剪切模量可通过对式（4）中的变量$ \left(\gamma -{\gamma }_{c}\right) $求导得到：

由骨架曲线和滞回曲线的时变剪切模量公式可知，等效剪应变的选取较为关键。采取应变偏量的第二不变量$ {J}_{2\varepsilon }\left({\varepsilon }_{ij}\right) $计算的等效剪应变可反映不同维度下的本构关系，但当直接采用其公式计算时，只能表达应变值大于0的状态，故将其改由增量的方式描述：

式中，$ {e}_{ij}^{0} $为应变偏量，$ {\gamma }_{{\mathrm{eq}}}^{t} $为t时刻等效剪应变；加载状态下，$ {\mathrm{sign}}=1 $，卸载时，$ {\mathrm{sign}}=-1 $；$ \Delta {\gamma }_{{\mathrm{incre}}}^{t+\Delta t}\left({e}_{ij}^{0}\right) $为等效剪应变增量；$ {\gamma }_{{\mathrm{gen}}}^{t}\left({e}_{ij}^{0}\right) $为t时刻广义剪应变，由第二不变量$ {J}_{2\varepsilon }\left({\varepsilon }_{ij}\right) $计算得到。

土体加卸载状态的判断方法采用等效剪应变优化算法（孙锐等，2004），可提高计算效率和精度。将当前状态下的应变值与转折点处的应变做差，即$ \gamma -{\gamma }_{{\mathrm{c}}} $，用新的应变偏张量$ {e}_{ij}^{0} $表示：

式中，$ {e}_{ij}^{0} $表示荷载每经历一次转向，应变张量将转折点$ {e}_{ij,{\mathrm{c}}} $作为零点重新开始记录。

1.4   Bryne孔压增量模型

本文应用了基于应变的Bryne动孔压模型，特征是考虑了土体的体积应变对动孔压发展的影响，在试验及模拟中仅需在不排水条件下进行，具有简单易操作的优点。同时，孔压发展也可直接与动应变幅值联系，有利于土体动孔压特性研究。

土中有效应力减小会导致土体体积变化，循环剪切作用引发的不可恢复体积压缩、孔隙水的排出和有效正应力的减小引起的土体体积变化，应满足饱和土体体积变化的相容条件。假定土的体积改变以压缩为正，孔隙水的体积改变排出时为正，则有：

式中，$ {\Delta \varepsilon }_{{\mathrm{v}},{\mathrm{d}}} $为往返剪切引起的体积压缩增量，$ {\Delta \varepsilon }_{{\mathrm{v}},{\mathrm{r}}} $为有效正应力减小导致的体积增量，$ {\Delta \varepsilon }_{{\mathrm{v}},{\mathrm{f}}} $为往返剪切引起的孔隙水变化量。假定地震动的瞬间饱和砂土不发生排水，则有：

饱和砂土回弹模量用$ {E}_{{\mathrm{r}}} $表示，每次往复剪切作用产生的孔压及减小的有效应力分别为$ \Delta u $和$ \Delta \sigma $，则有：

通过式（13）和式（14）可得到孔压增量与体积压缩增量的关系式为：

通过Martin和Finn试验结果分析，Bryne推导出了累积体应变$ {\varepsilon }_{{\mathrm{v}},{\mathrm{d}}} $和体应变增量$ {\Delta \varepsilon }_{{\mathrm{v}},{\mathrm{d}}} $的关系式为：

式中，$ {\gamma }_{a} $为第 N 次往复应力时的剪应变幅值；$ {\sigma }_{{\mathrm{c}}}' $为初始有效正应力，三维条件下采用八面体剪应变幅值$ {\gamma }_{{\mathrm{oct}}} $和初始八面体有效正应力$ {\sigma }_{{\mathrm{oct}}}' $；$ u $为超静孔隙水压力；$ {C}_{1}\mathrm{、}{C}_{2}\mathrm{、}m\mathrm{、}n\mathrm{、}k $为与土体性质有关的无量纲参数。

砂土系数$ {C}_{1}\mathrm{、}{C}_{2} $可通过Bryne推导出的相对密度$ {D}_{{\mathrm{r}}} $或修正标准贯入锤击数$ {N}_{1} $相关的经验表达式求得：

1.5   黏弹性边界

对地基边界做黏弹性处理并将地震作用加速度转换成等效节点力施加于边界处。不同时刻t下应力计算关系为：

式中，R为边界面处节点距震源的距离，$ {\mathrm{K}}_{{j}} $、$ {\mathrm{C}}_{{j}} $为常量。

边界面处的应力与位移、速度成正比，通过选择适当的弹簧阻尼器可模拟其应力状况，并且弹簧阻尼器的选择与边界面处的材料参数有关。通过在边界面处设置连续不断的弹簧阻尼器可以实现黏弹性边界的施加，弹簧刚度系数K与阻尼系数C的取值可通过下式确定：

式中，K为弹簧刚度；C为阻尼系数；R为波源距离边界的距离；c_p、c_s分别为介质压缩波、剪切波波速；G为剪切模量；ρ为介质密度；$ {\mathrm{\alpha }}_{\mathrm{N}} $、$ {\mathrm{\alpha }}_{\mathrm{T}} $分别为黏弹性边界的法向、切向修正系数，$ {\mathrm{\alpha }}_{\mathrm{T}} $取值区间为0.35～0.65，推荐值为0.5，$ {\mathrm{\alpha }}_{\mathrm{N}} $取值区间为0.8～1.2，推荐值为1。

节点控制长度选取节点两侧单元长度的各1/2，如图5所示。

图 5  节点控制边长（二维）

Figure 5.  Node control side length（Two-dimensional）

下载: 全尺寸图片 幻灯片

1.6   模型正确性验证

依据南京工业大学阮滨博士开展的饱和细砂不排水循环三轴试验及模拟结果，验证基于Davidenkov骨架曲线的孔压增量模型的正确性和适用性。数值仿真模型参考赵丁凤等（2017）建立，如图6所示。为使圆柱体土样上下部分受到均匀荷载，在其上下部分别设置刚性块。模型参数如表2所示。子程序中涉及参数如表3所示。

图 6  模型示意

Figure 6.  Schematic diagram of the model

下载: 全尺寸图片 幻灯片

表 2  模型参数

Table 2.  Model parameter

直径D/m	高度H/m	相对密度Dr	网格类型	单元数/个	天然孔隙比	初始有效围压/kPa	动剪应力比
0.0391	0.08	0.5	C3 D8 R	96	0.885	100	1.05

下载: 导出CSV 

| 显示表格

表 3  子程序参数

Table 3.  Subprogram parameters

$ \rho $/（kg·m−3）	${V}_{{\rm{s}}}$/（m·s−1）	A/B	$ {\mathrm{\vartheta }}_{\mathrm{d}} $/$ \mathrm{\vartheta } $	$ {\gamma }_{0} $	a	C1/C2	m/n	k/b
2 120	0.08	1.02/0.36	0.49/0.444	0.000 38	0.5	0.73/0.475	0.43/0.62	0.002 5/0.5

下载: 导出CSV 

| 显示表格

孔压比$ {\mathrm{\mu }}_{\mathrm{N}}/{\mathrm{\sigma }}_{\mathrm{c}}' $=1时，判定土样发生液化失稳，即试验停止。在该仿真的ABAQUS子程序中未能进行停止设定，故采用设定较长的时间因子进行计算，通过观察后处理的数据变化情况，选择停止模拟计算。

仿真所得的孔压时程曲线如图7所示，滞回曲线及应变时程曲线如图8、图9所示。由图7可知，孔压随着时间的变化而增加，但增长趋势不同，初始时刻增长较快，随后增长速率逐渐减慢，直到破坏前（$ {\mathrm{\mu }}_{\mathrm{N}}/{\mathrm{\sigma }}_{\mathrm{c}}' $=1）5 s左右，孔压增长速率会提高。轴向应变峰值和破坏时间如表4所示，轴向应变峰值分别为3.2 %（试验）和2.9 %（仿真），破坏时间分别为70 s（试验）和71 s（仿真）。综合对比不排水循环三轴试验与仿真结果可以发现，仿真与试验在前10 s的发展趋势一致；在64 s之后，试验试样在6 s内发生突变，仿真结果不会产生较大突变，将发生较缓慢的破坏；在10～64 s之间，仿真得到的孔压增长速率大于试验结果，这是由于目前的本构子程序无法实现滞回曲线捏缩，导致其应变路径在一定程度上发生变化，进而影响孔压的累积。仿真的试样破坏时间与试验测得的破坏时间误差在2 %左右。

图 7  试验与仿真孔压时程曲线对比

Figure 7.  Comparison of experimental and simulated pore pressure time-histories

下载: 全尺寸图片 幻灯片

图 8  滞回曲线

Figure 8.  Hysteresis curves

下载: 全尺寸图片 幻灯片

图 9  应变时程曲线

Figure 9.  Strain time-history curves

下载: 全尺寸图片 幻灯片

表 4  仿真与试验对比

Table 4.  Comparison between simulation and experiment

对比结果	试验	仿真
轴向应变峰/%	3.2	2.9
破坏时间/s	70	71.4

下载: 导出CSV 

| 显示表格

2.   结果与讨论

2.1   地震波强度对液化区分布的影响

盆地基岩面坡度为30 °时，不同幅值地震波作用下的孔压比分布如图10所示。随着地震波幅值的增大，液化区逐渐向深度方向发展，地表处土层塑性变形逐渐增大，主要表现为水平向的应变累积，产生了永久变形。液化区的发展存在一定规律，即盆地中心及近盆地边缘附近的孔压比发展较快，液化现象较严重。同时，观察阪神地震波0.1、0.2 g幅值作用下的孔压比发展云图，发现土层液化沿着深度方向不是逐层发生的，而是出现了非层状连续的液化分布。但这种现象并未表现在其他地震波作用的场地中，也说明了地震波频谱特性差异会影响场地液化分布情况。

图 10  基岩面坡度为30°时孔压比分布

Figure 10.  Distribution of pore pressure ratio with the slope of the bedrock surface of 30°

下载: 全尺寸图片 幻灯片

盆地基岩面坡度为45 °时，不同幅值地震波作用下的盆地内部孔压比分布如图11所示。基岩面坡度为45 °时，液化区的分布与30 °时基本一致，但仍存在一定区别。同等强度地震波作用下，坡度45 °时的液化区分布更广，塑性变形区更大，产生的永久变形更严重。阪神地震波作用下，场地沿深度方向的非连续液化分布减弱，仅在幅值为0.3 g时，盆地内部地表中心处以下20 m附近出现非液化区，但其孔压比已达到0.6～0.8。总体上看，在盆地内部表现出更显著的聚焦效应，液化区首先出现在盆地内部中心处附近，随着地震波强度的增大，逐渐向深度方向发展，但在0.1、0.2、0.3 g地震波幅值作用下，基岩面坡度处一部分土体均未发生液化现象。

图 11  基岩面坡度为45°时孔压比分布

Figure 11.  Distribution map of pore pressure ratio with the slope of the bedrock surface of 45°

下载: 全尺寸图片 幻灯片

盆地基岩面坡度为60 °时，不同幅值地震波作用下的盆地内部孔压比分布如图12所示。基岩面坡度为60 °时，盆地内部地表液化现象更显著，塑性区范围显著增大，且沿深度方向液化区域增长更迅速，盆地内部地表中心附近出现了不同程度的凹陷，盆地内边缘附近出现了凸起。同时，综合不同坡度的液化区分布云图可知，随着盆地基岩面坡度的增大，盆地内最不利位置有从边缘处向中心移动的态势。因为频谱有较大差异，盆地效应对入射波的频谱敏感，造成图10～图12中4条不同地震波作用下孔压分布出现差异。分析过程中模型内部存在反射应力波，黏弹性边界不能完全吸收反射应力波，网格划分不是完全对称，故网格没有细化到捕捉全频段地震波，导致对称部分节点处的结果存在差异。

图 12  基岩面坡度为60°时孔压比分布

Figure 12.  Distribution map of pore pressure ratio with the slope of the bedrock surface of 60°

下载: 全尺寸图片 幻灯片

2.2   盆地基岩面坡度对地表地震动特征的影响

2.2.1   基岩面坡度对加速度反应谱的影响

为对比基岩面坡度改变时加速度反应谱变化规律，将不同测点、不同基岩面坡度下的反应谱曲线进行了归纳分析，如图13、图14所示。

图 13  测点水平向反应谱（集集WGK0.2 g，5%阻尼比）

Figure 13.  Response spectrum of measurement points ①~⑦ in horizontal direction（Jiji WGK0.2 g, 5% damping ratio）

下载: 全尺寸图片 幻灯片

图 14  测点竖向加速度反应谱（集集WGK0.2 g，5%阻尼比）

Figure 14.  Measurement points ①~⑦ vertical direction acceleration response spectrum (Jiji WGK0.2 g, 5% damping ratio)

下载: 全尺寸图片 幻灯片

（1）水平向加速度反应谱

由图13可知，盆地基岩面坡度的变化会对地震动频谱特性产生一定影响。以测点①为例，基岩面坡度为30 °、45 °、60 °时，反应谱最大值分别为5.06、3.84、4.69 m/s²。不同周期情况下，反应谱存在一定规律，在周期0～0.6 s内，坡度30 °反应谱数值普遍大；在周期0.6～1.5 s内，坡度60 °反应谱数值普遍大；在周期1.5～4 s内，坡度45 °反应谱数值普遍大。盆地中心测点④加速度反应谱较好地反映了盆地聚焦效应与基岩面坡度的关系：当坡度为60 °时，反应谱在周期为3.5、4.5 s时再次出现了较大的加速度峰值。

将同一测点、不同工况的反应谱加速度峰值对应周期与边缘处对比发现，除测点④、45 °场地时的反应谱加速度峰值对应周期略微减小外，其余测点在不同角度下的反应谱加速度峰值对应周期均发生了不同程度的增大，特别是测点⑥，反应谱加速度峰值对应周期增大到了1.4 s。反应谱加速度峰值对应周期均集中在0.34～1.24 s。将不同测点的反应谱加速度峰值对应与边缘处对比发现，除测点①、30 °和测点②、45 °时的地震动加速度均放大1.04倍外，其余测点地震动加速度峰值均表现为缩小。这是因为土层的软化效应，阻尼增大，耗能增加，导致高频地震波被过滤，进而影响了峰值加速度。盆地内部所有测点反应谱曲线均表现出对长周期波显著的放大，且越靠近盆地中心处现象越明显，与孙悦等（2004）的研究结果基本一致。

（2）竖向加速度反应谱

由图14可知，除测点③外，其余各测点在不同角度下的地震动反应谱均表现出了一致规律，即随着盆地基岩面坡度的增大，反应谱加速度峰值随之增大。基岩面坡度为30 °、45 °、60 °时，测点①的反应谱加速度峰值分别为2.97、3.24、7.07 m/s²，测点②的反应谱加速度峰值分别为2.93、3.35、4.05 m/s²，测点④的反应谱加速度峰值分别为1.93、3.06、5.54 m/s²，测点⑤的反应谱加速度峰值分别为1.48、2.70、3.50 m/s²，测点⑥的反应谱加速度峰值分别为1.74、2.69、4.13 m/s²，测点⑦的反应谱加速度峰值分别为2.80、3.35、6.43 m/s²。通过反应谱加速度峰值的对比可知，除测点②外，其余测点的加速度峰值在45 °～60 °范围内变化最大，且最大增长达到2.18倍。综上所述，在地震动作用下，竖向加速度反应谱峰值随着基岩面坡度的增大而增大，且增大幅度逐级递增。

3.   结论

本文采用ABAQUS软件建立了3种不同基岩面坡度的二维T型沉积盆地模型，选取江油含增51 JYH、绵竹清平51 MZK、阪神KAK000-AT2、集集WGK-N近断层地震波，并经过0.1、0.2、0.3 g调幅后分别从基岩面输入，采用非线性有效应力法模拟了T型沉积盆地地震动效应，分析了地震波强度、盆地坡度对液化区分布及加速度反应谱、位移及速度时程曲线的变化规律。主要得出以下结论：

（1） 地表土层液化区的分布随着地震波强度的增大而增大，但因盆地的聚效效应，导致盆地内中心处的液化区发展最快。在阪神地震波作用下，盆地出现了非层状连续的液化分布。

（2） 地震波强度的增大导致地表土层发生液化的时间提前，且下层土的液化会减弱地震波的传播强度，导致上层地表土发生液化的时间延后。

（3） 地震动加速度反应谱、速度及位移时程曲线均能较好地反映盆地聚焦、放大效应。随着基岩面坡度的增大，盆地内边缘测点地震动反应减小，中心处的地震动反应增强。

（4） 反应谱加速度峰值对应周期集中在0.34～1.24 s。竖向加速度反应谱峰值随着基岩面坡度的增大而增大。