利用FLAC3D预测某滑坡体地震荷载作用下的稳定性

阚 露，张发明

(河海大学 地球科学与工程学院，南京 210098)

利用FLAC3D预测某滑坡体地震荷载作用下的稳定性

阚 露，张发明

(河海大学 地球科学与工程学院，南京 210098)

以澜沧江边某滑坡工程实例为背景，选用拟静力法对地震荷载进行模拟。利用ANSYS软件建立数值模型，并导入FLAC3D软件进行计算，预测滑坡稳定性。对比结果表明:天然工况条件下，滑坡体塑性区零星分布，剪应力集中带仅位于滑坡体表面，位移场分布合理;加载地震荷载后，滑坡体塑性区上下贯通，滑动面明显，剪应力集中带深入滑坡体内部，形成完整的潜在滑动面，而且位移大幅度增加;两种工况条件下滑坡体的破坏形式均以水平滑动为主，地震工况条件下的滑坡体会出现大面积失稳。该研究为抗震方案设计提供了有益参考。

地震;滑坡;稳定系数;拟静力法

收稿日期:2013－01－08
第一作者简介:阚 露(1988－)，男，安徽省滁州人，硕士，研究方向:岩体结构与工程稳定性，E-mail:klhhuedu@126．com。

我国西南部地区沟壑纵横、山高坡陡，是中国大陆内强震活动频度最高的地区。区域内，许多重大水利水电工程正在建设中或处于论证阶段。因此，研究地震对滑坡的影响不仅具有重要的理论意义，而且具有重要的实践意义。以5.12汶川地震为例，地震诱发了1 701处滑坡，1 844处岩石崩塌，以及1 093处边坡失稳[1]，并对很多水库造成了不同程度的危害。永久位移和安全系数是评价地震作用下边坡稳定性的两个主要指标。其中安全系数是国内采用的主要指标，计算方法有拟静力法、动力时程分析法、动力有限元强度折减法[2]等。

拟静力法是一种用静力学方法近似解决动力学问题的简易方法，具有简便、易操作的特点，在工程实践中得到了广泛的应用[3]。其基本原理是，将复杂的、不断变化的地震荷载简化成一个恒定的惯性力系，然后附加在研究对象上。如何确定地震的加速度是其核心问题[4]。动力分析法虽然计算精度高，但操作复杂，现有规范中亦没有明确的计算结果评价标准[5]，而且计算过程中需要大量的地震相关参数，例如地震波时程曲线等，然而，大部分中小型工程中缺乏观测资料，难以提供相关参数，甚至忽略了对地震工况的评估。基于此，笔者以澜沧江边某滑坡工程实例为背景，利用FLAC3D软件，采用拟静力法对地震荷载进行模拟，预测滑坡体的稳定性。

1 地震荷载模拟方法与有限差分原理

1．1 地震荷载模拟方法

根据拟静力法，利用大小和方向恒定的水平荷载来模拟地震工况下复杂多变的地震荷载。通过模拟地应力场将水平荷载加载到计算模型的每一个单元上。根据DL 5073—2000《水工建筑物抗震设计规范》得出滑坡某质点i的水平向地震惯性力:

Fi=0.11gGEi，

式中:GEi——集中在质点i的重力作用标准值;

g——重力加速度。

1.2 有限差分原理

1.2.1 FLAC有限差分公式

有限差分法就是在利用数值计算方法求解偏微分方程时，用有限差分近似公式代替每一处导数，从而将求解偏微分方程的问题转化为求解代数方程的问题。

由高斯散度定理有

式中:s——封闭曲面;

ni——曲面s的单位法向量;

f——标量、向量或张量;

xi——位置向量;

ds——增量弧长;

A——表面积。

定义f在面A上的梯度平均值，

将式(2)代入式(1)后得到

对于一个三角形子单元，

式中:Δs——三角形某一条边的边长，求和即在三角形三条边上加总;

＜f＞——相应边的平均值[6]。

1.2.2 应力应变计算公式

同理可以求出

由材料的本构方程和相应的边界条件，就可以求得应力增量。对于各向同性的材料，有

式中:σij——任意一点的应力张量，用来表示该点的应力状态;

λ、μ——拉梅常数;

θ——体积应变;

δij——系数，当i=j时，δij为1，否则，δij为零。

这样，通过上述各式的迭代求解，就可以得出每一迭代时步对应的单元的应力应变值[7]。

2 滑坡概况及数值模型

2．1 滑坡概况

滑坡地处澜沧江右岸，坡体主滑方向为N40°W，岸坡为阴坡，主滑方向长度为550 m，沿河长度为300 m，表面积为16.5×104m2。坡体前缘高程为750 m，后缘高程为980 m，滑坡地貌标志明显，平面展布呈箕形，以左右冲沟为边界，总体坡角约为21°，自然边坡高度285 m，岸坡植被中等发育，滑动面平均埋深约为15 m，估算体积方量为75×104m3。滑坡体发育于窝拖寨逆断层下盘，忙糯断层上盘，下伏基岩为二叠纪中粗粒花岗闪长岩，其表层物质以碎块石土层为主，结构松散(图1)。

