基于FLAC3D软件的降雨型滑坡稳定性分析

张磊 刘岁海 顾彩玉 冉谨荥

［关键词］ 降雨型滑坡；安全系数；FLAC3D；数值模拟；雅安市

［摘 要］ 以雅安市名山区前进乡林泉村一斜坡为例，通过现场调查、室内外试验以及数值模拟，在查明滑坡背景基础上，对其发育特征进行分析；基于FLAC3D软件，采用强度折减法计算处于天然工况和饱水状态下的滑坡安全系数，通过对比分析两种工况下的塑性图和剪应变增量云图来揭示该滑坡变形机制以及预测滑坡后续变化情况，结果表明：天然工况下安全系数为1.131，处于基本稳定状态，饱水状态下滑坡安全系数为1.042，处于不稳定状态；若受到暴雨等极端天气或在其他外界因素的影响，极有可能发生二次滑坡。

［中图分类号］ TU43  ［文献标识码］ A  ［文章编号］ 1000－0941（2023）05－0053－05

近年来，受人类活动和极端天气影响，我国地质灾害发生频繁，其中滑坡灾害的危害最为严重。我国西南部滑坡分布广泛，尽管许多大型的滑坡堆积物可以长时间保持稳定或基本稳定状态[1]，但众所周知，在地震、暴雨等极端条件的影响下，滑坡很可能被重新激活。为保障人民群众的生命及财产安全，应积极开展以滑坡稳定性分析为基础的滑坡隐患防治工作。传统的极限平衡法主要从物理角度出发，分析滑坡内部应力的变化以及应力之间的相互关系，进行理论性分析和力学计算［2-3］。但滑坡灾害是一个极其复杂多变的非线性过程，不能靠单一的线性描述进行准确分析。目前，数值模拟已成为滑坡稳定性分析的一种重要方法，而运用FLAC3D软件强度折减法功能进行数值模拟的优势就是在探究土质边坡的稳定性同时，也能有效地表现出边坡潜在滑移面的动态变化及所在位置［4］。

本研究以2019年7月23日雅安市名山区前进乡林泉村特大暴雨造成的滑坡为例（滑坡位置見图1）。目前该滑坡整体处于基本稳定状态，但在暴雨或连续降雨等不利因素影响下，滑坡将有可能发生失稳破坏，并威胁滑坡前缘4户13人的生命财产安全和公路安全，估算直接经济损失360万元。笔者在充分考虑各种因素对滑坡的影响下，选取FLAC3D软件作为数值模拟平台，采取强度折减法对林泉村滑坡进行稳定性分析，以期更有效地揭示滑坡变形机制。

1 研究区概况

1.1 研究区地质环境

研究区位于雅安市名山区，地理坐标为东经103°12′04.50″、北纬30°00′47.20″，处于成都平原与盆周山地过渡带，属亚热带季风性湿润气候区，日照少、无霜期长、四季分明，夏无酷日、冬无严寒，降水量充沛而集中，境内降水集中在夏季，尤以7月降水量最大，多大到暴雨、连续降雨。年平均气温15.4 ℃，年平均相对湿度82%。地势呈西北高、东南低的特点，由北西向南东倾斜，地形地貌以构造侵蚀低山为主，包括台状丘陵和浅丘平坝。主要河流有4条，大小支流37条，河流总体上呈现流量较小、流程较短的特点。地质构造属于新华夏系西部尾段，处于龙门山断裂带与四川断陷盆地的结合部，区内跨3个二级构造带，即天台山雁行带、成都平原凹陷带和熊坡背斜雁行带，形成一系列北北东向的褶皱和断层。滑坡区表层覆盖第四系全更新统滑坡堆积层（Q4del），主要是粉质黏土夹碎石，含量为5%～10%，厚0.5～15.0 m。滑坡范围之外的斜坡地带和滑体土之下为第四系全新统残坡积层（Q4el+dl），厚0.5～15.0 m，局部有碎石，上覆白垩系上统灌口组（K2g），其裂缝发育一般，内部填充粉质黏土，岩层产状为275°∠6°。

1.2 滑坡基本特征

滑坡位于雅安市名山区林泉村五组。滑坡后缘高610 m，前缘高553 m，前后缘高差57 m，剖面形态呈阶梯形，后部地形较陡，中部平缓，下部又较陡，平面呈圈椅状（见图2）。滑坡纵向上无转折，整体坡度20°～25°，坡向120°。滑坡主轴长约135 m，宽约77 m，面积约1.04万m2，厚约14 m，滑坡方量约14.56万m3，区内植被覆盖率达70%，主要为林木和农作物等。其中坡体上主要为林木，岩土层主要为粉质黏土夹碎石，为中型土质滑坡。

1.3 滑坡变形破坏特征

