索风营水电站Dr2岩质边坡渐进破坏全过程数值模拟

Numerical Simulation of the Whole Process of Progressive Damage of Dr2 Rocky Slope in the Suofengying Hydropower Station

摘要

为解决滑坡运动过程中连续-非连续过程的模拟难点，将有限单元法（FEM）和可变形圆化多边形离散单元法（DSDEM）相耦合，提出了模拟滑坡运动过程的FEM-DSDEM数值模拟方法。该方法采用有限单元法模拟基岩，使用可变形圆化多边形离散单元法模拟滑坡体，引入虚拟裂缝模型和Mohr-Coulomb准则模拟滑坡体的破碎与剪切破坏等从连续到非连续状态的转化，在截断边界上施加黏弹性人工边界，实现地震作用下滑坡体连续-非连续渐进演化过程的数值模拟。利用巴西圆盘试验和四点弯曲梁冲击试验，验证FEM-DSDEM方法模拟材料连续-非连续演化过程的可行性和准确性。开展了索风营水电站Dr2岩质边坡渐进破坏全过程数值模拟，研究了座滑破坏、局部崩塌破坏和滑移破坏这3种不同破坏模式下的滑坡体的动力特性和堆积体空间分布特征。结果表明：Dr2岩质边坡失稳时，滑动模式均沿滑裂面产生滑动破坏，局部发生崩塌破坏；滑坡失稳运动的时间短、速度快、对河床对岸的山体冲击力强；失稳破坏后形成的堆积物会冲击河道，破坏和淤堵进水口。研究成果可为同类岩质边坡渐进破坏全过程的研究提供数值模拟方案。

Abstract

Objective The destabilization process of high and steep rock slopes and its influence range are crucial for the safe operation of engineering projects. The FEM-DSDEM numerical simulation method is proposed to simulate the landslide motion process by coupling the finite element me-thod (FEM) with the deformable spheropolygon-based discrete element method (DSDEM) to achieve the numerical simulation of the continuou-discontinuous progressive evolution process of landslides under earthquake action. Methods The method uses the finite element method to simulate the bedrock and the deformable spheropolygon-based discrete element method to simulate the landslide body. The fictitious crack model (FCM) and the Mohr-Coulomb criterion were introduced to simulate the transformation from continuous to discontinuous states, including the crushing and shear damage of the landslide body. In addition, a viscous-spring boundary was imposed on the truncated boundary to realize the numerical simulation of the continuous-discontinuous gradual evolution process of the landslide body under earthquake effects. First, the feasibility and accuracy of the FEM-DSDEM numerical simulation method for reflecting the continuous-discontinuous evolution process of materials were verified using the Brazilian disc test and the four-point bending beam impact test. Second, the corresponding calculation program was developed to realize the numerical simulation of the continuous-discontinuous progressive damage evolution process of landslides. On this basis, the numerical simulation of the entire progressive failure process of the Dr2 rock slope at Suofengying Hydropower Station was conducted, and the dynamic characteristics of the landslide body and the spatial distribution characteristics of the accumulation body under three different failure modes, namely shear sliding failure, local collapse failure, and slip failure, were investigated. Results and Discussions In the Brazilian disc test, the peak load at the upper plate was 1 870 kN, and the tensile strength value obtained using the formula was 5.90 MPa, which was 5.36% lower than the input value. The result was within the allowable error range, verifying that the method could simulate the fracture process of brittle materials. In the four-point bending beam impact test, the numerical simulation results of the maximum deflection curve of the beam during the fracture process were compared to the experimental results, and the two sets of results were in good agreement. This finding verified that the method could simulate the fracture process of rock subjected to impact loading. The numerical simulation results of the entire progressive failure process of the Dr2 rock slope at Suofengying Hydropower Station showed that when the Dr2 rock slope became unstable, all sliding modes produced sliding failure along the sliding cleavage surface, accompanied by local collapse failure. The destabilizing movement of the landslide was short in duration, rapid in velocity, and generated a strong impact on the mountains located on the opposite bank of the riverbed. In addition, the accumulation of body formed after a destabilizing failure impacted the river channel, causing destruction and sediment blockage at the water intake. Among the three failure modes, the landslide occurring under the shear sliding failure mode blocked the river channel at the fastest rate. Owing to the insufficient strength of the shear surface, the sediments generated from the primary landslide could also induce secondary landslides within the shear surface. The average displacement of each gage point under the local collapse failure mode was the largest, and the peak velocity at each gage point was also the highest. In contrast, the average displacement at each measurement point was the smallest under the slip failure mode. Conclusions This study verified the capability of the FEM-DSDEM method to realize the numerical simulation of the continuous-discontinuous progressive evolution process of landslides under seismic action through the Brazilian disc test, four-point bending beam impact test, and numerical simulation of the entire progressive failure process of the Dr2 rock slope at the Suofengying Hydropower Station. In addition, the research results obtained for the Dr2 rock slope at the Suofengying Hydropower Station can provide numerical simulation solutions for investigating the entire progressive damage process of similar rock slopes. This study investigates the motion process and accumulation patterns of hazardous rock bodies after sliding down the slope under different failure modes. However, the interaction between the hazardous rock body and reservoir water after destabilization into the water was not considered. Therefore, the strong coupling effect and energy transfer between the landslide mass and the fluid can serve as important directions for future research. In addition, this study is based on a two-dimensional model, and some discrepancies can exist between the simulation results and actual conditions. Therefore, future studies can further explore the numerical simulation capability of this method under three-dimensional rock slope models.

Graphical abstract

关键词

离散单元法 / 有限单元法 / 连续-非连续 / 黏弹性人工边界 / 岩质边坡

Key words

discrete element method / finite element method / continuous-discontinuous / viscous-spring boundary / rock slope

本刊网刊

