引言

我国是一个地震多发的国家，尤其西部地区地震不断，且我国西部地形起伏多变，多高山峡谷。大量的强震动记录和震害调查表明，相较于平坦地形，复杂不规则地形对地震动影响更大（张建毅等，2012；郭明珠等，2013）。研究地形效应最直观、最有效的方法是实际场地强震动观测及进行实地震害调查和实测数据的分析整理。但强震动观测受制于小概率地震事件，在科学研究和应用上存在一定的局限性。基于此，学者们开展了解析研究和数值模拟来研究地形效应。解析法主要针对简单模型，如圆弧形（Trifunac，1972）、椭圆形（Wong，1982）、Ｕ形（Gao et al，2012）及Ｖ形（Zhang et al，2012）地形，通过建立简化力学模型并利用波函数展开法进行求解，从本质上解释不规则地形对地震动的散射问题，同时也能为数值模拟精确度的验证提供依据；而数值模拟对复杂地形的研究更具普适性，可弥补解析法在求解复杂地形时的不足，是目前研究地震动在局部场地中的传播及散射的主要技术手段，以有限差分法和有限元法为代表。Ning等（2018）基于有限差分法提出波场分离方法，模拟瑞雷波在二维河谷地形上的传播；李平等（2018）采用二维显式有限差分法研究梯形河谷场地深宽比、河谷坡角、覆盖层厚度及输入地震动强度对场地地震动的影响，并探究了各因素的影响程度；Jahromi和Karkhaneh （2019）用ABAQUS有限元软件研究了不规则地形之间的相互作用对地震响应的影响；孙纬宇等（2019）基于黏弹性边界结合有限元法，研究了SV波入射角度和河谷斜坡坡度对河谷地形地震动放大系数分布特征的影响。为了使数值模拟的结果更精确可靠，要准确地模拟真实场地，科学合理地对其设置人工边界。目前使用的局部人工边界主要有透射边界、黏性边界和黏弹性边界。其中黏弹性边界克服了黏性边界引起的低频漂移问题，能吸收边界上散射波的能量，较好地模拟无限地基的弹性恢复能力和辐射阻尼效应，在一些大型有限元软件中广泛应用（孙纬宇等，2016；马笙杰等，2020）。梁建文等（2014）采用黏弹性边界和等效结点力输入方法，对层状场地中三维凹陷地形进行地震响应分析；Chen等（2015）通过黏弹性人工边界条件，进行二维精细化非线性有限元建模，研究二维复合地形场地的地震效应特征；章小龙等（2017）运用结合黏弹性人工边界的显式有限元法，通过三维转二维的简化方法实现了SV波垂直入射下复杂三维地形地震反应的数值模拟。

地震动在河谷地形中的传播是一个极其复杂的过程，会对跨河谷结构产生较大的不利影响，所以除了从时域角度分析，还应就其频谱特性作进一步研究。王国新等（2013）通过对实际台阵记录的分析，发现地震动相干性随深度的增加而降低；Zerva和Stephenson （2011）以功率谱密度和相干性为评估指标探讨了实际冲积河谷场地的地震动特性；Imtiaz等（2018）通过滞后相干性，研究了场地几何不规则性对地震动的影响。目前，关于从时域和频域角度综合分析三维河谷地形对地震动影响的研究尚不充分，仍需作进一步研究。

本文基于西部某桥梁工程场地，利用ABAQUS有限元软件建立三维河谷地形，拟采用黏弹性人工边界结合隐式动力有限元法，进行动力响应分析，继而通过频谱分析，从时域和频域角度分别探讨河谷地形对地震动的影响，以期为开展河谷场地上桥梁三维地震反应分析提供非一致地震动输入，并进行三维河谷地形分析，以便确定因河谷地形产生的横桥向地震动。

1.   计算模型及计算方法

1.1   河谷场地模型及入射地震动

本文研究对象为西部地区一河谷地形，该地形属于半无限空间数值求解问题，通过人工边界的引入将其转化为有限域求解问题，形成几何尺寸为280 m×50 m×120 m的有限计算区域，剖面几何尺寸如图1所示。在此基础上，基于有限元计算精度要求并结合数值分析截止频率，确定六面体有限单元的几何尺寸为2 m×2 m×2 m，继而对分析区域进行有限单元离散，并建立三维有限元分析模型。之后通过动力有限元数值分析方法建立关于内节点的常微分方程组，并利用ABAQUS分析软件求解器实现其数值积分，其数值积分为隐式、非条件稳定的。假定该河谷场地为一均质且各向同性的弹性体，其密度为2 100 kg/m³，泊松比为0.35，弹性模量为860 MPa，剪切波速为390 m/s²。在河谷地形表面设置15个测点（图1），其中测点3与4，测点4与8，测点8与12，测点12与13之间的水平距离皆为30 m。

