引言

地震往往会触发大量山体滑坡等地质灾害. 2008年5月12日发生在我国四川省汶川县的M_S8.0大地震引发了3万5000多处滑坡(黄润秋，2009)，给当地人们的生命和财产造成了严重的损害.地震作用下的边坡动力稳定性问题一直备受国内外相关专家的关注.近些年来对众多地震滑坡震害现象的实际调查和理论试验研究(祁生文等，2003；徐光兴等，2008；黄润秋，2009; 殷跃平，2009)表明，地震动的幅值、频谱和持时这3个要素对地震边坡稳定性均有着不同程度的影响作用，地震动中不同幅值、频谱和持时的组合决定了地震作用下边坡等结构物的受灾程度(胡聿贤，1981).

峰值加速度是被较早接受的表征地震动强弱的参数，是控制地震动荷载对边坡稳定性影响的主要参数之一，也是边坡地震稳定性研究中使用最多的一个参数.目前已有较多研究通过现场调查或数值模拟方法进行了关于地震动峰值加速度对边坡动力响应的影响特性研究，并取得了很多成果(Wilson，Keefer，1985；Harp，Wilson，1995；Griffiths, Lane, 1999；祁生文等，2003；徐光兴等，2008；王秀英等，2009；言兴志等，2011).由于地震动频谱和持时较峰值加速度获取要相对复杂一些，因此针对地震动频谱和持时对边坡稳定性影响特性的研究较少，尤其在震后快速评估边坡稳定性时，选用峰值加速度这一指标具有一定的优势，但事实和已有研究均表明，具有不同频谱和持时的地震动对边坡稳定性的影响作用存在显著差异，因此，亟需对这些参数对边坡动力响应的影响规律加以研究.鉴于此, 本文拟基于有限元模拟方法选取不同幅值、频谱和持时的地震动加速度时程作为荷载输入，计算分析地震动参数对土质边坡响应的影响规律，以期为边坡地震稳定性的评价和地震动参数指标选取提供参考.

1.   计算模型与分析方法

1.1   边坡模型

为研究边坡在地震作用下的动力响应特性，建立如图 1所示的边坡模型.该模型参考徐光兴等(2008)，并在其基础上，将模型尺寸扩大，以便进一步消除边界影响.有限元模型如图 1所示，长为170 m，高为70 m，边坡坡角为34°，坡高30 m，边坡坡顶后缘长度为80 m.模型材料为均质材料；边坡土体的密度为2 070 kg/m³，弹性模量为90.8 MPa，泊松比为0.3，黏聚力为13.99 kPa，内摩擦角为25°.

图  1  边坡有限元模型示意图(修改自徐光兴等，2008)

Figure  1.  Schematic diagram of soil slope finite element model (modified from Xu et al, 2008)

下载: 全尺寸图片 幻灯片

1.2   动力方程与求解方法

对于边坡地震动力稳定性问题可以建立以下基本力学运动方程：

式中：M为质量矩阵，C为阻尼矩阵，K为刚度矩阵；，，u分别为响应的加速度、速度和位移列矩阵，为输入地震动加速度列矩阵.

地震是一种完全随时间变化的复杂荷载，边坡岩土体在地震作用下往往会进入弹塑性状态，这时便无法得到解析解，解析方法也不再适用，但通过数值计算可以得到结构动力反应的数值解.常用的数值计算方法有分段解析法、中心差分法、平均常加速度法、线性加速度法、纽马克(Newmark)-β法和威尔逊(Wilson)-θ法，这些方法的核心是假定结构在每一个微小时间步长内呈线弹性反应，通过逐步积分来获取该结构在整个时域内的最终响应.在有限元模拟软件ABAQUS中有隐式和显式两种算法，其中隐式算法是以纽马克-β法为基础，而显式算法是采用中心差分法来解决动力学问题.本文采用隐式算法.

1.3   土体本构模型

对于边坡土体材料，本文的模拟分析选用理想弹塑性本构模型，屈服准则采用摩尔-库仑破坏准则.该破坏准则广泛应用于岩土材料的理论与实践中.

1.4   边界条件

在模拟分析边坡等半无限体的动力问题时，为了防止地震波传播至模型边界处反射回模型，通常需要设置人工边界，允许必要的能量发散.目前用于解决动力问题边界条件设置的方法主要有简单截断边界、黏弹性边界、透射边界等，这些方法各有利弊, 其中黏弹性边界较为简单且适用性强，可较好地用于边坡地震稳定性问题中，因此本文采用黏弹性边界(刘晶波等，2006)，其中二维黏弹性人工边界等效物理系统的弹簧系数K和阻尼系数C可分别由下式求得：

