陶 骞，刘 超，朱志铭，杨 勇，邹祖银，2

(1.四川农业大学城乡建设学院，四川 都江堰611830;2.西南交通大学土木工程学院，四川 成都610031)

0 引言

自汶川地震以来，四川灾区最大的威胁来自于地震后的次生地质灾害——滑坡与泥石流［1～3］。2010年7月27日，雅安市汉源县万工乡集镇境内的二蛮山因受汶川大地震的波及，强降雨后岩体土壤吸收大量水分，导致大面积滑坡，致21人死亡，城镇建筑损毁搬迁，造成了重大人员及财产损失。二蛮山滑坡区域为一条高频泥石流沟［4］，在不良条件下极易再次滑坡。本文以二蛮山山体滑坡为例，在滑坡特征分析基础上，开展地震和暴雨工况下边坡稳定性分析研究，为震后高边坡安全性评价提供基础数据及相关技术依据。

近年来，数值模拟技术在边坡稳定性分析中有长足发展［5～10］，但其局限性也较明显。例如，基于刚体极限平衡理论建立的程序SLOPE/W过于依赖假设条件的真实性;采用有限元法的ANSYS软件有小变形假设，对于一般边坡工程，划分的单元网格数量庞大，计算量非常大。而新兴的快速拉格朗日法在边坡稳定性分析中可较好地取长补短，本文即采用基于该方法的FLAC3D软件对二蛮山滑坡的边坡进行分析计算。

1 研究区地质条件

1.1 地形地貌

二蛮山位于四川省雅安市汉源县万工乡集镇境内，北纬29°18'，东经102°42'，高山峡谷地貌。二蛮山西南侧槽状沟谷呈“V”字型展布，总体走向256°。其山顶海拔1967 m，万工乡集镇所在堆积扇海拔800～1000 m，滑坡滑源区顶部高程约1640 m，滑坡堆物前缘高程约960 m。二蛮山滑坡在沟谷右侧斜坡上留下的三角壁面如图1所示。

1.2 地质构造及地震烈度

二蛮山山体主要由二叠系峨眉山组 (P e)玄武岩和二叠系梁山组—阳新组 (P l—P y)灰岩构成。此外，在沟谷谷口存在第四系泥石流及崩坡积堆积扇，由泥石流堆积层和崩坡积层构成，物质组成为碎石土、块石土及粉质黏土。该区地震基本烈度Ⅶ度，地震分组第三组，加速度0.15 g。

1.3 地层岩性

二叠系峨眉山组玄武岩主要为暗紫色豆状玄武质凝灰岩或铁灰色致密块状玄武岩，岩体内柱状节理发育，斜坡浅表层风化严重，呈碎块状。梁山组以钙质、炭质黏土岩为主;阳新组以生物碎屑灰岩、白云岩为主，本研究区内主要是生物碎屑灰岩，灰岩表层风化溶蚀强烈，凹凸不平，多发育溶蚀沟槽。

1.4 水文地质条件

二蛮山三面环水，北西侧有白岩河环绕，南侧为大渡河，汉源县万工乡集镇位于大渡河中游左岸二级阶地上 (见图2)。

坡体水文条件具体如下:①孔隙潜水赋存于堆积层 (Q4sef)和崩坡积层 (Q4

c+dl)土壤，有一定保水性;碎石土土层地下水位受季节影响大，保水性较差，水量大。②基岩风化裂隙水赋存于岩石斜坡浅表层风化带中，其保水性和透水性好;大气降水是地下水的主要来源。③沟谷平均宽700 m，从源头至沟口堆积扇顶部长约1476 m，流域汇水面积约1.3 km2。

2 数值模型建立

2.1 模型基本信息

边坡采用三维模型，模型长300 m，宽300 m，平均高度170 m。坡面用31×31个关键点描述地表起伏，每个关键点距离为10 m。

图2 二蛮山地形地貌及水文条件Fig.2 Topography and hydrological conditions of Ermanshan Mountain

