0.   引言

边坡稳定性一直都是岩土地质研究中经久不衰的课题，地震及降雨引发边坡坍塌、滑坡以及泥石流等边坡失稳问题尤其如此。目前，单工况下边坡稳定性的研究层出不穷，多工况耦合作用下边坡稳定性的研究则报道较少。蒋中明等^[1]通过自编FISH函数实现了边坡三维饱和渗流计算分析；刘春玲等^[2]利用FLAC3D对某边坡地震稳定性进行了研究；齐信等^[3]基于GIS技术对汶川地震强震区地震诱发滑坡与后期降雨诱发滑坡控制因子进行了耦合分析；龚文俊等^[4]采用GEOSTUDIO软件以甘肃西和Ⅲ号滑坡为例研究了降雨和地震耦合作用对滑坡稳定性的影响。边坡稳定性与人类生产生活息息相关，因此，对边坡稳定性的研究显得尤为重要。

本文以成兰铁路岷江双线特大桥高边坡为例，结合多种数值分析方法的特点，开展天然工况、地震工况、降雨工况以及降雨-地震耦合工况下边坡稳定性研究和工程安全性评价，为单工况及耦合工况下边坡稳定性研究提供参考数据及技术依据。

1.   研究区概况

成兰铁路岷江双线特大桥高边坡位于四川省九寨沟县大录乡八郎沟玉瓦寨，边坡所在的2号岩堆大致呈扇形展布于斜坡上，岩堆体纵坡较陡，坡度约在30°－38°之间，坡面长满杂树，树林茂密。该区域处于我国地貌阶梯起伏最剧烈的深切陡峻高山峡谷。铁路线路沿岷江两岸而上，在该区域穿越了4条活动断裂带，共涉及规模不等的十几条活动断裂，平均滑动速率均大于0.5 mm/a以上^[5]。边坡区域地层岩性纷杂，不良地质分布广，在天然条件下会引起坡内岩土侧移导致轻微滑坡^[6]。岷江地震带属高烈度地震强震频发区，地震基本烈度为Ⅷ度，地震峰值加速度为0.2 g^[7]。该区域内气候多变，降水量较多，易形成暴雨。

2.   模型建立

2.1   数值模拟理论基础

本文采用分阶段弹塑性求解法、容重增加法与强度折减法的结合以及非线性动力反应分析法对边坡天然工况、降雨工况、地震工况及降雨-地震耦合工况下的应力应变特征进行研究。

2.1.1   分阶段弹塑性求解法

设置摩尔-库伦模型为材料的本构模型，先设置粘聚力和抗拉强度为大值，计算至平衡后，再设置粘聚力和抗拉强度为分析所采用的值求解，生成初始的应力场^[8]。

2.1.2   强度折减法

首先将边坡强度参数粘聚力c和内摩擦角$\phi$同时除以折减系数F，得到一组新的粘聚力和内摩擦角值(c′和$\phi$′)，然后，将c′和$\phi$′值作为一组新的材料参数输入计算；不断增加F值，直到边坡达到极限状态，发生失稳破坏，这时所对应的F值即为边坡的稳定性系数^[9]。用公式表示如下：

式中：F为折减系数；c为岩土边坡最初的粘聚力，kPa；$\phi$为岩土边坡初始内摩擦角，(°)；c′、$\phi$′分别为经过折减后的粘聚力和内摩擦角。

2.1.3   容重增加法

保持岩土的粘聚力及内摩擦角不变，逐步增加重力加速度G反复进行计算(逐步增加重力加速度等同于增加岩土体的容重)，直至边坡破坏，此时采用的重力加速度G_limit与实际重力加速度(G₀，通常取值为9.8 m/s²)之比即为该边坡的安全系数^[10]，即：

2.2   模型尺寸及结构面划分

研究所采用的边坡模型(见图 1)长180 m，宽80 m，高110 m。平均坡度降约0.48。坡体中上部易滑坡区域的坡度约36°。为便于分析，沿易滑坡区域对称轴取剖面1-1，易滑坡区对称轴与正北方向呈51°角。模型共由4个结构面组成。结构面Ⅰ是砂岩、板岩夹灰岩与上层的块石土以及碎石土的分界面，倾角40°；结构面Ⅱ是碎石土与块石土的分界面，倾角16°；结构面Ⅲ是块石土与粉质黏土的分界面，倾角34°；结构面Ⅳ是粉质黏土与上层覆土的分界面，倾角38°。边坡岩土物理力学参数(见表 1)根据现场勘测所获取资料确定。