式中：K_BN和K_BT分别为弹簧法向、切向刚度；ρ为介质质量密度；α_T和α_N分别为切向、法向黏弹性人工边界参数，二者的取值范围通常分别为0.35—0.65和0.8—1.2(刘晶波等，2006)，本文取其中值，即α_T=0.5，α_N=1.0；R为波源至人工边界点的距离；c_S和c_P分别为S波和P波波速，G为介质剪切模量，可通过下式分别求得：

式中，E和υ分别为结构材料的弹性模量和泊松比.

1.5   阻尼参数的确定

阻尼的设置是模拟计算中的一个关键环节，工程中常常采用瑞雷阻尼，即通过质量矩阵和刚度矩阵的线性组合来表述：C=αM+βK，其中α为质量阻尼系数，β为刚度阻尼系数.当计算选用的两阶阻尼比与结构阻尼比相等时，α和β可由下式确定：

式中，ω₁和ω₂为结构前两阶振型的圆频率.经计算，本文中边坡模型的前两阶圆频率分别为4.96 rad/s和5.45 rad/s，阻尼比ξ一般取0.05即可满足计算要求.

2.   输入地震动

地震荷载采用输入实际地震动记录加速度时程来实现，输入至模型底部的水平方向.选取6条特征周期T_g分别集中在0.25，0.35，0.45，0.55，0.65和0.90 s的地震动记录，过滤掉20 Hz以上的高频成分，并将加速度幅值统一调整为1 m/s²，截取时长为19.98 s. 图 2给出了所记录的6条地震动加速度时程及其相应的反应谱曲线.

图  2  输入地震动加速度时程(左)和对应反应谱(右)

Figure  2.  Input acceleration time-histories (left) and corresponding response spectra (right)

(a) T_g=0.25 s; (b) T_g=0.35 s; (c) T_g=0.45 s; (d) T_g=0.55 s; (e) Tg=0.65 s; (f) T_g=0.90 s

下载: 全尺寸图片 幻灯片

为模拟分析加速度幅值大小对边坡动力响应的影响，本文选取特征周期T_g=0.45 s的地震动加速度时程作为输入，将峰值加速度分别调整至1, 2, 3和4 m/s²，持时均设定为19.98 s.

为研究持时对边坡稳定性的影响，对于一条峰值加速度绝对值为a_max的地震动时程而言，将其持时定义为加速度的绝对值a₀首次与末次超过a_max/3的持续时间(胡聿贤，1981).据此可以得到输入T_g=0.45 s的地震动加速度时程的持时为T_d＝7.5 s，如图 3所示.

图  3  T_g=0.45 s的地震动持时

Figure  3.  Duration of the ground motion when T_g=0.45 s

下载: 全尺寸图片 幻灯片

3.   计算结果分析

3.1   地震动输入以及人工边界施加效果验证

地震荷载以及合理地施加人工边界是确保计算模拟正确性的前提.为了更全面地验证地震荷载和所施加人工边界的效果，在程序完成计算后，提取图 1中模型底部观测点A, B, C的加速度时程并进行对比分析，结果如图 4所示.可知，输出与输入的加速度时程曲线吻合程度较高，较好地验证了地震动时程输入的正确性和所施加人工边界的良好效果.

图  4  图 1所示边坡底部输出点A, B, C的加速度时程曲线

Figure  4.  Acceleration time-histories calculated at the floor points A, B, C of slope as shown in Fig. 1

下载: 全尺寸图片 幻灯片

3.2   地震动峰值加速度的影响

提取各地震动输入下各观测点的加速度、变形位移响应曲线和PGA放大系数，以分析峰值加速度对边坡响应的影响特性.

分别提取变形最为显著的坡肩处P₁点的竖向变形位移时程和坡脚处P₅点的水平向变形位移曲线，结果如图 5a和图 5b所示.可以看出，在P₁和P₅点，坡体反应的变形位移均随着输入地震动峰值加速度的增大而显著增大, 且输入的峰值加速度每增加1 m/s²，坡脚处P₅点和坡肩处P₁点的变形均增加0.2 m左右，这表明峰值加速度对边坡变形具有显著影响.

图  5  不同峰值加速度地震动输入下坡肩处P₁点的竖向位移时程(a)和坡脚处P₅点的水平位移时程(b). P₁和P₅点的位置见图 1

Figure  5.  Time-histories of vertical displacement of the point P₁ at the slope shoulder (a) and horizontal displacement of the point P₅ at the slope foot (b) under different input PGAs The locations of the points P₁ and P₅ are shown in Fig. 1

下载: 全尺寸图片 幻灯片

