一种模拟节理边坡破坏全程的改进SPH方法

An Improved SPH Method Based on Strength Reduction to Simulate Entire Process of Joint Slope Failure

摘要

裂隙岩质边坡失稳是重大地质灾害，其破坏涉及裂纹萌生、扩展与滑移耦合过程.传统数值方法难以兼顾连续破裂与非连续接触的模拟，且存在网格畸变、参数标定复杂等局限.开发了一种基于强度折减以及核函数改进的三维SPH算法，用于模拟三维裂隙岩体边坡破裂和裂纹扩展以及接触滑移全过程.在改进三维SPH方法中，通过强度折减以及带拉伸截断的摩尔库伦破坏准则实现了裂纹萌芽的判断，通过在改进的核函数内引入损伤标志，实现岩体的裂纹扩展，随后引入了损伤粒子点点三维接触准则，构建完整粒子和断裂粒子三维间接触力.首先采用三维单轴压缩试验来验证算法的可行性，并测定了不同倾角的单裂隙岩体的脆性断裂特征.随后将改进的SPH法以及强度折减理论运用于考虑不同节理倾角的多节理岩质边坡，模拟三维裂隙岩质滑坡破坏过程并评估其稳定性.研究结果表明基于强度折减改进的SPH方法在模拟三维裂隙岩坡破坏以及稳定性问题上具有计算效率高、参数标定少以及准确率高的优点，且该方法可被用于评估其他带结构面的岩质边坡的稳定性.

Abstract

The instability of fractured rock slopes is a major geological disaster, and its failure involves the coupling process of crack initiation, propagation and slip. Traditional numerical methods are unable to simultaneously simulate continuous fracture and discontinuous contact, and they have limitations such as grid distortion and complex parameter calibration. A three-dimensional SPH algorithm based on strength reduction and kernel function improvement is developed to simulate the entire process of fracture, crack propagation and contact slip of three-dimensional fractured rock mass slopes. In the improved three-dimensional SPH method, the determination of crack germination was achieved through strength reduction and the Mohr-Coulomb violation criterion with tensile truncation. The crack propagation of rock mass was realized by introducing damage marks in the improved kernel function. Subsequently, the three-dimensional contact criterion of damage particle dots was introduced to construct the three-dimensional contact force between intact particles and fractured particles. Firstly, three-dimensional uniaxial compression tests were adopted to verify the feasibility of the algorithm, and the brittle fracture characteristics of single-fracture rock masses with different inclinations were determined. Subsequently, the improved SPH method and the strength reduction theory were applied to multi-joint rock slopes considering different joint inclination angles to simulate the failure process of three-dimensional fractured rock landslides and evaluate their stability. The research results show that the improved SPH method based on strength reduction has the advantages of high computational efficiency, less parameter calibration and high accuracy in simulating the failure and stability of three-dimensional fractured rock slopes. Moreover, this method can be used to evaluate the stability of other rock slopes with structural planes.

Graphical abstract

关键词

光滑粒子流法 / 强度折减 / 岩质边坡 / 裂纹扩展 / 接触 / 工程地质学.

Key words

smooth particle hydrodynamics / strength reduction / rock slope / crack propagation / contact / engineering geology

0 引言

自然形成的岩石材料中包含由地质作用产生的各种裂纹，包括节理、层理、裂隙、孔洞等.这些裂纹的存在导致岩石材料呈现出非均质性和各向异性（Lee and Leon， 2011）.当岩石被施加载荷时，新的裂纹会在原先就存在于岩石中的原生裂纹上或者附近开始生长，逐步扩展，最终主导整块岩石的破裂行为（Fakhimi et al.， 2014）.由此可见，这些裂纹在很大程度上控制着岩石的力学性质，并且影响着岩质边坡这样的宏观岩石结构体的稳定性.因此，研究裂纹岩体的破裂机制对于工程建设的风险评估工作具有非常重大的意义.

目前研究裂隙岩体的破裂机制主要以试验和数值模拟方法为主.单轴压缩试验是指对岩石块体逐步增加所施加的单向载荷，使得岩石破裂，进而对岩石的裂纹扩展模式和其力学性质进行研究和测定的试验.它为研究在外部载荷作用下，岩石内部原生裂纹的生长和扩展对岩石破碎形式多产生的影响提供了有效途径.例如，在受到载荷的岩石内部，拉伸裂纹在节理倾角小于60°时，相邻裂纹尖端发生聚合，并且裂纹的聚合点随角度的增大从节理面上端移动到下端（Lee and Leon， 2011）.但是试验存在人为误差.造价昂贵等问题（Wang et al.， 2024），难以获取完整的破坏过程.

随着计算机技术和数值方法的不断进步，数值模拟已成为研究岩土材料破裂过程的重要工具（Wu et al.， 2023）.目前岩土领域的数值模拟方法主要分为连续性方法、非连续性方法和连续-非连续耦合方法.连续性方法基于连续介质假设，适用于材料未发生明显破裂或位移较小的场景.其中较为传统的方法为有限元法（FEM）（Zhang et al.， 2023）和有限差分法（FDM）（Do et al.， 2019），这两种方法依赖网格离散，虽计算高效，却在大变形场景下面临网格畸变困境.相比之下近场动力学（PD）方法基于非局部作用理论，基于粒子键断裂描述损伤（Wang et al.， 2018），却因缺乏接触准则难以模拟原生裂纹碎裂.而非连续性方法可以解决裂纹扩展和接触的问题，适用于材料已破裂或由离散块体/颗粒组成的场景，目前的主流研究为两种方法：离散元法（DEM）（Jia et al.， 2024）和不连续变形分析（DDA）（Hatzor et al.， 2004），但非连续性方法的缺陷在于其参数标定困难，计算所需时间较长.连续-非连续耦合方法结合了连续与非连续方法的优势，能够模拟材料从连续变形到破裂的全过程.在交界面处人为的假设太多且在交界面无法避免能量的损耗.较为常用的方法为有限离散元法（FDEM）和有限差分法-有限元耦合方法（FDM-DEM）.FDEM能够较好地反映岩石试样在单轴抗压强度测试、巴西抗拉强度测试和SHPB试验中的破裂情况，但对于裂缝处的局部应力和粒子速度，此方法的计算结果整体偏高（Fukuda et al.， 2020）；FDM-DEM方法能够很好地模拟岩石边坡的变形和破坏演变以及坍塌发展（刘毛毛等，2024），但大质量块体坠落模拟产生的动能较大，不同介质间的能量耗散较高（Wang et al.， 2022）.

