引言

肯尼亚内罗毕至马拉巴新建标轨铁路一期工程（以下简称“内马铁路一期工程”）位于非洲东部，全线位于东非高原区，线路横跨埃塞俄比亚-肯尼亚裂谷肯尼亚段，是世界最大的断层陷落带，沿线地势高低起伏较大，第四系覆盖层较厚。本工程所在区域地震地质条件复杂，研究基础薄弱，铁路跨越多个不规则局部地形，需开展专题研究确定场地地震动参数。

局部不规则地形、土层分布不均匀性是影响地表地震动的重要因素，对该问题的研究一直是地震工程领域备受关注的课题。工程上通常采用一维等效线性化方法确定水平成层场地地震动参数，但对于地形起伏较大的场地，仅进行一维地震反应分析是远远不够的，需采用二维或三维模型（Kamiyama等，2002；车伟等，2008）进行研究。战吉艳等（2013）采用ABAQUS有限元软件，利用修正Martin-Seed-Davidenkov黏弹性动力本构模型，研究了横向不均匀深厚软土场地非线性地震效应；赖强林等（2014）研究了复杂V型谷地一维模型与二维模型地震动参数的差异；王永幸等（2015）和杨笑梅等（2017）分别给出了二维土层地震反应分析频域、时域等效线性化解法。上述研究均认为，对于复杂不规则或横向不均匀场地，采用二维非线性模型能得到更合理的地震动参数。基于此，针对内马铁路一期工程某典型局部不规则场地，利用ABAQUS有限元软件建立二维场地计算模型，利用Python语言进行二次开发，采用等效线性化方法进行土层地震反应分析计算，得到场地地震响应，并探讨二维不规则地形地震动分布规律，为内马铁路一期工程抗震设防提供科学依据。

1.   工程概况与地震动输入

1.1   工程概况

拟分析的局部不规则地形工程场地位于东非大裂谷西翼山区，第四系覆盖层较厚，达80—100m，基岩以第四系玄武岩、凝灰岩为主。根据已有工程地质资料，确定工程地质剖面图如图 1所示，水平向长90m，竖向深100m。

图 1  工程地质剖面图

Figure 1.  Engineering geological profile

下载: 全尺寸图片 幻灯片

将工程场地各钻孔土层资料进行整理，可知该场地土体类别主要为黏土和细砂，钻孔ZK1与ZK2土体类型和各项参数（陆鸣等，2016）分别见表 1和表 2，其中土层沿浅部至深部依次进行编号，动三轴序号为10的土体为计算基底，其剪切模量比通常取1.0，阻尼比取0.05。

表 1  ZK1钻孔土层分布资料

Table 1.  Soil layer distribution data of ZK1 borehole

土层编号	土性描述	动三轴序号	层厚/m	波速Vs/m·s-1	密度/kg·m-3
1	黏土	1	3.5	191	2100
2	细砂	6	4.5	171	2050
3	细砂	6	4.5	182	2050
4	细砂	6	4.8	223	2050
5	黏土	2	4.2	249	2100
6	黏土	2	4.2	277	2100
7	细砂	7	4.6	284	2050
8	细砂	7	5.1	345	2050
9	黏土	3	6.0	341	2100
10	黏土	3	6.0	373	2100
11	黏土	3	6.0	409	2100
12	细砂	8	6.8	434	2050
13	黏土	4	2.5	445	2100
14	细砂	9	5.5	455	2050
15	细砂	9	5.5	475	2050
16	黏土	5	4.7	508	2100
17	黏土	5	5.0	534	2100
18	黏土	5	5.0	568	2100
19	计算基底	10	—	584	2200

下载: 导出CSV 

| 显示表格

表 2  ZK2钻孔土层分布资料

Table 2.  Soil layer distribution data of ZK2 borehole

土层编号	土性描述	动三轴序号	层厚/m	波速Vs/m·s-1	密度/kg·m-3
1	黏土	1	4.0	192	2100
2	黏土	1	5.0	209	2100
3	黏土	1	5.0	232	2100
4	细砂	6	3.5	253	2050
5	细砂	6	4.0	274	2050
6	黏土	2	4.0	275	2100
7	黏土	2	4.5	289	2100
8	细砂	7	5.5	345	2050
9	黏土	3	4.5	337	2100
10	黏土	3	4.8	362	2100
11	细砂	8	4.0	382	2050
12	细砂	8	4.0	418	2050
13	黏土	4	5.0	441	2100
14	黏土	4	5.0	477	2100
15	细砂	9	5.5	507	2050
16	细砂	9	5.0	523	2050
17	黏土	5	5.5	547	2100
18	黏土	5	6.0	569	2100
19	计算基底	10	—	590	2200