坡面槽状沟谷呈“V”字形展布，谷线约从模型右边界2/3处沿直线延伸到模型的左下角，平均纵坡坡降约0.55。沟谷左侧斜坡是滑体所在，坡度约为40°。为便于分析，沿滑体对称轴取剖面1-1，滑体对称轴与沟谷谷线呈西偏北14°角 (见图3)。

图3 1-1剖面位置图Fig.3 The location of 1-1 profile

2.2 内部结构面划分

模型内有Ⅰ、Ⅱ两个结构面。结构面Ⅰ是梁山组—阳新组灰岩与其上的峨眉山组玄武岩及上覆土层分界面，倾向305°，倾角40°;结构面Ⅱ是峨眉山组玄武岩与上覆土层分界面，倾向150°，倾角80°。结构面Ⅰ约从谷底切入模型，而结构面Ⅱ则通过模型底面左上角点切入模型。两个结构面将模型分为3部分，最底层的部分 (group3)是梁山组—阳新组灰岩，其上三菱台部分 (group2)是峨眉山组玄武岩，顶部大体积部分 (group1)则是上覆碎石土。各部分的物理力学参数根据现场勘查数据设置 (见表1)，模型各岩层关系如图4所示。

表1 岩层物理力学参数Table 1 Physical and mechanical parameters of rocks

图4 模型各岩层关系图Fig.4 Relationship of different parts of rock model

2.3 网格划分及边界设置

模型采用正四面体自由划分，并对近地表部分单元加密约10倍，地表某些畸形区域采用人工编辑，共划分了约62000个单元。

2.4 工况条件

考虑3种工况，即天然工况、降雨工况、地震工况，采用的计算方法是弹性求解法、强度折减法、非线性动力反应分析。

2.4.1 弹性求解法

将材料的本构模型设置为弹性模型，并将体积模量与剪切模量设置为大值，然后求解生成初始地应力场。

2.4.2 强度折减法

强度折减法中边坡稳定性安全系数 (Fs)定义为:

式中:R1——岩土体抗滑力，kN;R2——岩土体下滑力，k N。

强度折减法的要点是利用下述的公式调整岩土体的强度指标 (c和φF)，然后不断地迭加折减系数，直至其达到临界破坏，此时的折减系数即为安全系数Fs。

式中:cF——折减后的粘结力，kPa;φF——折减后的摩擦角，(°);Ftrial——折减系数。

2.4.3 非线性动力反应分析

利用周围区域真实密度得出网格节点集中质量，求解全部动力方程。在进行动力分析时，阻尼选择尤其重要，为减小计算量并且得到逼真的模拟结果，采用局部阻尼。

3 计算结果分析

3.1 控制条件设置

3.1.1 降雨工况

降雨是本案例主要工况，故在对位移、应变等图形分析的基础上，加以安全系数分析。计算时先采用弹塑性模型生成初始地应力场，再将模型所有非空区域定义为摩尔库伦模型，使用Slove Fos命令求解安全系数。

3.1.2 地震工况

采用弹塑性模型，为避免变形过大导致数据溢出，设定体系最大不平衡力为100，设置时间步长为1-6。采用局部阻尼，阻尼系数为0.314，模型导入汶川地震波，峰值加速度0.125 g，其动力加速度时程如图5所示。

图5 地震动力加速度时程图Fig.5 Earthquake acceleration velocity-time graph

3.2 各工况坡体变形特征分析

图6 为各工况边坡整体位移云图及1-1剖面位移云图。天然工况下，模型西北角变形显著，其位移高值中心在西北角坡脊上，有顺坡向下偏移趋势，天然工况总位移最大值约23.99 cm;在该脊周围，总位移量均匀降低 (见图6a、6 b)。

降雨工况下，滑体集中在北坡中部 (见图6 c、6 d)，与二蛮山滑坡情况吻合，其位移相对高值区分布于高程1549～1650 m范围内，降雨工况总位移最大值约17.15 m。与天然工况相比，高值区由西北角坡脊转移到西北坡面上，高值区范围缩小，但降雨工况滑坡位移增量1.60×101～1.71×101较天然工况位移增量2.25×10-1～2.40×10-1增大明显，与实际灾害情况比较吻合。

