基于有限元分析的缓倾滑坡模拟

朱 涛

(福建省地质工程勘察院，福州，350002)

缓倾滑坡[1-4]是众多滑坡中较为特殊的一类滑坡，其中的“缓倾”指的是滑动面的倾角较为平缓(为5 °～ 15°)，有时甚至近于水平。该类滑坡的滑动面土体内摩擦角往往大于滑面倾角，按照传统的极限平衡理论，没有特殊的原因难以发生滑坡，但实际上这类滑坡却屡见不鲜。我国的缓倾滑坡主要分布在三峡库区中段万州、重庆一带，四川盆地及周边地区，以及黄土高原等地。其中，事发于四川省巴中市南江县的一系列滑坡则是我国缓倾滑坡事故中最为典型的案例[5]。该县于2011年9月16日期间发生了千余处滑坡，这些滑坡沿着呈现“光面”形式的界面而滑动。该类边坡上覆土层较薄，岩层倾斜角度较缓，使得重力导致的下滑力较小，因而根据传统的极限平衡理论，应该具有较高的稳定性，一般不易发生失稳破坏。然而，数以千计的缓倾滑坡在降雨期间不约而同地陆续发生，有力地证明了其发生是由于十分特殊的成因机理和形成条件。虽然少量的缓倾滑坡发生并不会造成过大的危害，但这种大规模群集性的土质滑坡往往都发生在人员居住的密集区域，存在的大片良田和村落聚集地，因而必须引起足够的重视。

虽然目前国内外学者已经对常规滑坡进行了较为深入而全面的研究[6]，但是，对于近年来频繁发生，分布在广大农村的中小型缓倾滑坡的研究仍然相对较少。由于这类中小型缓倾滑坡给人不具备发生滑坡条件的错觉，因此在前期不会引起足够的重视，甚至还出现人为加剧滑坡形成与发展的行为。一旦滑坡形成，必然造成严重的经济损失或人员伤亡。而这些滑坡的治理经费会给经济原本就欠发达的农村地区和山区带来更大的压力。因此，笔者充分利用数值计算手段的优越性，选取三明明溪烽林场施工工地缓倾滑坡作为工程实例进行数值模拟分析，以期能够建立起该类滑坡的地质概念模型、变形破坏力学模式和数值计算模型，以便为将来该类滑坡的治理提供依据。

1 工程实例分析

所研究的缓倾滑坡位于三明明溪烽林场在建工地的开挖边坡，是福建省内较为罕见的土质缓倾滑坡。该滑坡位于侵入岩与沉积碎屑岩的接触带上，特殊的岩土体组合和独特的水文地质条件， 形成了该滑坡独特的变形破坏特征(图1)。

图1 明溪烽林场工地滑动面及滑床特征图Fig.1 Sliding surface and slide bed characteristic diagram in Mingxifeng forest site1—素填土；2—粉土；3—有机质土；4—泥岩残坡积土；5—花岗岩残坡积土；6—全风化花岗岩；7—强风化泥岩；8—强风化花岗岩；9—中风化花岗岩；10—滑动面；；12—钻孔

该滑坡主轴长230 m，前缘宽260 m，前后缘高差35 m，滑坡面积接近6×104m2，滑坡体的平均厚度为13 m，潜在体积76×104m3，主滑方向近正西，沿主滑方向坡面总体坡度8° ～ 10°，主滑段滑面倾角约6°。滑坡前缘为整平的建设用地，地面出现隆起现象(照片1)，在建厂房墙体开裂。滑坡前缘位于盆地中、前部，为地下水排泄区，出露地层为一套厚度数米至十数米的有机质土(照片2)，呈软塑或流塑状，具有流变性质。场地东侧为人工填土边坡(滑坡后缘)，宽度约210 m，高度7～9 m，坡度65°。

照片1 滑坡前缘的隆起和地面开裂Photo.1 Uplift of landslide front and ground cracking

照片2 滑坡前缘有机质土出露Photo.2 Organic soil exposure in landslide front

场地周边属剥蚀丘陵地貌， 总体地势为东高西低， 山顶标高一般 350～450 m， 个别山头超过500 m，山顶浑圆，山脊平缓圆润，相对高差60～90 m，最大相对高差 280 m，自然坡度一般 15°～25°。主要沟谷呈北东-南西向，为“U”型谷，沟底平缓开阔，形似一个口朝西的唢呐状。滑坡区微地貌为平缓的山间小盆地，南北宽约 370 m，东西长约450 m，三面环山，开口朝西，左右两侧为小溪沟，右侧溪沟为常年流水溪沟，溪沟中上部冲沟发育。