下载: 导出CSV 

| 显示表格

场地9种不同层位土体剪切模量和阻尼比随剪应变变化曲线（陆鸣等，2016）如图 2所示，图中土样编号1—9分别对应表 1和表 2中动三轴序号1—9，为不同土类不同深度土样动三轴试验结果。

图 2  9种土层剪切模量和阻尼比随剪应变变化曲线

Figure 2.  G/Gmax~γ and D~γ curves for 9 representative soil specimens

下载: 全尺寸图片 幻灯片

1.2   地震动输入

假定地震动为垂直入射，取值参考《肯尼亚内罗毕-马拉巴铁路CK0-CK120段地震危险性分析报告》（陆鸣等，2016），工程场地50年超越概率63%、10%和2%的水平向基岩加速度反应谱如图 3所示，每个概率水准下拟合3条人工地震动时程样本，拟合谱和目标谱周期控制点上的误差不超过5%，50年超越概率63%水准下的3条人工地震动时程样本曲线如图 4所示。

图 3  50年不同超越概率下水平向基岩加速度反应谱（阻尼比为5%）

Figure 3.  Acceleration response spectrum for rock site with different exceedance probabilities in 50 years (damping ratio 5%)

下载: 全尺寸图片 幻灯片

图 4  50年超越概率63%水准下的3条人工地震动时程样本曲线

Figure 4.  Acceleration response history curves for rock site with 63% probability of exceedance in 50 years

下载: 全尺寸图片 幻灯片

2.   二维等效线性化模型

2.1   二维等效线性化方法

本文采用等效线性化模型描述场地非线性特性，二维模型等效剪应变取值采用王笃国等（2016）给出的方法。二维计算模型示意如图 5所示，计算区域包括水平成层介质区和不规则地形区，计算区域两侧设置黏弹性人工边界，基岩基底为水平向。采用时域等效线性化方法进行计算，阻尼矩阵的选取采用Park等（2004）给出的完整Rayleigh阻尼公式，即[C]=α[M]+β[K]，目标频率f的选取参考马俊玲等（2013）的研究结果，采用第1阶和第3阶自振频率，阻尼系数α和β取值公式为：

图 5  二维模型计算示意图

Figure 5.  Schematic diagram for two dimensional simulation model

下载: 全尺寸图片 幻灯片

式中，ζ为土体阻尼比；ω₁=2πf₁；ω₃=2πf₃；f₁和f₃可由下式求得：

式中，n分别取值1和3；h_i为第i层土体厚度；v_si为第i层土体剪切波速；$H = \sum\limits_{i = 1}^N {{h_i}} $。

计算时对每个单元剪切模量和阻尼比与上次迭代计算结果进行对比，所有单元误差均小于5%时停止迭代，进行二维线弹性场地地震反应分析计算，得到场地地震响应。

2.2   人工边界地震动输入

计算区域两侧采用黏弹性人工边界，地震动输入时需将波动输入转换为等效荷载施加在人工边界上，并保证人工边界上的位移和应力与原自由场相同。本文采用刘晶波等（1998）给出的边界节点等效荷载计算公式：

式中，$p(t)$为边界点处等效荷载；$\sigma (x, y, z, t)$、$\dot u(x, y, z, t)$和$u(x, y, z, t)$分别为边界点处自由场求解得到的应力、速度和位移；c和k分别为黏弹性人工边界法向和切向弹簧阻尼系数，求解方法参见刘晶波等（1998）。

自由场求解采用一维等效线性化方法，一维计算模型可由计算区域两侧水平成层介质区得到。求得计算区域两侧自由场地震响应后，根据式（3）可得到等效荷载，将等效荷载分别施加于两侧人工边界，实现黏弹性人工边界地震动输入。

2.3   场地有限元模型建立