图  1  河谷地形剖面图

Figure  1.  Section of valley topography

下载: 全尺寸图片 幻灯片

入射地震动选用一宽0.2 s、持时1 s的近似脉冲，其位移时程曲线及傅里叶谱如图2所示，截止频率约为12 Hz，有限元计算时间步长取0.002 s。地震动从底部基岩处垂直入射，通过自由场从人工边界面输入，沿河谷地形水平x向振动。

图  2  输入脉冲位移时程（a）及其傅里叶振幅谱（b）

Figure  2.  Displacement time history （a） of input pulse and its Fourier amplitude spectrum （b）

下载: 全尺寸图片 幻灯片

1.2   三维黏弹性人工边界

在三维动力有限元分析中，为了消除人工边界的影响，在底边界和四周侧边界采用黏弹性人工边界条件模拟散射场自人工边界向计算区域外的辐射效应，实质是在人工截断边界上设置一系列并联的线性弹簧和黏滞阻尼器，如图3所示，其参数由围岩材料决定，具体计算公式（马笙杰等，2020）为：

图  3  黏弹性人工边界等效弹簧-阻尼器系统

Figure  3.  Viscous-elastic artificial boundary equivalent spring-damper system

下载: 全尺寸图片 幻灯片

式中：K_BN和K_BT分别为法向和切向弹簧的刚度系数；C_BN和C_BT分别为法向和切向阻尼器的阻尼系数；$ {\alpha _{\text{N}}} $和$ {\alpha _{\text{T}}} $为修正系数，本文考虑三维黏弹性边界情况，取$ {\alpha _{\text{N}}} $＝1.33，$ {\alpha _{\text{T}}} $＝0.67；R为波源至人工边界点的距离；ρ为介质密度，G为介质的剪切模量，C_P和C_S分别为P波和S波波速。

1.3   地震动输入

刘晶波和吕彦东（1998）将地震动转化为等效节点力来实现地震动输入，根据力的平衡方程可得到在模型边界点上需要施加的等效荷载，计算公式为：

式中：$ {F}_{\mathrm{B}} （ t ） $为边界节点处的等效节点力；${\sigma ^0} （ {{x_{\text{B}}}， {y_{\text{B}}}， t} ） $，${u^0} （ {{x_{\text{B}}}， {y_{\text{B}}}， t} ） $和${\dot u^0} （ {{x_{\text{B}}}， {y_{\text{B}}}， t} ） $分别为原自由场在边界处的应力、位移和速度，其中速度和应力可以根据边界处位移求得；K_B和C_B分别为黏弹性边界上的弹簧刚度系数和阻尼器的阻尼系数；A为节点控制面积。

在此基础上，刘晶波等（2018）提出新的地震动输入法，新方法可避免计算每一个边界节点处的自由场应力及由边界单元引起的附加荷载，也无需根据人工边界面的外法线确定等效节点力的作用方向，可简化等效节点力的计算过程。原有波动法在输入等效节点力后，有限元节点的动力方程为：

式中，M，C和K分别为模型质量矩阵、阻尼矩阵和刚度矩阵；$ {\boldsymbol{u}} $，$ \dot {\boldsymbol{u}} $和$ \ddot {\boldsymbol{u}} $分别为模型节点位移、速度和加速度向量；$ {{\boldsymbol{F}}} $为输入的等效节点力向量。输入等效节点力后，模型的位移场与自由波场一致，根据自由波场运动方程可得：

式中，上标0表示自由波场，若已知自由场，可根据式（5）计算等效节点力。原波动法需要对整个自由场模型进行计算，现仅需对子结构进行计算，即

式中，下标B，C分别表示边界单元节点和相邻土体节点。式（6）表明，只要保证边界单元节点和相邻土体节点的运动与自由场一致，则人工边界节点上的外荷载即为所需的等效节点力。

等效地震荷载输入的具体实施步骤为：① 根据三维河谷地形模型的空间尺寸建立三维人工边界子结构模型，如图4所示，对子结构模型最外侧节点设置弹簧—阻尼器；② 利用一维计算程序算得自由场位移时程u⁰；③ 对子结构模型所有节点输入相应的u⁰，进行子结构动力反应分析，计算求得人工边界节点处的反力F_B即为等效地震荷载；④ 用代码自动读取和写入大批量节点时程数据，施加到三维河谷地形模型中，实现地震动输入。