为研究坡体对PGA的放大规律，选择坡顶V₁₁点和坡面H₁₁点(其位置见图 1)两处来分析，图 6给出了这两点在不同峰值加速度地震动输入下PGA放大作用的变化.可见，坡体对输入PGA的放大作用随着输入地震动PGA的增大呈减小的趋势.当输入峰值为1 m/s²时，坡顶V₁₁点和坡面H₁₁点两处对输入PGA分别放大了1.52倍和1.68倍；随着地震动PGA增大至4 m/s²，两观测点处的PGA放大系数分别减小至1.10倍和1.15倍.这可能是由于随着输入地震动PGA的增大，土体内部剪应变增大而导致阻尼增大，滤波作用增强，抑制了PGA的放大效果.

图  6  不同峰值加速度地震动输入下观测点V₁₁和H₁₁(位置见图 1)处的PGA放大系数

Figure  6.  PGA amplification coefficient of acceleration response at the points V₁₁ and H₁₁, whose locations are shown in Fig. 1

下载: 全尺寸图片 幻灯片

3.3   地震动频谱的影响

同样选取边坡变形最明显的坡肩处P₁点的竖向变形位移时程和坡脚处P₅点的水平向变形位移和响应作为分析对象，其变形位移曲线分别如图 7a和7b所示.

图  7  不同特征周期T_g地震动输入时坡肩处P₁点竖向位移时程(a)和坡脚处P₅点水平向位移时程(b)

Figure  7.  Time histories of vertical displacement of the point P₁ at the slope shoulder (a) and horizontal displacement of the point P₅ at the slope foot (b) under ground motion with different characteristic period T_g

下载: 全尺寸图片 幻灯片

由图 7可知，随着输入地震动特征周期T_g的增大，边坡坡脚处和坡肩处的变形位移均逐渐增大.在特征周期T_g为0.45，0.55和0.65 s时的地震动作用下的边坡变形位移响应之间的变化相对于特征周期T_g为0.25和0.35 s时的小；但当输入地震动特征周期约为0.90 s时，观测点变形位移均明显增大.

为研究地震动特征周期对边坡PGA放大作用的影响，提取坡面各观测点在不同特征周期地震动作用下响应的峰值加速度，如图 8所示.可以看出，随着输入地震动特征周期的增大，坡面各点的PGA放大系数总体呈增大的趋势.

图  8  不同特征周期T_g的地震动输入时坡面各观测点PGA放大系数的变化

Figure  8.  Variation of PGA amplification coefficient of the points at the in slope surface under different input T_g

下载: 全尺寸图片 幻灯片

以上研究结果表明，输入地震动的特征周期对边坡响应的影响也十分显著，但规律较峰值加速度的要复杂.

3.4   地震动持时的影响

由图 5可以看出:在地震刚刚开始时，地震动加速度幅值较小，坡体尚处于弹性阶段，坡体观测点变形位移在零点附近波动；随着输入地震动加速度幅值的增大，动力作用增强，坡体变形位移迅速增大；直到第12 s左右，坡体观测点的变形位移在随后的地震动持续时间里基本维持不变，此时边坡产生了永久变形.整个地震动过程中坡体的变形在地震动持时T_d阶段内变化较为显著，不断地增大，这说明坡体变形位移受地震动持时的影响也十分显著.

4.   讨论与结论

分析边坡地震稳定性时需要考虑地震动特性对其的影响作用，地震动特性主要表现为加速度时程的幅值、频谱和持时等3个要素.本文采用单因素变量法，分别选取不同单一特性的地震动作为输入来分析其对边坡响应的影响作用，结果显示地震动的各个特性对边坡稳定性均有不同程度的影响.

PGA对边坡稳定性有显著的影响.边坡坡脚处和坡肩处的变形随地震动PGA的增大而单调增大，坡顶处PGA放大系数随着输入地震动PGA的增大呈现减小的规律，这可能是由于随着输入地震动PGA的增大，土体内部剪应变增大导致阻尼增大，坡体内部阻尼呈现明显的非线性状态，滤波作用明显增强，抑制了PGA的放大效应.

另一方面，地震动的频谱特性对边坡稳定性也有明显的影响.随着输入地震动特征周期的增大，边坡坡脚处和坡肩处的位移逐渐增大，特别是当特征周期大于0.65 s时，坡体位移响应明显增大，同时坡顶的PGA放大系数也基本上呈增大趋势.然而地震动频谱特性的影响规律较峰值加速度的复杂.此外，持时对边坡稳定性也有一定的影响.

综上，地震动特性参数中PGA幅值和频谱对边坡稳定性的影响效果较为明显，在边坡地震稳定性分析评价中应作为重点考虑因素.