0 引言

中国西部山区构造运动、风化作用等地质营力极为强烈，复杂地质环境和频繁人工活动作用孕育形成大量层状碎裂结构边坡^[1]。层状碎裂结构边坡具有节理面密集分布、岩体完整性差、稳定性差等显著特点，其滑动灾变过程更具强隐蔽性、大危害性、难防控性^[2]，层状碎裂结构边坡的滑动灾变过程特征是岩石力学工程和地质灾害防治领域研究的方向之一，可为边坡稳定性评价、灾害早期预警和防治措施提供理论支撑。

数值模拟包括连续介质方法和非连续介质方法，是解决重大科学问题和复杂工程难题必不可少的手段，以其重复性强、数据丰富等优势成为研究边坡变形破坏过程的有效工具之一。连续介质方法以均匀性假设为基础, 在描述节理岩体及其构成边坡的变形破坏过程和微观机制方面存在诸多缺陷^[3]；而离散单元法(Distinct Element Method，DEM)旨在探索岩体等非连续介质材料力学行为，力图从本质上探究边坡的变形破坏机制，不涉及复杂且尚未得到很好解决的节理和岩体本构模型，擅于模拟边坡的非连续、大变形破坏过程，能够更好地揭示节理岩体和岩质边坡微细观变形破坏机制^[4-5]。

国内外学者通过离散单元法模拟了岩质边坡的破坏过程，并探讨了岩质边坡的变形失稳破坏机制。Ivars等^[6]建立了节理岩体的离散单元模拟方法，为构建节理岩质边坡数值模型奠定了良好的基础。周喻等^[7]以赞比亚穆利亚希北露天铜矿边坡为工程背景，通过离散单元法分析了含2条顺倾断续节理的岩质边坡变形破坏细观力学机制。蒋明镜等^{[4, 8-10]}将基于室内试验得到的胶结接触模型植入了离散元商业软件，探讨了含2组节理的岩质边坡变形失稳破坏机理；Bhaduri等^[11]采用离散单元法分析了静动荷载作用下含不同倾角节理的岩质边坡破坏模式；江魏等^[12]立足于提升离散元强度折减法的计算效率，提出了基于BP神经网络的颗粒细观参数标定策略；黄达等^[13]以离心模型试验为原型，通过离散元数值模拟研究了反倾层岩质边坡的变形机理与稳定性影响因素；Scholtes等^[14]通过将改进节理面模型导入离散元商业软件，揭示了含多组节理面的岩质边坡阶梯式破坏机制；Zhou等^[15]利用含断续节理面的岩质边坡离散元数值计算探讨了岩质边坡的能量破坏准则；范昊天等^[16]采用离散单元法PFC2D软件模拟了苏州清明山滑坡，分析了滑坡体滑动演化过程与破坏机理。

国内外学者构建了节理岩质边坡离散元数值计算方法，集中研究了含一组或多组不同预设节理面岩质边坡的变形破坏过程及其机制。而工程实践中，由多组密集节理面切割形成的层状碎裂结构边坡分布尤其普遍，其危害性更为严重，迫切需要研究层状碎裂结构边坡灾变破坏过程特征和失稳破坏模式，进而提升层状碎裂结构边坡灾变过程与失稳破坏机制水平。然而，利用离散单元法研究层状碎裂结构边坡失稳破坏过程方面的研究却鲜见报道。

本研究通过假定粒间微观胶结参数服从Weibull分布，构建了非均匀性岩体的离散元数值计算模型，分别通过6组顺倾和反倾节理面切割均质边坡建立了顺倾层和反倾层碎裂结构边坡模型。探讨了不同均质度条件下单轴压缩试验结果，从胶结破坏、速度、动能等方面揭示了层状碎裂结构边坡的灾变破坏过程，并将破坏模式与室内试验结果进行了对比分析。

1 离散单元法数值计算模型

大量研究结果表明^{[4, 14]}，岩质边坡的变形破坏涉及到节理面滑动张开与岩体剪切破坏。因此，离散单元法层状碎裂结构岩质边坡的数值计算不仅要建立节理面数值计算模型，而且也要考虑非均匀性岩体模型。