初步的勘察分析结果表明，该滑坡主要受到强降雨影响，雨水大量渗入岩土体，增加了岩土体的容重，同时降低了土体的抗剪强度，增大了下滑力，同时，由于雨水的持续软化作用，边坡中形成了一条坡度十分缓的软弱滑带，从而导致了滑坡的发生。

2 滑坡工程实例的有限元模拟

2.1 有限元模型

采用ANSYS/LS-DYNA[7]软件对滑坡进行模拟，经过一定的简化，建立了单个单元厚度的滑坡有限元模型(图2)。边坡模型长度为370 m，左侧高度为110 m，右侧高度为50 m，单元划分从上至下(颜色由浅到深)的土层分别是粉土、有机质土、花岗岩残积土、全风化花岗岩和强风化花岗岩。其中，黑色单元为单层厚度的潜在滑动面单元，该滑动面是根据现场勘察、室内试验成果，并结合反演分析和工程类比的方法所获得的。

土体采用Mohr-Coulomb模型来描述，主要的本构参数为土体密度、黏聚力、内摩擦角和弹性模量，计算参数值采用试验得到的土体材料参数(表1)。

图2 滑坡的有限元模型Fig.2 Finite element model of landslide

2.2 计算工况

根据文献得知[8]，缓倾滑坡的主要诱发因素为地下水产生的静水压力和浮托力，以及由于地下水和降雨作用导致的土体强度参数值的降低(即生成了软弱滑动面)。再考虑到不同工况下的地下水位的特点，各工况之间的差别在于，地下水位越高，对滑坡后缘产生的静水压力以及滑坡体底部产生类似三角形分布的浮托力也就越大，因而对滑坡稳定性就越不利。同时，水位线以上的土体采用天然土体材料参数，而水位线以下的土体则需要采用饱和土体材料参数。

有鉴于此，通过软弱滑动面的生成、滑坡体后缘的静水压力和滑坡体底部的浮托力的施加来考虑降雨量大小，以及地下水对滑坡的影响，笔者采用以下3种不同地下水位工况进行综合分析。

表1 试验得到的土体材料参数

天然工况：地下水位取勘察钻孔的实际水位，水位线位于滑坡后缘裂隙以上10 m处，也即h1=10 m，因而最大静水压力和浮托力为γwh1=100 kPa。地下水位以上岩土体采用天然状态参数，地下水位以下采用饱和状态参数(图3-a)。

雨季工况：地下水位取高水位(根据雨季期间监测孔揭露的地下水水位)，水位线位于滑坡后缘裂隙以上20 m处，也即h2=20 m，最大静水压力和浮托力为γwh2=200 kPa。地下水位以上岩土体采用天然状态参数，地下水位以下采用饱和状态参数(图3-b)。

极端暴雨工况：考虑滑坡体处于全饱水状态，水位线位于坡体表面，也即h3=30 m，最大静水压力和浮托力为γwh3=300 kPa。滑坡岩土体全部采用饱和状态参数(图3-c)。

图3 不同地下水位计算工况的力学概化模型Fig.3 The generalized mechanical models for calculating under different groundwater level conditions

2.3 结果分析

2.3.1 安全系数计算结果

赵尚毅等[9]认为，对于一般的均质土坡平面应变问题，在ANSYS程序中将迭代次数设定为1 000次，将力和位移的收敛系数设定为0.000 01完全可以保证足够的计算精度。因此，在 ANSYS 程序的D-P强度准则中采用非关联流动法则(膨胀角为0)进行计算，力和位移的收敛标准系数均取为0.000 01，最大迭代次数为1 000次，一次性施加重力荷载，也即荷载增量步设置为1步，有限元求解器选用ANSYS程序提供的稀疏矩阵求解器选用全牛顿-拉普森迭代方法。

各工况的安全系数计算结果可以发现，天然工况的安全系数为0.95，雨季工况的安全系数为0.92，极端暴雨工况的安全系数为0.86，边坡基本都处于不稳定或者欠稳定状态，尤其是在雨季和极端暴雨情况下，该边坡稳定性较差。因此，亟需采取合适的治理措施，对边坡进行加固处理，以免发生滑塌。

此外，边坡被破坏时滑动面上节点位移和应变将产生突变，滑动面位于水平位移和应变突变的地方，因此可在ANSYS程序中通过绘制边坡等效应变云图来确定滑动面。有限元强度折减法滑动面展示(图4)，位移云图通过大小为0.1(即10%)的应变尺度来显示。显而易见，由于事先生成了强度参数较低的一层圆弧形的软弱滑动面单元，因而计算得到的最危险滑动面都发生于这层单元之中，整层单元都呈现大应变状态。

图4 有限元强度折减法滑动面展示Fig.4 Sliding surface display of strength reduction finite element method

2.3.2 滑坡变形