光滑粒子流体动力学（SPH）是一种基于拉格朗日框架的无网格粒子方法，通过核函数插值实现连续介质的离散化，由于不受网格约束，能够很好地模拟岩土材料的大变形行为.目前研究者们针对SPH在岩土工程领域的应用多数在于泥石流（Han et al.， 2019）和多场耦合（Yu et al.， 2023）等.SPH在岩石力学领域的研究处于起步阶段.目前SPH方法在模拟岩石破裂所运用的准则主要包括粘结单元、应力跌落法、虚拟键.粘结单元的原理是在相邻粒子间预设粘结界面单元，赋予其独立的强度准则，当界面应力超过阈值时，粘结单元失效，粒子作用力归零.此方法物理含义明确，但难以模拟随即裂纹的扩展（Bui et al.， 2008），例如，一种带有嵌入式框架过程区的新型SPH方法，可以用于岩石断裂建模，并且通过三点弯曲测试、巴西盘测试和半圆弯曲测试验证了方法的可行性（Wang et al.， 2019）；应力跌落法是当粒子受力条件满足破坏准则时，瞬间将其应力跌落至残余应力.此法适用于脆性材料，可以在计算中直接嵌入本构方程，但由于应力突变导致能量不守恒，相应的计算精度也受到一定限制（Xia et al.， 2024），利用基于应力跌落法的SPH方法，可以模拟在应力下降过程中的球应力不变假设下的非均质岩石的破坏过程（Shu et al.， 2020）；虚拟键是指在粒子间建立虚拟弹簧，通过刚度传递作用力，当键的应变超过阈值时，键发生断裂.此法支持随即裂纹扩展，但没有考虑键断裂后粒子之间仍然有继续发生接触的可能性（Yu et al.， 2023），例如，利用一种基于虚拟结合法的SPH方法对具有V形相交双裂纹特征的红砂岩进行剪切载荷测试，并分析其破坏模式和裂纹演变（Zhou et al.， 2023）.目前模型多以破裂为主，缺乏合理的破裂接触耦合准则，这限制了在单向荷载下模型捕捉材料合理的碎裂化过程的能力.目前利用SPH方法针对裂隙岩体进行单轴压缩实验模拟的研究还鲜有人涉足.因此，迫切需要一种破裂接触耦合的SPH方法来更精确地进行单轴压缩实验模拟以及捕捉岩体的碎裂和内部裂纹的生长扩展行为.

强度折减法是模拟边坡稳定性时常用的一种方法，其原理是通过逐步降低材料的强度参数如粘聚力、内摩擦角和抗拉强度或抗剪强度，直到系统达到临界破坏状态，此时的折减系数即为安全系数.强度折减法虽非真实破坏过程，但通过渐进失稳获取安全系数是边坡稳定性评价的常用标准（Liu et al.， 2015； Sun et al.， 2016； Yang et al.， 2019）.尽管折减过程是人为的，但岩体从弹性变形，到裂纹萌生，再到贯通滑移的渐进失稳过程与实际滑坡的力学机制一致.由于土体只存在抗剪强度而没有抗拉强度，故在分析土体边坡的稳定性时常常用到抗剪强度折减法.在FDM和FEM中，通常是同时折减粘聚力、内摩擦角和抗剪强度等指标来获取安全系数（Huang et al.， 2009；Dehghan et al.， 2022）.强度折减的SPH方法常用于分析土质滑坡的稳定性（Zhang et al.， 2019； Li and Wang， 2020），其有效性主要依赖于土体以剪切破坏为主导的特性.然而，对于岩质边坡或岩质滑坡，其破坏机制往往更为复杂，节理、裂隙的存在使得岩体表现出明显的拉裂破坏特征，此时抗拉强度成为控制稳定性的关键因素之一.因此，在SPH中利用抗拉强度折减法分析岩质滑坡稳定性具有重要的理论意义和实用价值，但目前尚未有人研究.本研究正是聚焦于这一关键空白，系统性地提出并应用SPH框架下的抗拉强度折减法来模拟和评估岩质滑坡的稳定性.这显著拓展SPH强度折减法在岩土工程复杂破坏分析中的应用，为理解岩质边坡的稳定性提供了全新的视角.

本研究提出了一种改进的三维SPH数值模拟方法，通过引入强度折减理论和核函数改进技术，成功构建了能够精确模拟三维裂隙岩体边坡从裂纹萌生、扩展到最终滑移破坏全过程的计算框架并获得该边坡的安全系数.该方法的核心创新体现在3个方面：首先，通过结合强度折减理论和带拉伸截断的摩尔库伦破坏准则，实现了对裂纹萌生阶段的准确判断；其次，在改进的核函数中引入损伤标志，有效模拟了裂纹的动态扩展以及滑面贯通过程；最后，提出了损伤粒子间的三维点点接触准则，完善了岩体破坏后断裂粒子与完整粒子间的力学相互作用模型.为验证方法的可靠性，研究首先开展了三维单轴压缩试验模拟，系统分析了不同倾角单裂隙岩体的脆性断裂特征.在此基础上，将改进的SPH方法与强度折减理论应用于多节理岩质边坡的稳定性分析，成功再现了不同节理倾角条件下三维裂隙岩质边坡的渐进破坏过程.并在讨论中综合比较了该算法的优缺点.

1 理论基础

1.1　离散方法

SPH方法的原理是将整个计算区域离散为一系列相互作用的粒子，每个粒子携带质量、密度、速度等物理量.这些粒子按照拉格朗日框架运动，意味着它们会跟随材料一起移动变形.与传统的网格方法不同，SPH不需要预先划分计算网格，这使得它特别适合处理大变形和非连续问题.

光滑粒子流体动力学（SPH）其场函数的积分近似可以表达为：