1.1 非均匀性岩体模型

2004年，Potyondy等^[17]通过理论推导建立了岩石胶结模型。该模型通过球颗粒及其之间平行胶结模拟由矿物晶粒几何体组成的岩石，可以较好地模拟岩石的弹性、断裂、破坏、强度、膨胀等物理力学特性，已被广泛运用于岩石介质的数值计算分析中。岩石胶结模型的力学响应和强度破坏准则见文献[17]。

岩石作为一种典型天然介质，由矿物晶粒、晶粒间接触和成岩过程造就的微细观缺陷构成，微细观缺陷主要包括各种形状和尺寸的孔洞、裂纹等。缺陷是岩石裂隙扩展贯通并发生破坏的诱发因素^[18]，导致岩石力学性质在空间上表现出非均匀性。然而，由岩石胶结模型构建的离散单元法岩石试样并未考虑决定岩石力学性能和破坏特性的主控因素-非均匀性^[19]，进而造成岩石胶结试样的应力应变曲线表现出强脆性特征，即峰前应力线性增大，峰后曲线急剧跌落，承载力瞬间下降。因而，基于岩石胶结模型的离散单元法岩石试样的应力应变曲线不符合泥岩、砂岩等易滑层状岩体的力学性质。

研究表明岩体脆性与其微细观结构的非均匀性相关^[20]，而岩体的非均匀缺陷服从于Weibull概率分布^[18]。因此，本研究假定离散单元法岩石试样微观参数服从Weibull分布，通过蒙特卡洛方法赋予颗粒之间胶结接触微观参数。平行胶结模型中的变形参数和强度参数均按式(1)所示的分布函数进行随机赋值。

(1)

式中，φ为概率密度；σ表征胶结的变形和强度参数；σ₀表征胶结的变形和强度参数均值；m为岩石试样的均质度，不同均质度条件下Weibull分布如图 1所示。由图可见，随着均质度的增加，胶结变形参数和强度微观参数的变化范围越小，试样的整体性质越均匀。

1.2 节理面模型

通常情况下，离散单元法采用降低微观强度和变形参数，删除颗粒等方式模拟节理面。这种做法会造成节理面颗粒影响起伏度，忽略节理面摩擦性，无法传递正应力等缺陷。为克服此问题，Pierce等^[21]提出了光滑节理接触模型，该模型可以更好地模拟节理面的力学行为并已得到广泛应用^[22]。

光滑节理接触模型如图 2所示。通过赋予位于光滑节理接触模型节理面两侧颗粒间节理面微观参数模拟结构面的力学行为。位于光滑节理接触模型节理面两侧的颗粒可以沿着其发生相对运动，并发生重叠。光滑节理接触模型节理面接触力包括法向力和切向力，分别通过粒间相对法向位移和相对切向位移进行累加更新并乘以相应刚度计算得到，见式(2)和式(3)。

(2)

(3)

式中，F_n为法向力；F_s为切向力；F_n0为初始法向力；F_s0为初始切向力；k_n为法向刚度；k_s为切向刚度；A为颗粒间节理面积；Δδ_n为颗粒间相对法向位移；Δδ_s为颗粒间相对切向位移。

非胶结光滑节理接触模型力学响应如图 3所示。由图可知，节理面法向仅能承受压力，压力随着法向压缩位移的增大而线性增加；切向力随着切向位移增加首先呈线性增大，当达到切向强度(法向压力与节理面摩擦系数乘积)后保持不变。

1.3 边坡模型微观参数

边坡模型微观参数包括颗粒微观参数、胶结微观参数、节理面微观参数3个部分(见表 1)。颗粒和胶结的弹性模量均值保持相同，法相刚度与切向刚度的比值为2.5，颗粒的摩擦系数选取0.5^[16]。通过模拟多组单轴压缩试验，并将应力应变曲线与室内试验相对比分析，选取切向强度均值为拉伸强度均值的3倍。节理面的法向刚度和切向刚度与胶结保持相同。

2 层状碎裂结构边坡模型 2.1 均质边坡模型

