金沙江乌东德库区必油照滑坡稳定性分析

张旭 周绍武 龚维强 潘金华 蒋树 宛良朋

摘要：在详细分析金沙江水电站库岸滑坡地质结构和区域地质环境的基础上，结合地层产状、水文气象等因素，利用Ansys和Flac 3D 建立必油照滑坡地质模型，利用强度折减法与有限差分程序进行多工况下的全过程动态数值模拟。计算坡体由初始化自重平衡作用到临界滑移失稳状态全过程的位移、速率变化趋势，分析滑坡破坏模式与滑移特征，并由剪切应变增量云图确定滑移面。通过数值模拟，阐释了滑坡破坏演化特征和失稳机理。坡体塑性变形与位移分析结果表明：滑坡在天然状态下基本稳定;在连续高强度降雨状态下欠稳定，局部失稳滑塌的可能性很大，这与现场调查分析的结果一致。利用Flac 3D 进行滑坡稳定性分析和评价简单、可行，且更加直观、方便，具有传统极限平衡法所无法比拟的优势。

关 键 词：库岸滑坡; 稳定性分析; 数值模拟; Flac 3D ; 金沙江

中图法分类号： P642 文献标志码： ADOI：10.16232/j.cnki.1001-4179.2019.01.023

在水电站开发过程中，库岸滑坡稳定性问题越来越突出，尤其水库蓄水后的库岸滑坡稳定性对大坝的安全和水电站的运行至关重要。随着我国水利水电事业的快速发展，为确保水电工程建设的顺利进行，研究库水位变化与库岸滑坡稳定性的关系已成为一项重要课题[1]。

目前，滑坡稳定性分析常使用的有极限平衡法和数值计算方法[2-3]，刚体极限平衡理论的缺点是难以动态反映滑坡堆积体在外界诱发因素作用下应力场的改变以及由其引起的变形情况，数值分析法克服了边坡形状不规则和材料不均匀的限制，在求解过程中，既考虑了单元体的变形协调，又考虑了岩体的破坏准则，故其计算结果更加精确合理，结果分析更加方便直观。关于滑坡变形失稳机制的研究，目前多采用野外监测数据分析、物理模型试验以及数值模拟分析等方法[4-6]。然而利用物理模型试验研究滑坡变形破坏所需成本较高、周期较长，并且难以模拟滑坡现场的复杂环境条件以及原始应力特征;采用监测数据分析，有时会因局部变形过大而造成对滑坡整体变形规律的判断误差。数值模拟方法以其高效、成本低廉且能模拟滑坡复杂的地形地貌及应力分布特征而广受青睐。目前常用的数值分析方法主要有离散元法、有限元法和有限差分法等。陆治斌[7]、周新等[8]采用离散元法分别对滑坡过程与危岩体变形演化进行了有效模拟;苗朝采用离散元法分析了地下水的渗流规律以及静水压力和物理软化作用对岩质滑坡变形破坏的影响[9];张旭[10]、李新卫等[11]采用有限元法对滑坡渗流场演化特征及其稳定性变化规律进行了饱和–非饱和数值分析;王世梅[12]、张永昌等[3]考虑库水位升降以及不同升降速率，采用有限元法分别分析了涉水边坡、古滑坡的稳定性变化规律;闫云明考虑变形参数，采用改进的有限元强度折减法进行了边坡稳定性评价[14];安海堂[15]、陈胜伟[16]等采用有限差分法分析了开挖对坡体的变形稳定影响，张旭基于有限差分的Fish语言编写能量程序[17]，进行了大型露天矿山及其开挖的边坡稳定性与能量演化机理研究。其中基于Flac 3D 的有限差分程序有效综合了离散元法和有限元法的优点，运用动态松弛方法和有限差分方法，计算过程简单、快速，且能较好地模拟不同介质在达到其强度极限或屈服极限时，所发生的变形破坏或塑性流动等力学行为，故广泛应用于边坡稳定性分析、支护工程设计、地下洞室开挖、路基填筑以及采矿区地面沉陷的预测等多个工程地质领域。如刘晓[18]、王双等[19]采用Flac 3D 分别就高速远程滑坡与黄土滑坡的动力稳定性进行了模拟分析;王忠福[20]、杨金旺等[21]采用Flac 3D 对滑坡、高边坡的变形破坏进行了稳定性模拟分析，均取得了较好的效果。

本文在前人研究成果的基础上，以金沙江必油照库岸滑坡为例，利用强度折减法的思想，以Flac 3D 数值模拟为平台，根据有限差分原理，确定滑坡体安全系数及相应的潜在破坏模式。利用Flac 3D 内置Fish语言，基于强度折减法编制滑坡体安全系数计算程序，计算坡体在天然状况的位移、应力应变以及在降雨、水位升降等不同工况下的安全系数，模拟滑坡实际滑移面与主滑方向，综合分析变形演化规律，预测分析滑坡稳定性。