库岸滑坡对大坝、水库等周边水工建筑物的运行安全有重要影响^[1‒3]，不仅会对自然生态环境产生严重的破坏，还对人民的生命财产安全产生了极大威胁^[4‒6]。

与连续介质方法相比，非连续介质法能真实反映块体碰撞与破坏等非连续本质，被广泛应用于滑坡体运动的研究，其中，一种较常用的方法是离散单元法^[7‒9]。离散单元^[10]根据单元几何特征可以分为圆盘（球）颗粒离散元^[11]和块体离散元两类。圆盘（球）颗粒离散元模型计算效率高^[12]，但外形表征性差^[13]。块体离散单元则可以更真实地描述块体形状，合理反映块体间的相互作用，但存在接触力计算效率低的问题^[14]。

圆化多边形离散单元法^[15‒17]将块体尖锐的角点处进行圆化处理，避免角点处法向奇异等问题，简化了接触判断，大大提高了接触检测的计算效率^[18‒19]，同时，具有很好的几何外形表征能力。然而，原始的圆化多边形离散单元基于刚体假设，不能描述单个离散单元的变形。

可变形圆化多边形离散单元法（DSDEM）^[20]将有限单元法与圆化多边形离散单元法相耦合，突破了圆化多边形单元刚体限制，可对任意形状单元的运动变形进行刻画。但是，在研究滑坡运动过程中，缺乏合适的模型模拟滑坡体连续-非连续渐进演化过程，无法准确反映滑坡体力学行为。

为了解决上述问题，本文提出FEM-DSDEM数值模拟方法，并开发了相应的计算程序。采用可变形圆化多边形离散单元法来分析岩石材料的变形与断裂演化特征，使用虚拟裂缝模型（FCM）^[21‒22]来确定拉伸破坏状态，基于拉伸截断的Mohr-Coulomb准则确定剪切破坏状态，同时，引入黏弹性人工边界模拟地震波输入^[23‒24]，用于模拟滑坡体连续-非连续渐进破坏演化过程。以索风营水电站Dr2岩质边坡滑坡为例，将本文方法应用到滑坡动态演化全过程分析和滑坡动力特性、堆积体空间分布特征等的预测中，为预测岩质边坡运动过程提供了一种思路和方法。

1 FEM-DSDEM数值模拟方法

1.1 FEM-DSDEM耦合分析方法

圆化多边形SP是基本多边形P与圆盘S的闵可夫斯基和，即由圆盘S扫过基本多边形的外轮廓所构成，如图1所示。圆化半径为\mathrm{0.05}r，其中，r为最大内切圆半径。

在FEM-DSDEM方法中，通常将计算域划分为上部边坡、下部基岩及交界部分3个区域，并分别记为A、B和C，FEM-DSDEM的耦合模型如图2所示。区域A为上部边坡，在滑动过程中经历了显著的变形和运动，因此，用DSDEM进行模拟；区域B为下部基岩，只会发生较小的变形，且不与滑坡体发生接触碰撞，采用FEM进行模拟；区域C同样是小变形区域，但会与区域A中的单元发生接触碰撞，因此，需要对区域A和C中的单元进行接触检测和接触力的计算。

利用FEM-DSDEM方法模拟地震作用下滑坡体连续-非连续渐进演化过程的步骤如下：

步骤1，建立几何模型，对几何模型进行分组，并用有限元网格对构件进行离散化。

步骤2，对于可能发生破碎的单元组，在相邻的有限单元之间嵌入0厚度的节理单元来模拟断裂过程^[25]，如图3所示。

步骤3，设置黏弹性人工边界，考虑外部荷载，确定材料参数和约束条件。

步骤4，预测t+\mathrm{\Delta }t时刻单元在每个节点上的位移。

步骤5，进行接触检测，并先使用式（1）～（4）计算接触力，再由式（5）得到t+\mathrm{\Delta }t时刻的等效节点力矢量。

当两个圆化多边形S{P}_{\mathrm{1}}和S{P}_{\mathrm{2}}接触时，接触形式可分为点-边接触和点-点接触，其中，边-边接触可简化为两次点-边接触，具体如图4所示。

单元间的接触力可分为法向接触力{\mathit{P}}_{\mathrm{n}}^{}和t+\mathrm{\Delta }t时刻的切向接触力{\mathit{P}}_{\mathrm{s}}^{t+\mathrm{\Delta }t}，分别表示为：

{\mathit{P}}_{\mathrm{n}}^{}={\displaystyle \sum _{i=\mathrm{1},j=\mathrm{1}}^{i=N_{V\mathrm{1}},j=N_{V\mathrm{2}}}}{\mathit{P}}_n(V_{\mathrm{1}}^i,V_{\mathrm{2}}^j)+{\displaystyle \sum _{i=\mathrm{1},j=\mathrm{1}}^{i=N_{V\mathrm{1}},j=N_{E\mathrm{2}}}}{\mathit{P}}_n(V_{\mathrm{1}}^i,E_{\mathrm{2}}^j)+{\displaystyle \sum _{i=\mathrm{1},j=\mathrm{1}}^{i=N_{E\mathrm{1}},j=N_{V\mathrm{2}}}}{\mathit{P}}_n(E_{\mathrm{1}}^i,V_{\mathrm{2}}^j)（1）

{\mathit{P}}_{\mathrm{s}}^{t+\mathrm{\Delta }t}={\mathit{P}}_{\mathrm{s}}^t+K_{\mathrm{s}}\cdot \mathrm{\Delta }{\mathit{\delta }}_{\mathrm{s}}^{t+\mathrm{\Delta }t}（2）