地震工况下，边坡整体在沟谷下端的位移较大，从剖面图 (见图6e、6f)可以看出，土层 (深蓝部分)较稳定，由土层向外，位移呈弧形扩散不断均匀递增，1-1剖面位移相对高值区主要分布于高程1590～1635 m范围内，在沟谷下端边沿达到最大值4.40×10-2。

3.3 各工况坡体稳定性分析

大量理论及工程实例分析结果表明，判断滑坡潜在滑动面或变形破坏部位，可根据其剪应变增量大小来判断，剪应变增量集中带通常是滑坡发生深部变形破坏的地带。

图6 各工况下边坡整体及1-1剖面位移云图Fig.6 Displacement of the whole slope and 1-1 profile under the various conditions

天然工况下，坡体剪应变增量峰值为1.40×101～1.45×101，集中在坡体顶部，沿基覆界面扩展延伸 (见图7a、7 b)，这是因为顶部基覆界面平缓，造成该处出现剪应变增量带。剪应变增量集中带相对较大的部分不足基覆界面20%，而锁固段占绝对的优势，控制了滑坡整体稳定性。

降雨工况下，坡体剪应变增量峰值为1.20×101～1.31×101，沿基覆界面延伸，中下部增量集中 (见图7c、7d)。与天然工况相比，剪应变增量有明显增大，且位置由顶部变为中下部。锁固段与天然工况相比明显缩小，表明在降雨工况下，边坡稳定性下降，滑动趋势明显。

地震工况下，坡体剪应变增量为2.25×10-4～2.48×10-4，其分布沿土层与峨眉山组玄武岩的结构面贯通 (见图7 e、7 f)。锁固段与天然工况及降雨工况下相比虽有明显放大，但其值降低为原值的10-4，表明在地震工况下，坡体整体松散，稳定性下降，滑动趋势愈加明显。

3.4 降雨工况下坡体稳定性

二蛮山属亚热带气候，干热少雨，雨季 (6—9月)降水量占全年70.6%。2010年7月23—25日，汉源县连续遭遇两次强降雨天气，累计降雨量达163 mm，致二蛮山沟谷内发生冲出量约5×104m3规模的泥石流灾害。强降雨过程中，泥石流对滑源区坡脚的掏蚀效应、地表水入渗饱水加载效应、坡体结构软化效应的迭合，诱发了滑坡发生。根据强降雨2 d后才发生滑坡的现象推测，雨水渗入坡体结构，致其物理性质软化可能是滑坡主要因素。

图7 各工况下边坡整体及1-1剖面剪切应变增量云图Fig.7 Shear strain increment of the whole slope and 1-1 profile under the various conditions

在模拟中，考虑了10 mm/d、25 mm/d、50 mm/d等3种降雨量情况;对剪切模量、体积模量等参数考虑了软化作用，进行了0.75、0.57、0.42倍折减;对容重考虑入渗饱水加载作用，进行了扩大，模拟计算后得到的安全系数见表2。

表2 不同强度折减安全系数Table 2 Safety coefficients under the different strength reductions

不折减时，边坡的安全系数Fs0值为1.88，安全储备很高，坡体安全，不会产生滑坡;各物理力学参数按第一组取值时，安全系数Fs1值为1.24，较折减前降低了34.04%，坡体安全储备降低，但还不至产生滑坡;各物理力学参数按第二组取值时，安全系数Fs2值为1.06，较折减前降低了43.62%，安全储备已所剩无几，坡体达到滑坡启动边缘，降雨量微增或者轻微震动可启动滑坡;各物理力学参数按第三组取值时，安全系数Fs3值为0.95，较折减前降低了49.47%，安全储备为0，坡体产生滑坡。分析可知，随着降雨量的增加，各参数软化，得到的坡体安全系数Fs1、Fs2、Fs3递减，边坡渐趋滑动。

3.5 地震工况下坡体监测点特征

2010年7月27日，二蛮山曾在第一次滑坡约30 min后，发生二次滑坡。二蛮山沟谷原为高频泥石流沟，滑坡区极易发生二次滑坡，结合地震灾区至今仍有余震发生的情况，对该坡地震工况下的动力特性分析尤为重要。