图  4  黏弹性人工边界三维子结构模型示意图

Figure  4.  The diagram of three-dimensional substructure model of viscous-elastic artificial boundary

下载: 全尺寸图片 幻灯片

1.4   地震动输入方法的精度验证

以三维弹性均匀半空间模型为例，验证黏弹性人工边界子结构地震动输入方法的可靠性。首先通过引入人工边界截取长、宽、高分别为60 m，60 m和30 m的有限计算区域，其密度为2 000 kg/m³，剪切波速为200 m/s，泊松比为0.3。结合有限元计算精度要求，确定有限离散单元大小为2 m×2 m×2 m，以此有限离散单元对计算区域进行离散，形成如图5a所示的弹性均匀半空间有限元离散模型。针对离散后的有限元模型建立动力运动方程，并利用ABAQUS软件中的ABAQUS/Standard模块求解计算区域内节点的地震反应。对于人工边界上节点的地震反应则采用黏弹性人工边界子结构法计算得到，子结构模型如图5b所示。入射地震动选用宽0.3 s、持时1 s的单峰脉冲地震动，其位移幅值为1 m，位移时程及其傅里叶振幅谱如图6所示。

图  5  弹性半空间三维有限元模型示意图

（a） 三维有限元模型；（b） 黏弹性人工边界子结构三维有限元模型

Figure  5.  Schematic diagram of elastic half-space three-dimensional finite element model

（a） Three-dimensional finite element model；（b） Three-dimensional finite element model of viscous-elastic artificial boundary substructure

下载: 全尺寸图片 幻灯片

图  6  输入脉冲位移时程（a）及其傅里叶振幅谱（b）

Figure  6.  Displacement time history of input pulse （a） and its Fourier amplitude spectrum （b）

下载: 全尺寸图片 幻灯片

计算得到垂直入射脉冲地震动作用下，三维弹性均匀半空间有限元模型的位移反应，图7给出了地表测点和底边界测点的位移时程。由图7可知，地表位移峰值为2 m，底边界位移峰值为1 m。由弹性均匀半空间波动理论可知，自由表面位移幅值应为入射地震动的两倍。同时在底边界处的上行波和下行波的幅值均为1 m，表明人工边界具有足够的精度，未对位移反应产生影响。这两方面验证了人工边界子结构地震动输入方法的可靠性。

图  7  模型顶面和底面位移时程

Figure  7.  Displacement time history of top and bottom surfaces of the model

下载: 全尺寸图片 幻灯片

2.   有限元数值模拟结果分析

2.1   地形效应时域特性分析

本文采用有限元法结合黏弹性边界，通过时域动力有限元方法研究河谷场地地震动.数值模拟得到了各测点的位移反应，以此为基础讨论了各测点地震动反应峰值和质点运动轨迹。

2.1.1   峰值结果分析

由于河谷场地水平x向即地面运动主方向有更大的地震响应，此处仅针对水平x向的响应进行讨论。图8为河谷地形表面各测点处的水平x向位移时程曲线，图9为各测点水平x向的位移峰值较输入地震动的放大系数图。由图8和图9可知，由于地形的不对称，河谷两边地震响应存在差异；河谷两岸坡顶处位移峰值放大显著，向两岸逐渐衰减至输入地震动的两倍；随坡高的增加，位移峰值呈增大趋势；测点4和测点12距测点8的水平距离相同，但放大系数存在明显差异，结合图1可知，这是由左坡下坡角大于右坡下坡角所致，即坡度越大，位移峰值放大越显著。这表明河谷地形效应引起了显著的地震动强度的差动效应。

图  8  河谷场地各测点水平x向位移时程曲线

（a） 左岸；（b） 左坡；（c） 右岸；（d） 右坡

Figure  8.  Displacement time history curve of horizontal x direction for each measuring point in valley site

（a） Left bank；（b） Left slope；（c） Right bank；（d） Right slope

下载: 全尺寸图片 幻灯片

图  9  河谷场地各测点水平x向位移峰值放大系数

Figure  9.  Amplification factor of peak ground displacement of horizontal x direction for each measuring point in valley site

下载: 全尺寸图片 幻灯片

2.1.2   质点运动轨迹