式（1）、（2）中，N_{V\mathrm{1}}、N_{V\mathrm{2}}和N_{E\mathrm{1}}、N_{E\mathrm{2}}分别为顶点数和边数，V_{\mathrm{1}}^i和V_{\mathrm{2}}^j为不同的顶点，E_{\mathrm{1}}^i和E_{\mathrm{2}}^j为不同的边，i、j为求和编号，上标t、\mathrm{\Delta }t均表示时间，K_{\mathrm{s}}为切向接触刚度，\mathrm{\Delta }{\mathit{\delta }}_{\mathrm{s}}^{t+\mathrm{\Delta }t}为t+\mathrm{\Delta }t时刻的切向增量位移，{\mathit{P}}_{\mathrm{s}}^t为t时刻的切向接触力。

同时，考虑库仑摩擦准则有：

{\mathit{P}}_{\mathrm{s}}^{t+\mathrm{\Delta }t}={\mathit{n}}_{\mathrm{s}}\mathrm{m}\mathrm{i}\mathrm{n}(\left|{\mathit{P}}_{\mathrm{s}}^{t+\mathrm{\Delta }t}\right|,\mu \left|{\mathit{P}}_{\mathrm{n}}^{t+\mathrm{\Delta }t}\right|)（3）

式中，\mu为摩擦系数，{\mathit{n}}_{\mathrm{s}}为t+\mathrm{\Delta }t时刻切向的单位矢量，{\mathit{P}}_{\mathrm{n}}^{t+\mathrm{\Delta }t}为t+\mathrm{\Delta }t时刻的法向接触力。

作用在接触单元上的净接触力{\mathit{P}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{c}\mathrm{t}}为：

{\mathit{P}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{c}\mathrm{t}}={\mathit{P}}_{\mathrm{n}}+{\mathit{P}}_{\mathrm{s}}（4）

式中，{\mathit{P}}_{\mathrm{s}}为切向接触力。

根据虚功原理，将S{P}_{\mathrm{1}}上的接触力转化为节点k的等效节点力矢量{\mathit{f}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{c}\mathrm{t}}^k：

{\mathit{f}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{c}\mathrm{t}}^k={\mathit{N}}^k{\mathit{P}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{c}\mathrm{t}}（5）

式中，{\mathit{N}}^k为单元S{P}_{\mathrm{1}}中接触点的形函数。

步骤6，使用式（6）和（7）计算节理单元应力，利用式（8）和（9）判断节理单元是否失效。如果节理单元不失效，则节理单元上的应力应转化为节点力是矢量{\mathit{f}}_{\mathrm{j}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}}^{}，可以用式（10）计算{\mathit{f}}_{\mathrm{j}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}}^{}。如果节理单元失效，则节点力为0。

正应力{\sigma }_{\mathrm{n}}和剪应力{\sigma }_{\mathrm{s}}分别表示为：

（6）

（7）

式（6）、（7）中，k_{\mathrm{n}}和k_{\mathrm{s}}分别为节理单元的法向刚度和切向刚度，f_{\mathrm{t}}和f_{\mathrm{s}}分别为抗拉强度和抗剪强度，{\delta }_{\mathrm{n}}和{\delta }_{\mathrm{s}}分别为节理单元的法向相对位移和切向相对位移，w_{f\mathrm{t}}和w_{f\mathrm{s}}分别为受力达到抗拉强度和抗剪强度时对应的法向和剪切位移，w_{f\mathrm{t}}+w_{Gf\mathrm{Ⅰ}}和w_{f\mathrm{s}}+w_{Gf\mathrm{Ⅱ}}分别为两个有限单元间完全破坏时最大的法向和切向位移，\phi为材料内摩擦角。

用虚拟裂缝模型和拉伸截断的Mohr-Coulomb准则来判断节理单元是否破坏。抗剪强度f_{\mathrm{s}}为：

（8）

式中，c为材料黏聚力。

法向位移和切向位移之间的关系若符合式（9）的条件，节理单元就会发生断裂：

{\left(\frac{{\delta }_{\mathrm{n}}-w_{f\mathrm{t}}}{w_{Gf\mathrm{Ⅰ}}}\right)}^{\mathrm{2}}+{\left(\frac{{\delta }_{\mathrm{s}}-w_{f\mathrm{s}}}{w_{Gf\mathrm{Ⅱ}}}\right)}^{\mathrm{2}}\ge \mathrm{1}（9）

节点力矢量{\mathit{f}}_{\mathrm{j}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}}^{}的计算式为：

{\mathit{f}}_{\mathrm{j}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}}^{}=\int N^{\mathrm{T}}\sigma \mathrm{d}\Omega（10）

式中，\mathit{N}为节理单元的形函数，\mathit{\sigma }为节理单元内的应力向量，\Omega为节理单元的积分区域。

步骤7，采用式（11）～（14）计算每个有限元在t+\mathrm{\Delta }t时刻的应力{\mathit{\tau }}_t+\mathrm{\Delta }\mathit{S}、应变{\mathit{\epsilon }}^{t+\mathrm{\Delta }t}，以及对应的内力矢量{\mathit{f}}_{\mathrm{i}\mathrm{n}}^{t+\mathrm{\Delta }t}。

节点内力矢量{\mathit{f}}_{\mathrm{i}\mathrm{n}}^{}^[26‒27]和t+\mathrm{\Delta }t时刻应变{\mathit{\epsilon }}^{t+\mathrm{\Delta }t}可分别表示为：

{\mathit{f}}_{\mathrm{i}\mathrm{n}}^{}=\int {\mathit{B}}_t^{\mathrm{T}}({\mathit{\tau }}_t+\mathrm{\Delta }\mathit{S})\mathrm{d}\Omega（11）

{\mathit{\epsilon }}^{t+\mathrm{\Delta }t}={\mathit{\epsilon }}^t+\mathrm{\Delta }\mathit{\epsilon }（12）