3种不同工况的滑坡水平向位移云图对比，其中黑色边框表示未变形的边坡轮廓，用于与变形后的边坡进行对照。为了更清楚地显示对比结果，将边坡变形都放大了3倍。

(1)由于软弱滑动面的存在，天然工况下的边坡也出现了较大的水平位移(图5-a)，滑坡体后缘区域出现了一定的位移集中现象(滑坡体后缘的位移云图整体呈深色)，而坡脚则出现一定的隆起。在滑坡现场也观测到了类似的滑坡前缘地面隆起、开裂的现象，因而，数值计算结果的准确性在一定程度上得以证明。

(2)由于地下水位的抬升，雨季工况下导致滑坡后缘的静水压力和滑坡底部的浮托力有所增大，使得边坡的整体变形也有所增加，前缘隆起现象也出现了加剧趋势，滑坡体的位移集中现象也明显扩大(图5-b)。

(3)在极端暴雨情况，由于边坡处于最不利的受力状态和土质条件下，因而边坡的位移也急剧增大，滑坡前缘的隆起明显远远超出原始地面(图5-c)，因此，亟需采取合适的治理措施，对边坡进行加固处理，以免发生进一步滑塌。

图5 3种不同工况下边坡的水平向位移云图(单位：m)(变形放大3倍)Fig.5 Horizontal displacement nephogram of the slopes under 3 different working conditions

通过3种不同工况对比也可以发现，对于这类坡度较缓、滑面倾角较小的边坡，一般情况下是难以发生明显滑坡的，即使坡体中生成了软弱滑带，但是由于滑动面的土体内摩擦角大于滑面倾角，滑坡变形相对也不会太剧烈。然而，对于暴雨工况，由于滑坡体底部浮托力的作用，加上后缘裂隙充盈水分使得滑坡推力大大增加，从而使得边坡出现剧烈的变形，前缘隆起急剧增大，导致边坡在倾角很小的情况下发生滑塌。

3 结语

(1)通过数值模拟计算，形象地还原了缓倾滑坡的变形与破坏，并且得到了不同工况和不同地下水条件下的边坡位移，能够为将来该类滑坡的研究和预判提供一条可行之路。

(2)所研究的边坡坡度较缓，因此，在良好土质条件下，边坡稳定性很强，安全系数达1.54，边坡变形也较小，说明这类坡度很缓的边坡，一般情况下是不容易发生滑坡的，缓倾滑坡的形成是受到特殊的地下水因素影响。

(3)通过调整作用在滑坡后缘的静水压力和滑坡底部的浮托力来表征降雨和地下水对缓倾滑坡的影响，结果表明，该缓倾滑坡在天然工况、雨季工况和极端暴雨工况下，安全系数分别为0.95，0.92和0.86。因此，从滑动面倾角上分析，该边坡一般应处于稳定状态，但实际上由于地下水所导致的力学作用，使得该边坡处于不稳定或者欠稳定状态，尤其是在雨季和极端暴雨情况下，稳定性较差，亟需采取合适的治理措施，对边坡进行加固处理，以免发生滑塌的危险。

(4)在天然工况、雨季工况和极端暴雨工况下，滑坡的位移逐渐增大，前缘的隆起逐渐变得更剧烈，滑坡体的整体位移集中现象的越发明显。

1 冯文凯，石豫川，柴贺军，等. 缓倾角层状高边坡变形破坏机制物理模拟研究. 中国公路学报, 2004, 17(2).

2 张群，许强，李江，等. 南江“9·16”群发性缓倾浅层土质滑坡特征与成因机制研究. 自然灾害学报, 2015, 24(3).

3 范宣梅，许强，黄润秋，等. 四川宣汉天台特大滑坡的成因机理及排水工程措施研究. 成都理工大学学报(自然科学版), 2006, 33(5).

4 黄国平，柳侃，叶龙珍. 明溪烽林场滑坡形成机理及防治措施. 江西理工大学学报, 2014, 35(5).

5 朱亮. 巴中黑山坡平推式滑坡成因机制及变形破坏演化过程分析.北京：中国铁道科学研究院, 2014.

6 许强，范宣梅，李园，等. 板梁状滑坡形成条件、成因机制与防治措施. 岩石力学与工程学报, 2010, 29(2).

7 Hallquist, J.O. . LS-DYNA Keyword User's Manual (971 R6.0.0). Livermore:Livermore Software Technology Corporation, 2012.

8 杜潇翔. 贵州省习水县程寨滑坡形成机制及稳定性研究.成都：成都理工大学, 2013.

9 赵尚毅，郑颖人，张玉芳. 极限分析有限元法讲座——Ⅱ有限元强度折减法中边坡失稳的判据探讨. 岩土力学, 2005, 26(2).