由于分析模型实为二维问题，故在数值分析时仅考虑了输入地震动水平方向（x向）及入射方向（y向）的变化，并讨论了地形引起的地震动转换效应，图10给出了河谷表面一些测点的质点运动轨迹。由图可知，由于地形的复杂性，地震动在该区域的散射也较为复杂，各测点除有水平x向位移外，皆有竖直y向运动分量，表明各测点均存在地震动转换；两岸靠近边界的测点主要以水平往复运动为主，而由于两岸坡顶处存在着面的相交及转化，出现较强烈的地震动转换，其运动轨迹呈近似逆进椭圆；坡上各测点侧向运动及转向运动较大；对比测点4与测点8及测点8与测点12可以看出，前者地震动转换较后者更为明显，这是由左坡下坡角较大导致的。上述分析结果表明，河谷地形效应能产生明显的地震动转换，其特征为逆进椭圆面波。

图  10  河谷场地部分测点质点运动轨迹

Figure  10.  Particle motion trajectory of some measuring points in valley site

下载: 全尺寸图片 幻灯片

2.2   地形效应频域特性分析

基于时域数值模拟结果，从频域角度进行了地形效应的相干性及相位特性的分析，即以测点1为参考点，计算了其余各测点的地震反应与测点1的相干系数及相位差。

2.2.1   相干性分析

对于两个计算点i和j，相干系数的计算公式（万珂羽，孙晓丹，2022）为：

式中，$ {{\overline S_{ii}} （ \omega ） }$，$ {{\overline S_{jj}} （ \omega ） }$，$ {{\overline S_{ij}} （ \omega ） }$分别为平滑后i，j点的自功率谱及两点间的互功率谱，本文采用汉明（Hamming）窗进行谱的平滑化，并据此计算得到了测点2至测点15相对于测点1的相干系数。相干系数是描述不同点间地震动相关程度的量，其值介于0至1之间，值越大表明两点间的相干性越强，0表示两点互不相干，1表示两点地震动过程完全一致。

图11展示了测点2至测点15与测点1之间地震动的相干系数随频率的变化关系。由图可知，两岸各测点间地震动的相干性较好，基本趋于1；坡上各测点地震动由于其复杂的散射效应，与测点1的相干性较差；随着坡高的降低，相干性逐渐降低；相干性最差的区域集中在测点7至测点10，即谷底及其周围小范围影响区域。由此可知，河谷地形效应将减小地震动的相干性。

图  11  河谷场地测点2至测点15与测点1的相干系数

（a） 左岸；（b） 右岸；（c） 左坡；（d） 右坡

Figure  11.  Coherence coefficients between measuring points 2 to 15 and measuring point 1 in valley site

（a） Left bank；（b） Right bank；（c） Left slope；（d） Right slope

下载: 全尺寸图片 幻灯片

2.2.2   相位差分析

图12为测点2至测点15与测点1之间地震动的相位差图。由图可知，两岸各测点间的地震动几乎不存在相位差；坡上各测点地震动与测点1之间的地震动存在显著的相位差，且相位差与场点位置相关，即各测点间皆存在不同的相位差，但相邻测点地震动相位差较相近。由此可以认为，河谷地形效应会导致地震动的相位差。

图  12  河谷场地测点2至测点15与测点1的相位差

（a） 左岸；（b） 右岸；（c） 左坡；（d） 右坡

Figure  12.  Phase difference between measuring points 2 to 15 and measuring point 1 in valley site

（a） Left bank；（b） Right bank；（c） Left slope；（d） Right slope

下载: 全尺寸图片 幻灯片

3.   讨论与结论

本文针对西部某桥梁工程三维河谷地形场地，采用与黏弹性人工边界相结合的时域隐式动力有限元方法，在垂直入射下分析了河谷地形对地震动的影响，得到以下主要结论：

1） 河谷地形对地震动强度有显著影响，且与河谷特征及场点位置相关，河谷底部地震动明显小于顶部和斜坡，河谷两岸坡顶处位移峰值放大最明显；随坡度的增大，位移峰值呈增长趋势。

2） 河谷地形能产生特征为近似逆进椭圆的地震动转换，且两岸坡顶处地震动转换显著；随坡度的增加，地震动转换更明显。

3） 河谷地形效应会减小地震动的相干性；河谷两岸地震动之间相干性较好，而与坡上地震动相干性较差；随边坡下降，地震动相干性逐渐降低，谷底地震动相干性最小。

4） 河谷地形效应会引发地震动的相位差，且与边坡上场点位置相关，但相邻场点相位差相近；河谷两岸地震动之间几乎不存在相位差，而边坡及谷底相位差较大。

5） 河谷地形会导致地震动明显的差动效应，即峰值差动、相位差，同时降低地震动的相干性。

本研究可为开展非一致激励下河谷场地的桥梁三维地震反应分析提供参考依据，后续可进一步进行土结相互作用研究。