利用ABAQUS有限元软件建立模型，采用Springs/Dashpots单元对黏弹性边界进行模拟，Springs/Dashpots单元为弹簧阻尼器单元，切向刚度K_bT=0.5G/R，法向刚度K_bN=G/R，切向阻尼系数C_bT=ρv_s，法向阻尼系数C_bN=ρv_p，其中G为自由场不同土层等效线性化计算得到的等效剪切模量，R近似取不规则地形中心点到人工边界点的距离，v_s和v_p分别为自由场不同土层等效线性化计算得到的剪切波波速和压缩波波速。土体单元大部分采用四边形平面应变单元CPE4R进行模拟，极少部分采用三角形平面应变单元CPE3进行模拟，本次计算剖面单元总数为33096个，节点总数为36369，长度、宽度方向单元尺寸基本小于1m。利用Python语言进行二次开发，实现了二维等效线性化方法的迭代算法，本次计算采用隐式H.H.T法求解动力方程，有限元模型如图 6所示，图中L1、L2、L3、R1、R2、R3为不规则地形不同位置处计算点，模型尺寸为375m×88m（长×宽）。

图 6  二维场地有限元模型

Figure 6.  Two dimensional FEM site model

下载: 全尺寸图片 幻灯片

3.   计算结果与分析

在50年超越概率63%、10%和2%水准下，取每个概率水准3条时程样本开展场地土层地震反应计算（共9种工况），分析峰值加速度和频谱特性，探讨局部不规则地形地震动分布规律。

3.1   峰值加速度

6个计算点50年超越概率63%、10%和2%水准下地震动峰值加速度计算结果见表 3。

表 3  6个计算点50年超越概率63%、10%和2%水准下地震动峰值加速度（单位：gal）

Table 3.  Peak ground motions for 6 points with 63%、10% and 2% probabilities of exceedance in 50 years (Unit: gal)

计算点		工况
		50年超越概率63%			50年超越概率10%			50年超越概率2%
		样本1	样本2	样本3	样本1	样本2	样本3	样本1	样本2	样本3
左侧	L1	-58.5	57.4	-61.1	144.5	154.0	151.1	-243.0	237.9	-245.3
	L2	-61.1	-62.6	-65.6	-158.2	-165.3	164.0	255.4	-249.2	-253.8
	L3	68.3	-67.9	76.1	-172.2	-176.1	174.9	-278.6	264.7	-272.0
右侧	R1	-63.7	-64.6	-62.7	158.2	-168.4	164.6	-250.2	255.6	254.0
	R2	-67.4	72.6	-66.1	165.4	-173.7	176.3	-255.6	-268.4	-263.6
	R3	81.5	82.9	-80.7	181.7	198.6	192.4	286.4	-283.5	-288.8

下载: 导出CSV 

| 显示表格

由表 3可知，同一工况下不同位置处地震动峰值加速度存在一定变化，底部峰值加速度较小，顶部峰值加速度较大，总体呈越靠近顶部峰值加速度越大的趋势；以底部峰值加速度为基准点，50年超越概率63%、10%和2%水准下的顶部峰值加速度放大系数分别为1.24—1.30、1.14—1.19和1.11—1.15，表明随着地震动输入强度的增加，放大系数呈略减小的趋势。

3.2   频谱特性

50年超越概率63%、10%和2%水准下不同位置处地震动加速度放大系数谱如图 7所示，限于篇幅问题，本文仅给出不规则地形左侧L1、L2和L3计算点每个水准下1条样本时程计算结果。由图 7可知：（1）3个工况下的地震动加速度放大系数谱均呈双峰或多峰特性；（2）同一概率水准下不同位置处的地震动加速度放大系数谱谱形存在一定差异，总体上看，顶部计算点长周期段地震动加速度放大系数较大。

图 7  不规则地形左侧不同位置处地震动加速度放大系数谱（阻尼比为5%）（一）

Figure 7.  Acceleration amplification facor spectrum for different positions in left side of irregular topography (damping ratio 5%)

下载: 全尺寸图片 幻灯片

4.   结论

基于内马铁路一期工程某局部不规则地形工程地质资料，利用ABAQUS有限元软件建立二维场地有限元模型，利用Python语言进行二次开发，实现二维等效线性化方法的迭代算法，共进行3个概率水准、9种工况下的土层反应计算，得出以下结论：

（1）不规则地形不同位置处地震动峰值加速度存在一定差异，通常情况下顶部峰值加速度高于底部，且随着地震动输入强度的增大，加速度放大系数呈略减小的趋势，50年超越概率63%、10%和2%水准下的顶部峰值加速度放大系数分别为1.24—1.30、1.14—1.19和1.11—1.15。

（2）不规则地形不同位置处地震动放大系数谱呈双峰或多峰特性，但谱形有所不同，主要体现在长周期段放大效应不同，顶部计算点长周期段地震动加速度放大系数较大。

本文给出的结果为内马铁路一期工程确定抗震设防要求提供了科学依据，同时可为类似跨越局部不规则地形工程结构地震动参数的确定提供参考。