图  1  边坡模型

Figure  1.  The model of slope

下载: 全尺寸图片 幻灯片

表  1  边坡岩石力学参数

Table  1.  Natural and rainfall condition parameters

工况	岩土体类型	D/(kg·m-3)	c/kPa	$\phi$/(°)	E/GPa	μ	σt/MPa	安全系数
天然	砂岩、板岩夹灰岩	2400	3.30E+07	55.0	2.30E+10	0.26	5.00E+05	1.31875
	块石土	2200	1.05E+04	40.0	8.20E+08	0.28	0.00E+00
	碎石土	2100	1.11E+04	40.0	8.10E+08	0.23	0.00E+00
	粉质黏土	1900	2.20E+04	20.0	1.00E+08	0.35	0.00E+00
25 mm/d(降雨量)	砂岩、板岩夹灰岩	2450	2.81E+07	46.7	2.30E+10	0.26	4.25E+05	1.00125
	块石土	2320	8.93E+03	34.0	8.20E+08	0.28	0.00E+00
	碎石土	2210	9.44E+03	34.0	8.10E+08	0.23	0.00E+00
	粉质黏土	2030	1.87E+04	17.0	1.00E+08	0.35	0.00E+00
35 mm/d(降雨量)	砂岩、板岩夹灰岩	2470	2.60E+07	43.8	2.30E+10	0.26	3.93E+05	0.93245
	块石土	2340	8.29E+03	31.4	8.20E+08	0.28	0.00E+00
	碎石土	2230	8.62E+03	31.4	8.10E+08	0.23	0.00E+00
	粉质黏土	2100	1.73E+04	15.7	1.00E+08	0.35	0.00E+00
45 mm/d(降雨量)	砂岩、板岩夹灰岩	2485	2.39E+07	39.9	2.30E+10	0.26	3.73E+05	0.91225
	块石土	2357	7.58E+03	29.2	8.20E+08	0.28	0.00E+00
	碎石土	2247	8.11E+03	29.2	8.10E+08	0.23	0.00E+00
	粉质黏土	2148	1.62E+04	15.1	1.00E+08	0.35	0.00E+00
注：D－容重；c－粘聚力；$\phi$－内摩擦角；E－弹性模量；μ－泊松比；σt－抗拉强度

下载: 导出CSV 

| 显示表格

2.3   控制条件设置

2.3.1   天然工况

天然工况即自重工况。赋予摩尔-库仑破坏准则，采用分阶段弹塑性求解法生成初始地应力场，利用自编强度折减法计算得出天然工况下边坡的安全系数。

2.3.2   降雨工况

降雨为该地区的主要工况。在天然工况的基础上，采用强度折减法的思想对岩土体的抗剪强度参数进行折减，其后利用自编强度折减法计算得出不同降雨条件下的边坡安全系数。

2.3.3   地震工况

研究边坡位于岷江地震带，地震频发。采用弹塑性模型以及摩尔-库仑破坏准则。为减少变形过大对结果的影响，体系最大不平衡力设定为100，时间步长设置为1~5，采用最小临界阻尼比0.005、最小中心频率12.8的瑞利阻力。在模型底部导入峰值加速度为3 m/s²的汶川地震波，其达到峰值的时间为9.97 s，持续时间为20 s。

2.3.4   降雨-地震耦合工况

考虑研究区特殊的环境条件，综合强度折减法和容重增加法，在降雨工况的基础上施加地震动力荷载，比较分析在不同降雨量和地震共同作用下边坡的应力应变特征及其稳定性。各耦合工况设置情况见表 2。

表  2  耦合工况设置

Table  2.  Multi load-case set

工况	降雨条件/(mm·d-1)	地震条件
耦合工况一	25	汶川地震波
耦合工况二	35	汶川地震波
耦合工况三	45	汶川地震波

下载: 导出CSV 

| 显示表格

3.   各工况计算及结果分析

3.1   天然工况及降雨工况