式（11）、（12）中：{\mathit{B}}_t为应变-位移矩阵；{\mathit{\tau }}_t为t时刻的柯西应力向量；{\mathit{\epsilon }}^t为在t时刻的应变；\mathrm{\Delta }\mathit{S}为第二类Piola-Kirchhoff应力增量，如式（13）所示；\mathrm{\Delta }\mathit{\epsilon }为t时刻到t+\mathrm{\Delta }t时刻的应变增量，如式（14）所示。

\mathrm{\Delta }\mathit{S}=\mathit{D}\mathrm{\Delta }\mathit{\epsilon }（13）

\mathrm{\Delta }\mathit{\epsilon }={\mathit{B}}_t{\mathit{u}}_{t+\mathrm{\Delta }t}（14）

式（13）、（14）中，\mathit{D}和{\mathit{u}}_{t+\mathrm{\Delta }t}分别为t时刻到t+\mathrm{\Delta }t时刻的材料本构增量矩阵和单元节点位移增量。

步骤8，采用Velocity Verlet算法求解运动控制方程式（15），计算t+\mathrm{\Delta }t时刻的加速度{\mathit{a}}_{}^{t+\mathrm{\Delta }t}、速度{\mathit{v}}_{}^{t+\mathrm{\Delta }t}和位移{\mathit{u}}_{}^{t+\mathrm{\Delta }t}。

\mathit{M}\mathit{a}={\mathit{f}}_{\mathrm{e}\mathrm{x}}^{}-{\mathit{f}}_{\mathrm{i}\mathrm{n}}^{}-\mathit{C}\mathit{v}（15）

式中：\mathit{M}和\mathit{C}分别为质量矩阵和瑞利阻尼矩阵；\mathit{v}和\mathit{a}分别为速度和加速度；{\mathit{f}}_{\mathrm{e}\mathrm{x}}^{}为总外力矢量，包括接触力的等效节点力矢量{\mathit{f}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{c}\mathrm{t}}^{}、外部总载荷矢量{\mathit{f}}_{\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}}^{}和节理单元的节点力矢量{\mathit{f}}_{\mathrm{j}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}}^{}。

1.2 黏弹性人工边界

为了准确模拟地震波的输入，引入最早由Deeks等^[28]提出的黏弹性人工边界。

阻尼器和弹簧的力学参数由基岩材料决定^[29]，人工边界节点上的等效刚度系数K_{\mathrm{b}}^{}和阻尼系数C_{\mathrm{b}}^{}的计算式分别为：

K_{\mathrm{b}}=\frac{G}{\mathrm{2}{r}_{\mathrm{b}}}（16）

C_{\mathrm{b}}=\rho C_{\mathrm{s}}（17）

式（16）、（17）中，\rho和G分别为介质密度和剪切模量，C_{\mathrm{s}}为剪切波速，r_{\mathrm{b}}^{}为波源到人工边界点的距离。

1.3 数值模拟方法验证

1.3.1 巴西圆盘试验

巴西圆盘试验模型如图5所示。网格尺寸为3.7～8.7 mm，圆盘半径为0.1 m，在板的顶部施加恒定的垂直速度v_y=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{2}}\mathrm{m}/\mathrm{s}。力学参数如下：杨氏模量E=\mathrm{21.00}\mathrm{G}\mathrm{P}\mathrm{a}，泊松比\upsilon =\mathrm{0.2}，抗拉强度f_{\mathrm{t}}=\mathrm{2.50}\mathrm{M}\mathrm{P}\mathrm{a}，Ⅰ型断裂能G_fⅠ=\mathrm{120}\mathrm{N}/\mathrm{m}，黏聚力c=\mathrm{20.00}\mathrm{M}\mathrm{P}\mathrm{a}，内摩擦角\phi =45°，Ⅱ型断裂能G_fⅡ=\mathrm{1}\mathrm{500}\mathrm{N}/\mathrm{m}。

圆盘的裂缝发展过程及x方向的正应力{\sigma }_{xx}分布变化过程如图6所示。最初在圆盘与加载板接触点处出现应力集中，之后应力高值区向边缘扩散，见图6（a）、（b）。t=0.064 s时，圆盘中心出现裂缝，此时正应力达到抗拉强度，见图6（c）。t=0.080 s时，出现贯穿裂缝，圆盘被破坏，见图6（d）。

当施加696 kN荷载力时，圆盘中心的正应力约为2.20 MPa。除两个加载板附近应力集中处外，{\sigma }_{xx}和{\sigma }_{yy}（y方向的正应力）的本文数值解与文献[30]的解析解一致，如图7所示。上板处峰值荷载为1 870 kN，利用正应力公式可以求得圆盘抗拉强度为5.90 MPa，与给定的抗拉强度相差5.36%，说明了本文方法模拟连续-非连续演化过程的准确性。

1.3.2 四点弯曲梁冲击试验

为了验证本文方法的准确性，选择由Murray等^[31]开展的四点弯曲梁冲击试验进行模拟，模型如图8所示。力学参数如下：杨氏模量E=\mathrm{25.80}\mathrm{G}\mathrm{P}\mathrm{a}，泊松比\upsilon =\mathrm{0.15}，密度\rho =\mathrm{2}\mathrm{320}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}，抗拉强度f_{\mathrm{t}}=\mathrm{2.70}\mathrm{M}\mathrm{P}\mathrm{a}，Ⅰ型断裂能G_fⅠ=\mathrm{98}\mathrm{N}/\mathrm{m}。

模拟混凝土梁的破坏过程结果如图9所示。首先，裂缝出现在冲击材料的正下方梁下表面；然后，向上延伸，在梁端1/4处也相继出现裂缝，直至梁完全断裂；最后，梁被分解成5段，与文献[31]实验结果一致。

2 索风营水电站Dr2岩质边坡渐进破坏过程数值模拟