通过调查得知，该滑坡失稳前已发生局部变形，并形成连续贯通的地面裂缝。最早于2019年7月20日下午出现迹象，在这期间，由于居民的房屋老化，加之屋后坡体大多数为泥质沙土，因此在连续降雨的影响下，坡体早已出现张拉裂缝并持续扩大，滑坡两侧则有少许变形。随后几天降雨持续，斜坡上部土壤吸收过多水分，质量增加，向下部的压力也在缓慢增加，但由于植被较多，因此暂未形成滑动；而中部与下部多为粉质黏土，土体较为松散，随着雨水渗入，土体有效应力降低，容重增加，发生滑动变形。这一变形的时间差，造成临空面的形成，在中部和前缘的张拉力下，滑坡后部于23日晚上发生变形破坏。通过调研发现，林泉村滑坡为中距碎石土质滑坡，应力释放较为分散，滑坡发生后，后缘滑床出露，滑土体堆积于中部平台，裂缝较为分散，分布于中后部1处、中部4处及中下部2处，前缘明显鼓胀，且由于滑坡发生后应力释放，因此在侧缘形成一处土体溜滑和一处地面下沉，屋后墙受损且裂缝贯通发育。通过无人机拍摄和现场实地勘察，得到该滑坡各部位的典型特征，见图3。

2 滑坡失稳影响因素分析

2.1 构造活动

滑坡区位于龙门山强烈抬升区前缘的凹陷区，自中三叠世末期，印支运动促使龙门山构造带及宝兴背斜古海槽逐渐隆起，形成前缘凹陷之后，经历了晚三叠世印支运动强烈抬升和下陷， 张磊等：基于FLAC3D软件的降雨型滑坡稳定性分析燕山运动的反复抬升、下陷，后受喜山运动抬升、下陷及后期挤压，褶皱成山。由于区内凹陷与附近多隆起，因此山体前缘和下部质量增加，推力较大，且坡度较陡，形成许多细小的张拉裂缝，剪切破坏明显，裂缝内的松散物质在长期的风化侵蚀和流水侵蚀共同作用下，特别容易造成滑动危害。

2.2 降雨因素

雅安市气象局资料显示，名山区2019年5—7月的累计降雨量达到986.5 mm，对比数据可知，这几个月的降雨量相当于该区全年降雨量的62%，特别是7月的降雨量更是达到了640.2 mm，是往年的2～3倍。据资料，名山区在7月11日发布了暴雨蓝色预警，当日预计3 h内降雨量将达到50 mm以上。降雨在滑坡破坏前持续了17 d，总降雨量为430.1 mm，7月13—21日，滑坡区持续降雨，特别是在滑坡发生前的7月21日20时至22日8时，发生本月最大降雨，为172.5 mm。对于这次滑坡来说，长时间的降雨一方面使得坡体前缘下切侵蚀速度加快，另一方面使得地下水位上升，减小了滑体的法向应力，在张拉应力及中上部的推力作用下，斜坡进一步发生失稳。

3 基于FLAC3D的滑坡变形分析

图4为滑坡剖面示意。采用ANSYS软件建立滑坡模型，并导入FLAC3D，运用其中的强度折减法进行计算，在得到滑坡的安全系数后，通过滑坡的最大剪切应变增量以及不同工况下的塑性区分布图来揭示滑坡变形机制。

该滑坡模型沿x轴方向长为210 m，沿y轴方向宽为1 m，沿z轴方向最高点为85 m。为了计算精确，堆积体、滑动体的网格尺寸为0.5 m，基岩部分的网格尺寸为1.5 m，模型总共划分23 274个节点、68 915个单元，各边界点均耦合。边坡模型顶部边界为自由边界，模型四周约束法向位移，底部刚性约束，所得到的模型见图5。采用FLAC3D内置的mohr-coulomb本构模型进行模拟，所得参数见表1。

在天然工况下，滑坡的安全系数为1.131，处于基本稳定状态，而最大剪应力位于边坡上缘。图6为天然工况下最大剪切应变云图，由图6可知，滑动面的大致位置在边坡上缘及中部，最大剪切应变增量为0.351，初步判断内部存在滑动趋势，但剪切带内部不连续，仅前缘及中部台阶处有着较明显的剪切变形，表明该处滑坡抗滑力是由滑带前缘提供的，可以判断出未形成连续的滑动面，整体仍处于稳定状态。

图7为天然工况下滑坡塑性区分布，由图7可知，滑坡的前缘及中部大部分未在张拉屈服状态，曾进入过张拉屈服状态，但最终未发生塑性破坏；仅有坡体上缘表面出现张拉塑性區域，内部整体稳定，未出现明显塑性贯通区。这与滑坡发生前，其中部出现的拉张裂缝相对应。由此可知，在天然状态下，该滑坡的中间部分堆积体虽具备变形能力，但滑坡前缘堆积体稳定不变形的状态，为其提供了充足的抗滑力，保障坡体保持稳定。但是，若遭暴雨或持续降雨，则滑带含水率将上升，滑坡堆积体质量增加，塑性区沿滑带将进一步扩展拉伸至坡体前缘，坡体表面受拉伸作用影响，破坏范围也将逐渐变大。