图1 滑坡剖面示意Fig．1 Geological profile of landslide

2.2 几何模型

FLAC3D软件采用显示差分法求解微分方程，不形成刚度矩阵，每一步计算仅需要很小的内存。在求解过程中通过叠加每一时步的小变形获得大变形求解，从而仅占用计算机很少的内存即可，并且在模拟材料的塑性破坏和塑性流动方面具有明显优势。而地震工况下滑坡的破坏失稳常常伴随着巨大的坡体变形，所以文中采用FLAC3D软件预测滑坡在地震工况下的稳定性。由于FLAC3D前处理功能较弱，建立复杂的三维模型非常困难，故选用ANSYS软件进行前期建模，划分计算网格，再导入FLAC3D程序中计算，以优化分析过程。

模型以主滑方向的反方向为y轴正方向，竖直向上为z轴正方向，指向河流上游方向为x轴正方向。竖直方向坐标采用实际高程坐标，建立直角坐标系，如图2所示。模型分为两个部分，共有11 518个单元，2 777个节点，其底面高程为650 m，沿平行河流方向长度为869 m，沿主滑方向长度为1 077 m，采用Solid45单元计算。计算模型边界约束形式为，侧边界只对水平方向进行约束，底边界在水平和竖直方向均进行约束，模型的上部边界取为自由面。FLAC3D程序中，节点速度是主要变量，所以选取模型的边界条件是通过约束模型边界的节点速度实现的[8]，即模型底部边界的水平、竖直方向的速度约束和四周边界水平方向的速度约束。在计算程序中，将边界上xvel或yvel的值设置为0，以达到约束边界的目的。利用弹性求解法生成初始地应力场的步骤是，将体积模量和剪切模量设置为大值并将材料的本构模型设置为弹性模型，然后进行求解。初始地应力场生成以后，赋予材料相关参数，设置重力，加载水平荷载，并采用摩尔－库伦模型进行计算。

图2 三维模型网格Fig．2 Three-dimensional block group model

根据实际地勘报告，计算采用的各层物理力学参数如表1所示，其中，基岩为微风化和新鲜岩体。

表1 岩、土体物理力学指标Table 1 Physico-mechanical parameters of rock and soil

3 结果分析

3．1 塑性区分布规律

天然工况条件下塑性区零星分布于滑坡体上游边界的冲沟附近，如图3a所示。塑性区分布高程(z坐标)范围为790～830 m，平行河流方向(x坐标)范围约为550 m。而地震工况条件下塑性区已连成大片，基本覆盖了整个滑坡表面，滑坡体上部以张拉屈服为主，中部及下部以剪切屈服为主，如图3b所示。通过对比得知，在地震工况条件下，滑坡体塑性区面积急剧增加。

从天然工况条件下x=560 m的塑性区剖面(图4a)分布区域来看，滑坡体的屈服区域很小，仅在滑坡体表面零星出现，塑性区分布于y坐标(主滑方向)270～480 m、z坐标790～830 m之间，厚度为5 ～10 m。这说明滑坡体发生剪切和张拉破坏的可能性很小，即使发生，也仅是局部区域，不会出现整体性滑动。而从地震工况条件下同一剖面的塑性区分布图(图4b)来看，塑性区已经上下贯通，距离滑坡体表面深度为10～20 m，可以明显看出滑动面位置，很可能发生整体性滑动。

图3 整体塑性区分布Fig．3 Plastic zones distribution on whole

图4 x=560 m剖面塑性区Fig．4 Plastic zones distribution on X=560 m plane

3.2 剪应力(剪应变增量)分布规律

判断堆积体的(潜在)滑动面(带)，可根据其剪应力(应变增量)来判断。剪应力较为集中或剪应变增量较大(绝对值)的部位，则为其(潜在)滑动面(带)，变形破坏也都沿此处发生;剪应力分布较为分散(均匀)或者剪应变增量较小、基本上未发生变化的部位，一般不会有潜在滑动面产生，因此，这些部位也不会发生较大的变形和破坏。天然工况条件下剪应力集中带主要出现在滑坡体中上部、上游边界冲沟附近，高程为880 m左右(图5a);而地震工况条件下剪应力集中带的表面积约为天然工况条件下表面积的两倍(图5b)。

图5 剪应变增量立体展示Fig．5 Shear strain increment on whole

天然工况条件下x=450 m的剪应变增量剖面见图6a，剪应力集中带仅出现在滑坡表面以下约5 m范围内，坡体内部没有明显的剪应力集中。地震工况条件下同一剖面的剪应变增量见图6b，可以看出，潜在滑动面的位置已经深入坡体内部，很可能会发生整体性滑动。

3.3 位移场分布规律