1 工程概况

必油照滑坡位于金沙江右岸必油照村西侧，下距乌东德坝址73.3 km，上距龙川江口江边乡约33.5 km。滑坡前缘至金沙江，高程约900 m，后缘至必油村平台西侧山坡一带，分布高程1 210 m左右，两侧以冲沟为界，宽约350～750 m，纵向长约1 080 m，平面面积约0.58 km2。根据钻孔揭露，滑坡前缘厚约50 m， 1 160m高程平台位置厚约156 m，后缘厚约110 m，滑坡平均厚度約60～70 m，滑坡体积约4 640万m3。滑坡区冲沟发育，切割深度一般为60～100 m，地表地形凌乱，滑坡体两侧冲沟切割形成坡度60 °的陡坎。在滑坡纵向上发育有两级平台，1 160 m高程为滑坡体平台，平台顺河向长约300 m，宽约170 m;1 360～ 1 390m高程之间为必油照村平台，宽约170～200 m，长约500 m。滑坡前缘至1 360 m高程平台之间地形坡度为20 °～30 °，1 160 m高程滑坡体平台至必油照村平台之间地形坡度约20 °～25 °，局部地形较陡，坡度约30 °～40 °，地形概貌如图1所示。

滑体主要由碎石土，碎、块石土及砂岩似基岩块体组成。碎石土主要分布在滑坡体平台一带的表层部位;块石、碎石土为滑坡体的主要物质成份;按砂岩块石中是否含石膏条带，砂岩似基岩块体可分为砂岩块体层及含石膏条带的砂岩块体层，含石膏条带的砂岩块体层主要分布在滑坡体前缘高程1 030 m以上;砂岩块体层主要分布在滑坡体的前缘及上游侧冲沟，970～1 020 m高程之间，块石主要为紫红色砂岩。滑坡体似基岩块体岩体极为破碎，大部分呈散体状，层理依稀可辨认，地表岩层产状变化大（见图1）。根据钻孔揭露的岩芯，基/覆界面上未见明显的滑带特征，但在前缘竖井井深28 m处揭露到浅层滑带，滑带土为粉质黏土夹砾。

图1 滑坡地形概貌 Fig.1 Landform of landslide

2 稳定性分析

2.1 计算方法与本构模型

采用Flac 3D 中Fish语言自编强度折减程序，对稳定性进行有效计算，其中的参数强度折减为

c′=c/F（1）

φ=arctantanφ F（2）

式中，F为强度参数折减系数;c，φ分别为有效凝聚力和内摩擦角;c′，φ′分别为折减后的有效凝聚力和内摩擦角。

本构模型为各向同性弹塑性模型，采用Mohr-Coulomb （M-C）准则，其力学模型为

f=σ1- 1+sinφ1-sinφ σ3- 2ccosφ1-sinφ （3）

式中，σ1，σ3分别为最大、最小主应力;f为屈服函数，且f>0时，材料处于塑性流动状态，f<0时，材料处于弹性变形阶段，f=0时，处于弹、塑性的临界状态。

以上即为剪切破坏判据，拉伸破坏判据为

f t=σt-σ3=0（4）

式中，σt为岩体抗拉强度。

有限差分程序中，岩体的体积模量、剪切模量分别由下式计算：

K= E3（1-2v） （5）

G= E2（1+v） （6）

式中， K、G分别为体积模量、剪切模量;E为弹性模量;ν 为泊松比。

2.2 数值模型建立

根据该区工程地质条件，建立滑坡三维地形如图2所示，其中典型的中轴线剖面二维地质图如图3所

示，剖面厚度取0.1 m，采用Ansys建立中轴剖面有限元模型，导入Flac 3D 中，进行类似二维平面应变模式的有限差分计算。滑坡形成背景为红层顺向坡，钻孔揭露基/覆面为碎石土，滑带特征不明显，前缘竖井揭露的浅层滑带土为粉质黏土夹砾，综合选定的滑坡体岩土力学参数见表1。

2.3 灾变演化机理分析

在坡体岩层分界处及坡脚处，分别布置特征测点与特征单元体进行滑坡灾变过程的有效监控，局部布置如图4所示。

天然自重工况下（折减系数 F =1.0），由总位移云图5（a）可以看出，主变形区发生在坡体上顶部，主要以自重下的沉降变形为主，中部陡坡段出现局部小变形，总体较稳定，无滑移趋势。当折减系数 F =1.14时，变形计算恰好不收敛，局部坡段塑性区贯通，位移变化明显，其中坡体中部陡坎段变形量最大且最为集中，为潜在表层失稳风险区，如图5（b）所示。故稳定性系数为 Fs=1.13。