f(r)\approx {\int }_{D}f\left(r^{\text{'}}\right)W\left(\left|r-r^{\text{'}}\right|,h\right)\mathrm{d}{r}^{\text{'}},(1)

式中：f(r\text{'})表示宏观变量；D表示计算域；W\left(\left|r-r^{\text{'}}\right|,h\right)是一个核函数，其值取决于\left|r-r^{\text{'}}\right|和光滑长度h；f(r)表示积分近似的结果.

核函数是积分近似当中的重要组成部分，所以对于核函数的选择也十分关键.核函数需要满足归一化、delta函数属性和紧缩条件.

{\int }_{D}W\left(r-r^{\text{'}},h\right)\mathrm{d}{r}^{\text{'}}=\mathrm{1}，(2)

W\left(r-r^{\text{'}},h\right)=\mathrm{0},\left|r-r^{\text{'}}\right|>\kappa h .(3)

方程（2）为核函数的归一化过程；方程（3）为紧缩条件，其中\kappa为光滑因子，与光滑长度h一同定义了函数的作用条件.本研究所选用的核函数为受到广泛使用的3次样条核函数，其表达式为：

.(4)

1.2　SPH方程

1.2.1　控制方程

在数值模拟过程当中，每个粒子需要满足连续性、动量守恒、运动方程，其表达式如下：

\left\{\begin{array}{c}\frac{\mathrm{d}{\rho }_i}{\mathrm{d}t}={\rho }_i{\sum }_{j=\mathrm{1}}^N \frac{m_j}{{\rho }_j}\left({\mathit{v}}_i^{\alpha }-{\mathit{v}}_j^{\alpha }\right)\frac{\partial W_{ij}}{\partial {\mathit{x}}_i^{\alpha }}\\ \frac{\mathrm{d}{v}_i^{\alpha }}{\mathrm{d}t}={\sum }_{j=\mathrm{1}}^N m_j\left(\frac{{\sigma }_i^{\alpha \omega }}{{\rho }_i^{\mathrm{2}}}+\frac{{\sigma }_j^{\alpha \omega }}{{\rho }_j^{\mathrm{2}}}+\delta {\Pi }_{ij}+T_{ij}\right)\frac{\partial W_{ij}}{\partial {\mathit{x}}_i^{\alpha }}+\frac{{\mathit{f}}_i^E}{m_i}\\ \frac{\mathrm{d}{\mathit{x}}_i^{\alpha }}{\mathrm{d}t}={\mathit{v}}_i^{\alpha }\end{array}\right.,(5)

式中：{\rho }_i为粒子密度，\alpha与\omega为笛卡尔分量，T_{ij}为人工应力，{\Pi }_{ij}为人工粘度.

1.2.2　本构方程

对于岩石这一类固体材料，柯西应力张量{\sigma }^{\alpha \omega }被分解为静水压力P和偏应力张量\tau，其关系式如下：

{\sigma }^{\alpha \omega }=-P{\delta }^{\alpha \omega }+{\tau }^{\alpha \omega }，(6)

式中：P为静水压力，由体积变形引起，通过状态方程与密度关联.{\delta }^{\alpha \omega }为克罗内克符号，{\tau }^{\alpha \omega }是偏应力张量，表征材料的剪切响应，由应变率与旋转率共同决定.偏应力张量由以下本构方程描述：

\stackrel{.}{{\tau }^{\alpha \omega }}=\mathrm{2}G({\epsilon }^{\alpha \omega }-\frac{\mathrm{1}}{\mathrm{3}}{\delta }^{\alpha \omega }{\epsilon }^{\gamma \gamma })+{\tau }^{\alpha \gamma }{R}^{\omega \gamma }+{\tau }^{\omega \gamma }{R}^{\alpha \gamma }，（7）式中：G为剪切模量，{\epsilon }^{\alpha \omega }为应变率张量，R为旋转率张量，其表达式如下：

{\epsilon }_i^{\alpha \omega }=\frac{\mathrm{1}}{\mathrm{2}}{\sum }_{j=\mathrm{1}}^N\frac{m_j}{{\rho }_j}\left(\left(v_i^{\alpha }-v_j^{\alpha }\right)\frac{\partial W_{ij}}{\partial x_{ij}^{\omega }}+\left(v_i^{\omega }-v_j^{\omega }\right)\frac{\partial W_{ij}}{\partial x_i^{\alpha }}\right)，(8)

R_i^{\alpha \omega }=\frac{\mathrm{1}}{\mathrm{2}}{\sum }_{j=\mathrm{1}}^N\frac{m_j}{{\rho }_j}\left(\left(v_i^{\alpha }-v_j^{\alpha }\right)\frac{\partial W_{ij}}{\partial x_i^{\omega }}-\left(v_i^{\omega }-v_j^{\omega }\right)\frac{\partial W_{ij}}{\partial x_i^{\alpha }}\right) .(9)

1.2.3　时间集成

蛙跳算法（leap-frog algorithm）作为SPH方法常用的时间积分算法，其所用步数较少，效率较高（Liu et al.， 2008），粒子的各项性质随时间推进计算，如下式：

\left\{\begin{array}{ll}t& =t_{\mathrm{0}}+\mathrm{\Delta }t\\ {\rho }_i\left(t_{\mathrm{0}}+\mathrm{\Delta }t/\mathrm{2}\right)& ={\rho }_i\left(t_{\mathrm{0}}\right)+(\mathrm{\Delta }t/\mathrm{2})\mathrm{\Delta }{\rho }_i\left(t_{\mathrm{0}}\right)\\ v_i\left(t_{\mathrm{0}}+\mathrm{\Delta }t/\mathrm{2}\right)& =v_i\left(t_{\mathrm{0}}\right)+(\mathrm{\Delta }t/\mathrm{2})\mathrm{\Delta }{v}_i\left(t_{\mathrm{0}}\right)\\ x_i\left(t_{\mathrm{0}}+\mathrm{\Delta }t/\mathrm{2}\right)& =x_i\left(t_{\mathrm{0}}\right)+(\mathrm{\Delta }t/\mathrm{2})\mathrm{\Delta }{v}_i\left(t_{\mathrm{0}}\right)\end{array}\right.,(10)

式中： t_{\mathrm{0}}为初始时间，\mathrm{\Delta }t为时间步长.

1.3　改进策略

1.3.1　损伤处理

现有SPH损伤模型（如虚拟键法）仅通过键断裂模拟裂纹，忽略断裂后接触.本文改进核函数引入损伤标志，其物理意义为局部刚度衰减因子：当损伤变量达到阈值时粒子完全失效，强制粒子退出连续场并激活接触准则，实现破裂-接触耦合机制.

模拟过程中会出现材料断裂或者失效的情况，并且需要根据摩尔库伦准则针对材料内部的拉伸裂隙和剪切裂隙进行评估，故此需要引入损伤模型.摩尔库伦准则判断材料的剪切与拉伸破坏形式如下：

\left\{\begin{array}{c}{\sigma }_{\mathrm{f}}\ge {\sigma }_{\mathrm{t}}\\ {\tau }_{\mathrm{f}}\ge c+{\sigma }_{\mathrm{t}}\mathrm{t}\mathrm{a}\mathrm{n}\varphi \end{array}\right.,(11)

式中：{\sigma }_{\mathrm{f}}和{\tau }_{\mathrm{f}}为岩石临界失效时的最大拉伸应力和剪切应力，c为粘聚力，\varphi为内摩擦角.

通常引入损伤变量D_i，建立累计损伤模型对材料的断裂进行模拟，其关系式如下：

D_i={\int }_{\mathrm{0}}^t{d}_{i}dt,d_i=f({\sigma }_i,{ϵ}_i)，(12)

式中:D_i为累计损伤变量，d_i代表由应力与应变对材料造成的损伤，上述积分反映为损伤变量D_i随应变或应力变化所造成的历史累积.当D_i=1时，代表岩石完全损伤，粒子失效并且停止传递力的作用；当时，代表岩石部分损伤，刚度逐渐退化.

为了实现将损伤粒子与其他粒子之间的联系切断，本研究引入了损伤标志\mathrm{\zeta }来改进核函数，改进方式如下式：

W_{\mathrm{d}\mathrm{a}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{d}}=\zeta W\left(r_{ij},h\right) .(13)

当粒子损伤变量D_i达到阈值时，损伤标志\zeta =\mathrm{0}，该粒子脱离与其他粒子的相互作用；而对于未受损的粒子，则\zeta =\mathrm{1}.

1.3.2　接触处理

为了实现岩体损伤破坏后沿破坏面的变形、滑移行为，在完整粒子和损伤粒子间加入了接触力.当一损伤粒子与完整粒子之间的距离达成一定的条件时加入接触，条件为：

I^n=\frac{h_i+h_d}{\mathrm{2}}-l_{id}， (14)

式中：l_{id}为完整粒子i与损伤粒子d之间的距离，h_i和h_d分别为两粒子的光滑核半径.当条件达到I^n>\mathrm{0}时，即引入接触力.此时粒子所受到的合力由接触力和其他外力组成，接触力可以分解为法向接触力和切向接触力，其计算方程如下：

F_i^{\mathrm{c}}=F_i^{\mathrm{c}\mathrm{n}}+F_i^{\mathrm{c}\mathrm{s}}，(15)

\left\{\begin{array}{l}{F}_i^{\mathrm{c}\mathrm{n}}=k_n{I}^{\mathrm{n}}n\\ F_i^{\mathrm{c}\mathrm{s}}={\left\{F_i^{\mathrm{c}\mathrm{s}}\right\}}_{\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}}+k_{\mathrm{s}}\mathrm{\Delta }{I}^{\mathrm{s}}\end{array}\right.，(16)

\left\{\begin{array}{c}\frac{\mathrm{1}}{k_{\mathrm{n}}}=\frac{\mathrm{1}}{E_i{r}_i}+\frac{\mathrm{1}}{E_d{r}_d}\\ \frac{\mathrm{1}}{k_{\mathrm{s}}}=\frac{\mathrm{1}}{E_i{\kappa }_i{r}_i}+\frac{\mathrm{1}}{E_d{\kappa }_i{r}_d}\end{array}\right.，(17)

式中：F_i^{\mathrm{c}\mathrm{n}}和F_i^{\mathrm{c}\mathrm{s}}分别为法向和切向接触力，\mathrm{\Delta }{I}^{\mathrm{s}}是剪切增量，n为单位法向量，k_{\mathrm{n}}和k_{\mathrm{s}}分别代表法向/剪切刚度，它们由剪切模量和粒子的光滑核半径共同定义，E是剪切模量，\kappa是剪切和法向刚度之比.

当粒子间的滑动力超过最大摩擦力而产生滑移时，其滑动状态下的剪切力的表达式如下：

F_i^{\mathrm{c}\mathrm{s}}=\mu ‖F_i^{\mathrm{c}\mathrm{n}}‖\frac{F_i^{\mathrm{c}\mathrm{s}}}{‖F_i^{\mathrm{c}\mathrm{s}}‖}‖F_i^{\mathrm{c}\mathrm{s}}‖ .(18)

1.3.3　强度折减法

此时，按照强度折减方法，在0.2 s时施加岩石体的实际强度，经过多次强度折减迭代后，最终折减强度被定义为边坡在滑动面形成时的安全系数（FOS）.基于SPH的边坡强度折减理论如下所示：在此，基于具有拉伸断裂阈值准则的强度折减法被纳入到SPH中.剪切强度和拉伸强度均根据公式（16）进行变化：

\left\{\begin{array}{l}{c}_{\mathrm{f}}=c/K\\ R_{\mathrm{t}\mathrm{f}}=R_{\mathrm{t}}/K\\ {\varphi }_{\mathrm{f}}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}(\mathrm{t}\mathrm{a}\mathrm{n}\varphi /K)\end{array}\right.，（19）

式中：K为强度折减系数，c、R_{\mathrm{t}}和\varphi分别为实际的粘聚力、抗拉强度以及内摩擦角，c_{\mathrm{f}}和R_{\mathrm{t}\mathrm{f}}以及{\varphi }_{\mathrm{f}}分别为折减后的粘聚力、抗拉强度以及内摩擦角.K值定义为当坡体滑动面贯通时的FOS.粘聚力c、内摩擦角\varphi、抗拉强度R_{\mathrm{t}}同步折减，折减终止条件为：当坡体内部的滑动面贯通时，立即停止折减，并将此时的折减系数确定为FOS.

1.3.4　粒子搜索方法

在模拟节理裂隙的结构特征时，可采用基于粒子搜索的几何构建方法.以节理裂隙的生成为例，其具体实施过程如图3.

明确裂隙的几何参数，包括裂隙的起始坐标和终止坐标.在三维空间中，根据这两个端点坐标，可以确定裂隙的空间位置，并在连接线上按固定间隔布置一系列搜索点.每个搜索点配置特定的作用范围，该范围以搜索半径来界定.这里需要控制两个关键尺寸参数：搜索点之间的布置间距需小于模型中基础粒子的平均分布间隔，而单个搜索点的作用范围则应大于基础粒子间的典型距离.当搜索点布置完成后，系统会检测每个搜索点作用范围内包含的基础粒子.这些被搜索范围覆盖的基础粒子即被判定为构成节理裂隙几何边界的组成部分.根据模拟需求，可对这些标识粒子执行两种操作：直接从计算模型中移除，或者将其转移到特定的边界粒子集合中参与后续计算.该方法的核心优势在于通过搜索点间距和搜索半径两个参数可以实现精准调控裂隙的空间位置.搜索点间距越小，生成的裂隙边界越光滑；搜索半径越大，则裂隙的影响带宽度越大.整个过程的计算效率取决于搜索算法的优化程度，特别是大规模模型中的邻近粒子快速检索技术.

2 数值模拟与分析

本文分别对已建模的三维单裂隙岩体进行单轴压缩试验以及对节理岩质边坡进行稳定性分析，并将结果与实际实验结果进行对比，以验证改进SPH方法模拟脆性断裂和边坡失稳的有效性.

2.1　裂隙岩体单轴压缩模拟

建立了与实验（Yang and Huang， 2017）相对应的单裂隙岩体的SPH三维模型（图4）.三维岩体模型尺寸为50 mm × 22 mm× 100 mm，粒子数量为250 000个，粒子初始间距为1 mm.考虑1 mm宽，10 mm长的预制裂隙，其角度设定为0°、30°、45°、60°以及90°.为实现稳定静载，施加于上边界虚粒子的加载速率为0.2 m/s，下边界虚粒子固定.计算的时间步长为1.5×10^-7 s/step.SPH模拟岩石所需的关键参数，如密度、弹性模量、泊松比和内摩擦角，可以直接从实际实验中获取，如表1所示.相比之下，粘聚力和抗拉强度需要通过标定确定.即通过拟合模拟的应力-应变曲线与研究岩石的单轴压缩实验的应力-应变曲线，采用试错法估计其值.人工粘度系数β₁和β₂分别为1和2.

图4为SPH方法计算得到的裂隙倾角θ = 45°的渐进过程失效.根据Griffit初始断裂准则，当最大主应力达到抗拉强度时，预制裂缝两端产生翼型裂纹.在连续载荷作用下，翼型裂纹从预制裂缝尖端开始，沿曲线路径扩展，逐渐接近裂缝中点并与垂直载荷方向平行.翼型裂纹由于拉伸应力引起，表面光滑平整，这与试验高度一致.通过比较实验结果中典型翼型裂纹以及交合裂纹（图4c），SPH的模拟结果与Yang and Huang（2017）的实验结果非常一致.因此，SPH方法被认为可有效地模拟三维非连续岩体的裂纹扩展和脆性破坏特征.图5为不同预制裂隙角度条件下裂隙岩体的单轴压缩破坏特征，可以发现角度变化影响了最终剪切-拉伸裂纹的分布，其中90°以及0°的节理裂隙岩体在高应力作用下的侧面水平剪切裂纹分布加剧，导致了更强的侧向失稳破坏.

2.2　裂隙岩质边坡失稳模拟

在本节中，将SPH方法应用于一个具有代表性的间断多裂隙岩质边坡（图6），该边坡具有不同倾角θ（0°、30°、45°、60°、-30°、-45°和 -60°）的一组平行裂隙.由完整粒子构成的边坡高度为32 m，长度为38 m，宽度为12 m，不同倾角的一组间歇性裂隙的尺寸和分布情况如图7所示.此外，这些裂隙沿z方向穿透岩质边坡.损伤粒子用于表示预制裂隙.所建立的数值模型由114 285个粒子组成，初始粒子间距dx = 0.4 m，光滑核半径设定为1倍初始粒子间距.计算的时间步长为3×10^-5 s.通过使用边坡的底部边界粒子进行x、y、z全方向固定，同时在x方向上固定前后边界，y方向上固定左右边界.岩石基质和层理面的物理力学性质也在表2中有所详细说明.

SPH模拟岩石所需的关键参数，如密度、弹性模量、泊松比和内摩擦角，可以直接从实际实验中获取，如表2和表3所示.考虑到层理的机械参数获取十分困难，这里进行了简化，只考虑了基于带拉伸截断的摩尔库伦准则的强度参数弱化，其余机械参数与层间基质保持一致.相比之下，粘聚力和抗拉强度的值需要通过标定确定.即，通过拟合模拟的应力-应变曲线与研究岩石的单轴压缩实验的应力-应变曲线，采用试错法估计其值.它们的值表示粒子之间的结合强度，被视为微观参数，类似于DEM法中的微观参数（武东阳等，2021）.因此，它们的值通常大于岩石样本的实际值（邓树新等，2019）.

在现场尺度上对节理岩质边坡进行的SPH模拟分为两个步骤.首先通过在模型中应用真实的弹性模量、泊松比、密度和无限强度来平衡边坡模型的初始应力（竖向应力）.当竖向应力在0.2 s（即2 000 step）时逐渐达到平衡状态时，监测点的自重应力达到0.54 MPa（根据ρgh求得）.

图8展示了现场尺度边坡在裂纹扩展与渐进式破坏过程以及相应的速度分布，以及其与边坡向下滑动模式所报道结果的对比.计算始于一个小变形阶段，在此阶段模型中的应力达到平衡，弹性变形占主导地位.随后，由于应力集中，沿着边坡趾部附近的裂隙开始出现剪切裂纹.当这些裂纹穿透边坡时，形成了可以发生滑动的阶梯状破坏面，这在图7c~7e中通过FDEM、DEM和扩展有限元模拟得到了验证.通过SPH得到的阶梯状破坏模式与文献中报道的现场观察（Brideau et al.， 2009； Zhou and Chen， 2019）一致，随后边坡进入大变形阶段.值得注意的是，基于SPH颗粒模型（图8a~8b在2.0~2.5 s时）有助于完全向地面滑动.在滑动过程中，岩块撞击基底台阶，伴随着剧烈的振动，导致大量拉伸和剪切裂纹的产生，形成了大量的碎屑和小块岩石并发生飞溅.最终的岩块沉积情况与台阶路径破坏后的实地观察结果一致（Brideau et al.， 2009）.总体而言，上述模拟结果表明，SPH方法能够精确捕捉整个岩质边坡三维破坏过程，包括裂纹的初始和扩展、岩坡的摩擦滑动.与在没有接触处理的情况下进行裂纹模拟相比，可以实现裂纹扩展，但难以实现断裂面相互分离以及岩块滑动，这是由于未能处理岩块之间的摩擦接触所致.

图9比较了不同节理倾角下节理岩坡的滑动面和相应的速度场.根据强度折减理论，岩体的裂纹扩展和速度受预制倾角的影响.可以发现，与边坡表面朝向边坡面倾斜的层面相比，边坡内部的层面滑动面通常要更深一些，而与边坡表面背离倾斜的层面相比，这一点与文献Huang and Gu（2017）中观察到的滑动面情况是一致的.

图10对比了不同节理倾角下岩质边坡的裂纹分布、速度场以及最终沉积变形量.与二维情况不同，在三维岩质边坡的最终沉积中，在z方向上观察到显著的速度差异.相对于其他倾角，节理岩质边坡的崩解程度在θ =0°时最大，在θ =60°时最小.这种差异主要是由于破坏模式所致，即高倾角的节理岩质边坡容易形成局部剥落破坏，而θ=0°则容易形成整体破坏.这种岩质边坡的破坏模式直接导致了θ=0°节理岩体的边坡最大安全系数.通过安全系数对比（图11），此外θ=0°的安全系数最大，θ=-30°安全系数最小.三维节理岩质边坡的稳定性安全系数大于反倾层，这与二维情况不同，主要是因为反倾层边坡的破坏是滑动破坏.

图12和图13为裂纹数量以及边坡临坡面顶部位移随着时间的变化规律，从图中可看出以剪切裂纹在滑移破坏中起主导作用，其中θ =0°节理边坡的剪切裂纹最甚.对比不同节理倾角下的FOS，发现节理倾角θ=0°时的滑动破坏产生的位移最大，θ=60°的位移最小.

3 讨论

改进SPH方法在模拟岩石脆性破坏方面展现出独特优势.该方法采用拉格朗日框架和无网格技术，能够有效处理材料的大变形问题.通过优化设计的链表搜索算法和蛙跳算法，显著提升了相邻粒子相互作用的计算效率，这使得边坡尺度下裂纹扩展和破坏过程的模拟既快速又可靠.该方法的一个关键特点是其微观参数具有明确的物理意义，特别是可以直接关联到材料的黏聚力和抗拉强度等基本力学性质.与DEM方法相比，改进SPH方法在利用强度折减法的边坡模拟中表现出明显的效率优势.

大尺度节理边坡数值模拟的计算效率体现在计算成本上.为了比较SPH与传统方法的计算成本，选取了块状边坡滑动倾倒破坏的数值模拟进行研究.在三维和二维尺度下，比较了SPH和DEM（即PFC3D）达到相同破坏状态所需的计算时间.对于SPH以及DEM，采用多维同时比较，三维与二维模型对应的粒子数量分别为47 830和4 783.本研究中SPH运行程序的硬件环境为联想R9000P.处理器为AMD Ryzen 7-5800H，配备Radeon图形8核.基础频率为3.20 GHz，最大加速频率为4.40 GHz.在本计算中，采用单核单线程.在三维尺度下，SPH模型和DEM模型达到相同破坏状态（即相同的最大水平位移）所需的计算时间（图14a~14b）分别为约22 min和252 min.此外，在二维尺度下，SPH模型和DEM模型模拟相同破坏状态所需的计算时间（图14c~14d）分别为35.2 s和440.0 s，SPH耗时降低到DEM耗时的10%以下.鉴于此，SPH可为评估大尺度岩体非连续运动高效计算提供帮助.

SPH方法在计算架构上展现出突出的并行计算优势.其开源代码框架通过多线程优化和动态负载均衡技术，能够在单核单线程模式下保持高效运算.例如，在进行多工况敏感性分析时（如16组参数组合），通过Code-Block集成开发环境的多窗口并行模式，可将计算任务分配到多个独立线程中同步执行，理论上可实现线性加速比，即计算效率提升16倍.这一特性尤其适用于高性能计算平台，当CPU主频较高且支持超频时，SPH代码的计算效率可进一步提升.相比之下，传统DEM通常依赖于商业软件（如PFC3D）的内置并行化模块，尽管其能够充分利用多核CPU资源，但由于接触力链计算的固有复杂性（如接触检测与力链更新的高计算开销），即使在CPU利用率达到100%的情况下，其计算效率仍显著低于SPH方法.因此，在相同硬件配置条件下，SPH方法能够更高效地利用计算资源，为大规模节理边坡的数值模拟提供更强大的技术支持.这一优势使得SPH方法在工程实际应用中具有更好的可扩展性和计算效率.

相较于DEM，SPH在工程尺度岩质边坡模拟中表现出更强的工程适用性.DEM虽能精细刻画颗粒间接触力学行为，但其接触参数标定（如滑动摩擦系数、法向/切向刚度比）与岩体宏-细观力学响应的映射关系存在不确定性，参数标定过程也较为繁琐.例如，需要通过与充分约束的实验数据集进行校准，来模拟岩石的典型力学行为，通常需要一组微观和宏观参数（通常超过10个）（Kazerani et al.， 2010）.而SPH通过本征状态方程与人工粘度模型直接关联粒子应力-应变响应与岩石宏观力学参数，其参数体系与带拉伸截断的Mohr-Coulomb准则理论等岩体破坏准则兼容性更高，并且SPH模型所需的大部分微观参数由于是偏微分方程求解，都可以直接采用宏观物理真实值，无需进行微观参数到宏观参数的大量校准.在计算效率方面，SPH因无需显式求解接触力链网络，其单时间步长内的计算复杂度为O（N），显著低于DEM的O（N²）量级（Cleary et al.， 2016）.结合CPU并行计算架构，SPH可高效实现百万级粒子规模的边坡滑移全过程模拟，为边坡时效稳定性分析提供高分辨率数值工具.

在SPH中，对于较为复杂的边界，传统的粒子边界条件有两种（Monaghan， 2005），第1种类型的粒子（型号I）设置在固定边界上，第2种类型的粒子（型号II）分布在边界领域外，类型I的虚粒子被用于对内部粒子施加边界排斥力，防止内部粒子穿透边界.当类型I的虚粒子成为邻近边界处实粒子的相邻粒子时，则会在沿两粒子的中心线处对实粒子产生一个作用力.类型II虚粒子没有固定的参数，它们是在每个计算步中由对应的实粒子对称产生的，用于提高SPH在边界处的计算精度（Monaghan， 2005）.未来可以将对核函数进行进一步的改进和校正，优化复杂边界处理，以实现更为真实的滑坡扩展过程模拟.

针对不同的工况以及所用方法的差异，一般会采用不同的FOS.传统数值模拟方法中，土质边坡FOS<1.3时会被认定为是不稳定的，而岩质边坡通常FOS<1.5（Liu et al.， 2015）.而关于节理密度与FOS之间的关联性问题，会在未来的研究中进一步探索.

当然，改进SPH方法也存在一些需要改进的方面.在拉应力作用逐渐减小的情况下可能出现数值不稳定，边界条件的处理相对复杂，这在一定程度上限制了其在固体力学中的应用.相比之下，该方法在流体力学领域表现更为出色，特别是在涉及水、热、力和化学过程的多场耦合问题中.虽然SPH方法在实际工程应用中仍面临挑战，但通过改进SPH对岩质边坡破坏的模拟以及对抗拉强度折减法的推广，为该方法在工程实践中的推广积累了重要经验.

4 结论

本研究提出了一种基于强度折减与核函数改进的三维SPH算法，通过裂隙岩体单轴压缩试验模拟和裂隙岩质边坡失稳模拟，验证了该方法的准确性与适用性，并揭示了不同裂隙倾角下的破坏机制.主要结论如下：

（1）通过引入抗拉强度折减、带拉伸截断的摩尔库伦破坏准则以及改进的核函数，该算法成功模拟了三维岩体边坡从裂纹萌生、扩展到最终滑坡的全过程，无需依赖预制裂隙即可实现自然破裂演化.相较于离散元法（DEM），本方法在计算效率、精度及参数标定简便性方面具有显著优势.

（2）在单轴压缩试验中，模拟结果与物理实验一致，成功再现了翼型裂纹和交合裂纹的典型破坏模式，验证了算法对裂隙扩展机理的捕捉能力.多节理裂隙边坡的失稳模拟表明，节理倾角对安全系数具有决定性影响：当倾角为0°时，边坡呈现整体滑移破坏模式，安全系数最大；随着倾角增大，局部剥落破坏趋势增强，稳定性显著降低.这一规律为岩质边坡工程风险评估提供了理论依据.

（3）尽管本方法在三维岩体破裂模拟中表现优异，但计算效率仍可进一步提升，例如通过GPU并行计算优化大规模粒子系统的求解速度.进一步提高SPH粒子分辨率，以刻画更精细的裂纹网络及微观损伤演化过程，增强对复杂地质条件的适用性.

本研究为岩质边坡破坏机理分析及灾害预警提供了新的数值模拟手段，后续可结合实际工程监测数据进一步验证并拓展其应用场景.

参考文献

[1] Brideau, M. A., Yan, M., Stead, D., 2009. The Role of Tectonic Damage and Brittle Rock Fracture in the Development of Large Rock Slope Failures. Geomorphology, 103(1): 30-49. https://doi.org/10.1016/j.geomorph.2008.04.010

[2] Bui, H. H., Fukagawa, R., Sako, K., et al., 2008. Lagrangian Meshfree Particles Method (SPH) for Large Deformation and Failure Flows of Geomaterial Using Elastic-Plastic Soil Constitutive Model. International Journal for Numerical and Analytical Methods in Geomechanics, 32(12): 1537-1570. https://doi.org/10.1002/nag.688

[3] Cleary, P. W., Pereira, G. G., Lemiale, V., et al., 2016. Multiscale Model for Predicting Shear Zone Structure and Permeability in Deforming Rock. Computational Particle Mechanics, 3(2): 179-199. https://doi.org/10.1007/s40571-015-0073-4

[4] Dehghan, A. N., Khodaei, M., 2022. Stability Analysis and Optimal Design of Ultimate Slope of an Open Pit Mine: A Case Study. Geotechnical and Geological Engineering, 40(4): 1789-1808. https://doi.org/10.1007/s10706-021-01993-8

[5] Deng, S. X., Zheng, Y. L., Feng, L. P., et al., 2019. Application of Design of Experiments in Microscopic Parameter Calibration for Hard Rocks of PFC^3D Model. Chinese Journal of Geotechnical Engineering, 41(4): 655-664 (in Chinese with English abstract).

[6] Do, N. A., Dias, D., Dinh, V. D., et al., 2019. Behavior of Noncircular Tunnels Excavated in Stratified Rock Masses—Case of Underground Coal Mines. Journal of Rock Mechanics and Geotechnical Engineering, 11(1): 99-110. https://doi.org/10.1016/j.jrmge.2018.05.005

[7] Fakhimi, A., Lanari, M., 2014. DEM-SPH Simulation of Rock Blasting. Computers and Geotechnics, 55: 158-164. https://doi.org/10.1016/j.compgeo.2013.08.008

[8] Fukuda, D., Mohammadnejad, M., Liu, H. Y., et al., 2020. Development of a 3D Hybrid Finite-Discrete Element Simulator Based on GPGPU-Parallelized Computation for Modelling Rock Fracturing under Quasi-Static and Dynamic Loading Conditions. Rock Mechanics and Rock Engineering, 53(3): 1079-1112. https://doi.org/10.1007/s00603-019-01960-z

[9] Han, Z., Su, B., Li, Y. G., et al., 2019. Numerical Simulation of Debris-Flow Behavior Based on the SPH Method Incorporating the Herschel-Bulkley-Papanastasiou Rheology Model. Engineering Geology, 255: 26-36. https://doi.org/10.1016/j.enggeo.2019.04.013

[10] Hatzor, Y. H., Arzi, A. A., Zaslavsky, Y., et al., 2004. Dynamic Stability Analysis of Jointed Rock Slopes Using the DDA Method: King Herod’s Palace, Masada, Israel. International Journal of Rock Mechanics and Mining Sciences, 41(5): 813-832. https://doi.org/10.1016/j.ijrmms.2004.02.002

[11] Huang, D., Cen, D. F., Ma, G. W., et al., 2015. Step-Path Failure of Rock Slopes with Intermittent Joints. Landslides, 12(5): 911-926. https://doi.org/10.1007/s10346-014-0517-6

[12] Huang, D., Gu, D. M., 2017. Influence of Filling-Drawdown Cycles of the Three Gorges Reservoir on Deformation and Failure Behaviors of Anaclinal Rock Slopes in the Wu Gorge. Geomorphology, 295: 489-506. https://doi.org/10.1016/j.geomorph.2017.07.028

[13] Huang, M. S., Jia, C. Q., 2009. Strength Reduction FEM in Stability Analysis of Soil Slopes Subjected to Transient Unsaturated Seepage. Computers and Geotechnics, 36(1/2): 93-101.

[14] Jia, Q. J., Liu, X. M., Tan, X., 2024. Crack Evolution Mechanism of Stratified Rock Mass under Different Strength Ratios and Soft Layer Thickness: Insights from DEM Modeling. Soils and Foundations, 64(6): 101534. https://doi.org/10.1016/j.sandf.2024.101534

[15] Kazerani, T., Zhao, J., 2010. Micromechanical Parameters in Bonded Particle Method for Modelling of Brittle Material Failure. International Journal for Numerical and Analytical Methods in Geomechanics, 34(18): 1877-1895. https://doi.org/10.1002/nag.884

[16] Lee, H., Jeon, S., 2011. An Experimental and Numerical Study of Fracture Coalescence in Pre-Cracked Specimens under Uniaxial Compression. International Journal of Solids and Structures, 48(6): 979-999. https://doi.org/10.1016/j.ijsolstr.2010.12.001

[17] Li, L., Wang, Y., 2020. Identification of Failure Slip Surfaces for Landslide Risk Assessment Using Smoothed Particle Hydrodynamics. Georisk: Assessment and Management of Risk for Engineered Systems and Geohazards, 14(2): 91-111. https://doi.org/10.1080/17499518.2019.1602877.[LinkOut]

[18] Liu, M. B., Liu, G. R., Zong, Z., 2008. An Overview on Smoothed Particle Hydrodynamics. International Journal of Computational Methods, 5(1): 135-188. https://doi.org/10.1142/s021987620800142x

[19] Liu, M. M., Shi, Z. M., Li, B., 2024. Analysis of Dynamic Response and Failure Mode of Bedding Rock SlopesSubject to Strong Earthquakes Based on DEM⁃FDM Coupling. Earth Science, 49(08):2799-2812.

[20] Liu, S. Y., Shao, L. T., Li, H. J., 2015. Slope Stability Analysis Using the Limit Equilibrium Method and Two Finite Element Methods. Computers and Geotechnics, 63: 291-298. https://doi.org/10.1016/j.compgeo.2014.10.008

[21] Monaghan, J. J., 2005. Smoothed Particle Hydrodynamics. Reports on Progress in Physics, 68(8): 1703-1759. https://doi.org/10.1088/0034-4885/68/8/r01

[22] Shu, Q., Wang, X. B., Zhao, Y. F., et al., 2020. Numerical Simulation of Failure Processes of Heterogeneous Rock Specimens under Assumption of Invariant Spherical Stress during Stress Drop. Rock and Soil Mechanics, 41(10): 3465-3472. https://doi.org/10.16285/j.rsm.2020.0077

[23] Sun, C. W., Chai, J. R., Xu, Z. G., et al., 2016. Stability Charts for Rock Mass Slopes Based on the Hoek-Brown Strength Reduction Technique. Engineering Geology, 214: 94-106. https://doi.org/10.1016/j.enggeo.2016.09.017

[24] Wang, X. M., Xiao, Y. J., Shi, W. B., et al., 2022. Forensic Analysis and Numerical Simulation of a Catastrophic Landslide of Dissolved and Fractured Rock Slope Subject to Underground Mining. Landslides, 19(5): 1045-1067. https://doi.org/10.1007/s10346-021-01842-y

[25] Wang, Y. N., Bui, H. H., Nguyen, G. D., et al., 2019. A New SPH-Based Continuum Framework with an Embedded Fracture Process Zone for Modelling Rock Fracture. International Journal of Solids and Structures, 159: 40-57. https://doi.org/10.1016/j.ijsolstr.2018.09.019

[26] Wang, Y. T., Zhou, X. P., Kou, M. M., 2018. Peridynamic Investigation on Thermal Fracturing Behavior of Ceramic Nuclear Fuel Pellets under Power Cycles. Ceramics International, 44(10): 11512-11542. https://doi.org/10.1016/j.ceramint.2018.03.214

[27] Wang, Z. H., Zhu, Z. M., Zhou, L., et al., 2024. Dynamic Mechanical Properties and Failure Characteristics of Layered Composite Rock Containing a Tunnel-Shaped Hole. Theoretical and Applied Fracture Mechanics, 129: 104217. https://doi.org/10.1016/j.tafmec.2023.104217

[28] Wei, W. B., Cheng, Y. M., Li, L., 2009. Three-Dimensional Slope Failure Analysis by the Strength Reduction and Limit Equilibrium Methods. Computers and Geotechnics, 36(1-2): 70-80. https://doi.org/10.1016/j.compgeo.2008.03.003

[29] Wu, D., Li, H. B., Fukuda, D., et al., 2023. Development of a Finite-Discrete Element Method with Finite-Strain Elasto-Plasticity and Cohesive Zone Models for Simulating the Dynamic Fracture of Rocks. Computers and Geotechnics, 156: 105271. https://doi.org/10.1016/j.compgeo.2023.105271

[30] Wu, D. Y., Yu, L. Y., Su, H. J., et al., 2021. Experimental Study and PFC^3D Simulation on Crack Propagation of Fractured Rock-Like Specimens with Bolts under Uniaxial Compression. Rock and Soil Mechanics, 42(6): 1681-1692 (in Chinese with English abstract).

[31] Xia, C. Z., Shi, Z. M., Kou, H. J., 2024. A Modified Smoothed Particle Hydrodynamics Method Considering Residual Stress for Simulating Failure and Its Application in Layered Rock Mass. Journal of Mountain Science, 21(6): 2091-2112. https://doi.org/10.1007/s11629-023-8362-5

[32] Yang, S. Q., Huang, Y. H., 2017. An Experimental Study on Deformation and Failure Mechanical Behavior of Granite Containing a Single Fissure under Different Confining Pressures. Environmental Earth Sciences, 76(10): 364. https://doi.org/10.1007/s12665-017-6696-4

[33] Yang, Y. T., Sun, G. H., Zheng, H., et al., 2019. Investigation of the Sequential Excavation of a Soil-Rock-Mixture Slope Using the Numerical Manifold Method. Engineering Geology, 256: 93-109. https://doi.org/10.1016/j.enggeo.2019.05.005

[34] Yu, S. Y., Ren, X. H., Zhang, J. X., 2023. Using an Improved SPH Algorithm to Simulate Thermo-Hydro-Mechanical-Damage Coupling Problems in Rock Masses. Case Studies in Thermal Engineering, 47: 103085. https://doi.org/10.1016/j.csite.2023.103085

[35] Zhang, H. J., Wu, S. C., Zhang, Z. X., et al., 2023. Reliability Analysis of Rock Slopes Considering the Uncertainty of Joint Spatial Distributions. Computers and Geotechnics, 161: 105566. https://doi.org/10.1016/j.compgeo.2023.105566

[36] Zhang, W. J., Zheng, H., Jiang, F. Y., et al., 2019. Stability Analysis of Soil Slope Based on a Water-Soil-Coupled and Parallelized Smoothed Particle Hydrodynamics Model. Computers and Geotechnics, 108: 212-225. https://doi.org/10.1016/j.compgeo.2018.12.025

[37] Zhou, X. P., Chen, J. W., 2019. Extended Finite Element Simulation of Step-Path Brittle Failure in Rock Slopes with Non-Persistent En-Echelon Joints. Engineering Geology, 250: 65-88. https://doi.org/10.1016/j.enggeo.2019.01.012

[38] Zhou, Z. Q., Zhao, Y., Bi, J., et al., 2023. Shear Mechanical Properties and Failure Modes of Rock with V-Shaped Intersecting Double-Cracks. Theoretical and Applied Fracture Mechanics, 124: 103755. https://doi.org/10.1016/j.tafmec.2023.103755

基金资助

湖南省科技创新计划资助(2024AQ2044)

湖南省地质院科技计划项目(HNGSTP202516)

国家自然科学基金-博士生基础研究项目(424B2055)

深地国家科技重大专项(2024ZD1001400)