离散元模型生成方法包括定点成样法、等向压密法、粒径放大法、分层欠压法。其中定点成样法主要用于离散单元法验证研究，等向压密法适用于生成密样，分层欠压法主要用于生成离散单元法三轴试样。为更好地控制边坡应力环境，本研究采用粒径放大法生成均匀边坡的数值计算模型，步骤如下。

(1) 采用膨胀法生成一个立方体试样，由69 640个颗粒组成，颗粒最小半径为0.60 mm，最大半径为1.050 mm，平均半径为0.825 mm，孔隙率为0.35。

(2) 固定上、下端和侧向约束，采用离散元伺服系统，均匀地增大或减小所有颗粒的半径，直至达到700 kPa的竖向应力。

(3) 为增加密实度，同时剔除接触数目较少颗粒的同时，增加颗粒半径，经过不断的调整，最终所达到的平均接触数目为4.81个，满足岩石试样密实度的要求。

(4) 离散单元法岩质边坡试样及其力链分布如图 4所示。按图 4(a)所示要求进行削坡，边坡倾角为α，边坡长度为120 mm，高度为70 mm，宽度为30 mm。通过图 4(b)所示力链分布可知，所生成的离散单元法边坡试样粒间力分布较为均匀。

2.2 层状碎裂结构边坡模型

原位勘测结果表明，层状碎裂结构边坡包含密集节理面，节理面方位很多，其稳定性主要取决于岩块间的镶嵌情况和咬合力，不再沿着某个或某组节理面发生滑动破坏。为生成碎裂结构边坡，按照《工程岩体分级标准》(GB/T 0218—2014)要求，利用多组节理面切割均质边坡模型构建层状碎裂结构岩体，进而形成层状碎裂结构边坡。根据节理面与边坡倾角之间的关系，层状碎裂结构边坡可分为顺倾层边坡和反倾层边坡2种类型^[23]。本研究所生成的节理面主要用于切割均质岩体，进而形成碎裂状结构岩体边坡。

每个节理面组由2个全连通无限长的节理面组成。为了均匀切割均质边坡体，切割边坡岩体的节理面组中心位置在均质边坡模型中均匀设置，如图 5所示。每个中心点生成6个节理面组，节理面组倾向与边坡倾向相同；6个顺倾层边坡的节理面组倾角分别为15°，30°，45°；6个反倾层边坡节理组倾角分别为-15°，-30°，-45°。

为探索节理面摩擦性质、节理面间距、边坡倾角对顺倾层和反倾层边坡灾变过程的影响规律，选取节理面摩擦系数为0.1，0.3，0.5，0.7，0.9，节理面间距为2.5，7.5，12.5，17.5 mm，边坡倾角为30°，45°，60°，75°，开展10余组离散单元法数值计算，方案如表 2所示。

3 离散单元法数值计算结果

鉴于层状碎裂结构边坡变形破坏涉及到非均匀性岩体破坏和节理面滑动，为解决岩石胶结模型模拟岩石时遇到的峰后应力陡降缺陷，首先通过系列均质度条件下离散单元法单轴压缩试验探讨岩体脆性力学行为特征，然后从胶结破坏数、速度、动能等方面分析了顺倾层和反倾层边坡灾变过程特征。

3.1 离散单元法岩体单轴压缩试验

不同均质度条件下由离散单元法单轴压缩试验得到的应力和胶结破坏数与应变之间的关系如图 6所示，其中胶结破坏数可用于表征岩体裂纹扩展。由图可知，随轴向应变增加，轴向应力首先线性增加，达到峰值强度后逐渐减小，峰后应力呈现渐进性跌落，并未发生瞬间跌落，与非坚硬岩石的物理力学行为相一致。随应变增加，胶结破坏数首先保持不变，直至0.15%，此阶段没有新裂纹萌生，仅表现为原生缺陷的闭合，之后胶结破坏数快速增大，最后稳定于2.5万左右；随着均质度增加，单轴压缩强度先增加，当均质度m大于7之后轴向压缩强度基本保持不变。

