1、DOI：10.19645/j.issn2095-0144.2022.12.005收稿日期：2022-10-10作者简介：张检召（1997-），男，重庆人，硕士研究生，主要从事地质灾害评价与预测，E-mail：。基于SPH的溃屈式滑坡失稳过程分析张检召1，瞿华南2（1.成都理工大学 地质灾害防治与地质环境保护国家重点实验室，四川 成都 610059；2.四川省地质矿产勘查开发局 九九水文地质工程地质队，四川 绵阳 621000)摘要：光滑粒子流体动力学（Smoothed Particle Hydrodynamics，SPH）方法已经广泛应用于滑坡地质灾害的研究，然而对大多数滑坡失稳过程的研究仍然
2、是在二维模型或相同岩性基础上开展的。为了研究三维模型下软硬互层岩质顺层滑坡失稳的运动状态，基于光滑粒子流体动力学方法，采用修正Mohr-Coulomb准则并耦合强度折减法对滑坡失稳运动开展数值研究，模拟滑坡失稳后的运动过程。结果表明：结合强度折减法，光滑粒子流体动力学方法可以较准确地模拟滑坡的渐进破坏，并获取潜在滑移面；采用该模型模拟南国大道顺层岩质滑坡失稳变形后，得到的结果与现场结果十分吻合，结果精度相比二维模型提升约20%，同时模型还可以较为清晰地反演软硬互层滑坡轻微弯曲变形、弯曲隆起、滑动面贯通的溃屈式破坏过程；随着折减系数的增加，滑坡危险范围逐渐增加，但增加速率逐渐减缓。关键词：软硬互
3、层滑坡；光滑粒子流体动力学；稳定性分析；强度折减法；滑坡大变形中图分类号：P642.22文献标志码：A文章编号：2095-0144（2022）12-0019-06第 58 卷 第 12 期2022 年 12 月GANSU WATER RESOURCES AND HYDROPOWER TECHNOLOGY甘 肃 水 利 水 电 技 术Vol.58，No.12Dec.，20221前言我国的地质灾害爆发十分频繁，其中又以滑坡灾害最为严重，分布最广，数量最多1。据统计，历史上滑坡灾害至少对400多个乡镇和10 000多个村庄造成严重破坏，滑坡灾害点数量高达41万，总面积1.74106km2，占国土总面
4、积的18.10%。自1995年以来，连续多年有不少于1 000人因滑坡灾害而丧命2。因此，对滑坡稳定性进行研究，并分析其滑动过程及危险范围具有重要意义。目前，针对滑坡稳定性的研究方法主要有极限平衡法和有限元法3-4。极限平衡法中对整体滑坡的假设限制了该方法的应用5；有限元法相对成熟，但由于其网格的特性，当滑坡失稳发生变形后会产生较大的畸变，从而导致计算结果不准确甚至计算中止。特别是对一些高速远程滑坡的模拟预测，有限元法将不再适用。鉴于网格模型存在的这种局限性，非连续介质的无网格模型开始广泛应用于滑坡等地质灾害的模拟、预测及评价。其中光滑粒子流体动力学在大变形问题和追踪流体自由液面等方面相比传统
5、数值模型取得了更好的结果。国内外专家学者在一定程度上开展了滑坡失稳运动方面的数值研究，但大多是在二维模型的基础上展开分析，或是只考虑同一种岩性下的滑坡失稳过程，综合考虑软硬互层的三维滑坡模型仍需深入研究。为此，作者以光滑粒子流体动力学（SPH）方法为基础，开展了三维模型下软硬互层岩质顺层滑坡失稳的运动状态研究，探讨了滑坡失稳的破坏过程、位移变化和危险范围。2基于SPH方法的滑坡稳定性计算模型2.1滑坡控制方程光滑粒子流体动力学（SPH）是拉格朗日无网格方法，该方法将连续的流体或固体用一系列离散的点来表示，每个点承载了物质的质量、密度和速度等特性，通过求解每个点的运动，最终得到整个系统的力学特征
6、。其定义如下6：f(x)=f()x W()x-x,h dx（1）式中：W(x-x，h)为光滑核函数。对于滑坡的模拟，采用连续性方程和动量方程作为控制方程，它们的拉格朗日描述如下：didt=ij=1NmjjvijWijxj（2）19dvdt=j=1Nmji+ji+jWijxi+f（3）式中：v为颗粒速度；为应力张量；f为外力。岩土体采用的是修正Mohr-Coulomb（摩尔-库仑）准则。标准的Mohr-Coulomb屈服表面F可以表示为：F=-Psin+K()J2-ccos（4）式中：P为压力；为内摩擦角；K（）为应力罗德角函数；J2为应力偏张量的第二不变量；c为黏聚力。标准的Mohr-Coul
7、omb准则中，当剪切应力为0时，土体所受到的压力为负，这显然与实际情况不符。因此，ABBO等7对Mohr-Coulomb准则提出修正，修正后的准则与坐标轴平滑相交。修正后的Mohr-Coulomb屈服面表达式为：F=-Psin+J2K()2+a2sin2-ccos（5）式中：a表示修正后的屈服面与标准Mohr-Coulomb屈服面之间的拟合度。当a=0时，式（5）转换为标准的Mohr-Coulomb准则。2.2滑坡强度折减法为分析滑坡体的稳定性，引入剪切强度折减法（SRM）8中的折减系数（Fs），通过不断降低土体抗剪强度来最终自动形成滑面，此时滑坡破坏的临界折减系数被认为是该滑坡体的安全系数，
8、公式如下：c=cFs（6）=arctan(tan Fs)（7）式中：c和分别为滑坡体的黏聚力和内摩擦角。3案例分析案例中的岩质滑坡处位于达州市达川区南国大道二段，2020年8月中旬，该地区经历一次强降雨，在暴雨过后斜坡发生滑动。根据现场调查，滑坡体主要岩性为侏罗纪中统沙溪庙组的紫红色泥岩、泥质粉砂岩互层（图1），整体较完整，局部有大块石崩落，滑坡体左侧有2条纵向裂缝，间距6 m，延伸长度40 m，岩层产状30242，斜坡坡向287，坡度4045，岩层产状与斜坡近乎一致，滑坡为顺层岩质滑坡。斜坡纵向长150 m，横向宽约200 m，坡高约80 m。滑坡后缘有明显的滑坡后缘壁，滑坡后缘下错37 m
9、。滑坡处由于无明显临空面，导致滑坡前缘隆起变形，挤压南国大道二段道路往西偏移，阻断了200 m长的道路，严重威胁到南国大道的交通运行及斜坡对面在建小区的安全。根据现场调查，研究区岩性为泥岩和砂岩互层，且前缘无临空面，砂岩为含水层，泥岩为不透水层，在强降雨作用下岩体内的水不易排出，导致岩体软硬结构面易饱和，当层面倾角（42）明显大于结构面内摩擦角峰值时，该滑坡体发生溃屈式破坏。具体的演变过程可以分为三个阶段。轻微弯曲变形阶段：在连续强降雨条件下，雨水沿着卸荷裂隙和节理进入岩体，当滑面处下滑力大于其抗剪强度时，滑坡体开始缓慢下滑，而滑坡体前端受阻，导致滑坡体坡脚轻微隆起变形，局部岩体被压碎；弯曲隆
10、起阶段：随着变形加剧，滑坡体前缘持续隆起，岩体后缘出现拉张裂缝，局部出现崩落或滑落现象；滑动面贯通失稳阶段：滑坡体滑动面贯通，滑坡体整体失稳下挫。基于无人机拍摄所得影像建立研究区数字高程模型（DEM），再将数字高程模型导入ArcGIS 软件，并提取研究区地形点云文件，然后将地形导入Rhino三维建模软件建立研究区三维地形，并根据实测剖面（图2）建立软硬互层滑坡模型，之后将建好的模型导入HyperMesh 软件划分网格并转换成SPH粒子，最终得到了研究区域的地形计算模型（图3）。图3中深色部分为砂岩，浅色部分为泥岩。粒子距离设置为1.5 m，最终生成176 300个滑坡体粒子。滑坡参数通过现场取
11、样并根据 GB/T50266-2013 工程岩体试验方法 开展室内试验获得，研究采用的主要参数见表1。图1研究区滑坡全貌后缘壁J2s前缘隆起图例南国大道堆积区后缘壁滑坡边界2022年第12期甘肃水利水电技术第58卷20图4显示了南国大道二段滑坡坡顶检测点S1的位移变化，二维（2D）与三维（3D）模拟结果均显示坡顶位移具有先增大后趋于稳定的特点，但通过与现场实测数据的对比发现，3D 模拟结果误差为8.8%，2D模拟结果误差为38.8%。图5为南国大道二段滑坡最终坍塌稳定后的剖面图。图2研究区典型实测剖面34534033533032532031531030530029529028528027527
12、0265260255250245高程/m107南国大道二段道路范围前缘隆起范围后缘壁高7 mQel+dl4J2s30242S9S7S1S4S8S2S5S3S60102030405060708090100110120130140150160170180190水平距离/m345340335330325320315310305300295290285280275270265260255250245高程/m图例沥青路面Qel+dl4第四系残坡积层J2s侏罗系中统沙溪庙组砂岩泥岩粉质黏土检测点滑面图3研究区模型表1岩体参数取值岩体名称泥岩砂岩软弱结构面密度/（g/cm3）2.562.552.56杨氏模量
13、/GPa20.050.01.5泊松比0.300.250.35黏聚力/MPa0.9802.2000.045内摩擦角/（）27.044.015.52D结果3D结果121086420位移/m05101520时间/s现场调查图4坡顶位移2D、3D模拟结果对比图5滑后地形与2D、3D模拟结果对比375350325300275250Z/m地形线实测剖面线3D结果2D结果050100150200X/m第12期张检召，等：基于SPH的溃屈式滑坡失稳过程分析第58卷21为验证数值模型的可靠性，参考CHANG等9使用验证数值模型的准确性，可表示为模拟结果和实际发生的滑坡堆积体的重叠程度（-22），当=2时表示完美
14、的重叠，当=-2时表示没有重叠。其中的计算公式为：=AxAobserved-AyAobserved-AzAobserved+VxVobserved（8）式中：Ax为预测正确的面积；Ay为预测错误的面积；Az为预测缺失的面积；Aobserved为实际滑动的面积；Vx为预测正确的体积；Vobserved为实际滑动的体积。经计算，3D模拟结果=1.86，精度为93%，2D模拟结果=1.42，精度为71%。与二维SPH模型比较而言，三维模型能够更加真实地反映自然状态下滑坡体的受力条件，反演滑坡从启动到堆积的过程，并具有较好的仿真结果，计算结果见表2。滑坡体塑性应变云图见图6。受地形地貌及结构面的影响，
15、坡脚首先出现破坏，发生轻微弯曲变形 图6（a）；随着变形的加剧，由于滑坡坡脚无临空面，下滑受阻，滑坡处于弯曲隆起阶段，滑坡中部出现剪切裂缝，后缘发生拉张破坏 图6（b）；经过一段时间的变形，滑坡塑性应变不断增大，坡顶拉张开裂加剧，发生整体失稳下滑 图6（c）。滑坡体塑性变形特征说明了滑坡受地形及软弱结果面的影响作用较为明显，软弱结构面为滑坡发生溃曲式破坏提供了发展空间。表2数值精确性参数确定参数Ax/m2Ay/m2Ay/m2Vx/（104m3）Aobserved/m3Vobserved/（104m3）3D1 838.1940.65104.8634.751 898.6935.901.862D1
16、534.260.00366.820.151 898.690.191.42图6滑坡动力学过程模拟应变云图（c）t=7.9 s（d）t=20.0 sEffective Plastic Strain0.0060.0050.0040.0030.0020.0010.0000.1380.1150.0920.0690.0460.0230.000Effective Plastic StrainEffective Plastic Strain1.2331.0280.8220.6170.4110.2060.000为研究滑坡体在失稳过程中的运动特征，在滑坡体上布置了监测点S1S9（图2）。比较了滑坡体表层监测点S1S3、滑坡体内监测点S4S6以及滑面监测点S7S9的最大位移（图7）。从图7中可以看出，随着滑坡体内深度的增加，其各监测点的最大位移总体呈减小的趋势，曹琰波等10在岩质（a）t=0.1 s（b）t=1.2 sEffective Plastic Strain0.7640.6370.5090.3820.2550.1270.0002022年第12期甘肃水利水电技术第58卷22滑坡的失稳运动模拟中也发现表