基于动态强度折减-改进矢量和法的边坡稳定性分析

王思源; 童小东

doi:10.11779/CJGE20220548

基于动态强度折减-改进矢量和法的边坡稳定性分析 English Version

王思源^1,,
童小东^{1, 2, ,}

1.
东南大学土木工程学院, 江苏南京 211189
2.
东南大学混凝土及预应力混凝土结构教育部重点实验室, 江苏南京 211189

基金项目:

国家自然科学基金项目 51078082

详细信息

作者简介:
王思源(1998—)，男，硕士研究生，主要从事岩土工程数值模拟方面的研究工作。E-mail: wsyseuswjtu@163.com

通讯作者:
童小东, E-mail: tongxiaodong88@163.com

中图分类号: TU435
计量
- 文章访问数: 379
- HTML全文浏览量: 73
- PDF下载量: 149
出版历程
- 收稿日期: 2022-05-02
- 网络出版日期: 2023-02-23
- 刊出日期: 2023-06-30

Stability analysis of slopes based on dynamic strength reduction- improved vector sum method

WANG Siyuan^1,,
TONG Xiaodong^{1, 2, ,}

1.
School of Civil Engineering, Southeast University, Nanjing 211189, China
2.
Key Laboratory of Concrete and Prestressed Concrete Structures, Ministry of Education (Southeast University), Nanjing 211189, China

摘要

摘要: 在边坡稳定性分析中，极限平衡法与矢量和法均假定边坡失稳是由剪切破坏造成的，因此难以正确求解坡体上部存在拉裂区这类复杂边坡的安全系数，应用整体强度折减法虽能求得前述复杂边坡的安全系数，却难以搜索得到包含坡体上部拉裂区的完整滑动面。针对上述问题，采用动态强度折减法来搜索边坡滑动面，并提出了可兼顾边坡张拉破坏与剪切破坏的改进矢量和法，从而可更好地求解边坡稳定性问题。算例表明，应用所提出的方法不仅能搜索得到完整的边坡滑动面，而且求得的安全系数与整体强度折减法基本相同。所提出的方法适用于土质、岩质以及土岩组合等各类型边坡的稳定性分析，且受主滑方向、网格尺寸、模型尺寸的影响较小，有助于推动边坡稳定性分析方法的进一步完善。
- 边坡稳定性 /
- 张拉破坏 /
- 动态强度折减法 /
- 改进矢量和法
Abstract: In the slope stability analysis, both the limit equilibrium method and the vector sum method assume that the slope instability is caused by shear failure, accordingly, it is difficult to correctly solve the safety factor of the complex slope with tensile failure at the upper part, while the global strength reduction method can obtain the safety factor, it is difficult to search for the complete sliding surface containing the tensile failure zone in the upper part of the slope. In view of the above problems, the dynamic strength reduction method is adopted to search for the sliding surface of slope, and an improved vector sum method which can take into account the tension failure and shear failure is proposed, so that the slope stability problem can be solved better. The results of three cases show that the dynamic strength reduction method can completely reflect the sliding surface of slopes, and the safety factor obtained by the improved vector sum method is basically the same as that by the global strength reduction method. The proposed method is applicable to the stability analysis of various types of slopes, such as soil, rock and soil-rock combinations, and is less affected by the main sliding direction, mesh size and model size, which helps to promote the further improvement of slope stability analysis methods.
- slope stability /
- tensile failure /
- dynamic strength reduction method /
- improved vector sum method

HTML全文

0. 引言

岩石、混凝土等脆性材料的内部常常含有大小不一的微裂隙,这些裂纹在一定的条件下会发生扩展继而贯通最后发展为宏观裂缝,这将直接导致结构体强度的降低,严重时会发生破坏。岩石等脆性材料中裂纹萌生、相互发展贯通的机理始终是岩石力学中需要探究的问题,多年来,国内外相关学者已经取得丰硕的成果^[1-4]。由于岩石类材料内部含有大小不一的节理面,同时这些节理面的方向并不固定,因此在利用理论手段探究裂纹扩展机理时十分复杂,往往需要进行大量的简化。利用试验手段研究岩石破裂机理时存在周期长、成本高、试验条件不易控制等问题。因此学者们常采用数值模拟的手段探究岩石破坏机理。

有限元法在模拟结构变形方面取得巨大成功^[5-7],但在模拟裂纹扩展时会存在尖端奇异性等问题,并且往往需要引入外部准则,在计算裂纹扩展过程中需要多次划分网格,虽然后续发展了扩展有限元理论^[8],但在计算不连续位移时仍需其他准则,模拟岩石中复杂的裂纹交叉问题时仍存在较大的困难。近年来新发展的颗粒离散元,在模拟裂纹扩展方面取得了不错的效果^[9],但在模拟岩石材料连续区域时不能表明连续变化的特征,因为岩石的连续区域还是属于连续介质。分子动力学可以较好的反应结构微观特点,但在模拟岩石类材料裂纹扩展问题时存在效率低等问题^[10]。

近年来一种利用非局部特性来描述物质力学行的方法出现^[11],称之为近场动力学。在该理论中不存在位移的微分方程,而是通过积分方程来描述材料的力学特征^[12],而且利用积分方程可以很好的避免传统方法中尖端奇异性的问题。近场动力学理论中已包含相应的破坏准则,因此在模拟岩石裂纹扩展时不需要额外的判别准则。近场动力学理论最早由silling提出,黄丹等^[13-14]将该理论引入国内,针对混凝土破坏等结构破裂问题进行了模拟。周小平等^[15]、谷新保等^[16]利用近场动力学方法解决岩石类材料的破坏问题,模拟了岩石中裂纹在不同条件下的扩展。王超等^[17]在解决冰桨接触及潜艇破冰上浮问题中应用了近场动力学。Zhang等^[18]提出了正交异性的弹性模型并进行断裂问题分析。Cheng等^[19]利用近场动力学方法对功能梯度材料的破坏进行了分析。朱其志等^[20]运用该理论模拟了不同倾角的单裂缝及三裂缝的岩石试件破坏。