脆性指数被广泛用于表征岩石的力学行为特征，国内外学者先后从强度、应变、能量、硬度等角度定义脆性指标，并进行了工程应用。本研究利用基于岩石峰前起裂及峰后特征的脆性评价方法^[19]分析不同均质度条件下离散单元法岩样的脆性特征。应力应变曲线峰前、峰后、综合脆性指数与均质度之间的关系如图 7所示。当均质度小于7时，峰前和综合脆性指标随均质度增加减小，而峰后脆性指数近乎保持不变；当均质度大于7时，峰前、峰后、综合脆性指标基本保持不变。

由图 6和图 7可知，运用本研究提出的基于Weibull分布的胶结接触微观参数赋值方法可以解决峰后应力应变曲线陡降的离散单元法模拟缺陷，通过改变均质度可模拟出不同峰前和综合脆性特征。本研究所提方法可以更好地模拟坚硬岩石和软弱岩石的力学响应行为。

3.2 边坡变形破坏过程

考虑到大尺寸岩质边坡数值模型计算难度，结合本研究模拟灾变过程需求，加之通常采用的强度折减法和重力增加法数值计算结果具有相似性^[24]，因而参照离心机模型试验原理，通过重力增加法模拟边坡灾变过程。为得到边坡模型滑动破坏时所对应的重力加速度，分别将赋予表 2所示参数的顺倾层边坡和反倾层边坡置于不同的重力环境中进行反复计算分析，得到顺倾层边坡和反倾层边坡滑动破坏时分别对应的3 500倍和4 000倍重力加速度。

为使离散单元法边坡岩体的运动与实际情况保持一致，需对离散单元法边坡数值模型赋予阻尼。离散元软件中存在局部阻尼和接触点阻尼。局部阻尼主要通过附加给颗粒一个阻尼力来实现阻尼的功能，该类阻尼会减小颗粒的运动速度，增加计算时间，但是对胶结强度没有任何影响。通过与室内试验结果^[25]对比，选取0.03作为岩体运动时的局部阻尼系数。接触点阻尼用于模拟岩块碰撞后动能转化为热能的消耗程度，这是一个不可逆转的过程，参照相关研究成果^[26-27]选取法向和切向阻尼比均为0.25。

3.2.1 胶结破坏数

顺倾层边坡胶结破坏数与时步之间的关系如图 8所示，其中胶结破坏可表征裂隙扩展。由图可知，施加3 500倍重力加速度后，胶结破坏数迅速增大，达到峰值后保持不变；胶结破坏数峰值随节理面摩擦系数减小，随边坡倾角增加而增大；除节理面间距为17.5 mm之外，胶结破坏峰值数随节理面间距增加而减小，其中17.5 mm的节理面间距导致边坡模型主要沿着节理面滑动。这说明顺倾层边坡节理面越光滑，边坡越陡，节理面间距越小，边坡变形失稳范围越大，裂隙延伸扩展越多。

反倾层边坡胶结破坏数与时步之间的关系如图 9所示。由图可知，除了节理面摩擦系数较大边坡之外，施加4 000倍重力加速度后，胶结破坏数迅速增大，而后近似保持不变；节理面摩擦系数≥0.3时，胶结破坏数增速缓慢，且峰值较小；相对顺倾层边坡，节理面摩擦系数对胶结破坏数影响更大；胶结破坏数峰值随节理面摩擦系数减小，节理面间距减小，边坡倾角增大而增大。这说明反倾层边坡节理面越光滑，节理面间距越小，边坡越陡，边坡变形破坏范围越大, 裂隙延伸扩展越强烈。

3.2.2 速度 3.2.2.1 速度场

速度是评估滑坡灾害危害性的重要指标，可用于揭示层状碎裂结构边坡灾变过程特征。本研究将颗粒速度量值按照不同颜色划分为7等分，如图 10所示。速度场中箭头的方向代表速度的方向，箭头的长短表示速度的大小。本研究选取节理面摩擦系数为0.1、节理面间距为2.5 mm、边坡倾角为60°的顺倾层和反倾层边坡模型的速度场变化过程阐述其灾变过程特征。