Dr2岩质边坡位于索风营水电站右坝肩上方、引水发电系统进水口处，是坝址区最大、最严重的危岩体。Dr2岩质边坡的全景图、边坡位置图及上游侧视图分别如图10（a）、（b）和（c）所示。

2.1 Dr2岩质边坡的形态与结构特征

Dr2岩质边坡的不同剖面形式相近，具有二维化特征，因此，本文基于图11所示的典型剖面进行二维计算分析。

边坡顶高程约为1 080 m，坡脚高程约为900 m。上部T₁m地层岩层倾角一般在12°～17°之间，下部T₁y³地层岩层倾角相对较陡，岩层倾角一般在15°～25°之间。T₁m地层倾角由上部地层和中部地层（T₁m¹和T₁m²）组成。海拔1 070 m以上有一个坡度为5°～10°的由白云岩组成（T₁m²）的缓坡台地，在海拔960～1 070 m之间有一个坡度近70°的由石灰岩（T₁m¹）组成的陡壁，海拔960 m以上的体积约为4.66×10⁵ m³。Dr2岩体沿河道被3条明显的卸载裂缝L₁、L₂和L₃切割，如图12（a）所示。其中，后缘连接L₁的拉伸断裂裂缝将岩石边坡与后方完整岩层分离，如图12（b）所示。

Dr2岩质边坡的基底为T₁y³，为灰绿色、紫红色泥岩夹泥灰岩及灰岩。基底下部与埋藏深度为20～50 m的Ⅲ号堆积体相连。边坡和基底分布在J₁、J₂、J₃、J₄ 4个软弱夹层中。

2.2 地震波的输入

索风营Dr2岩质边坡坝址区基本地震烈度为6度，设计地震动加速度峰值为0.05g，地震动加速度反应谱特征周期为0.35 s。通过对现有的地震记录进行人工处理，依据《工程场地地震安全性评价技术规范》（GB 17741—2005）建立人工地震波。地震波时间步长为0.02 s，持续时间为10.00 s。在人工边界节点上施加等效荷载，在地基的边界上采用黏弹性人工边界输入地震波。

2.3 计算模型及参数设置

对Dr2岩质边坡地质条件和边界条件进行分析，确定了3种可能发生的破坏模式：

1）座滑破坏。危岩体底部T₁y³地层内发育多条夹层，可能形成危岩体向外剪切滑出的外缘不利结构面，危岩体后缘拉裂面已基本贯通。T₁y³泥岩层强度相对较弱，上百米高的陡立岩柱直接作用在下部软岩基座，由于Dr2岩质边坡的自重和其他外载荷作用，泥岩层被压缩，可能形成间歇破碎剪切面，剪切面可能与前缘夹层和后缘L₁拉伸裂缝连接形成滑动面。

2）局部崩塌破坏。Dr2岩质边坡高度为180 m，宽度仅为29～37 m。当危岩体产生一定变形或因结构面参数降低时，危岩体顶部的局部块体有可能引发局部崩塌破坏。其岩质边坡的突出岩体被L₃沿河切割，形成薄板岩体。随岩体的进一步风化作用，裂缝进一步延伸，可能导致局部崩塌破坏。

3）滑移破坏。主要表现为沿夹层的顺层滑动。由于危岩体岩层产状倾向上游，在近直立结构面为侧向面切割上游分离形成的块体，存在沿层间夹层向坡外滑移的可能。其中，最大可能的滑动方向为顺岩层真倾角方向，偏移真倾角方向越大，其岩层视倾角越小，滑动的可能性越小。由于Dr2岩质边坡后缘L₁裂隙的高贯通率，在破坏过程中，分布区岩层往往沿J₁软弱夹层上游向外滑移，最终会使滑移破坏。

根据3种破坏模式，综合考虑危岩体所处地理位置后分析认为，危岩体的破坏应以座滑破坏模式为主，其次为崩塌、滑移破坏模式，其中，崩塌破坏又以局部的崩塌破坏为主。

利用FEM-DSDEM数值模拟方法对Dr2岩质边坡滑坡的运动过程进行数值模拟，可得到3种破坏模式下滑坡的动力特性和堆积体空间分布特征。根据工程情况，建立数值计算模型。模型水平方向跨度1 298 m，竖直方向最大跨度689 m，3种破坏模式的数值模型分别如图13、14和15所示，每种破坏模式均选择4个监测点。其中：座滑破坏数值模型的单元总数为8 951，节点总数为11 648；局部崩塌破坏数值模型的单元总数为4 594，节点总数为4 948；滑移破坏数值模型的单元总数为9 288，节点总数为12 161。

采用钻探、坑探、地球物理勘探等综合方法对边坡地质条件进行调查，通过测试确定材料性能参数。根据勘察报告得知，Dr2岩质边坡底部的摩擦系数为0.30，T₁m的抗拉强度为3.0 MPa。所有节理单元的法向刚度值和切向刚度值均为\mathrm{5.0}\times {\mathrm{10}}^{\mathrm{7}}\mathrm{N}/\mathrm{m}，其余计算参数^[32]如表1所示。

2.4 Dr2岩质边坡3种破坏模式下的模拟结果

2.4.1 座滑破坏

座滑破坏的模拟结果如图16所示。在重力作用下，滑坡体最初沿滑动面下滑；Ⅲ号堆积体在上部岩体的推动下沿滑动面加速时，会出现空中抛射现象，如图16（b）所示；由于另一侧山体的阻挡作用，滑坡会直接冲击河床，迅速堵塞河道，如图16（c）所示；山体滑坡产生的堆石在24.30 s时已经完全堵塞了河道，如图16（e）所示。

座滑破坏的滑坡的位移、速度模拟结果分别如图17（a）、（b）所示。

