基于零厚度粘聚力单元的降雨入渗下岩质边坡稳定性分析
(1. 武汉大学土木建筑工程学院, 武汉 430072
2. 中国电建集团中南勘测设计研究院有限公司, 长沙 410014)
摘要: 岩质边坡往往存在着大量节理、裂隙等过水通道, 短时强降雨或持续降雨条件下雨水大量进入坡体, 极易诱发岩质边坡发生滑坡灾害, 因此研究降雨入渗作用下岩质边坡的稳定性具有重要意义。本文以ABAQUS数值仿真平台为基础, 通过嵌入零厚度粘聚力单元实现降雨条件下岩质边坡变形破坏过程的连续-非连续模拟, 以矢量和法安全系数为衡量指标, 重点研究降雨强度、降雨时长和降雨入渗位置对边坡稳定性的影响。结果表明: 随着降雨强度和降雨时长的增大, 边坡体内裂隙深度逐渐变大, 边坡稳定性逐渐降低; 降雨入渗的初始位置距离坡肩越远, 降雨诱发的裂隙延伸长度就越短, 边坡稳定性越好。研究结果为降雨条件下岩质边坡工程治理提供参考。
关键词: 粘聚力单元, 矢量和法安全系数, 岩质边坡, 降雨, 数值模拟
DOI: 10.48014/bcce.20220830001
引用格式: 翁磊, 吴志军, 马亮亮. 基于零厚度粘聚力单元的降雨入渗下岩质边坡稳定性分析[J]. 中国土木工程通报, 2023, 1(1): 1-14.
文章类型: 研究性论文
收稿日期: 2022-08-31
接收日期: 2022-10-02
出版日期: 2023-06-28
0 引言
我国幅员辽阔，67%以上的面积被山区覆盖，滑坡灾害广泛存在且频发，严重威胁着国民生命及财产安全。降雨是诱发边坡失稳最普遍和最主要的自然因素。岩质边坡往往存在大量结构面、节理、裂隙以及断层等导水带，地表雨水很容易通过这些通道进入坡体[1]。坡体雨水渗流作用导致岩体抗剪强度大大降低，并且大量雨水进入岩体后不仅增加了坡体自身重量，还增大了坡体所受荷载，导致坡体稳定性降低，甚至引发严重的大规模滑坡灾害[2]。因此，研究降雨作用对岩质边坡稳定性的影响具有重要意义。
早期由于缺乏理论基础与科学工具，人们主要依靠实践经验，对边坡稳定性评价以长期定性观测为主，并配合地质分析等方法对其进行定性描述。20世纪60年代初，随着理论基础的日渐成熟，边坡稳定性分析中融入了力学机制分析方法，这段时期形成了极限平衡法。近十几年来，随着计算机技术的日渐成熟和数值算法的不断发展，数值分析法在边坡稳定性分析中得到了应用。至此，边坡稳定性研究由实践经验到理论研究、由定性到定量逐渐丰富和深入。
常见的边坡稳定性数值分析方法包含两大类，一类是基于连续介质力学的数值分析方法，包括有限元法、边界元法和有限差分法等，如付建新等[3]基于交替隐式有限差分法，研究了不同降雨条件、边坡内部不同位置的瞬态体积含水率分布。M.Mukhlisin等[4]采用有限元强度折减法，对非饱和土质边坡稳定性进行了分析。刘金龙等[5]利用有限元法探讨了降雨条件下边坡的稳定性，认为裂隙的存在为雨水入渗提供了通道，降低了边坡稳定性。林姗等[6]验证了虚单元法在分析边坡稳定性分析时的可靠性。这类连续介质方法可较好反映边坡在降雨作用下的应力应变演化过程，模拟边坡破坏前弹、塑性变形阶段的演化情况，但是这些连续方法难以模拟边坡裂隙扩展及破坏这一非连续灾变过程[7]。为此，Goodman等[8]采用节理单元及改进后的扩展有限元[9]或广义有限元[10]方法实现了裂纹扩展的模拟，但模拟裂纹的数量有限。另一类方法为基于非连续介质力学的数值分析方法，包括离散元法、非连续变形分析法和数值流形法等，如吴志坚等[11]采用颗粒流方法研究了裂隙型和非裂隙型土坡的动力响应和变形演化过程。张郭新等[12]基于DDA方法，指出水降低了边坡坡脚的稳定性，从而在一定程度上引发了边坡失稳。蒋明镜等[13]基于离散元方法模拟岩质边坡的稳定性，结果表明离散元方法能更真实地反映岩体间的相互作用和边坡破坏的渐进过程。这些基于非连续介质的数值分析方法可以较好地处理滑坡过程中的大变形问题，模拟边坡变形发展以及后期破坏的全过程，可以直观呈现边坡体的滑落、崩解过程[14-16]，获得了较多研究成果，但是非连续方法计算效率较低，不适用于复杂边坡的计算[17,18]。
近年来，结合连续数值分析方法和非连续数值分析方法各自优势而发展起来的连续-非连续数值分析方法在岩土工程领域得到了较多应用[19-22]。鉴于边坡失稳过程是一个从变形到破坏、从连续到非连续演变的过程，本文基于ABAQUS平台及二次开发端口，实现了一种基于零厚度粘聚力单元的连续-非连续数值计算方法，采用矢量和法安全系数并实现其程序化、智能化求解，研究了不同降雨入渗对岩质边坡稳定性的影响规律，成果为降雨条件下岩质边坡工程治理提供理论参考。
1 零厚度粘聚力单元本构模型改进
1.1 零厚度粘聚力单元
本文采用粘聚力单元实现模拟岩体的连续-非连续变形破坏过程，基本思路为:首先在模型中将粘聚力单元沿实体单元边界全局插入，岩体变形阶段的连续性行为由实体单元本构关系表征，岩体发生破裂的不连续行为则由粘聚力单元本构关系表征，从而实现岩体变形破坏的连续-非连续全过程仿真模拟。为此，基于ABAQUS平台，通过嵌入零厚度粘聚力单元从而生成连续-非连续边坡模型。ABAQUS程序中，粘聚力单元的本构模型为双线性本构模型
，如图1所示，模型中的损伤变量D由下列公式计算:
(1)
其中，
为整个加载历史期间有效位移的极值，
为损伤起始时所对应的有效位移，
为材料完全失效时相应的有效位移;独立模式(只有法向或切向发生变形)下粘聚力单元的有效位移为各个方向上相应的相对位移;混合模式下(法向和切向组合变形)下粘聚力单元的有效位移为[24]:
(2)
图1 双线性本构模型粘聚力单元
(a)法向牵引力vs法向位移;(b)切向牵引力vs切向位移
Fig.1 Bilinear constitutive model of the cohesive element
(a)normal traction vs normal displacement; (b)tangential traction vs tangential displacement
图1(a)表示了粘聚力单元的法向本构关系，粘聚力单元的法向牵引力
随着法向相对位移
的增加而线性增大。当法向牵引力达到其峰值
，线性损伤演化法则生效，粘聚力单元的法向刚度
开始不断退化，直到粘聚力单元完全破坏失效(对应失效位移
处)，其法向本构关系可用下式表达:
(3)
式中，
为损伤起始时的法向相对位移。
图1(b)表示了粘聚力单元的切向本构关系，粘聚力单元的两个切向牵引力
和
随着切向相对位移
和
的增加而线性增长。当切向牵引力达到峰值时，线性损伤演化法则生效，粘聚力单元的切向刚度
和
开始逐渐退化，直至发生破坏失效(分别对应失效位移
和
)。粘聚力单元的两个切向本构关系可用下式表达:
(4)
(5)
1.2 本构模型改进及验证
由于ABAQUS程序中粘聚力单元自带的损伤起始判据中没有采用Mohr-Coulomb强度准则来确定最大切向牵引力。因此，粘聚力单元本构无法反映边坡岩体的粘聚力和内摩擦角两个特征参数，无法真实描述边坡岩体的力学性质。因此，本文在粘聚力单元原有双线性本构模型的基础上，引入Mohr-Coulomb强度准则来确定损伤起始时的切向牵引力[35]，即把峰值切向牵引力表征为法向牵引力tn、粘聚力c和内摩擦角φ的函数，如式(6)所示。图2表示了改进后的粘聚力单元模型中峰值切向牵引力、粘聚力、内摩擦角间的关系:
(6)
为验证改进的粘聚力单元本构模型的有效性，进行了5组不同条件下的直剪模拟。图3为直剪模型示意图，一个粘聚力单元连接两个实体单元。实体单元的边界条件为:单元1的左边施加向右的水平推力，上部施加向下的法向压力;单元2四周全部为位移约束。通过模拟，得到不同法向压力下粘聚力单元的峰值切向牵引力，然后判断其是否满足Mohr-Coulomb准则。根据文献[25]，粘聚力单元的主要输入参数如表1所示。
图2 基于Mohr-Coulomb准则的损伤起始判据
Fig.2 Damage initiation criterion based on the Mohr-Coulomb criterion
图3 改进的粘聚力单元本构模型验证
Fig.3 Validation of the modified constitutive model of cohesive element
表1 粘聚力单元参数取值
Table 1 Values of the cohesive element parameters
|
单元类型 |
参数名称 |
取值 |
单位 |
|
实体单元 |
弹性模量E |
1×1010 |
Pa |
|
|
泊松比v |
0.25 |
|
|
|
密度ρ |
2700 |
kg/m3 |
|
粘聚力单元 |
法向刚度E33 |
1×109 |
N/m3 |
|
|
切向刚度G13 |
1×108 |
N/m3 |
|
|
切向刚度G23 |
1×108 |
N/m3 |
|
|
内聚力c |
3×104 |
Pa |
|
|
内摩擦角φ |
20 |
° |
设置5组法向压力分别为0.05MPa、0.10MPa、0.15MPa、0.20MPa和0.25MPa。模拟时，首先在单元1上边界以0.2kN/s的加载速率施加法向压力直至目标值，然后采用位移控制方式在单元1左边界以0.01mm/s的速度向右进行加载直到粘聚力单元发生破坏。
另外，对改进的粘聚力单元进行了法向拉伸模拟，得到了单元拉伸破坏的峰值拉应力。图4表示了不同法向压力下峰值切向牵引力模拟结果，图中红色直线为线性拟合结果。同时，图4给出了法向拉伸作用下的峰值拉伸力。由模拟结果可知，内摩擦角φ为20.6°、内聚力为0.034MPa，与表1中输入的参数非常接近，证实了本文基于Mohr-Coulomb强度准则改进的粘聚力单元模型的有效性。
图4 算例结果验证
Fig.4 Validation of the results of the algorithm
2 矢量和法安全系数及算法实现
2.1 矢量和法安全系数计算方法
常见的计算边坡安全系数的方法主要包括极限平衡法和强度折减法。前者是沿滑动面对力积分的代数和求解，对于直线型及圆弧型边坡滑面，安全系数的计算根据力矩平衡而得到，具有明确的力学和物理含义，但是对于这两种特殊滑面以外的其他滑面，安全系数计算则缺少科学依据[23]。对于后者，主要是按同一系数同等程度地对岩土体的抗剪强度参数进行折减，这与岩土体的实际强度特性存在较大的差异[26]。针对上述两种方法的不足[27,28]，郑文博等[29]从边坡的应力状态出发，首先采用数值流形方法获得边坡的应力场分布，然后通过网格节点信息构造有向加权图，将边坡稳定性分析问题转化为图论问题，从而寻找边坡最危险滑裂面和安全系数，该方法较好地统一了岩/土质边坡的稳定性分析。另一方面，葛修润院士从边坡稳定问题的力学角度出发，认为边坡破坏时坡体的滑动具有方向性，提出了矢量和法安全系数计算方法[30]。本节将结合本文所使用的连续-非连续边坡模型，对边坡的矢量和法安全系数求解方法进行推导，并开发相应的求解程序。
根据文献[31]，矢量和法安全系数定义为边坡潜在滑动面上提供抗滑力的各力沿潜在滑动趋势方向
投影的代数和
与提供下滑力的各力沿此方向上投影代数和
之比，即:
(7)
图5为矢量和法安全系数的求解示意图，滑动面上任意一点的应力均可分解成垂直于该点切线的正应力σni和平行于该点切线的切应力στi;滑动面上任意一点切线与x轴正向的夹角用αi表示。其中，滑动面上任意一点的下滑力Ti由该点处的正应力σni和切应力στi分别在投影方向θ上分力的合力提供，即:
(8)
另一方面，滑动面上任意一点的抗滑力Ri由该点处的抗剪强度τfi和法向反力σ'ni分别在投影方向θ上分力的合力提供，即:
(9)
因此，根据式(7)~(9)，可以求得矢量和法安全系数。
图5 矢量和法安全系数求解示意图
Fig.5 Schematic diagram of the vector method safety factor solution
2.2 算法实现
根据矢量和法安全系数的计算原理，基于Matlab平台开发了边坡矢量和法安全系数的求解程序，如图6所示。具体来说，主要可分为以下5个步骤:
图6 矢量和法安全系数计算流程图
Fig.6 Flow chart of the vector method safety factor calculation
(1)获取边坡在各种内外载荷(包括自重、降水等)作用下的应力场，本文通过采取连续-非连续数值方法计算得到，获取滑动面上单元节点的位置坐标，用以确定下滑方向;
(2)获取基本信息后，以滑动面上每一个单元网格为积分单位，通过单元上的节点坐标信息，计算出每个微分段的长度Δli及其下滑方向αi;
(3)获得每个微分段下滑方向后，计算边坡体的整体下滑方向;
(4)将每个微分段上的下滑力和抗滑力分别沿整体下滑方向进行投影，求和得到整体下滑力T(θ)和整体抗滑力R(θ);
(5)最后利用式(7)，即可求解边坡矢量和法安全系数K(θ)。
2.3 对比验证
图7 数值求解矢量和法安全系数验证
(a)模型尺寸;(b)节理边坡;(c)边界条件
Fig.7 Validation of numerical solution of the vector method safety factor
(a)model dimensions;(b)joined slope;(c)boundary conditions
为验证基于连续-非连续数值模型得到矢量和法安全系数的有效性，在ABAQUS平台中建立与文献[32]中边坡模型相同的数值模型，该模型含有一条直线型节理，模型的几何尺寸如图7(a)所示。坡面上有一不规则的未脱落体通过节理与边坡主体粘结成一个整体。因此，该边坡的潜在滑动面就是沿着该节理的一条直线。ABAQUS平台中建立的连续-非连续模型如图7(b)所示，在边坡节理位置处嵌入零厚度粘聚力单元，通过定义粘聚力单元的参数模拟节理的力学行为。模型的边界条件如图7(c)所示，边坡左侧边界Ⅰ及右侧边界Ⅲ均为水平位移约束，边坡底部边界Ⅱ为固定位移约束。
模型的主要计算参数如表2所示，其中实体单元的参数与该节理边坡的真实参数保持一致，具体数值来源于文献[32]，粘聚力单元的参数通过试错法反复校核[33,34]得到。校核时对比的理论解为节理内摩擦角φ=20°、粘聚力c=50kPa对应的矢量和法安全系数。表3列出了数值模拟结果、理论解以及误差，可以看到理论值和模拟值误差仅为1.54%，验证了本文基于连续-非连续模型的矢量和安全系数求解方法的有效性。
表2 节理边坡模型材料参数
Table 2 Material parameters of the joined slope model
|
单元类型 |
参数名称 |
取值 |
单位 |
|
实体单元 |
弹性模量E |
1×1010 |
Pa |
|
|
泊松比v |
0.25 |
|
|
|
密度ρ |
2700 |
kg/m3 |
|
粘聚力单元 |
法向刚度E33 |
1.5×109 |
N/m3 |
|
|
切向刚度G13 |
1.5×108 |
N/m3 |
|
|
切向刚度G23 |
1.5×108 |
N/m3 |
|
|
内聚力c |
3×104 |
Pa |
|
|
内摩擦角φ |
20 |
° |
表3 模拟结果与理论值对比表
Table 3 Comparison table of the simulation results with the theoretical values
|
名称 |
安全系数 |
|
理论值 |
2.343 |
|
模拟值 |
2.379 |
|
误差 |
1.54% |
为进一步验证基于连续-非连续数值方法求解矢量和安全系数的普适性，保持内摩擦角20°不变，分别设置6组不同的节理粘聚力，即c=15kPa、25kPa、35kPa、45kPa、55kPa、65kPa开展数值模拟分析。图8汇总了由数值模拟和理论计算得到的边坡安全系数结果。可以看到，不同节理粘聚力工况下，由数值模拟与理论得到的安全系数值非常吻合，进一步验证了基于连续-非数值方法求解矢量和安全系数的正确性，同时验证了粘聚力单元参数的合理性。
图8 不同节理粘聚力下安全系数模拟值与理论值对比
Fig.8 Comparison between the simulated and theoretical values of safety factor for different joint cohesion
3 降雨入渗下岩质边坡稳定性分析
岩质边坡往往存在大量地质结构面、节理、裂隙以及断层等，地表降水很容易通过这些通道进入坡体中，使边坡岩体抗剪强度大大降低，影响边坡整体稳定性，甚至造成大规模滑坡灾害。鉴于粘聚力单元能很直接地模拟液体在裂隙间的渗流运移过程这一优势，本节将借助基于粘聚力单元的连续-非连续方法模拟降雨入渗作用下边坡的稳定性，主要考虑降雨强度、降雨时长和入渗位置对边坡稳定性的影响。
3.1 边坡数值模型
计算模型如图9所示，该边坡坡顶处存在一条与坡肩距离为5m、深度为3m的原生裂隙。输入的岩体物理力学参数与表2保持一致，按照1.1节中的方法在实体单元边界插入零厚度的粘聚力单元，构建连续-非连续数值模型。
图9 边坡模型几何尺寸
Fig.9 Geometric dimensions of slope model
由于受降雨作用之前，边坡体为有初始应力而无初始应变的状态，因此只有得到与实际情况相符合的初始应力场，才能为后续的降雨分析打好基础[36-38]。一般认为，岩体地应力平衡后产生的位移数量级相对较小，达到10-4，可认为对后续数值模拟的影响极小。图10为ABAQUS模拟地应力平衡后竖直方向的位移云图，可以看到边坡竖直方向的位移变化数量级在10-5~10-6。因此，该岩质边坡模型地应力平衡结果满足要求，可用于开展后续数值模拟研究。
通过在坡顶原生裂隙处施加垂直向下的流速模拟降雨入渗边界条件，随着降雨时长和降雨强度的变化，节理单元在孔隙压力和重力作用下发生损伤劣化，甚至发生开裂失效。同时，降雨入渗条件下边坡的潜在滑动面为开裂单元所组成的裂隙面，可通过潜在滑动面上单元和节点的应力状态求解矢量和法安全系数。
图10 地应力平衡后竖直方向位移云图
Fig.10 Vertical displacement distribution after ground stress equilibrium
3.2 降雨强度对边坡稳定性的影响
根据《降水量等级》(GB/T 28592—2012)标准，大雨到特大暴雨的降雨强度分别为R=60mm/d，120mm/d，240mm/d，480mm/d，本模型中坡顶面积为27.4m2。因此，可以计算得到对应上述降雨强度下原生裂隙处雨水入渗速率分别为1.14L/min，2.28L/min，4.56L/min，9.12L/min。
控制降雨时长为48h，分别对模型施加设计好的四种降雨强度R=60mm/d，120mm/d，240mm/d和480mm/d。边坡内裂隙深度演化模拟结果如图11所示。图12汇总了降雨强度和裂隙深度之间的关系。从图中可以发现，降雨强度越大，边坡体内裂隙逐渐向坡脚延长贯穿。这是由于在降雨作用下，裂隙内的水压增大，且雨水在裂隙和岩体进行渗流扩散，在一定程度上软化了岩体和裂隙间的强度并增大了岩体自重，最终导致裂隙向边坡深部发生扩散，且扩散方向逐渐向边坡面偏移，使边坡的稳定性下降[40]。特别是当降雨强度达到480mm/d时，裂隙的宽度和深度都明显增加，深度达到了14.6m，即将形成一条延伸到坡脚的完整滑移带。
图11 不同降雨强度下裂隙深度变化
Fig.11 Changes in the fracture depth under different rainfall intensities
边坡受到不同降雨强度降雨48h后，其水平方向的位移变化如图13所示。从图中可以发现裂隙到坡面之间的坡体在降雨作用下向坡面方向发生了明显的滑移，并且滑移程度与降雨强度呈正相关。随着降雨强度从60mm/d增加到480mm/d，坡体产生的最大水平位移依次为1.9cm，2.6cm，4.5cm，8.0cm，表明降雨强度越大，坡体产生的滑移量越大，其稳定性也随之降低。
图12 降雨强度与裂隙深度的关系
Fig.12 Relationship between the rainfall intensity and fracture depth
图13 不同降雨强度下的坡体水平位移图
Fig.13 Horizontal displacement of slope under different rainfall intensities
为了更准确地表征降雨强度与边坡稳定性的关系，利用2.1节矢量和法安全系数分别计算出了不同降雨强度持续作用48h后边坡的安全系数，结果如图14所示。从图中更加直观地发现，降雨时长相同时，随着降雨强度的逐渐增强，边坡的安全系数不断减小，说明边坡的稳定性随着降雨强度的不断增强而逐渐降低。这是因为降雨强度越大时，短时间内进入裂隙中的雨水就越多，地下水位迅速上升，孔隙水压力随之增大，不仅会在水压劈裂作用下使裂隙进一步扩展延伸，还会对边坡岩体强度造成更大程度的弱化效应，使得边坡稳定性减小。
图14 降雨强度与边坡安全系数的关系
Fig.14 Relationship between the rainfall intensity and slope safety factor
3.3 降雨时长对边坡稳定性的影响
本节在控制降雨强度不变的情况下，探究不同的降雨时长与边坡稳定性之间的关系。降雨强度保持为240mm/d不变，设计了如表4所示的四种降雨时长，入渗位置与图9保持一致。
表4 降雨时长设计
Table 4 Rainfall duration design
|
降雨强度R(mm/d) |
降雨时长t(h) |
降雨量P(mm) |
|
240 |
12 |
120 |
|
|
24 |
240 |
|
|
48 |
480 |
|
|
72 |
720 |
图15给出了不同降雨历时下的边坡体内裂隙开度及长度的变化情况。从图中可以发现，随着降雨时长的增加，边坡体内的裂隙长度在原生裂隙的基础上逐渐扩展、不断发育变长。在降雨时长t=12h时，裂隙的长度扩展延伸到7.0m深，裂隙延伸的方向几乎没有改变;然后随着降雨时长的继续增加，在t=24h时，裂隙深度延伸到9.8m，并且裂隙延伸的方向由原来的竖直向下开始向坡脚方向发生偏移;随着降雨时间的继续增加，到t=48h时，发现裂隙扩展的方向已经显著向坡脚发生了偏转，此时裂隙在竖直方向上的长度已经达到了12.8m;随着降雨持续到t=72h，裂隙的深度达到了16.9m，已经形成了一条明显延伸到坡脚的滑移带，随着降雨的继续增加，此滑移带极有可能在坡脚处贯通，形成完整的滑裂面，边坡由此而发生失稳破坏。
图15 不同降雨时长下裂隙深度变化
Fig.15 Changes in the fracture depth under different rainfall duration
图16给出了不同降雨时长下，边坡体内裂隙长度的发育情况。从图中可以发现裂隙长度除了随降雨时长的增加而增加外，其增加幅度还有明显的不同;在t=12~24h期间，边坡体内的裂隙发展较快，增长幅度较大;随着降雨的持续，在t=24~48h期间，裂隙扩展延伸的程度趋于平缓;然后随着降雨时长的继续增加，在t=48~72h期间，裂隙的增长幅度又变得加快，但是没有降雨时间为t=12~24h期间的增长幅度大。说明降雨在前期12~24h和后期48~72h对边坡的稳定性影响较大。
图16 降雨时长与裂隙深度的关系
Fig.16 Relationship between the rainfall duration and fracture depth
不同降雨时长条件下水平方向位移变化如图17所示，可以发现原生裂隙到坡面之间的坡体在降雨作用下产生了滑移，并且滑移程度与降雨强度呈正相关，间接反映了边坡稳定性的降低。以上模拟结果是因为随着降雨时间增加，雨水入渗量越大，裂隙水压力和岩体软化程度越高，岩质边坡所受重力也增加，边坡裂隙扩展就越深。利用2.1节的矢量和法安全系数，计算出了降雨时长分别为12h、24h、48h、72h时对应的边坡安全系数，结果如图18所示。结果表明，随着降雨时长的增加，边坡的安全系数逐渐减小，边坡稳定性逐渐降低。
3.4 降雨入渗位置对边坡稳定性的影响
本节在控制降雨强度和降雨时长不变的情况下，探究不同的入渗位置对边坡稳定性的影响规律。保持降雨强度为240mm/d，降雨时长为48h，设计了三种不同的入渗位置来探究降雨入渗位置对边坡稳定性的影响。以坡肩为基准分别取距坡肩5m，10m和15m的三个初始入渗位置，裂隙深度均为3m，如图19所示。
图17 不同降雨时长下坡体水平位移图
Fig.17 Horizontal displacementdistribution of the slope under different rainfall duration
图18 降雨时长与边坡安全系数的关系
Fig.18 Relationship between the rainfall duration and slope safety factor
图19 不同降雨入渗位置示意图
Fig.19 Schematic diagram of different rainfall infiltration locations
不同入渗位置和裂隙长度的关系的模拟结果如图20所示。从图中发现相同降雨作用下，入渗的初始位置不同，原生裂隙扩展延伸的深度会有很大差异。当降雨入渗的初始位置为距坡肩a1=5m处的位置Ⅰ时，边坡体内原生裂隙扩展延伸的深度达到了12.8m，并且裂隙延伸的方向沿坡面发生了明显的偏转;当降雨初始入渗位置为距坡肩a2=10m处的位置Ⅱ时，原生裂隙延伸的深度明显减少，只有7.5m长;当降雨初始入渗位置为距坡肩a3=5m处的位置Ⅲ时，原生裂隙扩展延伸的深度继续减少，只达到5.9m。对比三种不同的降雨入渗位置，发现降雨入渗的位置对边坡体内原生裂隙延伸的深度有很大影响;入渗的位置距离坡肩越远，坡体内原生裂隙延伸长度就越短，说明对边坡稳定性的影响也就越小。
图20 降雨入渗位置与裂隙深度的关系
Fig.20 Relationship between the rainfall infiltration location and fracture depth
边坡体内原生裂隙延伸深度与降雨入渗位置间的关系如图21所示。从图中很直观地发现，随着降雨入渗初始位置到坡肩的距离越来越远，裂隙延伸的深度变得越来越短;从侧面反映出降雨入渗的初始位置距离坡肩越远，对边坡稳定性的影响也就越小。
图21 降雨入渗位置与裂隙长度的关系
Fig.21 Relationship between the rainfall infiltration location and fracture length
相同强度、时长的降雨作用在不同入渗位置处时边坡体水平位移变化如图22所示。从图中发现不同的降雨入渗位置会导致坡体产生不同的水平位移，位置Ⅰ、Ⅱ、Ⅲ所对应的最大水平位移分别为4.5cm、2.0cm、1.7cm;当降雨入渗的位置距离坡肩越近时，产生的水平位移越大;随着降雨入渗的位置逐渐远离坡肩，坡体产生的水平位移随之逐渐降低。说明降雨入渗位置越远离坡肩，对边坡稳定性的影响越小。
采用边坡矢量和法安全系数，分别计算出相同条件的降雨从三个不同位置入渗后所对应的边坡安全系数，计算结果如图23所示。图中的结果表明随着降雨入渗位置到坡肩的距离由5m变为10m，再变到更远的15m，边坡的安全系数逐渐增大，边坡变得越来越稳定。说明降雨入渗的位置能影响边坡的稳定性，降雨入渗位置距离坡肩越远，降雨对边坡稳定性的影响就越小。
图22 不同降雨入渗位置下的坡体水平位移图
Fig.22 Horizontal displacement of the slope under different rainfall infiltration locations
图23 降雨入渗位置与边坡安全系数的关系
Fig.23 Relationship between the rainfall infiltration location and slope safety factor
4 结论
(1)开发了基于ABAQUS平台的连续-非连续算法，对边坡矢量和法安全系数进行求解，实现了矢量和法安全系数的程序化、智能化求解，极大提高了求解效率。
(2)随着降雨强度的逐渐增强，边坡体内的裂隙深度逐渐扩展加深，边坡的安全系数逐渐减小，稳定性降低。
(3)随着降雨时长的增加边坡体内原生裂隙不断发育变长，边坡的安全系数随之逐渐减小，稳定性逐渐降低。
(4)降雨入渗的初始位置距离坡肩越远，原生裂隙延伸的长度越短，边坡的安全系数随之而增大，稳定性越高。
利益冲突: 作者声明无利益冲突。
[①] *通讯作者 Corresponding author:吴志军，zhijunwu@whu.edu.cn
收稿日期:2022-08-31; 录用日期:2022-10-02; 发表日期:2023-06-28
基金项目:国家自然科学基金项目(42077246，52004182)资助
参考文献(References)
[1] Senthilkumar V, Chandrasekaran S S, Maji V B. Rainfallinduced landslides: Case study of the Marappalam Landslide, Nilgiris District, Tamil Nadu, India[J]. International Journal of Geomechanics, 2018, 18(9): 05018006.
https://doi.org/10.1061/(ASCE)GM.1943-5622.0001218
[2] 王述红, 何坚, 杨天娇. 考虑降雨入渗的边坡稳定性数值分析[J]. 东北大学学报(自然科学版), 2018, 39(8): 1196-1200.
https://doi.org/10.12068/j.issn.1005-3026.2018.08.026
[3] 付建新, 宋卫东, 杜建华. 考虑二维降雨入渗的非饱和土边坡瞬态体积含水率分析[J]. 工程科学学报, 2015, 37(4): 407-413.
https://doi.org/10.13374/j.issn2095-9389.2015.04.002
[4] Mukhlisin M, Baidillah M R, Ibrahim A, et al. Effect of soil hydraulic properties model on slope stability analysis based on strength reduction method[J]. Journal of the Geological Society of India, 2014, 83(5): 586-594.
https://doi.org/10.1007/s12594-014-0087-1
[5] 刘金龙, 栾茂田, 王吉利, 等. 降雨条件下土坡饱和-非饱和渗流及稳定性分析[J]. 岩土力学, 2006, 27(S1): 103-107.
https://doi.org/10.16285/j.rsm.2006.s1.058
[6] 林姗, 郭昱葵, 孙冠华, 等. 边坡稳定性分析的虚单元强度折减法[J]. 岩石力学与工程学报, 2019, 38(S02): 3429-3438.
https://doi.org/10.13722/j.cnki.jrme.2019.0030
[7] 王学滨, 芦伟男, 钱帅帅, 等. 静水压力条件下开挖直径及卸荷时间对巷道围岩变形-开裂影响的连续-非连续方法模拟[J]. 应用力学学报, 2020, 37(04): 1841- 1848, 1880.
https://doi.org/10.11776/cjam.37.04.B023
[8] Goodman R E, Taylor R L, Brekke T L. A model for the mechanics of jointed rock[J]. Journal of The Soil Mechanics and Foundations Division, 1968, 94(SM3): 637-659.
https://doi.org/10.1061/JSFEAQ.0001133
[9] Moёs N, Belytschko T. Extended finite element method for cohesive crack growth[J]. Engineering Fracture Mechanics, 2002, 69(7): 813-833.
https://doi.org/10.1016/S0013-7944(01)00128-X
[10] Strouboulis T, Babuška I, Copps K. The design and analysis of the generalized finite element method[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 181(1): 43-69.
https://doi.org/10.1016/S0045-7825(99)00072-9
[11] Wu Z, Zhang D, Wang S, et al. Dynamic-response characteristics and deformation evolution of loess slopes under seismic loads[J]. Engineering Geology, 2020, 267: 105507.
https://doi.org/10.1016/j.enggeo.2020.105507
[12] 张国新, 雷峥琦, 程恒. 水对岩质边坡倾倒变形影响的DDA模拟[J]. 中国水利水电科学研究院学报, 2016, 14(03): 161-170.
https://doi.org/10.13244/j.cnki.jiwhr.2016.03.001
[13] 蒋明镜, 江华利, 廖优斌, 等. 不同形式节理的岩质边坡失稳演化离散元分析[J]. 同济大学学报(自然科学版), 2019, 47(02): 167-174.
https://doi.org/10.11908/j.issn.0253-374x.2019.02.002
[14] 李祥龙, 唐辉明, 胡巍. 层面参数对顺层岩质边坡地震动力破坏过程影响研究[J]. 岩土工程学报, 2014, 36(03): 466-473.
https://doi.org/10.11779/CJGE201403009
[15] Fakhimi A, Wan F. Discrete element modeling of the process zone shape in mode I fracture at peak load and in post-peak regime[J]. International Journal of Rock Mechanics and Mining Sciences, 2016, 85: 119-128.
https://doi.org/10.1016/j.ijrmms.2016.03.014
[16] 沈华章, 郭明伟, 王水林, 等. 基于离散元的边坡矢量和稳定分析方法研究[J]. 岩土力学, 2016, 37(02): 592-600.
https://doi.org/10.16285/j.rsm.2016.02.033
[17] 刘泉声, 邓鹏海, 毕晨, 等. 深部巷道软弱围岩破裂碎胀过程及锚喷-注浆加固FDEM数值模拟[J]. 岩土力学, 2019, 40(10): 4065-4083.
https://doi.org/10.16285/j.rsm.2018.1032
[18] 严成增. 模拟水压致裂的另一种二维FDEM-flow 方法[J]. 岩土力学, 2017, 38(06): 1789-1796.
https://doi.org/10.16285/j.rsm.2017.06.029
[19] Bicanic N, Munjiza A, Owen D R J, et al. From Continua to Discontinua-A Combined Finite Element/Discrete Element Modelling in Civil Engineering[C]. Interna-tional Conference on Computing in Civil and Structural Engineering, 1995: 1-13.
[20] Mahabadi O K, Lisjak A, Munjiza A, et al. New combined finite-discrete element numerical code for geomechanical applications[J]. International Journal of Geomechanics, 2012, 12(6): 676-688.
https://doi.org/10.1061/(ASCE)GM.1943-5622.0000216
[21] Zhang H H. Analysis of Crack Interaction Problem by the Numerical Manifold Method[J]. Advanced Materials Research, 2012, 446-449: 797-801.
https://doi.org/10.4028/www.scientific.net/AMR.446-449.797
[22] Yang Y, Zheng H. Direct Approach to treatment of contact in numerical manifold method[J]. International Journal of Geomechanics, 2016, 17: 4016012.
https://doi.org/10.1061/(ASCE)GM.1943-5622.0000714
[23] Alfano G. On the influence of the shape of the interface law on the application of cohesive zone models[J]. Composites Science & Technology, 2004, 66(6): 723-30.
https://doi.org/10.1016/j.compscitech.2004.12.024
[24] Camanho P P, Davila C G. Mixed-Mode Decohesion Finite Elements for the Simulation of Delamination in Composite Materials[R]. NASA/TM-2002-211737 2002.
[25] 甘亮. EGS-E水平巷道水力压裂与热交换数值模拟初步研究[D]. 武汉: 武汉大学, 2019.
[26] 郭明伟, 葛修润, 李春光, 等. 基于矢量和方法的边坡稳定性分析中整体下滑趋势方向的探讨[J]. 岩土工程学报, 2009, 31(4): 577-583.
https://doi.org/1000-4548(2009)04-0577-07
[27] 郑颖人, 赵尚毅. 有限元强度折减法在土坡与岩坡中的应用[J]. 岩石力学与工程学报, 2004, 23(19): 3381-3388.
https://doi.org/1000-6915(2004)19-3381-08
[29] 郑文博, 庄晓莹, 李耀基, 等. 基于流形方法和图论算法的岩/土质边坡稳定性分析[J]. 岩土工程学报, 2013, 35(11): 2045-2052.
https://dx.doi.org/10.3969/j.issn.1000-4548(2013)11-2045-08
[30] 丰定祥, 吴家秀, 葛修润. 边坡稳定性分析中几个问题的探讨[J]. 岩土工程学报, 1990, 12(3): 1-9.
https://dx.doi.org/10.19345/j.cnki.1671-0037.2013.06.057
[31] 葛修润. 岩石疲劳破坏的变形控制律、岩土力学试验的实时X射线CT扫描和边坡坝基抗滑稳定分析的新方法[J]. 岩土工程学报, 2008, 30(1): 1-20.
https://dx.doi.org/1000-4548(2008)01-0001-20
[32] 沈华章. 岩土应变软化分析及其在边坡工程中的应用 [D]. 北京: 中国科学院大学, 2016.
[33] Koyama T, Jing L. Effects of model scale and particle size on micro-mechanical properties and failure processes of rocks-A particle mechanics approach[J]. Engineering Analysis with Boundary Elements, 2007, 31(5): 458-472.
https://doi.org/10.1016/j.enganabound.2006.11.009
[34] Zheng W, Zhuang X, Tannant D D, et al. Unified continuum/ discontinuum modeling framework for slope stability assessment[J]. Engineering Geology, 2014, 179: 90-101.
https://doi.org/10.1016/j.enggeo.2014.06.014
[35] Ma G, Wei Z, Chang X L. Modeling the particle breakage of rockfill materials with the cohesive crack model[J]. Computers and Geotechnics, 2014, 61: 132-143.
https://doi.org/10.1016/j.compgeo.2014.05.006
[36] 马云峰. 基于ABAQUS 岩土工程中地应力平衡的探讨[J]. 科技与创新, 2014(8): 61-61.
https://doi.org/10.15913/j.cnki.kjycx.2014.08.006
[37] 刘爱华, 杨清, 吴均平. ANSYS三维地应力场数值模拟方法应用研究[J]. 地质力学学报, 2013, 19(2): 133-142.
https://doi.org/1006-6616(2013)02-0133-10
[38] 代汝林, 李忠芳, 王姣. 基于ABAQUS的初始地应力平衡方法研究[J]. 重庆工商大学学报(自然科学版), 2012, 29(9): 76-81.
https://doi.org/1672-058X(2012)09-0076-06
Stability Analysis of Rock Slope under Rainfall Infiltration Based on A Zero-Thickness Cohesive Element
(1. School of Civil Engineering, Wuhan University, Wuhan 430072, China
2. Powerchina Zhongnan Engineering Corporation Limited, Changsha 410041, China)
Abstract: There are often a large number of joints, fissures and other water channels in rock slopes. Under the condition of short-term heavy rainfall or continuous rainfall, a large amount of rainwater enters the slope body, which can easily induce landslide disaster on rock slopes. Therefore, it is of great significance to study the stability of rock slope under the infiltration of rainfall. Based on the ABAQUS numerical simulation platform, this paper implements continuous-discontinuous simulation of the deformation and damage of rock slopes under rainfall conditions by embedding zero-thickness cohesive elements and focuses on the effects of rainfall intensity, rainfall duration and rainfall infiltration location on the slope stability using vector method safety factors as indicators. The results show that as the rainfall intensity and rainfall duration increase, the depth of fissures in the slope gradually becomes larger, and the stability of the slope gradually decreases. The farther the initial location of rainfall infiltration is from the shoulder of the slope, the shorter the rain-induced fissure extension length is and the better the stability of the slope is. The research results provide reference for the engineering treatment of rock slope under rainfall conditions.
Keywords: Cohesive element, vector method safety factors, rock slope, rainfall, numerical simulation
DOI: 10.48014/bcce.20220830001
Citation: WENG Lei, WU Zhijun, MA Liangliang. Stability analysis of rock slope under rainfall infiltration based on a zero-thickness cohesive element[J]. Bulletin of Chinese Civil Engineering, 2023, 1(1): 1-14.