以键为基础的近场动力学概念清晰,应用较为广泛^[21],其物质点之间在一定范围内发生相互作用。PMB模型中一点处的破坏由临界断裂能来判定,若点之间变形储存的能量大于临界值则发生破坏,不再发生相互作用。在近场动力学中发生破坏点处的力学参数如弹性模量保持不变,而在以往的混凝土、岩石破裂的理论与试验中,人们发现裂纹尖端处有效弹性模量是不断降低的,这会促进裂纹的发展,该理论也被称为应变软化准则^[22]。

本文结合前人的研究,在近场动力学中引入应变软化的应变能密度（SED）准则,将岩石类材料裂纹尖端处的力学指标与该点的残余应变能建立联系,同时引入威布尔分布理论以满足岩石材料的异质性。利用所提出的方法对不同倾角单裂纹岩石试件及双平行裂纹岩石试件在单轴压缩条件下的裂纹扩展过程进行了模拟,并分析其裂纹扩展机理,最后将模拟结果与室内试验结果对比。结果表明,本文所提出的方法可以很好的模拟单裂纹及双平行裂纹岩石试件的裂纹扩展情况。

1. 近场动力学方法

1.1 基本理论

与连续介质力学相比,近场动力学理论最大的不同是利用积分方程代替微分方程,可以较好的处理不连续问题。如图1所示,该方法基本原理是在一定范围内点之间存在相互作用,点x在t时刻的控制方程为^[21]

图 1 物质点相互作用

Figure 1. Interaction of material points

下载: 全尺寸图片幻灯片

$ρ (x) \ddot{u} (x, t) = \int_{H_{x}} f (x^{'}, x, u, u^{'}, t) d V^{'} + b (x, t)$ ,

(1)

式中, $x^{'}$ 为距离 $x$ 半径 $δ$ 内相互作用点,在 $t$ 时刻位移分别为 $u (x, t)$ 与 $u (x^{'}, t)$ ,初始时的相对位置向量为 $x^{'} - x = ξ$ ,发生变形时的相对位移向量为 $u (x^{'}, t) -$ $u (x, t) = η$ ,因此变形后的相对的位置向量为 $ξ + η$ 。方程中 $ρ (x)$ 与 $\ddot{u} (x, t)$ 分别点 $x$ 的密度与加速度,而 $b (x, t)$ 为外力荷载, $f$ 为相互作用力函数。

点之间的相互作用,称为“键作用”,可以比作弹簧中的弹性力。在连续介质力学中,相互作用是通过应力直接接触,在有限距离内相互作用的概念是近场动力学理论与传统理论的基本区别。近场动力学目前已经发展了键型、常规态型与非常规态型等分支,本文在键为基础的近场动力学中引入应变能密度（SED）软化准则,另外两种在这里暂时不做讨论。

1.2 微弹性模型

与本节以微弹性模型为例介绍近场动力学基本理论,在该模型中相互作用力大小相等,方向相反。作用力函数f可由微势能求得^[21]：

$f (ξ, η) = \frac{\partial ω (ξ, η)}{\partial η} (\forall ξ, η)$ 。

(2)

微势能 $ω$ 为每个键中单位体积平方的能量,代表物质点作用关系,因此一点处应变能密度W为