由图17（a）可见：各滑坡体的平均最大位移为160 m；滑坡前缘1号监测点的位移最大，达到了170 m，滑动时间在17 s左右；位于滑坡前缘底部的2号监测点位移较小，仅为60 m，这说明其堆积部位接近滑坡启动部位；位于滑坡上半部分的3、4号监测点的位移约为150 m。结合图16可以发现，3、4号监测点没有滑入山谷，而是聚集在Ⅲ号堆积体的剪切面上。由于剪切面强度不足，这些滑坡的沉积物可能会在剪切面内引起二次滑坡。

由图17（b）可见：滑坡的整个过程持续约35.00 s，在滑坡发生初期，由于受地震荷载作用，滑坡前缘监测点（1、2号监测点）的速度明显增大，这两点的峰值速度分别为22.0和17.5 m/s；位于滑坡后缘的3、4号监测点由于高程较高，峰值速度出现较晚；滑坡体在整个过程中受到重力和地震荷载的共同作用，经历了启动、加速、减速和堆积等阶段。

2.4.2 局部崩塌破坏

局部崩塌破坏的模拟结果如图18所示。最外层山体与Dr2岩石边坡分离，并沿整个L₃段滑动，滑坡在撞击Ⅲ号堆积体时破裂；一部分碎石堆积在Ⅲ号堆积体上，另一部分碎石流入乌江，最终沉积在河床底部，如图18（f）所示；岩体局部崩塌滑动时间很短，形成的岩块沿Ⅲ号堆积体顶部滑动和抛射；局部崩塌破坏的岩石量较小，滑坡并未堵塞电站进水口，但高速运动的岩体会严重威胁工程运行安全；1号监测点最终堆积在Ⅲ号堆积体上，其他3个监测点与Ⅲ号堆积体碰撞后流入山谷。

局部崩塌破坏滑坡的位移、速度模拟结果分别如图19（a）、（b）所示。

由图19（a）可见：各监测点的平均最大位移为210 m；2号监测点位于滑块体底部，其初始滑块速度较快，但持续运动时间很短；3号监测点位于滑坡中下部，位移较小，说明3号监测点没有落入山谷，而是最终堆积在Ⅲ号堆积体上。这些堆积体会对Ⅲ号堆积体产生一定的影响，可能导致Ⅲ号堆积体发生剪切滑动。

由图19（b）可见，除2号监测点外，1、3、4号监测点均有两个峰值，原因可能是滑坡在撞击Ⅲ号堆积体时发生破裂，2号监测点位堆积在滑坡启动部位，而1、3、4号监测点位继续发生滑动。

2.4.3 沿J₁软弱夹层滑移破坏

Dr2岩质边坡沿J₁软弱夹层滑移破坏的模拟结果如图20所示。在滑坡初期，Ⅲ号堆积体由于抗拉强度低，再加上重力和地震荷载的共同作用，最先发生滑移。随后，Dr2岩质边坡产生倾斜滑动现象，由于Dr2岩质边坡底部泥岩强度较低，在滑动过程中被压碎，部分破碎的泥岩碎屑在危险岩体滑动过程中被挤出，如图20（a）所示。接着，Dr2岩质边坡在滑动过程中与其他岩体碰撞发生破碎，如图20（b）～（d）所示。最终，破碎的岩石和泥岩一起滑入河道，并堵塞了河道。

地震荷载作用下，滑移破坏滑坡的位移、速度模拟结果分别如图21（a）、（b）所示。

由图21（a）可见：各滑动体的平均最大位移为140 m；前缘的滑动体相对较快开始发生滑移破坏，但运动的持续时间非常短，监测到的最大位移为220 m。

由图21（b）可见：整个滑坡持续了35.00 s；滑坡发生初期，滑坡前缘1、2号监测点速度明显增大，其峰值速度分别为17.0和22.0 m/s；位于滑坡后缘的3、4号监测点由于高程较高，速度峰值出现较晚。

3 结 论

本文提出了FEM‒DSDEM数值模拟方法。将有限单元法与可变形圆化多边形离散单元法相耦合，采用有限单元法模拟基岩，使用可变形圆化多边形离散单元法模拟滑坡体，并引入黏弹性人工边界来模拟地震波的输入，开发了相应的计算程序，实现了滑坡体连续-非连续渐进破坏演化过程数值模拟。

验证了本文方法模拟断裂全过程的可行性和准确性。通过巴西圆盘试验验证了该方法能够模拟脆性材料的断裂过程，通过四点弯曲梁冲击试验验证了该方法可以模拟岩石承受冲击载荷后断裂的过程。

本文以实际工程为例，研究了Dr2岩质边坡的连续-非连续渐进破坏演化过程及其动力特性、堆积体空间分布特征。3种破坏模式的模拟结果表明：Dr2岩质边坡失稳时，滑动模式均沿滑裂面产生滑动破坏，局部发生崩塌破坏；滑坡失稳运动的时间短，速度快，对河床对岸的山体冲击力强；失稳破坏后形成的堆积物冲击河道，破坏和淤堵进水口。3种破坏模式中，座滑破坏模式下发生的滑坡堵塞河道速度最快，由于剪切面强度不足，一次滑坡的沉积物可能在剪切面内引起二次滑坡。局部崩塌破坏模式下各测点的平均位移最大，各测点的峰值速度也最大。滑移破坏模式下各测点的平均位移最小。

本文研究了不同破坏模式时危岩体下滑后的运动过程和堆积形态，但未考虑危岩体失稳入水后与库水的相互作用，滑坡体与流体之间的强耦合和能量传递可以作为后续的研究方向。此外，本研究基于二维模型，模拟结果与实际结果可能存在一定的差距，未来可以探索该方法在三维岩质边坡模型下的数值模拟。

参考文献

[1] Zhou Jiawen, Chen Mingliang, Qu Jingkun,et al.Research and prospect on disaster-causing mechanism and pre-vention-control technology of reservoir landslides[J].Advanced Engineering Sciences,2023,55(1):110-128.