以竖向位移为例，对比3处测点在稳定和失稳状态下的变化如图6所示，位于陡坡段的测点1从初始的天然自重平衡到失稳滑移全过程， 位移均大于平缓台阶处测点2和坡脚测点3，测点2、3位移基本保持不变，所处坡段为稳定段，测点1所在的陡坡段位移呈现无限增长趋势，即表现出局部失稳滑移现象。

对比稳定、失稳两种状态下的总速率分布如图7所示，临滑失稳时主要表现为坡体中下部浅表层运动，且坡表速率大于深部，但总体均较小，在坡表陡坡段相对集中。此时坡体最大主应力分布如图8所示，坡体上部由于受横向位移影响，浅表层局部分布有一定的小范围受拉区，若岩体较脆弱，将可能会产生拉裂缝，下部表层受中、上部的推移作用表现为受压状态。

对比典型区域的3个特征单元体最大剪应变增量变化如图9所示。位于陡坡段的A、B单元剪应变增量在破坏时增加得较为明显，且陡坎处B单元大于A单元，表明单元B所在的部位作为剪出口剪出。C单元全过程值为0，表明坡脚无剪切破坏行为。

由破坏时的滑坡最大剪应变增量云图10可以看出，在滑坡中部的坡表小范围坡度较陡的浅层部位，即特征单元体A、B所在部位，出现塑性贯通区，将会产生浅表层失稳滑移。在坡体后缘深层部位有一较大的逆冲塑性区，但未与坡表连通，将保持稳定，不会产生深部断裂或错位。

对比全过程的坡体拉、剪破坏体积变化趋势如图11所示。

剪破坏体积明显大于拉破坏体积。详细进行量化对比发现， F=1.0时，剪破坏占总破坏体积的97.53%，F=1.14时为97.75%;拉破坏在F=1.0时为2.47%，F=1.14时为2.25%（见表2）。明显地，滑坡体无论是自重稳定状态还是失稳滑移， 均表现出以 受剪为主的屈服破壞特征。由于滑坡失稳表现的是局部陡坡段的浅表层滑移，失稳时剪破坏所占的体积比相对于天然自重状态增大的不明显，拉破坏体积不变，表明未出现拉裂破坏。

由此可见，局部坡度对滑坡稳定性影响较大，可造成滑坡灾害。因此，在前期可对局部较陡坡段进行有针对性的坡度放缓处理，或局部加固，可有效防止滑移风险。

2.4 多工况稳定性计算

联合强度折减法与有限差分程序，计算其它工况下的滑坡稳定性结果见表3。滑坡体整体在库水位蓄水短暂工况条件下稳定性系数将略减小;在水库蓄水期如遇暴雨，稳定性几乎处于临界状况，处于欠稳定，浅层滑坡体可能产生变形失稳，滑坡整体稳定性较差;如在蓄水期遇地震，将会发生失稳导致水库滑坡，由此还将产生涌浪。

从现场情况来看，滑坡呈前陡后缓的形态，整体坡角约20 °，基岩面呈前缓后陡形态，前缘倾角约15 °，后缘约36 °，经地表调查，滑坡体除冲沟内见局部变形滑动现象外，其它无明显变形迹象，目前天然现状为稳定。

从以上稳定性分析计算结果来看，滑坡体在极限工况的稳定性较差，由于滑坡滑面坡角较缓（平均坡度约20 °），滑坡整体不会产生高速滑动情况，产生的局部浅层滑移引起的涌浪作用也较小。

3 结 语

本文以金沙江乌东德库区必油照滑坡为例，论述了Flac 3D 在滑坡稳定性分析评价中的应用。表明天然状况下该坡体较稳定，在水库蓄水、降雨等条件下稳定性较差，极端地震条件将有失稳风险。同时阐释了滑坡灾害的致灾演变机理。实践表明，通过Flac 3D 对滑坡稳定性进行数值模拟和分析评价是合理、可行的。该方法计算简单，可迅速求得滑坡稳定性系数，直观性强，方便模拟任何地貌、任意形状滑动面滑坡的失稳破坏过程，且可考虑岩土体的非线性大变形，模拟分析不同工况，有效克服了传统极限平衡分析和有限元法的不足。

参考文献：

[1]Du J，Yin K，Lacasse S.Displacement prediction in colluvial landslides，Three Gorges Reservoir，China[J].Landslides，2013，10（ 2）：203 - 218.

[2]郑颖人，赵尚毅.边（滑）坡工程设计中安全系数的讨论[J].岩石力学与工程学报，2006，25（9）：1937-1940.

[3]赵尚毅，郑颖人，时卫民，等.用有限元强度折减法求边坡稳定安全系数[J].岩土工程学报，2002，24（3）：343-346.

[4]卢书强，易庆林，易武，等.三峡库区树坪滑坡变形失稳机制分析[J].岩土力学，2014，35（4）：1123-1130.