图8为饱水状态下最大剪切应变云图，由图8可知，饱水状态下滑动面剪切应变增量出现较为明显的连续区域，最大剪切应变增量0.055 6。在饱水状态下边坡上缘的表层以下大约15 m的位置为剪切应变最大值处，且有不断向下延伸的趋势，结合滑坡在饱水状态下的安全系数1.042，判断滑坡处于欠稳定状态。

图9为饱水状态下滑坡塑性区分布。由图9可知，在暴雨条件下，滑坡体的屈服状态与天然状态下有明显差别，滑坡体处于饱水状态时坡体下滑的规模比天然状态下大得多，且塑性破坏出现明显贯通区域直至边坡下缘，形成滑动带，通过塑性区云图初步判断滑坡即将失稳。由剪应力自滑坡上缘向中部乃至前缘递减可以看出滑坡此时虽并未形成滑动面，但降雨渗入作用使得滑带逐渐变得饱和，当中间堆积体提供的抗滑力不足时，便会在重力作用下造成滑坡体滑动。值得注意的是，沿着强风化带界面已出现明显剪切塑性破坏区域，若不及时采取边坡加固措施，则可能出现新的滑动下推力，造成边坡大面积滑塌。

根据一系列数值模拟结果，得到不同工况下的水平位移云图和最大主应力云图，见图10、图11。结合林泉村滑坡后的地形地貌特点以及周围环境特征分析得出，林泉村滑坡由于强降雨或连续降雨导致斜坡坡体后缘下错、中部拉裂、前缘鼓胀突出、房屋变形这一过程，受地质环境影响，其失稳过程可以分为4个阶段：

（1）前缘剪出。受长期降雨或强降雨影响，坡体前缘土体达到饱和，并且研究区属于凹陷区域，受长期下切侵蚀作用，加上中上部的岩土体吸收雨水产生更大的向滑坡区下部的作用力，使得坡体前缘率先剪出，但是由于整个滑坡区有一定的高差，因此未能整体剪出。

（2）牵引滑动。在前缘剪出后，发生前一级滑动，而中部的堆积物在雨水渗透导致质量变大的情况下，通过前一级滑动所提供的临空变形空间，在前缘牵引力以及底部扬压力的作用下发生牵引滑动，使第二、第三阶梯形成比之前更大的临空变形面。

（3）高速滑动。随着前缘剪出，中部也在张拉力作用下发生撕裂，造成中部临空。而长时间的降雨渗透导致原斜坡的堆积土逐渐变得饱和，除增加自身容重外，还浸湿软化了原来的土层，使得其后缘失稳，在其后部静水压力增大和底部扬压力、前缘牵引力共同作用下沿基覆界面高速滑动，并发生抛洒现象。

（4）复合滑动。在前面3个滑动后，沿着强风化带界面出现明显剪切塑性破坏区域，加上雨水渗透，增大了岩土体容重，坡体质量急速增长，大大降低了土体的抗剪强度，张拉裂缝不断扩大，在前缘牵引力和后部推移的作用下，滑体跟随堆积体一起滑动，但是又随着坡度逐渐趋于平缓，下滑力减小，速度降低。

4 结 语

对林泉村滑坡进行了野外实地勘察，并开展了相关力学试验，对比分析了天然工况和饱水状态下滑坡的抗剪强度变化情况，揭示了滑坡变形机制。对林泉村滑坡进行了滑坡稳定性评价，得到以下结论：

（1）林泉村滑坡平面形态大致呈椭圆形，主滑方向为120°～130°，滑坡面积约为1.04万m2，厚约14 m，滑坡方量约14.56万m3，为降雨型滑坡。

（2）由基于FLAC3D的数值模拟分析结果可知，天然工况下滑坡安全系数为1.131，处于基本稳定状态；饱和状态工况下的滑坡安全系数为1.042，处于失稳状态。滑坡前缘滑动具有牵引式特征，中部沿强风化带界面有明显剪切破坏，后缘具有推移临空特征。

（3）受地形特征影响，堆积层、滑动带为滑坡变形提供了基础，若该地发生持续降雨或者暴雨，则该滑坡将会继续发生失稳破坏。

[参考文献]

[1] 邱煜珩，舒中潘，张军，等.基于FLAC 3D的推移式滑坡变形破坏模式及稳定性分析[J].四川建筑，2020，40（4）：171-174.

[2] 罗忠行，雷宏权.基于FLAC3D的米贝复式滑坡稳定性分析[J].中国地质灾害与防治学报，2020，31（4）：52-62.

[3] 刘光华，熊超，赵鹏.滑坡抗剪强度参数反演数值模拟研究[J].重庆交通大学学报（自然科学版），2013，32（5）：969-973.

[4] 李勤光.基于强度折减理论的西南某库岸滑坡稳定性FLAC3D分析[J].铁道勘察，2012，38（1）：52-55.

［作者简介］ 张磊（1996—），男，贵州遵义人，硕士研究生，研究方向为工程地质勘察与评价；通信作者刘岁海（1971—），男，陕西咸阳人，副教授，硕士，研究方向为工程地质勘察与评价。

［收稿日期］ 2022-07-15

（责任编辑 杨傲秋）