为研究斜坡的动力反应规律，在模型1-1剖面高程1721 m和1623 m处分别设置1个监测点，在高程1589 m处设置2个监测点，其编号依次为:P1、P2、P3、P4(见图8)，用以监测各测点加速度、速度、位移，及单元体的最大主应力、最小主应力等特征值。

图8 监测点布置示意图Fig.8 Layout of monitoring points

3.5.1 速度时程分析

图9为模拟得到的各监测点动力特征图，由图9 a可知，P2、P3点速度时程相似，0.098 s时，速度同时达最小值，分别为-0.38 m/s、-0.21 m/s，随后正向增加，0.500 s时速度达0.90 m/s，其后缓和增长至1.20 m/s。

P1点速度时程与P2、P3点相似，但速度时程波动频率降低，0.62 s时，速度绝对值最大约1.63 m/s，其后发展趋于平缓。

P4点速度时程远异于其他监测点，2.00 s内处于负值区，曲线波动明显、离散性大，0.73 s时，其速度绝对值达到最大值，为2.62 m/s。

分析表明，在远离地表的前提下，速度时程随高程增加产生滞后效应;随着距地表距离减少，速度时程波幅增大，离散性增大。

3.5.2 位移时程分析

由图9b可知，P2、P3点竖直位移时程接近，在0.20 s内显下降趋势，随后呈近似一次函数上升，P3点竖直位移在整个持时中均略高于P2点。

P1点竖直位移时程与P2、P3点相似，但略滞后于P2、P3点约0.33 s，随持时发展，P2、P3与P1点呈渐远趋势。

P4点位移时程远异于其他监测点，2.00 s内处于负值区，呈近似一次函数下降。2.00 s时，P1、P2、P3、P4点位移分别达到0.13 m、0.18 m、0.8 m、-0.33 m。

图9 各监测点动力特征曲线Fig.9 Dynamic characteristics of the monitoring points

分析表明，位移时程随高程增加产生滞后效应;随着距地表距离的减少，位移时程趋于线性;坡脚处的P4位移为负增长，表明地震作用使其受压;而其余监测点位移为正增长，所在区域均出现拉应力，有滑动趋势。

3.5.3 剪切应变增量时程分析

由图9c可知，P3点剪切应变增量时程保持在P2点下方约0.30μm，在0.05～0.50 s持时内，曲线减速递增，0.50 s时曲线斜率接近于0，其后约0.26 s持时内有一上升台阶，在1.00 s后近似零增长趋势。

在0.48 s持时以内，P1点剪切应变增量显著上升，至0.48 s时达到5.79μm，在0.48～1.35 s内，剪切应变增量在4.42～5.79μm范围波动，其后围绕5.45μm上下波动。

P4点剪切应变增量时程与P1点相似，但变化频率更快，离散性更大，在0.12 s迅速上升至4.49μm，其后，波幅渐小，最终围绕4.62μm上下波动。

分析表明，在一定范围内，随着高程的增加，各监测点剪切应变增量随之增大且趋势相近;随着与地表距离的减小，各监测点剪切应变增量波动频率加快，波幅增加。

4 结论

本文运用FLAC3D比较详细地对汉源二蛮山滑坡进行建模分析，并模拟了天然、降雨、地震3种不同工况。为使模拟更真实，采用了不同本构模型，得到了相应应力云图、应变云图、位移云图、系统不平衡力曲线、地震加速度时程等图形，并对坡体的稳定性进行评价分析。

通过降雨工况下坡体稳定性分析，得到强度折减安全系数表，为量化分析降雨量与滑坡启动关系分析提供基础。

假设地震工况得到的各监测点动力特性，模拟了汉源二蛮山二次震发滑坡的发生。

本文直观地反应了汉源二蛮山滑坡产生的过程;详细地对于二蛮山坡体地震工况进行模拟，为二次震发滑坡提供理论依据和预警，研究成果也可为龙门山、大小凉山、龙泉山、大巴山、金沙江峡谷、岷江上游段峡谷地区等其他地震滑坡频发区的滑坡预警提供借鉴。