[5]赵代鹏，王世梅，谈云志，等.库水升降作用下浮托减重型滑坡稳定性研究[J].岩土力学，2013，34（4）：1017-1024.

[6]唐晓松，郑颖人，唐辉明.边坡变形破坏演化特征的数值分析[J].重庆大学学报，2013，36（10）：101-113.

[7]陆治斌，蔡俊超，郑达.基于离散元的小坝滑坡演化机理研究[J].人民长江，2017， 48（S2）：116-119.

[8]周新，巨能攀，张潇敏.水平岩层软弱基座危岩体演化过程研究[J].人民长江，2015， 46（22）：40-44.

[9]苗朝，石胜伟，谢忠胜，等.红层缓倾岩质斜坡地下水作用机制及稳定性分析[J].人民长江，2016，47（18）：50-55.

[10]张旭， 谭卓英， 周春梅.库水位变化下滑坡渗流机制与稳定性分析[J].岩石力学与工程学报，2016，35（4）：713-723.

[11]李新衛，晏鄂川，李辉.暴雨作用下铜仁地区典型滑坡变形机理研究[J].人民长江， 2016，47（7）：44-49.

[12]王世梅，刘佳龙，王力，等.三峡水库库水位升降对谭家河滑坡影响分析[J].人民长江，2015，46（8）：83-86.

[13]张永昌，钱辈贝，郭飞.三峡水库水位升降作用下树坪滑坡复活机理研究[J].人民长江，2015，46（17）：57-62.

[14]闫云明，崔树军，刘见宝.强度折减法改进及在边坡稳定性评价中的应用[J].人民长江，2017，48（7）：79-82.

[15]安海堂，余会明，张卫强.土质边坡开挖试验及其数值模拟研究[J].人民长江，2016，47（19）：48-52.

[16]陈胜伟，刘先林，陈宇.岑水高速六活口滑坡稳定性分析与治理研究[J].人民长江， 2016，47（2）：27-31.

[17]张旭.基于能量的露天矿边坡灾变时空演化与多模型综合评价[D].北京：北京科技大学，2017.

[18]刘晓，唐辉明，胡新丽，等.金鼓高速远程滑坡形成机制及动力稳定性[J].岩石力学与工程学报，2012，31（12）：2527-2537.

[19]王双，裴向军，张晓超，等.黄土滑坡地震液化数值模拟分析—以宁夏石碑塬滑坡为例[J].人民长江，2015，46（11）：25-29.

[20]王忠福，刘汉东，何思明.阿海水电站进场公路K62 滑坡动态变形监测及三维数值分析[J].岩土工程学报， 2013，35（S2）：94-100.

[21]杨金旺，陈媛，张林，等.基于地质力学模型试验综合法的顺层岩质高边坡稳定性研究[J].岩石力学与工程学报，2018，37（1）：131-140.

引用本文：张 旭，周绍武，龚维强，潘金华，蒋 树，宛良朋.金沙江乌东德库区必油照滑坡稳定性分析[J].人民长江，2019，50（1）：124-129.

Stability analysis of Biyouzhao landslide on bank of reservoir on Jinsha River

ZHANG Xu  ZHOU Shaowu3， GONG Weiqiang4， PAN Jinhua5， JIANG Shu1， WAN Liangpeng3

（1. Postdoctoral Workstation， China Three Gorges Corporation，Beijing 100038， China; 2. College of Civil & Hydraulic Engineering，Tsinghua University，Beijing 100084，China; 3. China Three Gorges Construction Management Co.， Ltd， Chengdu 610000， China; 4. Chengmenshan Copper Mine of Jiangxi Copper Corporation， Jiujiang 332100， China; 5. China Gezhouba Group Three Gorges Construction Engineering Co.， Ltd， Yichang 443000， China）

Abstract：Based on the detailed analysis on geological structure and regional geological environment of landslide on the bank of reservoir on the Jinsha River， stratum occurrence， hydrometeorology and other factors， a geological model for Biyouzhao landslide was established by using ANSYS and Flac3D. Then combined with strength reduction method and finite difference program， the displacement and displacement rate of the landslide in the full process from the initial deadweight balance to critical slipping instability under multiple working conditions were calculated and the failure mode and slipping characteristics of the landslide were analyzed， and the slipping plane was determined by the shear strain incremental cloud map. Through numerical simulation， the evolution characteristics and instability mechanism of landslide failure were explained. The results of plastic deformation and displacement analysis showed that the landslide was basically stable in natural state and unstable under the condition of continuous intensive rainfall， which was consistent with the results of field investigation and analysis. Using Flac3D to analyze and evaluate the landslide stability is simple， feasible， and more intuitive and convenient than the conventional method.

Key words： reservoirs bank landslide; stability analysis; ANSYS; Flac3D; Jinsha River