模拟过程中，考虑不同降雨情况对边坡稳定性的影响，设置25 mm/d、35 mm/d、45 mm/d等3种降雨情况，对摩擦角、粘聚力、体积模量、剪切模量等参数分别进行0.85、0.78、0.72倍折减。同时，考虑不同降雨情况下雨水入渗的加载作用，结合容重增加法，对土体容重进行不同情况的扩大。各降雨工况岩土体参数与其安全系数见表 1。

由表 1可知，自然工况下边坡的安全系数为1.318，有较高的安全储备，边坡处于稳定状态。降雨量为25 mm/d时，安全系数为1.001，较天然工况降低24.23%，安全储备所剩无几，边坡处于临界滑坡状态，微小的不利影响即可启动滑坡。降雨量为35 mm/d时，安全系数为0.932，较天然降雨工况降低29.29%，安全储备为0，此时边坡的抗滑力小于下滑力，边坡发生破坏。降雨量为45 mm/d时，安全系数为0.912，较天然降雨工况降低30.80%，安全储备为0，较35 mm/d降雨工况，滑体产生更大的位移。

图 2显示了降雨工况下安全系数随降雨强度的变化，由图 2可知，当降雨量为0时，边坡处于稳定状态；随着降雨量的增加，边坡稳定性逐渐降低，且当降雨量小于25 mm/d时降低幅度较大，降雨量大于25 mm/d时降低幅度较小。该曲线的拟合方程为：

图  2  安全系数与降雨强度关系曲线

Figure  2.  The relationship curve between safety factors and rainfall intensity

下载: 全尺寸图片 幻灯片

图 3为各降雨工况下边坡的位移分布图。从图 3可知，不同降雨工况下的位移高值区均分布于坡腰中上部，但范围随降雨量增加逐步扩大。降雨量为25 mm/d时，最大位移值为39.3 mm(见图 3a)，位移高值区分布于高程2683~2714 m范围内。降雨量为35 mm/d时，最大位移值为341.7 mm(见图 3b)，位移高值区分布于高程2682~2719 m范围内。降雨量为45 mm/d时，最大位移值为528.1 mm(见图 3c)，位移高值区分布于高程2681~2733 m范围内。

图  3  各降雨工况下位移

Figure  3.  The shift diagram of different rainfall conditions

下载: 全尺寸图片 幻灯片

分析可知，降雨工况与天然工况相比，位移高值区分布大致相同，随降雨量增加，雨水对土体的软化和加载作用增强，位移高值区范围沿坡面逐步扩大，且边坡稳定性逐渐降低。当降雨量达到35 mm/d时，滑体位移由39.3 mm突变至341.7 mm，表征其在雨水入渗和加载的综合作用下发生局部滑坡，这与实际工程情况较吻合。

3.2   地震工况

数值模拟模型采用弹塑性模型及摩尔-库仑准则。模型阻尼采用瑞利阻尼，其最小临界阻尼比取0.005，最小中心频率取12.8。施加汶川地震动力作用，汶川地震波的峰值加速度为3 m/s²，持时为20 s。汶川地震波加速度时程曲线见图 4。

图  4  汶川地震波加速度时程

Figure  4.  The acceleration diagramin in Wenchuan seismic condition

下载: 全尺寸图片 幻灯片

在模型的关键位置布置4个监测点，具体位置(见图 5)为：A点位于高程2740 m的坡肩处，B点位于高程2723 m的上坡腰处，C点位于高程2670 m的下坡腰处，D点位于高程2627 m的坡脚处，记录各监测点的加速度、位移(见图 6、图 7)。

图  5  监测点位置分布示意图

Figure  5.  The position diagram of monitoring sites

下载: 全尺寸图片 幻灯片

图  6  各监测点加速度

Figure  6.  The acceleration diagram of the monitoring sites

下载: 全尺寸图片 幻灯片

图  7  各监测点位移

Figure  7.  The shift diagram of the monitoring sites

下载: 全尺寸图片 幻灯片

由图 6可知，边坡坡面各监测点加速度响应差异较大。坡脚D点加速度峰值最小，几乎为0；坡肩A点其次，加速度在±0.6 m/s²范围内变动；坡腰C点的加速度在±1.23 m/s²范围内变动；坡腰B点的加速度峰值最大，加速度在-3.3 m/s²~+1.4 m/s²范围内变动。