不同时步下顺倾层边坡速度场如图 11所示(图中t为时步)。

(1) 重力环境施加初期边坡速度场呈杂乱无序状态，随时步增加速度场逐渐集中于中、上部岩体，中、上部岩体底部速度增大；中、上部岩体发生变形并在其底部形成剪切带，剪切带从坡顶向坡脚扩展，滑坡体下错；滑动面基本贯通，边坡滑动处于蠕动-挤压阶段。

(2) 当时步由5×10⁴步增加至10×10⁴步时，下部岩体抗滑段被完全剪断，剪切带贯通，圆弧滑动面形成，滑动面处颗粒速度逐渐减小，越靠近坡面滑坡体的速度越大，至此边坡变形处于滑动-剧滑阶段。

(3) 随着时步的继续增加，滑体内部速度差别逐渐缩小直至稳定，但仍有少量颗粒滚落，此时边坡变形处于稳定压密阶段。

不同时步下反倾层边坡速度场如图 12所示。具体描述如下。

(1) 重力环境施加之初，边坡内颗粒首先呈杂乱无序的状态，中部临空面岩体速度随时步增加而增大，中、上部发生开裂变形；当时步由2.5×10⁴步增至4.5×10⁴步时，中、下部岩体速度增大，张拉裂缝向坡脚扩展直至贯通，发生翻折倾倒，靠近坡面岩体速度较大。

(2) 随着坡面中部临空面岩体反倾滑落，坡顶岩体失去支撑，后缘岩体开裂，坡顶岩体的速度增加，初始速度方向近乎垂直坡面，岩体发生翻折倾倒。

(3) 当时步达到7.0×10⁴步时，坡顶岩体向下翻滚滑动，随时步继续增加，滚落下来的岩土体沿坡面滑动，此时越靠近表面，岩体的速度越大。最终坡体表面处的岩土体达到稳定，但仍有少量的颗粒滚落。

综上所述，顺倾层边坡的灾变过程主要表现为：中、上部岩体变形—坡顶至坡脚剪切带扩展—抗滑段剪断—圆弧滑动。反倾层边坡的变形失稳具有从中部临空岩体至坡顶岩体的阶段性，中部临空面岩体主要表现为：变形—中部开裂—裂隙向坡脚扩展—翻折倾倒。坡顶岩体主要表现为：后缘开裂—翻折倾倒—坡面滑动。本研究模拟得到的反倾层边坡灾变过程特征与模型试验结果基本一致^[28]。

3.2.2.2 最大速度

不同节理面摩擦系数、节理面间距、边坡倾角条件下顺倾层边坡最大速度与时步之间的关系如图 13所示。由图可知，随着时步增加, 最大速度总体上先迅速增加而后呈波动式减小，节理面摩擦系数和间距对最大速度峰值的影响较小，但对峰后曲线变化规律影响较大。边坡倾角对最大速度、峰前和峰后曲线变化规律影响均较为显著。

不同节理面摩擦系数、节理面间距、边坡倾角条件下反倾层边坡最大速度与时步之间的关系如图 14所示。由图可知，最大速度随时步增加先迅速增大而后呈波动式减小，边坡倾角大于30°时才发生失稳破坏，这与黄润秋等^[29]研究结果相一致。节理面摩擦系数和节理面间距对速度峰值近乎无影响，而边坡倾角对速度峰值影响较为显著。

对比图 13和图 14可知，由于反倾层边坡呈倾倒破坏，所以峰后最大速度减小较快；不论是顺倾层边坡还是反倾层边坡，节理面间距和摩擦系数对速度峰值的影响较小，边坡倾角对速度峰值影响较为显著；顺倾层边坡节理面摩擦系数和间距对峰后曲线变化规律影响相对反倾层边坡更加显著。

3.2.3 动能

动能的变化可以用来评估灾害的危害程度。考虑颗粒平动速度和转动的动能计算为：

(4)

式中，M为颗粒的数目；m_i为每个颗粒的质量；v_i为颗粒的平动速度；I_i为每个颗粒的转动惯性矩；ω_i为颗粒的转动速度。

