1. 引言
滑坡是指在降雨、地震、地下水、人类工程活动等作用下，岩土体沿一定软弱面向下滑动的自然现象[1]，我国作为一个多山国家，长期遭受地质灾害的侵扰[2]，其中滑坡造成的灾害损失尤为显著。其中降雨作为诱发滑坡的常见因素之一，众多学者通过对降雨型滑坡进行了不断的探索。
对于滑坡的研究，随着计算机技术的发展，数值模拟技术因其成本低、操作方便、可重复试验等优点，成为了反演分析滑坡运动过程和动力学机制的最为高效的手段之一[3]。因此现阶段大量专家学者多采用数值模拟手段对滑坡破坏过程进行研究。离散元数值方法为代表的非连续数值方法块体之间相互独立，运动方程随迭代时间步计算不断更新，计算过程不受变形量限制，可有效地模拟介质的开裂、分离等非连续现象，可以反映岩土工程中的变形机理、破坏过程、变形破坏造成的不良影响等[4]。李龙起等[5]采用接触黏结模型，通过颗粒离散元的方法研究某土质边坡在降雨作用下的破坏模式及坡内土体的力学响应，探究坡体破坏模式与重要物理量之间的内在联系。常晁瑜[6]等以下马达子滑坡为研究对象，基于颗粒流数值模拟软件对黄土地震滑坡和地震液化型滑坡进行对比研究，提出了地震滑坡的失稳变形模式及液化滑坡在地震作用下失稳破坏特征，地震滑坡与地震液化型滑坡防治以及后续类似工况滑坡研究提供参考。赵州等[7]运用颗粒流软件方法，针对滑坡的破坏运动过程展开了系统而深入的研究。在研究过程中，采用颗粒离散元数值模拟软件进行双轴压缩试验，对新府山滑坡模型的细观参数进行了标定，从破坏特征和滑坡强度两方面，探究了新府山滑坡的变形破坏机理与运动特征，为该滑坡的深入认识和后续治理提供了重要的理论依据。吴伟乐[3]等应用离散元方法，建立三维地质模型，对在强降雨作用下的碎屑岩滑坡远程运动成灾模式开展研究，为防灾减灾提供定量化科学依据。周健等[8]利用颗粒流的方法研究细观参数取值差异对于土质边坡破坏形式的影响，展示了颗粒流用于观察岩土体的破坏过程的优势。
本文以河南省伊川县某滑坡为研究对象，以前期的地质勘察等资料为基础，利用颗粒离散元软件对滑坡在降雨作用下的变形破坏过程进行模拟，监测滑坡内不同位置的位移、速度的变化情况，分析该滑坡在运动过程中不同阶段下的破坏变形特征，为该滑坡的预防治理提供参考依据。
2. 研究区工程概况
2.1. 滑坡规模及形态特征
研究区位于伊川县酒后镇南部，属低山丘陵地貌，低山区地形起伏较大，沟谷深切，坡体走向为近东西方向，坡体倾向为近似正北方向，地形呈南高北低，滑坡体前缘高程264 m，后缘最高程为350 m，坡体相对高差约为85 m。滑坡所处坡体自然坡度23~29˚，根据地质调查，此地为老滑坡，历史上有多次滑动，目前仍有少量滑动，坡体呈凹型，台阶状，坡体上为当地村民为防止坡体滑动种植的树木，植被较为发育。滑坡体纵长约590 m，平均宽度约400 m，面积0.243 km2。滑坡体滑体厚度2 m~21.5 m，体积约75 × 104 m3，主滑坡方向约342˚。滑坡形态规模如图1所示，根据滑体体积划分为中型滑坡。滑坡具体形态规模见图1。
滑坡变形与场区的地层岩性、地形地貌、水的作用和人类工程活动有极大的关系，滑坡在枯水季节处于基本稳定状态，在丰水季节局部处于不稳定状态，有变形发生。
Figure 1. Scale of landslides
图1. 滑坡形态规模
据该地区房屋是依靠山坡就势而建，呈明显的西北低，东南高。位于坡体前缘下方的部分房屋、院墙或地面发生不同程度的开裂现象，且基本上都垂直于滑坡方向见照片，开裂严重的房子无法居住的情况下已搬离被废弃，开裂不严重的房屋，对墙体裂隙进行修补后仍在其居住。
坡体呈台阶状，滑坡体部分后缘近似裸露，便于雨水通过后缘下渗，在滑坡体上种植有大量的树木，生长有大量杂草，并且滑坡体部分部位还存在凹陷不利于地表水排泄，从而增大了地表水下渗，增大了地下水的赋存量，降低了滑坡的抗滑力，加大了滑坡的下滑力。并且滑坡体前缘部分已被当地村民建房时开挖，造成坡体向不利的方向发展，从而坡体滑移变形，形成滑坡，由此可见，暴雨或持续降雨是引起滑坡的关键因素。
2.2. 滑坡结构特征
根据钻探揭露及野外现场鉴别，地层上部为第四系粉质黏土，下伏古近系洛阳组红褐色泥质砂岩，工程地质剖面图见图2，各地层的特征分述如下。
Figure 2. Landslide profile
图2. 滑坡剖面图
滑坡体下伏泥质砂岩层处于稳定状态，斜坡表层的粉质黏土在天然状态下处于基本稳定状态、在连续降雨、暴(大)雨的影响下局部处于不稳定状态。滑坡体下伏强风化岩层的倾向与坡体方向大体一致，为顺向坡，未见软弱结构面或软弱夹层，不会产生基岩层内的滑动。
3. 数值模拟计算
3.1. 数值模型建立
通过对现场勘察资料的深入分析，已经确定研究区滑坡的滑动面。在此基础上，本研究选用Ball-Wall法构建模型，利用颗粒单元对滑体进行填充，借助墙单元来模拟滑床，详细建模流程如下所示。
首先根据滑坡的工程地质剖面图，绘制1∶1的几何模型；将其导入颗粒流软件中，并生成相应的墙体边界；在墙体内生成的颗粒，并设置颗粒大小与孔隙率，采用颗粒最小半径为0.1 m，最大半径为0.15 m，共生产111387个颗粒，并经过迭代计算，使其平衡；为土体颗粒赋予颗粒流软件中所使用的细观参数，最后删除相应滑体的边界约束墙；对滑坡模型进行速度、位移、颗粒间接触力与接触力矩清零，使其达到最初状态，在滑坡体的不同部位设置5个监测点，监测其位置颗粒的位移与速度，分析滑坡的运动特征；最后，对坡体施加重力，使坡体在重力作用下开始缓慢变形。数值模型见图3。
Figure 3. Numerical model
图3. 数值模型
3.2. 微观参数选取
Table 1. Granular microparameters of slope rock and soil
表1. 滑坡体岩土体颗粒细观参数
密度/ (g/cm−3) |
有效模量 |
法向刚度/(N/m) |
切向刚度/(N/m) |
法向粘结刚度/(N/m) |
切向粘结刚度/(N/m) |
胶结抗拉强度/Pa |
胶结抗剪强度/Pa |
1780 |
2.8e6 |
1.3e7 |
1.3e7 |
6.1e7 |
6.1e7 |
6.7e4 |
6.7e4 |
颗粒离散元软件中参数采用的是颗粒之间的微观参数，而实际岩土体中的宏观参数很难与微观参数对应，颗粒不能直接实现宏观材料的性质，但二者之间存在一定的联系[9]。鉴于降雨是导致坡体滑动的关键诱发因素，所以对处于降雨工况下的坡体处展开模拟计算，采用线性粘结模型，根据前人的研究并结合现场实际情况，滑体在降雨工况下的微观力学参数见表1。
3.3. 滑坡破坏特征分析
为探索滑坡的破坏特征，对滑坡运动过程中不同时间步长下的位移以及速度状态的数据进行保存，通过观察不同关键时间节点下滑坡的位移以及速度状态，分析滑坡破坏过程、总结滑坡的破坏特征。
Figure 4. 10 000, 30000 and 50000 step velocity cloud map
图4. 1万、3万与5万时步的速度云图
通过观察速度云图，运行至1万时步时，坡体中前部的颗粒为速度最大区域，因为该部位坡体的坡度最大，在重力的作用下产生的影响最大，滑坡中部向上部位的近端出现少量裂纹，并且有逐渐扩大的趋势，意味着坡体正处于蠕滑变形阶段。滑坡后缘的速度其次，有向下滑动的趋势，此时坡体的速度分布以中部最大的坡度转折点为界大致分成了两段，表现形式较为一致，速度大小从高处到低处逐渐减小。运行至3万时步时，相较于1万时步，整体速度增加量较少，是因为建模时给颗粒赋值，原始的平衡状态被打破，颗粒为了重新达到平衡状态而快速运动，最大速度分布对于1万时步时的位置向前，滑坡中部位置裂隙数量增加，说明随着时间的推移，滑坡中前部的颗粒蠕滑变形加剧；滑坡后缘部位的颗粒速度增加缓慢。运行至5万时步时，速度分布较前阶段发生明显变化，说明随着时间的推移，最大速度分布逐渐从滑坡中部向坡脚过渡；伴随重力的作用下，对坡脚的挤压效应也愈发明显，坡脚处的位移增大，并有少量的颗粒从坡脚处被挤压出坡体；滑坡后缘部位的颗粒速度减小。此时坡体的速度分布逐渐转变，以中部最大的坡度转折点为界大致分成了速度从高处到低处逐渐增大的两部分。滑坡不同时间步长下的速度云图见图4。
Figure 5. 80000 and 130000 step displacement cloud map
图5. 8万及13万时步位移云图
通过观察位移云图，当运行至8万时步时，滑坡前缘部位的坡体在其后部坡体的推动下向滑床外部滑动。滑坡中部的裂纹逐渐增多并向坡内发展，并且滑坡后缘出现拉张裂缝，且有逐渐贯通的趋势。运行至13万时步，滑坡后缘出现滑塌现象，滑体呈现出整体下滑的趋势，坡脚颗粒滑出近10米，形成蠕滑–拉裂型滑坡。滑坡不同时间步长下的位移云图见图5。
3.4. 滑坡运动特征分析
如图3中圆点所示，分别选取坡体中的不同位置，其圆心处颗粒设为监测点，依次设置5个监测点，颗粒坐标分别为(34.9086, 101.4565)，(140.9658, 86.918)，(252.3173, 77.6999)，(346.5243, 57.2863)，(492.811, 29.6542)，各个位置监测点编号如图3所示，通过监测颗粒的位移与速度随时步的变化来分析整个坡体的运动过程，根据得到的曲线图与边坡整体变形破坏特征相结合进行分析，研究坡体破坏特征与位移、速度之间的关系。位移–时步曲线图与速度–时步曲线图见图6与图7。
Figure 6. Displacement change curve
图6. 位移变化曲线
Figure 7. Velocity change curve
图7. 速度变化曲线
根据图6中的位移–时步曲线图并结合速度与位移云图可知，前期5个监测点记录到的位移变化缓慢，滑坡此时土体存在微小蠕变，4万步后，监测点位移增速普遍加快，D监测点的位移要大于其他部位的监测数据，是因为D监测点所处的位置坡率最大，加速度较大，产生较大位移，E监测点处于坡脚位置，同时受到D监测点附近颗粒的推挤，位移变化较A、B、C监测点也较为明显，同时与位移云图的结果相对应。运行至8万步，A监测点的位移增速明显变缓，由于前期滑坡前缘速度相对其他部位较大，产生较大位移，由于后缘土层较薄，同时需向前推移土体，故位移增速逐渐减缓；由于斜坡中前缘滑动较快，使后缘产生张拉裂缝，随着裂缝逐渐扩张，后缘失去了前部的牵引作用，仅靠重力驱动向下滑塌，虽然A监测点位移增速有所加快，但其与前部的位移差越来越大；后续随着滑动面被贯通，所有监测点均位移呈现线性增大的趋势，有整体下滑的趋势。
通过观察图7速度–时步曲线图可知，模拟开始初期，不同部位的颗粒的运动速度变化都呈现出增加的趋势，其中D监测点的增速最快，其次为A监测点，与速度云图所对应。7万至8万时步，A、D监测点速度进度下降阶段，是由于A、D两监测点区域颗粒以较大速度向前滑移，率先开始向其前方颗粒挤压，导致速度降低，其中A监测点速度骤降，由于后缘土层较前方的土层较薄，相对于D点的下降趋势更为明显，与位移–时步曲线图相对应，后续A、D两监测点速度有所增快，表明后缘局部形成张拉裂隙，前缘颗粒开始滑出滑床。B、C、E受到各自后方速度较快的颗粒推移，表现出速度持续增加的状态，随着时间的推移，挤压作用愈发减小，速度在短时间内减小，但在滑动面逐渐贯通，速度变化较为稳定，有整体下滑的趋势。
4. 结论
本文以河南省伊川县某滑坡为研究对象，通过建立1∶1的数值模型并基于颗粒离散元的数值计算方法，研究该滑坡在降雨工况下的破坏模式及运动特征，得出结论如下。
(1) 通过观察不同时步下的速度与位移云图，该滑坡在降雨工况下整体破坏特征为：中部剪切–前缘牵引–后缘拉裂，表现为蠕滑–拉裂式破坏。
(2) 通过观察位移–时步曲线图与速度–时步曲线图，前期位移增速较慢，坡体处于蠕滑阶段，坡体中前部的位移及速度最大，该位置率先破坏，其次为后缘部位；中前部与后缘因速度差挤压前方颗粒，其区域的土体速度降低；随着时间推移，前缘部位颗粒开始剪出，滑出滑床，后缘出现张拉裂缝，出现滑塌现象，故速度开始增加；后续速度变化较为稳定，有整体下滑的趋势。