［1］ 张永双，雷伟志，石菊松，等.四川5·12地震次生地质灾害的基本特征初析 ［J］.地质力学学报，2008，14(2):109～116.ZHANG Yong-shuang，LEI Wei-zhi，SHI Ju-song，et al.General characteristics of 5.12 earthquake-induced geo-hazards in Sichuan［J］.Journal of Geomechanics，2008，14(2):109～116.

［2］ 苏生瑞，张永双，郝莉莉，等.汶川地震断裂带附近建筑物的灾害效应模拟研究 ［J］.地质力学学报，2011，17(4):381～387，409.SU Sheng-Rui，ZHANG Yong-sheng，HAO Li-li，et al.Numerical study on the hazard effect of a building near the fault produced by the Wenchuan earthquake［J］.Journal of Geomechanics，2011，17(4):381 ～387，409.

［3］ 甘建军，黄润秋，李前银，等.都江堰—汶川公路汶川地震次生地质灾害主要特征和形成机理 ［J］.地质力学学报，2010，16(2):146～158.GAN Jian-jun，HUANG Run-qiu，LI Qian-yin，et al.Formation mechanism of geo-hazards triggered by Wenchuan M s 8.0 earthquake along Dujiangyan-Wenchuan highway［J］.Journal of Geomechanics，2010，16(2):146 ～158.

［4］ 许强，董秀军，邓茂林，等.2010年7.27四川汉源二蛮山滑坡-碎屑流特征与成因机理研究 ［J］.工程地质学报，2010，18(5):609～622.XU Qing，DONG Xiu-jun，DENG Mao-lin，et al.The Ermanshan rock slide-debris flow of July 27，2010 in Hanyuan，Sichuan:Characteristics and failure mechanism ［J］.Journal of Engineering Geology，2010，18(5):609～622.

［5］ 刘春玲，祁生文，童立强，等.利用 FLAC3D分析某边坡地震稳定性 ［J］.岩石力学与工程学报，2004，23(16):2730～2733.LIU Chun-ling，QI Sheng-wen，TONG Li-qiang，et al.Stability analysis slope under earthquake with FLAC3D［J］.Chinese Journal of Rock Mechanics and Engineering，2004，23(16):2730～2733.

［6］ 何专，姚令侃.地震和暴雨工况下边 (滑)坡稳定性分析方法评价 ［J］.水土保持通报，2009，29(3):178～181.HE Zhuan，YAO Ling-kan.Assessment on methods of analyzing slope(landslide)stability under earthquake and rainstorm［J］.Bulletin of Soil and Water Conservation，2009，29(3):178 ～181.

［7］ 陈星，李建林，王家成.基于FLAC3D的边坡地震反应分析 ［J］.中国农村水利水电，2010，(1):58～61.CHEN Xing，LI Jian-lin，WANG Jia-cheng.Seismic response analysis of slope based on FLAC3D［J］.China Rural Water and Hydropower，2010，(1):58 ～61.

［8］ 陈昌凯，阮永芬，熊恩来.地震作用下边坡稳定的动力分析方法 ［J］.地下空间与工程学报，2005，1(6):1054～1057.CHEN Chang-kai，RUAN Yong-fen，XIONG En-lai.Dynamic analysis method of slope stability in earthquake function［J］.Chinese Journal of Underground Space and Engineering，2005，1(6):1054 ～1057.

［9］ 史丹，陈蕴生，韩信，等.岩质边坡地震动力稳定性研究进展 ［J］.西北农林科技大学学报:自然科学版，2006，34(2):147～150.SHI Dan，CHEN Yun-sheng，HAN Xin，et al.Seismic dynamic stability of rock slope research advance ［J］.Journal of Northwest Sci-Tech University of Agriculture and Forest:Natural Science Edition，2006，34(2):147～150.

［10］ 常士骠，张苏民.工程地质手册 (第四版)［M］.北京:中国建筑工业出版社，2006:167～169.CHANG Shi-biao，ZHANG Su-min.Engineering geological manual(fourth edition)［M］.Beijing:China Architecture＆Building Press，2006:167 ～169.