$W = \frac{1}{2} \int_{H_{x}} ω (ξ, η) d V^{'}$ ,

(3)

式中,1/2代表每个键储存的能量对 $x$ 贡献一半的应变能密度。在微弹性模型中,相互作用力的方向与变形方向相同,作用力 $f$ 可以表示为

$f (ξ, η) = f (ξ, η) \frac{ξ + η}{| ξ + η |}$ 。

(4)

键的变形通过伸长率 $s$ 表示,不仅能体现点之间的位置变化,而且可以控制键的破坏裂,

$s = \frac{| ξ + η | - | ξ |}{| η |}$ 。

(5)

$| ξ |$ , $| ξ + η |$ 分别为变形前后的长度,与应变定义类似。取值可为正或负,表示键的拉伸或者压缩,当变形超过临界值 $s_{0}$ 时,键会发生不可逆的破坏。微弹性模型具体可参考文献[21],微势能函数 $ω$ 的表达式为

$ω (ξ, η) = \frac{c s^{2} | ξ |}{2}$ ,

(6)

式中, $c$ 为微弹性模量,与弹性模量E类似,代表材料的变形特征。根据式（2）,（4）可得作用力函数：

$f (ξ, η) = c s μ \frac{ξ + η}{| ξ + η |}$ ,

(7)

式中, $μ$ 为控制键断裂的函数, $μ = {\begin{cases} 1 (s < s_{0}) \\ 2 (其他) \end{cases}$ 。

依照式（7）,作用点之间的伸长率大于临界值 $s_{0}$ 时,作用力会不可逆的消失。Silling^[11]推导了不同条件下的c值,并证明该模型在三维模型中泊松比为1/4,平面应力模型中的泊松比为1/3。在三维模型中：

$c = \frac{18 k}{π δ^{4}}$ ,

(8)

式中, $k$ 为体积模量, $δ$ 为作用半径。临界伸长率 $s_{0}$ 由单位面积开裂所需能量的临界释放率 $G_{0}$ 确定：

$s_{0} = \sqrt{\frac{5 G_{0}}{9 k δ}}$ 。

(9)

近场动力学方法中无需引入额外准则来判定破坏,而是通过研究点连接所有键的破坏情况,来定义该点处的损伤：

$D (x) = 1 - \frac{\int_{H_{x}} μ (η, ξ) d V}{\int_{H_{x}} d V}$ ,

(10)

式中, $D (x)$ 值域为0到1,0表明点 $x$ 处还未破坏,而1表明完全破坏。

1.3 离散化与时间积分

为求解平衡方程（1）,首先进行模型均匀离散化,某时刻 $x_{i}$ 与邻域内 $x_{j}$ 相互作用的离散方程为

$ρ {\ddot{u}}_{i}^{n} = \sum_{j} f (u_{j}^{n} - u_{i}^{n}, x_{j} - x_{i}) V_{j} + b_{i}^{n}$ ,

(11)

式中, $n$ 为时间步, $V_{j}$ 为域内相互作用 $x_{j}$ 的体积,离散时作用点有时并不是完全在处于邻域半径内,需按比例进行折减：

$V_{j} = {\begin{cases} Δ x^{3} (| ξ | \leq δ - 0.5 Δ x) \\ (\frac{δ + 0.5 Δ x - | ξ |}{Δ x}) Δ x^{3} (δ - 0.5 Δ x < | ξ | < δ + 0.5 Δ x) \\ 0 (| ξ | \geq δ + 0.5 Δ x) \end{cases} 。$

时间离散可采用Verlet-Velocity显示差分：

$\begin{array}{l} {\dot{u}}_{n + 1 / 2} = {\dot{u}}_{n} + \frac{Δ t}{2} {\ddot{u}}_{n}, \\ u_{n + 1} = u_{n} + Δ t u_{n + 1 / 2}, \\ {\dot{u}}_{n + 1} = {\dot{u}}_{n + 1 / 2} + \frac{Δ t}{2} {\ddot{u}}_{n + 1} 。 \end{array}}$

(12)

参考文献[21]可得到稳定积分的时间步：

$Δ t \leq \frac{β_{safe} {(| ξ |)}_{\min}}{{(c_{k})}_{\max}}$ 。

(13)

2. 应变能密度软化准则

2.1 应变能密度方程

应变能密度理论（strain energy density, SED）,最早由Sih^[23]在研究不同加载条件下材料破坏时提出。SED同时考虑了裂纹尖端处的多个应力分量,可较好地应用于岩石等脆性材料在复杂的边界条件下的破坏过程研究。变形而储存在单元中的能量称为应变能W,SED为单位体积储存的能量^[23]：

$\frac{d W}{d V} = \int_{0}^{ε_{i j}} σ_{i j} d ε_{i j}$ ,

(14)

式中, $σ_{i j}$ , $ε_{i j}$ 分别为应力与应变分量,近场动力学与连续介质力学通常通过应变能密度建立联系。岩石类材料裂纹扩展时通常伴随着不可逆的能量耗散,会引发岩石逐渐破坏,继而导致材料强度与裂纹尖端处力学参数的降低,可利用应变软化概念描述岩石的渐进破坏现象。

岩石应力到达峰值后若持续加载,材料强度会随着岩石变形增加不断降低,因变形引发强度降低的特点,称为应变软化^[22],如图2所示。

图 2 岩石应变软化应力应变曲线

Figure 2. Strain-softening stress-strain curves of rock

下载: 全尺寸图片幻灯片

岩石单元应力达到A点,由于该阶段未产生破坏因此如果此时卸载,其路径沿AO返回。若持续加载,应力到达D点时,岩石单元已发生一定破坏,面积OABDG为目前吸收的应变能W_a,如果在D点卸载,其路径沿DE返回（不考虑塑性变形则沿DO返回）,面积OABDE（不考虑塑性时为OABDO）表示已耗散应变能W,而面积DEG（不考虑塑性时为DOG）表示可恢复应变能W_e,面积GDFCN则代表残余应变能W_r,面积OABDFCN表示岩石单元全应变能W,以上关系为

$(\frac{d W}{d V}) = {(\frac{d W}{d V})}_{a} + {(\frac{d W}{d V})}_{r}$ ,

(15)

${(\frac{d W}{d V})}_{a} = {(\frac{d W}{d V})}_{d} + {(\frac{d W}{d V})}_{e}$ 。

(16)

2.2 SED软化准则的近场动力学

岩石单元抵达屈服点A后,开始逐渐发生破坏,此时可采用折减弹性模量E的方式降低单元的强度,也称有效弹性模量。如图3所示,随着加载持续进行,岩石单元逐渐发生破坏,从而产生能量耗散,不同耗散阶段岩石单元有效弹性模量分别为 $E_{1}$ , $E_{2}$ , $E_{3}$ , $E_{4}$ , $\dots$ , $E_{n}$ （仅考虑弹性时为虚线均经过原点O）。

图 3 岩石破坏后的弹性模量折减

Figure 3. Reduction of elastic modulus after failure of rock

下载: 全尺寸图片幻灯片

结合图2,3,岩石单元发生破坏之后,对应的单元有效弹性模量可以表示为残余应变能密度 $W_{r}$ 的函数：

$E_{n} = f (W_{r})$ 。

(17)

$E_{n}$ 与 $W_{r}$ 之间为复杂的非线性关系,为了简化计算,将有效弹性模量 $E_{n}$ 均分N多个不同的值：

$E_{n} = \frac{N + 1 - n}{N}$ 。

(18)

因此计算时每个单元的有效弹性模量 $E_{n}$ 随单元的破坏程度变化,残余应变能密度 $W_{r}$ 越小,有效弹性模量越小：

（1）当 $\frac{d W}{d V} < {(\frac{d W}{d V})}_{A} = \int_{0}^{ε_{A}} σ_{i j} d ε_{i j}$ 时

岩石单元处于弹性变化阶段,没有发生损伤,有效弹性模量：

$E_{n} = E$ 。

(19)

（2）当 ${(\frac{d W}{d V})}_{A} \leq \frac{d W}{d V} < {(\frac{d W}{d V})}_{C}$ 时

岩石单元处于硬变软化阶段,逐渐发生损伤,有效弹性模量：

$E_{n} = \frac{N + 1 - n}{N}$ 。

(20)

（3）当 ${(\frac{d W}{d V})}_{C} \leq \frac{d W}{d V}$ 时

岩石单元处于残余应力阶段,基本完全损伤,有效弹性模量可认为等于残余弹性模量：

$E_{n} = 0.05 E$ 。

(21)

基于软化准则的近场动力学微观机制可由图4解释,以点 $i$ 为研究对象,变形后点 $i$ 处的应变能储存在与 $j$ 连接的键中并由式（3）计算,结合式（6）,若点 $i$ 发生损伤破坏,应变能密度为

$W = \frac{1}{4} \int_{H_{x}} μ c s^{2} | ξ | d V^{'}$ 。

(22)

图 4 应变能微观机制

Figure 4. Microscopic mechanism of strain energy

下载: 全尺寸图片幻灯片

结合2.1节可得

$W = {(\frac{d W}{d V})}_{e} + {(\frac{d W}{d V})}_{r}$ 。

(23)

断裂键储存的能量为耗散应变能 $W_{d}$ ,结合3.2节可得点 $i$ 处有效弹性模量为 $E_{n}$ ,计算点 $i$ 与邻域内点作用力时微模量可取两点的平均值：

$f (ξ, η) = \frac{c_{i} + c_{j}}{2} s μ \frac{ξ + η}{| ξ + η |}$ 。

(24)

岩石类材料内部含有大量发育程度及开裂条件不同的微裂隙,每个点具有不同的裂纹扩展条件,因此本文引入服从威布尔分布的临界破坏条件来描述岩石的异质性,如下函数定义：

$f (u) = \frac{m}{u_{0}} {(\frac{u}{u_{0}})}^{m - 1} \exp [- {(\frac{u}{u_{0}})}^{m}]$ ,

(25)

式中, $u$ 为单元参数（如强度或弹性模量）,比例参数 $u_{0}$ 为参数 $u$ 的平均值, $m$ 为形状参数定义分布函数的形状,不同 $m$ 的分布函数如图5所示, $m$ 越小微观特性越不均匀,岩体质量及完整性越差,反之岩体质量越好,实际模拟时可根据不同岩石的完整程度进行调试以选择合适的 $m$ 值。这里 $u$ 取临界能量释放率 $G_{0}$ ,则每个单元为 $G_{i}$ ,伸长率为

$s_{i} = \sqrt{\frac{5 G_{i}}{9 k δ}}$ 。

(26)

图 5 概率密度图

Figure 5. Map of probability density

下载: 全尺寸图片幻灯片

计算点i与邻域内点作用时判定破坏的临界伸长率同样取两点的平均值。

3. 单裂纹岩石破坏模拟

3.1 数值模拟方案

为验证方法的有效性,本节对含不同倾角裂纹岩石在单轴压缩情况下的破坏过程进行了模拟,并且与前人试验结果进行对比。利用文献[24]的试验模型,如图6所示,岩石裂纹的角度分别为30°,45°,60°,在端施加持续位移边界条件,得到试件裂纹扩展过程。

图 6 预制单裂纹模型

Figure 6. Model for prefabricated single crack

下载: 全尺寸图片幻灯片

含裂纹岩石试件尺寸100 mm×50 mm×20 mm,中间含有长度28 mm不同倾角的预制裂纹,由于裂纹在厚度方向影响较小,因此计算时采用平面应力方式进行计算,试件离散为200×100计算节点,并且节点的间距 $Δ x = 0 .0005 m$ ,作用半径依照^[19]取 $δ = 3 Δ x$ 。计算中弹性模量取 $E = 4.71 GPa$ ,而密度取 $ρ = 1810$ ${kg/m}^{3}$ ,依照式（11）确定稳定时间步,临界伸长率可由岩石抗拉强度 $f_{t}$ 得

$s_{0} = \frac{f_{t}}{E}$ 。

(27)

试件抗拉强度 $1 .5 MPa$ ,由式（27）可得临界伸长率 $s_{0} = 0.00318$ 。

3.2 形状参数影响

由于在计算时引入了服从威布尔分布的临界破坏条件,因此本节探究形状参数m对计算结果的影响,并确定适合本次岩石破坏模拟的m值,模拟方案如3.1节所述,预制角度取30°,m值从非均匀到均匀依次取2,6,10,20,50,100。依照文献[24]取相同的加载条件,计算总时间步为1400步,最终计算结果见图7。

图 7 形状参数影响

Figure 7. Influences of shape parameter

下载: 全尺寸图片幻灯片

由图7可知,形状参数m对最终计算结果具有较大影响,在相同的加载条件下,形状参数m越小,不均匀性越强,试件破坏程度越强,呈现出一定的弱化效应。

形状参数m为2时,不均匀性非常强时,含裂纹试件的破坏形式不规则,但整体还是沿翼型裂纹、次生共面剪切裂纹及反翼型裂纹破坏。

形状参数m为6,10,20时,不均匀性较强时,含裂纹试件出现明显的翼型裂纹、次生共面剪切裂纹及反翼型裂纹破坏。

形状参数m为50,100时,不均匀性较弱时,含裂纹试件破坏比较规则,在相同的加载步下仅出现明显对称的翼型裂纹。

结合实际破坏情况与本节计算结果,本文计算时形状参数m取10,其余计算时,可根据实际岩石的不均匀程度进行调整。

3.3 不同倾角裂纹扩展模拟

确定合适的形状参数m后,本节将探究不同预制倾角裂隙在单轴压缩条件下裂纹扩展机理,并将最终计算结果与前人试验及模拟结果做出对比,以验证新方法的有效性,数值结果可以为试验现象提供一个清晰的解释,增强裂纹扩展机理的理解。

模拟方案如3.1节所述,预制裂纹的角度分别为30°,45°,60°,在边界施加与前文相同的持续位移边界条件,形状参数m由3.2节计算结果取10,计算总时间步为1400步。

预制角度45°裂纹随时间扩展如图8所示,模拟结果表明,随加载过程持续,新产生的裂纹（包括翼型裂纹、剪切裂纹及反翼型裂纹）均从裂纹尖端产生并扩展。

图 8 单裂隙裂纹扩展过程

Figure 8. Crack growth process of single crack

下载: 全尺寸图片幻灯片

在加载的初期阶段,产生垂直于预制裂纹的翼型裂纹,扩展一段距离后裂纹方向趋向于最大主应力方向,翼型裂纹扩展到与预制裂纹长度大致相等时停止扩展,此时次生共面剪切裂纹（以下简称剪切裂纹）出现并持续扩展,扩展一段距离后出现反翼型裂纹,反翼型裂纹与翼型裂纹位置相同方向相反,以上裂纹扩展机理与这与LI等^[25]试验观测及ZHU等^[20]模拟结果基本相同,本次模拟裂纹扩展机理如图9所示。

图 9 单裂隙裂纹扩展机理

Figure 9. Mechanism of single crack growth

下载: 全尺寸图片幻灯片

为验证本文模拟结果的有效性,将本文结果与文献[24]的试验及模拟结果做出对比,由于文献[22]仅展示了最终结果,因此将预制倾角30°,45°,60°最终模拟结果与原文做出对比,如图10所示。

图 10 不同倾角裂纹扩展结果对比

Figure 10. Comparison of crack growth results with different inclination angles

下载: 全尺寸图片幻灯片

当预制裂纹 $α =$ 30°时,裂纹的主要模式为翼型裂纹与剪切裂纹,翼型裂纹发展到一定程度后产生剪切裂纹最终产生宏观破坏,此外在试验中由于岩石非均质性较强,有可能发生裂纹的合并,而且在模拟中低角度裂纹对反翼型裂纹起到一定的抑制作用,与试验结果吻合。

当预制裂纹 $α =$ 45°时,翼型裂纹扩展方向首先与预制裂纹垂直,随后沿加载方向扩展,该现象本文模拟结果与试验结果非常一致,翼型裂纹扩展到一定程度后产生剪切裂纹与反翼型裂纹最终引发破坏。

当预制裂纹 $α =$ 60°时,裂纹扩展方式与前面类似,试验结果中下半部分翼型裂纹扩展明显模拟结果可以很好的捕捉到这一过程,此外反翼型裂纹扩展程度较大,并且可能与剪切裂纹合并导致宏观破坏。

通过对比不同的结果,可以看出本文模拟出的结果与前人试验、模拟结果具有良好的一致性,并且本文方法可以清晰的反映出裂纹扩展的路径及某一点的损伤程度,此外本文方法无需类似文献[24]引入区别拉伸与剪切破坏的准则,所有的破坏均是在同样的判别准则下发生。

与经典近场动力学类似模拟结果相比本文方法可以很好的捕捉到翼型裂纹先垂直预制裂纹随后按最大主应力方向的裂纹扩展过程进一步揭示裂纹扩展机理,另外基于威布尔分布的临界破坏准则能够较好地体现岩石非均质性,反映非均匀破坏机理。

4. 双平行裂纹岩石破坏模拟

4.1 数值模拟方案

为进一步验证方法的有效性,本节对含不同倾角预制双平行裂纹岩石在单轴压缩情况下的破坏过程进行了模拟,并且与文献[24]结果进行对比。

含裂纹岩石试件尺寸、离散方法、参数选取均与3.1节一致,不在此阐述,所不同的是中间预制45°的平行双裂纹如图11所示。

图 11 预制双平行裂纹模型

Figure 11. Model for prefabricated double parallel cracks

下载: 全尺寸图片幻灯片

4.2 双裂纹扩展过程

预制角度45°双平行裂纹随时间扩展如图12所示,模拟结果表明,随加载过程持续,新裂纹从裂纹尖端萌生、扩展、相互贯通。

图 12 双裂隙裂纹扩展过程

Figure 12. Crack growth process of double cracks

下载: 全尺寸图片幻灯片

如图12所示,加载初期阶段,预制双裂纹扩展机理与单裂纹时类似,分别产生独立的翼型裂纹,随着加载的持续进行,产生的翼型裂纹逐渐接近预制裂纹,继而发生裂纹的第一次贯通。

裂纹发生第一次贯通以后,预制双裂纹处萌生剪切裂纹,产生的位置为别为上预制裂纹的左侧,下预制裂纹的右侧,随着剪切裂纹持续扩展,逐渐与第一次贯通时产生的裂纹连接,从而引发第二次贯通,并且形成闭合的破坏环。

裂纹发生第二次贯通以后,裂纹的扩展形式仍以剪切裂纹为主,并未出现明显的反翼型裂纹,宏观裂纹的扩展沿着上预制裂纹的右侧剪切裂纹与下预制裂纹的左侧剪切裂纹方向,最终引发破坏。

4.3 结果对比

本文近场动力学模拟结果与文献[24]试验、数值模拟结果对比如图13所示。

图 13 双平行裂纹扩展结果对比

Figure 13. Comparison of double-parallel crack growth results

下载: 全尺寸图片幻灯片

通过对比,可以看出本文模拟结果与试验结果具有很好的一致性,不仅可以清晰模拟出试验中每一条裂纹,并且裂纹扩展的位置也与试验基本相同,试验结果也存在数值模拟出现的闭合破坏环,因此可以验证本文提出方法的有效性。

5. 结论

（1）基于SED准则的近场动力学在模拟岩石类材料裂纹扩展方面具有较大的优势,能够较直观地反应结构破裂的过程,与传统方法相比不存在尖端奇异性的问题,也无需要引入拉伸剪切等外部准则。

（2）基于威布尔分布的临界破坏条件可较好的描述岩石异质性,结合SED准则能够较真实的模拟裂隙岩石单轴压缩条件的裂纹扩展过程,反应单裂纹单轴压缩条件下的破裂机理。

（3）预制双平行裂纹岩石破坏过程可分为3个阶段：翼型裂纹上下贯通；剪切裂纹贯通形成闭合破坏环；剪切裂纹形成宏观裂纹引发破坏。

图 1 屈服接近度求解示意图

Figure 1. Schematic diagram of yield approach index

下载: 全尺寸图片幻灯片

图 2 动态强度折减法搜索滑动面流程图

Figure 2. Flow chart of sliding surface search by dynamic strength reduction method

下载: 全尺寸图片幻灯片

图 3 滑动面受力分析

Figure 3. Stress analysis diagram of sliding surface

下载: 全尺寸图片幻灯片

图 4 土质边坡计算模型

Figure 4. Calculation models of soil slopes

下载: 全尺寸图片幻灯片

图 5 土质边坡滑动面

Figure 5. Sliding surfaces of soil slopes

下载: 全尺寸图片幻灯片

图 6 岩质边坡计算模型

Figure 6. Calculation model of rock slope

下载: 全尺寸图片幻灯片

图 7 岩质边坡滑动面

Figure 7. Sliding surface of rock slope

下载: 全尺寸图片幻灯片

图 8 岩质边坡塑性区

Figure 8. Plastic zones of rock slope

下载: 全尺寸图片幻灯片

图 9 安全系数与网格尺寸的关系曲线

Figure 9. Relationship between safety factor and mesh size

下载: 全尺寸图片幻灯片

图 10 安全系数与模型尺寸的关系曲线

Figure 10. Relationship between safety factor and model size

下载: 全尺寸图片幻灯片

图 11 应力边界条件模型

Figure 11. Model of stress boundary condition

下载: 全尺寸图片幻灯片

表 1 土质边坡土体物理力学指标的取值

Table 1 Physico-mechanical indices of soil slopes

土休		黏聚力/kPa	内摩擦角/(°)	重度/ (kN·m^-3)
均质土坡		20.00	20.0	20.00
非均质土坡	^#1土层	49.00	29.0	20.38
	^#2土层	0	30.0	17.64
	^#3土层	7.84	20.0	20.38
	^#4土层	0	30.0	17.64

下载: 导出CSV

表 2 土质边坡计算结果对比

Table 2 Comparison of calculated results of soil slopes

计算方法	均质土坡		非均质土坡
计算方法	安全系数	相对误差/%	安全系数	相对误差/%
整体强度折减法	1.293	—	1.363	—
Morgenstern-Price	1.271	-1.70	1.440	+5.65
本文方法(1)	1.277	-1.24	1.410	+3.45
本文方法(2)	1.275	-1.39	1.410	+3.45
本文方法(3)	1.278	-1.16	1.410	+3.45
传统矢量和法(1)	1.291	-0.15	1.234	-9.46
传统矢量和法(2)	1.284	-0.70	1.234	-9.46
传统矢量和法(3)	1.293	0	1.233	-9.54
应力代数和法	1.275	-1.39	1.412	+3.60
注：本表以整体强度折减法的计算结果为基准进行比较；(1)、(2)、(3)分别表示以极限抗滑力方向、滑动力方向、滑入点指向滑出点方向作为边坡的主滑方向(其中，均质土坡主滑方向与水平方向的夹角分别为28.21°，26.13°，29.25°；非均质土坡主滑方向与水平方向的夹角分别为24.98°，25.10°，21.41°)。

下载: 导出CSV

表 3 岩质边坡岩土体物理力学指标

Table 3 Physico-mechanical indices of rock slope

介质	弹性模量/GPa	泊松比	黏聚力/MPa	内摩擦角/(°)	天然重度/(kN·m^-3)	抗拉强度/MPa
f231	2.0	0.28	0.9	22.8	25.8	0
卸荷裂隙	2.0	0.28	2.0	36.0	26.2	0
V2	2.0	0.27	1.8	21.8	22.1	0
V1	4.0	0.27	2.0	26.5	24.5	1.00
IV	6.0	0.26	7.0	38.6	25.8	3.35
III2	7.5	0.23	17.5	51.3	26.2	6.00
III1	8.0	0.24	15.0	50.2	26.2	5.50
III	9.0	0.22	20.0	52.5	26.5	7.00

下载: 导出CSV

表 4 岩质边坡计算结果对比

Table 4 Comparison of calculated results of rock slopes

计算方法	岩质边坡
计算方法	安全系数	相对误差/%
整体强度折减法	1.659	—
Morgenstern-Price	12.800	+671.55
本文方法(1)	1.657	-0.12
本文方法(2)	1.660	+0.06
本文方法(3)	1.652	-0.42
传统矢量和法(1)	1.854	+11.75
传统矢量和法(2)	1.841	+10.97
传统矢量和法(3)	1.773	+6.87
应力代数和法	1.980	+19.35
注：(1)、(2)、(3)对应的边坡主滑方向与水平方向的夹角分别为34.09°，36.96°，51.08°。

下载: 导出CSV

参考文献(19)

[1]	吴志轩, 张大峰, 孔郁斐, 等. 基-填界面开挖台阶对顺坡填筑高边坡稳定性影响研究[J]. 工程力学, 2019, 36(12): 90-97. https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX201912011.htm WU Zhixuan, ZHANG Dafeng, KONG Yufei, et al. Study on stability influence of high slope foundation-fill interfacial excavation steps[J]. Engineering Mechanics, 2019, 36(12): 90-97. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX201912011.htm
[2]	SARKAR S, CHAKRABORTY M. Stability analysis for two-layered slopes by using the strength reduction method[J]. International Journal of Geo-Engineering, 2021, 12(1): 1-22. doi: 10.1186/s40703-020-00130-3
[3]	杨文辉, 章定文, 闫茜, 等. 深浅根混种法加固边坡稳定性的数值分析[J]. 东南大学学报(自然科学版), 2020, 50(1): 161-168. https://www.cnki.com.cn/Article/CJFDTOTAL-DNDX202001021.htm YANG Wenhui, ZHANG Dingwen, YAN Qian, et al. Numerical analysis on stability of slope reinforced by combination of deep and shallow roots[J]. Journal of Southeast University (Natural Science Edition), 2020, 50(1): 161-168. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DNDX202001021.htm
[4]	郑颖人, 赵尚毅. 有限元强度折减法在土坡与岩坡中的应用[J]. 岩石力学与工程学报, 2004, 23(19): 3381-3388. doi: 10.3321/j.issn:1000-6915.2004.19.029 ZHENG Yingren, ZHAO Shangyi. Application of strength reduction fem in soil and rock slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(19): 3381-3388. (in Chinese) doi: 10.3321/j.issn:1000-6915.2004.19.029
[5]	陈国庆, 黄润秋, 周辉, 等. 边坡渐进破坏的动态强度折减法研究[J]. 岩土力学, 2013, 34(4): 1140-1146. doi: 10.16285/j.rsm.2013.04.040 CHEN Guoqing, HUANG Runqiu, ZHOU Hui, et al. Research on progressive failure for slope using dynamic strength reduction method[J]. Rock and Soil Mechanics, 2013, 34(4): 1140-1146. (in Chinese) doi: 10.16285/j.rsm.2013.04.040
[6]	陈国庆, 黄润秋, 石豫川, 等. 基于动态和整体强度折减法的边坡稳定性分析[J]. 岩石力学与工程学报, 2014, 33(2): 243-256. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201402006.htm CHEN Guoqing, HUANG Runqiu, SHI Yuchuan, et al. Stability analysis of slope based on dynamic and whole strength reduction methods[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(2): 243-256. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201402006.htm
[7]	侯世伟, 马士贺, 李宏男, 等. 基于局部强度阶梯折减法的边坡渐进破坏研究[J]. 防灾减灾工程学报, 2020, 40(1): 72-78. https://www.cnki.com.cn/Article/CJFDTOTAL-DZXK202001010.htm HOU Shiwei, MA Shihe, LI Hongnan, et al. Research on progressive slope failure based on stepwise reduction method of local strength[J]. Journal of Disaster Prevention and Mitigation Engineering, 2020, 40(1): 72-78. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DZXK202001010.htm
[8]	王述红, 朱承金, 张紫杉, 等. 基于动态强度折减DDA法的边坡多滑面稳定性分析[J]. 煤炭学报, 2019, 44(4): 1084-1091. https://www.cnki.com.cn/Article/CJFDTOTAL-MTXB201904013.htm WANG Shuhong, ZHU Chengjin, ZHANG Zishan, et al. Stability analysis of multi-slip surface of slope based on dynamic strength reduction DDA method[J]. Journal of China Coal Society, 2019, 44(4): 1084-1091. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-MTXB201904013.htm
[9]	葛修润. 岩石疲劳破坏的变形控制律、岩土力学试验的实时X射线CT扫描和边坡坝基抗滑稳定分析的新方法[J]. 岩土工程学报, 2008, 30(1): 1-20. http://www.cgejournal.com/cn/article/id/12718 GE Xiurun. Deformation control law of rock fatigue failure, real-time X-ray CT scan of geotechnical testing, and new method of stability analysis of slopes and dam foundations[J]. Chinese Journal of Geotechnical Engineering, 2008, 30(1): 1-20. (in Chinese) http://www.cgejournal.com/cn/article/id/12718
[10]	孙建生. "矢量和法"抗滑稳定计算模型的力学概念错误[J]. 岩土工程学报, 2021, 43(5): 975-980. doi: 10.11779/CJGE202105024 SUN Jiansheng. Mechanical concept errors in anti-sliding stability computational model of "vector sum method"[J]. Chinese Journal of Geotechnical Engineering, 2021, 43(5): 975-980. (in Chinese) doi: 10.11779/CJGE202105024
[11]	王伟, 陈国庆, 郑水全, 等. 考虑张拉-剪切渐进破坏的边坡矢量和法研究[J]. 岩土力学, 2019, 40(增刊1): 468-476. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2019S1070.htm WANG Wei, CHEN Guoqing, ZHENG Shuiquan, et al. Study on the vector sum method of slope considering tensile-shear progressive failure[J]. Rock and Soil Mechanics, 2019, 40(S1): 468-476. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2019S1070.htm
[12]	张海涛, 罗先启, 沈辉, 等. 基于矢量和的滑面应力抗滑稳定分析方法[J]. 岩土力学, 2018, 39(5): 1691-1698, 1708. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201805019.htm ZHANG Haitao, LUO Xianqi, SHEN Hui, et al. Vector-sum-based slip surface stress method for analysing slip mass stability[J]. Rock and Soil Mechanics, 2018, 39(5): 1691-1698, 1708. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201805019.htm
[13]	DIEDERICHS M S, KAISER P K, EBERHARDT E. Damage initiation and propagation in hard rock during tunnelling and the influence of near-face stress rotation[J]. International Journal of Rock Mechanics and Mining Sciences, 2004, 41(5): 785-812.
[14]	周辉, 张传庆, 冯夏庭, 等. 隧道及地下工程围岩的屈服接近度分析[J]. 岩石力学与工程学报, 2005, 24(17): 3083-3087. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200517012.htm ZHOU Hui, ZHANG Chuanqing, FENG Xiating, et al. Analysis of rock mass stability in tunnel and underground engineering based on yield approach index[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(17): 3083-3087. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200517012.htm
[15]	郭明伟, 葛修润, 李春光, 等. 基于矢量和方法的边坡稳定性分析中整体下滑趋势方向的探讨[J]. 岩土工程学报, 2009, 31(4): 577-583. http://www.cgejournal.com/cn/article/id/13220 GUO Mingwei, GE Xiurun, LI Chunguang, et al. Study on potential sliding direction in slope stability analysis based on vector sum method[J]. Chinese Journal of Geotechnical Engineering, 2009, 31(4): 577-583. (in Chinese) http://www.cgejournal.com/cn/article/id/13220
[16]	刘素锦, 郭明伟, 李春光. 三维边坡整体主滑方向的确定[J]. 岩土力学, 2018, 39(S2): 37-44. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2018S2006.htm LIU Sujin, GUO Mingwei, LI Chunguang. Determination of main sliding direction for three-dimensional slope[J]. Rock and Soil Mechanics, 2018, 39(S2): 37-44. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2018S2006.htm
[17]	陈栋, 李红军, 朱凯斌. 基于新主滑趋势方向的矢量和边坡稳定分析方法[J]. 岩土力学, 2021, 42(8): 2207-2214, 2238. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202108016.htm CHEN Dong, LI Hongjun, ZHU Kaibin. Vector sum analysis method for slope stability based on new main sliding trend direction[J]. Rock and Soil Mechanics, 2021, 42(8): 2207-2214, 2238. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202108016.htm
[18]	孙建生. 基于有限元应力载荷宏观刚性滑裂面极限平衡抗滑稳定计算理论[J]. 岩石力学与工程学报, 2018, 37(4): 862-875. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201804009.htm SUN Jiansheng. Theoretical model of stability calculation of macro-rigid sliding planes with FEM stress load based on limit equilibrium[J]. Chinese Journal of Rock Mechanics and Engineering, 2018, 37(4): 862-875. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201804009.htm
[19]	彭少文, 贺金明, 孙冠军, 等. 高陡边坡深部卸荷张拉裂缝的成因机制探讨[J]. 人民长江, 2015, 46(14): 32-35. https://www.cnki.com.cn/Article/CJFDTOTAL-RIVE201514010.htm PENG Shaowen, HE Jinming, SUN Guanjun, et al. Discussion on genetic mechanism of large-depth unloading tension fracture in high and steep slope[J]. Yangtze River, 2015, 46(14): 32-35. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-RIVE201514010.htm

施引文献(22)

期刊类型引用(11)

1.	王恒泰. 基于有限元强度双折减系数法的非均质边坡稳定性分析. 山西建筑. 2025(08): 63-68+72 . 百度学术
2.	刘红修. 基于有限元强度折减法的三维边坡安全系数研究. 铁道建筑技术. 2025(05): 184-187+197 . 百度学术
3.	孙士国，周鑫，彭是焱. 基于残余位移增量判据的土石围堰边坡稳定性分析. 水电能源科学. 2024(02): 81-84 . 百度学术
4.	周东波，梅源，杨彤. 基于响应面法的砂质黄土路堑高边坡稳定性研究. 中国安全生产科学技术. 2024(07): 139-147 . 百度学术
5.	高俊，韩小柯，韩文超，党发宁，廖昌隆. 大倾角坝基上混凝土面板堆石坝滑移变形特征研究. 应用力学学报. 2024(04): 812-821 . 百度学术
6.	郭慧敏，杨兴隆，韩晓东，王探宇，亢有为，余振，梁兴龙. 广西苍梧路堑边坡稳定性分析. 科技和产业. 2024(15): 253-257 . 百度学术
7.	吴罡，阮永芬，李鹏辉，槐以高，乔文件. 极限平衡法与强度折减法在边坡稳定分析中的问题讨论. 贵州大学学报(自然科学版). 2024(05): 70-77 . 百度学术
8.	黄乐源，田城航，仉文岗，孟轩宇，闫国强，王鲁琦. 降雨作用下三维矿山边坡的力学状态响应. 土木与环境工程学报(中英文). 2024(06): 70-79 . 百度学术
9.	刘赞. 基于FLAC 3D的边坡支护方案优化研究. 西部交通科技. 2023(02): 50-53 . 百度学术
10.	侯世伟，焦宏宇，韩俊艳，张皓，李艳凤. 考虑衰减滞后的边坡分层强度折减原理. 土木工程与管理学报. 2023(06): 32-38 . 百度学术
11.	丘永富. 新疆蒙库铁矿露天采场南帮边坡稳定性分析及控制. 现代矿业. 2022(09): 237-240 . 百度学术