[2] 周家文,陈明亮,瞿靖昆,等.水库滑坡灾害致灾机理及防控技术研究与展望[J].工程科学与技术,2023,55(1):110-128.

[3] Li Changdong, Long Jingjing, Jiang Xihui,et al.Advance and prospect of formation mechanism for reservoir landslides[J].Bulletin of Geological Science and Technology,2020,39(1):67-77.

[4] 李长冬,龙晶晶,姜茜慧,等.水库滑坡成因机制研究进展与展望[J].地质科技通报,2020,39(1):67-77.

[5] Zhan Qinghua, Wang Shimei, Wang Li,et al.Analysis of failure models and deformation evolution process of geological hazards in Ganzhou city,China[J].Frontiers in Ear-th Science,2021,9:731447. doi:10.3389/feart.2021.731447

[6] Xia Guoqing, Liu Chun, Xu Chong,et al.Dynamic analysis of the high-speed and long-runout landslide movement process based on the discrete element method:A case study of the Shuicheng landslide in Guizhou,China[J].Advances in Civil Engineering,2021,2021:8854194. doi:10.1155/2021/8854194

[7] Liu Guangyu, Xu Wenjie, Tong Bin,et al.Study on dynamics of high-speed and long run-out landslide hazards based on block discrete element method[J].Chinese Journal of Rock Mechanics and Engineering,2019,38(8):1557-1566. doi:10.13722/j.cnki.jrme.2019.0158

[8] 刘广煜,徐文杰,佟彬,等.基于块体离散元的高速远程滑坡灾害动力学研究[J].岩石力学与工程学报,2019,38(8):1557-1566. doi:10.13722/j.cnki.jrme.2019.0158

[9] Jiang Baode, Li Xiuchun, Luo Haiyan,et al.A comparative analysis of heterogeneous ensemble learning methods for landslide susceptibility assessment[J].China Civil Engineering Journal,2023,56(10):170-179.

[10] 江宝得,李秀春,罗海燕,等.异质集成学习在滑坡易发性评价中的对比研究[J].土木工程学报,2023,56(10):170-179.

[11] Liu Dingzhu, Cui Yifei, Wang Hao,et al.Assessment of local outburst flood risk from successive landslides:Case study of Baige landslide-dammed lake,upper Jinsha river,eastern Tibet[J].Journal of Hydrology,2021,599:126294. doi:10.1016/j.jhydrol.2021.126294

[12] Shen Haohan, Zhang Hai, Fan Junkai,et al.Influence of contact radius on rock mechanical property and its application in discrete element method software EDEM[J].Rock and Soil Mechanics,2022,43(Supp1):580-590.

[13] 申浩翰,张海,范俊锴,等.离散单元法软件EDEM中接触半径对岩石力学特性的影响及其应用[J].岩土力学,2022,43():580-590.

[14] Wang Huanling, Sha Cong, Xu Weiya,et al.Research on str-ength of soil-rock mixture based on particle discrete element method[J].China Civil Engineering Journal,2020,53(9):106-114.

[15] 王环玲,沙聪,徐卫亚,等.基于颗粒离散元的土石混合体强度影响研究[J].土木工程学报,2020,53(9):106-114.

[16] Xiao Haobo, Qi Tianqi, Yang Shuhan,et al.Macro and mi-cro-behaviors of ellipsoidal particle system using 3D DEM simulation[J].Advanced Engineering Sciences,2023,55(6):78-86.

[17] 肖浩波,漆天奇,杨舒涵,等.椭球颗粒体系宏、细观特性的3维离散元分析[J].工程科学与技术,2023,55(6):78-86.

[18] Jiang Wei, Yan Jinzhou, Ouyang Ye,et al.Calibration of micro parameters of particles in granular discrete element method to assess slope stability by strength reduction method[J].Advanced Engineering Sciences,2023,55(5):50-60. doi:10.15961/j.jsuese.202200185

[19] 江巍,闫金洲,欧阳晔,等.边坡稳定性强度折减颗粒离散元法分析的细观参数标定策略[J].工程科学与技术,2023,55(5):50-60. doi:10.15961/j.jsuese.202200185

[20] Tan Pan, Rao Qiuhua, Li Zhuo,et al.A new method for quantitative determination of PFC3D microscopic parameters considering fracture toughness[J].Journal of Central South University(Science and Technology),2021,52(8):2849-2866. doi:10.11817/j.issn.1672-7207.2021.08.030

[21] 谭攀,饶秋华,李卓,等.考虑断裂韧度的PFC-3D细观参数标定新方法[J].中南大学学报(自然科学版),2021,52(8):2849-2866. doi:10.11817/j.issn.1672-7207.2021.08.030

[22] Galindo-Torres S A, Pedroso D M, Williams D J,et al.Breaking processes in three-dimensional bonded granular materials with general shapes[J].Computer Physics Communications,2012,183(2):266-277. doi:10.1016/j.cpc.2011.10.001

[23] Zhang Yongshuang, Guo Changbao, Lan Hengxing,et al.Reactivation mechanism of ancient giant landslides in the tectonically active zone:A case study in Southwest China[J].Environmental Earth Sciences,2015,74(2):1719-1729. doi:10.1007/s12665-015-4180-6

[24] Wei M D, Dai F, Xu N W,et al.Discussion on “a calibration methodology to obtain material parameters for the representation of fracture mechanics based on discrete element simulations”[J].Computers and Geotechnics,2017,86:246-248. doi:10.1016/j.compgeo.2017.03.001

[25] Alonso-Marroquín F, Wang Yucang.An efficient algorithm for granular dynamics simulations with complex-shaped objects[J].Granular Matter,2009,11(5):317-329. doi:10.1007/s10035-009-0139-1