由图 7可知，坡脚D点位移几乎为0；坡肩A点在地震波持时3.5 s时开始移动，11 s时位移发生突变，14 s时位移达到最大值7.3 cm；坡腰C点在3.5 s时开始移动，12 s时位移发生突变，17 s时位移达到最大值-11.5 cm；坡腰B点在1 s时开始移动，12 s时位移发生突变，18 s时达到最大值-57 cm。

分析认为，在地震波持时18 s时各监测点位移均已达到稳定状态，表征各监测点在汶川地震作用下产生塑性变形。坡腰C点所在土层(粉质黏土)与坡腰B点所在土层(砂岩、板岩夹灰岩)在地震动力作用下的位移有极大差异，表征发生较大相对位移，即C点处土层发生局部破坏，边坡处于不安全状态。

3.3   降雨与地震耦合工况

由各耦合工况坡腰B点加速度图(见图 8)可知，从0 s开始各工况坡腰B点就有较大的加速度，其中耦合工况一最小，耦合工况二其次，耦合工况三最大。耦合工况一坡腰B点处加速度从0 s时的-1.27 m/s²锐减到0，然后在-1.1~1.08 m/s²范围内变动；耦合工况二坡腰B点处加速度从0 s时的4.71 m/s²锐减到0，然后在-2.01~3.15 m/s²范围内变动；耦合工况三坡腰B点处加速度从0 s时的-7.23 m/s²锐减到0，然后在-4.02~3.15 m/s²范围内变动。

图  8  各耦合工况坡腰B点加速度

Figure  8.  The acceleration diagram of slope waist point B in each coupling condition

下载: 全尺寸图片 幻灯片

图 9显示了各耦合工况下坡腰B点的位移。由图 9可知，耦合工况一、二、三的坡腰B点位移在0~17 s的范围内大致呈线性增加。耦合工况一坡腰B点处的位移在17.8 s时达到最大值1.9 m，之后趋于稳定；耦合工况二坡腰B点处的位移在18.3 s时达到最大值5.2 m，之后趋于稳定；耦合工况三坡腰B点处的位移在18.5 s时达到最大值13.6 m，之后趋于稳定。

图  9  各耦合工况坡腰B点位移

Figure  9.  The shift diagram of slope waist point B in each coupling condition

下载: 全尺寸图片 幻灯片

由降雨与耦合工况坡腰B点位移对比图(见图 10)可知，降雨-地震耦合工况下坡腰的位移明显比降雨工况大。降雨量为25 mm/d时，坡腰位移在降雨工况下为0.039 m，在耦合工况下为1.9 m；降雨量35 mm/d时，坡腰位移在降雨工况下为0.34 m，耦合工况下为5.2 m；降雨量45 mm/d时，坡腰位移在降雨工况下为0.53 m，在耦合工况下为13.6 m。

图  10  不同降雨与耦合工况坡腰B点位移对比

Figure  10.  The shift contrast figure of slope waist point B in different rainfall and coupling condition

下载: 全尺寸图片 幻灯片

分析认为，降雨工况使得边坡土体的容重和含水率增加，粘聚力和内摩擦角降低，地震动力作用效应得到了放大，导致在降雨-地震耦合工况下，边坡发生破坏，其位移量达到单一降雨工况和单一地震工况时的数十倍。

4.   结论

结合工程实例，运用FLAC3D有限差分软件详细地对成兰铁路岷江双线特大桥高边坡进行建模，并模拟、分析了天然、降雨、地震以及降雨-地震耦合等4种工况下边坡的位移、加速度等变化特征，对坡体的稳定性进行安全性评价。

通过天然工况和降雨工况坡体稳定性分析，得到天然工况和不同降雨工况的强度折减安全系数表以及降雨量与安全系数关系式，表明随降雨量增加，边坡稳定性逐步降低，为量化分析降雨量与滑坡的关系提供基础数据。

通过对比不同降雨工况与不同耦合工况，得到降雨工况和耦合工况的关系曲线，表明降雨和地震的共同作用，会大大削弱边坡的稳定性，使边坡稳定性较降雨和地震单独作用时要低，极其容易造成滑坡。故须对类似特殊环境的边坡采取合适的工程措施，增加其稳定性。

较详细地对成兰铁路岷江双线特大桥高边坡进行数值模拟分析，为单工况以及耦合工况下边坡的稳定性研究提供理论依据。同时，也对大巴山、龙门山、龙泉山等其他滑坡频发山区的研究和预警提供借鉴。