节理面摩擦系数为0.1、节理面间距为2.5 mm、边坡倾角为60°的顺倾层和反倾层边坡动能与时步的关系曲线如图 15所示。随着时步的增加，动能首先迅速增加，而后又迅速减小，但反倾层边坡动能增长较快，且所达到的峰值较大。这说明反倾层边坡冲击力更大，危害更强。

3.3 失稳破坏模式

重力荷载作用下顺倾层边坡沿圆弧滑动面发生滑动破坏^[30]，破坏模式为顺倾层剪切-圆弧滑动，大部分岩土体堆积在坡脚处。图 16为离散单元法边坡破坏失稳模式。

图 16(a)为模型试验得到的层状边坡模型失稳破坏模式，其为折线型滑动破坏，但滑坡体未完全散落至坡脚。由此表明，离散单元法数值计算与模型试验得到层状碎裂边坡失稳破坏模式近乎相同，但是由于离散单元法数值模型和试验模型节理面发育程度、岩体破碎程度等不同导致滑动面形状和滑坡体滑动范围不同。

图 16(b)所示的反倾层边坡失稳模式主要表现为中、下部岩体①拉裂, 翻折倾倒; 上部岩体②翻折倾倒, 坡面滑动，破坏面为折线型；而模型试验得到的失稳破坏模式为后缘拉裂倾倒, 前缘压裂倾倒，破坏面也为折线型。由此表明，离散单元法数值计算与模型试验得到的破坏面基本一致，但是坡体的破坏顺序相反，主要是由于模型试验采用基底摩擦试验装置开展，其受力作用与离散单元模型的受力作用不同造成的。

离散单元法数值计算结果与模型试验结果虽然在滑动面形状、破坏顺序等方面有所差异，但总体上二者失稳破坏模式具有一致性，由此表明了离散单元法边坡数值模拟结果的合理性与可靠性。

4 结论与展望

层状碎裂结构边坡灾害具有隐蔽性强、危害性大、防控难度大等特点，其稳定性评价和灾害防治均需探明边坡灾变过程特征。为利用离散单元法探讨层状碎裂结构边坡灾变过程，提出了可表征非均匀性的岩体模型，建立了顺倾层和反倾层碎裂结构边坡模型，开展了不同均质度条件下离散单元法单轴压缩试验和不同节理面摩擦系数、节理面间距、边坡倾角条件下10余组边坡数值模拟计算。从胶结破坏数、速度、动能等方面分析了层状碎裂结构边坡的灾变过程特征，并将破坏模式与室内模型试验结果进行了对比分析。

(1) 当均质度≤7时，随均质度增大，单轴压缩强度增大, 脆性指数减小，而均质度>7后应力应变曲线和脆性指数近乎保持不变；改变均质度可模拟不同峰前和综合脆性特征，克服了采用岩石胶结模型模拟岩石时遇到的峰后应力陡降缺陷。

(2) 节理面越光滑，边坡越陡，节理面间距越小，裂隙延伸扩展越多，破坏范围越大；相比节理面间距和摩擦系数，边坡倾角对速度峰值的影响更为显著；反倾层边坡动能增长更快，峰值更大。

(3) 顺倾层边坡表现出中、上部岩体变形—剪切带由坡顶至坡脚扩展—抗滑段剪断—圆弧滑动灾变过程；而反倾层边坡表现出中部临空面岩体变形—中部开裂—裂隙向坡脚扩展—翻折倾倒、坡顶岩体后缘开裂—翻折倾倒—坡面滑动灾变过程。

(4) 本研究数值计算结果与多个室内模型试验结果相吻合，验证了离散单元法数值计算模型和计算方法的正确性和合理性。

本研究围绕层状碎裂结构边坡稳定性评价和灾害防治实际需求，提出了一种非均匀性岩体和层状碎裂结构边坡灾变过程的离散单元法数值模拟新思路，所建立的数值模型可很好地模拟非均匀性岩体和层状碎裂结构边坡的力学行为和灾变过程。后续将会利用该数值模型探讨层状碎裂结构边坡的宏微观失稳破坏机制。