[26] Alonso-Marroquín F.Spheropolygons:A new method to simulate conservative and dissipative interactions between 2D complex-shaped rigid bodies[J].Europhysics Letters,2008,83(1):14001. doi:10.1209/0295-5075/83/14001

[27] Liu Lu, Ji Shunying.Bond and fracture model in dilated polyhedral DEM and its application to simulate breakage of brittle materials[J].Granular Matter,2019,21(3):41. doi:10.1007/s10035-019-0896-4

[28] Liu Lu.Dilated polyhedron based discrete element method and its applications on the analysis of ice loads on marine structures[D].Dalian:Dalian University of Technology,2019.

[29] 刘璐.扩展多面体离散元方法及其在海洋结构冰载荷分析中的应用[D].大连:大连理工大学,2019.

[30] Zhou Yangshiqi, Zhao Lanhao, Shao Linyu,et al.A deformable spheropolygon-based discrete element method[J].Rock and Soil Mechanics,2022,43(7):1961-1968.

[31] 周洋诗琦,赵兰浩,邵琳玉,等.可变形圆化多边形离散单元法[J].岩土力学,2022,43(7):1961-1968.

[32] Hou Yongkang.Research on fictitious crack analytic modelling and concrete fracture process zone under combined stress conditions[D].Shijiazhuang:Shijiazhuang Tiedao U-niversity,2020.

[33] 侯永康.复杂应力状态下虚裂纹模型解析与混凝土断裂过程区研究[D].石家庄:石家庄铁道大学,2020.

[34] Wang Xuebin, Tian Feng, Ma Bing,et al.A three-dimensional crackable Lagrangian element method for modelling tensile cracking processes of rock-like materials[J].Chinese Journal of Applied Mechanics,2023,40(5):1171-1179. doi:10.11776/j.issn.1000-4939.2023.05.023

[35] 王学滨,田锋,马冰,等.岩石类材料拉裂过程模拟的三维可开裂拉格朗日元方法[J].应用力学学报,2023,40(5):1171-1179. doi:10.11776/j.issn.1000-4939.2023.05.023

[36] Jing Pengxu, Yin Chao, Lijun Men,et al.Dynamic response analysis of rock slope based on viscoelastic artificial boundary condition[J].Journal of Water Resources and Architectural Engineering,2021,19(1):190-194. doi:10.3969/j.issn.1672-1144.2021.01.031

[37] 景鹏旭,尹超,门丽君,等.基于黏弹性人工边界条件的岩质边坡动力反应分析[J].水利与建筑工程学报,2021,19(1):190-194. doi:10.3969/j.issn.1672-1144.2021.01.031

[38] Zhao Boming, Xia Chen.A study on visco-elastic artificial boundary conditions under multiple source input[J].China Civil Engineering Journal,2015,48(Supp1):147-151.

[39] 赵伯明,夏晨.多源输入条件下的黏弹性人工边界研究[J].土木工程学报,2015,48():147-151.

[40] Yan Chengzeng.Simulating thermal cracking of rock using FDEM-TM method[J].Chinese Journal of Geotechnical Engineering,2018,40(7):1198-1204. doi:10.11779/CJGE201807005

[41] 严成增.FDEM-TM方法模拟岩石热破裂[J].岩土工程学报,2018,40(7):1198-1204. doi:10.11779/CJGE201807005

[42] Zhou Xiaoping, Cheng Hao.Multidimensional space method for geometrically nonlinear problems under total Lagrangian formulation based on the extended finite-element method[J].Journal of Engineering Mechanics,2017,143(7):04017036. doi:10.1061/(ASCE)EM.1943-7889.0001241

[43] Wang Lu, Xu Fei, Yang Yang.Improvement of the total Lagrangian SPH and its application in impact problems[J].Chinese Journal of Theoretical and Applied Mechanics,2022,54(12):3297-3309.

[44] 王璐,徐绯,杨扬.完全拉格朗日SPH在冲击问题中的改进和应用[J].力学学报,2022,54(12):3297-3309.

[45] Deeks A J, Randolph M F.Axisymmetric time-domain tr-ansmitting boundaries[J].Journal of Engineering Mechanics,1994,120(1):25-42. doi:10.1061/(asce)0733-9399(1994)120:1(25)

[46] Shu Cheng, Wu Yonghong, Du Mengxiang.Seismic response analysis of structure supported by foundations with different locations based on viscoelastic artificial boundary[J].Industrial Safety and Environmental Protection,2022,48(2):47-52. doi:10.3969/j.issn.1001-425X.2022.02.011

[47] 舒丞,吴永红,杜孟翔.基于粘弹性人工边界的掉层结构地震响应分析[J].工业安全与环保,2022,48(2):47-52. doi:10.3969/j.issn.1001-425X.2022.02.011

[48] Patel S, Martin C D.Evaluation of tensile Young’s modulus and Poisson’s ratio of a bi-modular rock from the displacement measurements in a Brazilian test[J].Rock Mechanics and Rock Engineering,2018,51(2):361-373. doi:10.1007/s00603-017-1345-5

[49] Murray Y, Abu-Odeh A, Bligh R.Evaluation of LS-DYNA Concrete Material Model 159[R].Washington, DC:Federal Highway Administration,2007.

[50] Xing Yangyang, Niu Zhiwei, Zeng Shuyuan.Analysis on dynamic stability of Dr2 dangerous rock mass at Suofengying Hydropower Station[J].Water Resources and Hydropower Engineering,2019,50(3):179-185.

[51] 邢洋阳,牛志伟,曾树元.索风营水电站Dr2危岩体动力稳定分析[J].水利水电技术,2019,50(3):179-185.

基金资助

国家重点研发计划项目(2022YFC3005402)

中国水利水电科学研究院水利部水工程建设与安全重点实验室开放研究基金项目(202209)

长江科学院开放研究基金项目(CKWV20231171/KY)