天然工况条件下的滑坡体变形主要出现在滑坡体中部和中上部区域，如图7、8所示。竖直方向位移总体表现为下沉，最大沉降区与水平方向最大位

图6 x=450 m剪应变增量剖面Fig．6 Shear strain increment on X=450 m plane

图7 滑坡体竖直方向位移云图Fig．7 Vertical displacement contour on whole

移区位置基本一致，局部沉降值最大为0.71 cm (图7a)。出现局部最大沉降的区域位于滑坡体中上部、上游边界冲沟附近，高程(z坐标)范围为890～930 m，平行河流方向(x坐标)范围为550～590 m，表面积小于100 m2;在滑坡体中部，高程为830～860 m出现隆起，最大值为0.29 cm。地震工况条件下滑坡体沉降区和隆起区的位置与天然工况条件下的基本一致，但位移量在数量级上有很大差别。局部沉降值最大为49.77 cm，在滑坡体中下部，最大隆起值为50.83 cm(图7b)。

水平方向位移主要表现为倾向河流方向。天然工况的水平方向局部最大位移值为2.10 cm，出现区域与竖直方向局部最大沉降区一致，滑坡体其他部位的水平方向位移值大都在0～1 cm之间(图8a)。地震工况下局部最大水平位移值为200.93 cm，位移值过大，判断为滑坡体失稳(图8b)。

图8 滑坡体水平方向位移云图Fig．8 Horizontal displacement contour on whole

根据强度折减法得出天然工况下滑坡体稳定系数为1.19，地震工况条件下滑坡体稳定系数为0.88。

4 结论

文中结合工程实例，选用拟静力分析法对一滑坡地震工况下的稳定性进行数值分析，主要得到以下结论:

(1)天然工况条件下，滑坡体塑性区分布零星，剪应力集中带仅位于滑坡体表面，位移场分布合理;加载地震荷载后，滑坡体塑性区上下贯通，滑动面明显，剪应力集中带深入滑坡体内部，形成完整的潜在滑动面而且位移大幅度增加。可以判断，地震工况条件下滑坡体会出现大面积失稳。

(2)天然工况和地震工况两种条件下，滑坡体的水平位移值均约为竖直位移值的三倍，可见滑坡体的破坏形式以水平滑动为主。

[1] 徐梦珍，王兆印，漆力健．汶川地震引发的次生灾害链[J]．山地学报，2012，30(4):502－512．

[2] 叶海林，黄润秋，郑颖芝，等．地震作用下边坡稳定性安全评价的研究[J]．地下空间与工程学报，2009，5(6):1248－1252．

[3] 李 亮，褚雪松，庞 峰，等．地震边坡稳定性分析的拟静力方法适用性探讨[J]．世界地震工程，2012，28(2):57－63．

[4] 包世华．结构动力学[M]．武汉:武汉理工大学出版社，2004．

[5] 宫经伟，严新军，夏新利，等．拟静力法在溢洪道闸室抗震稳定计算中的应用[J]．人民黄河，2010，32(1):107－108．

[6] 邹 栋，郑 宏．快速拉格朗日法及其在边坡稳定性分析中的应用[J]．矿业研究与开发，2005，25(5):84－87．

[7] 黄润秋，许 强．显式拉格朗日差分分析在岩石边坡工程中的应用[J]．岩石力学与工程学报，1995，14(4):346－354．

[8] 陈育民，徐鼎平．FLAC/FLAC3D基础与工程实例[M]．北京:中国水利水电出版社，2009．泥的钢纤维混凝土密实度得到提高，从而使得超声波声速、回弹值和混凝土抗压强度都得到提高。这与在普通混凝土中得到的结论是一致的。

(编辑 荀海鑫)

Stability prediction of landslide under seismic load using FLAC3Dsoftware

KAN Lu， ZHANG Faming
(School of Earth Sciences＆Engineering，Hohai University，Nanjing 210098，China)

This paper building on a landslide associated with a reservoir construction on Lancang River discusses the simulation of seismic load by selecting pseudo-static method，the development of numerical model using ANSYS software，the introduction of FLAC3Dsoftware for calculation，and the prediction of the landslide stability．The comparison shows that，under the natural condition，there occur the sporadic distribution of the plastic zones of the slope，the concentration of the shear stress only on the surface of the landslide，and the rational displacement of field;but under seismic load，there arise the complete connection of the plastic zones of the slope，the obvious slip plane，and the extension of the shear stress concentration area deep into the interior body of the slope，resulting in a complete potential sliding surface，marked by a significant increase in the displacement;the failure mode of the landslides under the two conditions is dominated by the horizontal slip while landslide due to the seismic condition is accompanied by large areas of instability．The study serves as a

for the anti-seismic designs．

earthquake;landslide;stability factor; pseudo-static method

10．3969/j．issn．1671－0118．2013．01．017

P642．22

1671－0118(2013)01－0077－06

A