三维土石混合体边坡大变形SPH-DEM模拟

Large Deformation Analysis of 3D Soil-Rock Mixture Slopes Using SPH-DEM Method

摘要

自然界中的大多数边坡由土体与块石混合构成，一旦发生大变形破坏，可能严重威胁周围人民生命安全及基础设施建设运维.针对土石混合体边坡研究中存在的精细化建模及块石与土体耦合计算难题，建立了一种三维土石混合体边坡高保真建模技术，并提出了基于SPH-DEM耦合的三维土石混合体边坡大变形模拟方法，进一步分析了块石对边坡大变形冲击过程的影响，并预测了青海省浪加滑坡体再次滑动后冲击大坝附属建筑物动力过程，结果表明：块石含量和位置分布显著影响滑坡冲击过程，块石含量越高，冲击力时程曲线的差异性越显著，若块石与建筑物直接碰撞，峰值冲击力较不考虑块石时提高约30%.浪加滑坡再次滑动后前缘最大运动距离达108 m，启闭房、施工营地所受峰值冲击力分别高达自身重力的20倍和4倍，极易被滑坡体冲毁.该研究成果可为土石混合体滑坡的灾害预测和危险性分析提供参考依据.

Abstract

Most natural slopes are composed of soil-rock mixtures, whose large deformation and failure pose severe threats to human safety and infrastructure. To overcome the challenges in high-fidelity modeling and coupled soil-rock mixture interaction analysis, this study develops a high-fidelity 3D modeling technique for soil-rock mixture slopes and proposes an SPH-DEM coupling method to simulate the large deformations. It further analyzes the impact of boulders on the deformation and failure process of the slope and predict the variation of the impact on dam appurtenant structures after the reactivation of the Langjia landslide in Qinghai Province. The results reveal that boulders within the landslide body significantly increasing the landslide velocity, and the content and position of the boulders affect the impact process. Higher boulder content leads to a more pronounced difference in the impact force time-history curve. When boulders collide directly with buildings, the peak impact force increases by about 30% compared to scenarios without boulders. After reactivation, the maximum movement distance of the front edge of the Langjia landslide in Qinghai Province reaches 108 m, with peak impact forces on the intake and construction camp structures reaching 20 times and 4 times their own weight, respectively, making them highly vulnerable to destruction by the landslide. The findings of this study provide valuable insights for disaster prediction and risk analysis of soil-rock mixture landslides.

Graphical abstract

关键词

土石混合体边坡 / 大变形 / 精细化建模 / 冲击过程 / SPH-DEM方法 / 工程地质学.

Key words

soil-rock mixture slope / large deformation / refined modeling / impact process / SPH-DEM method / engineering geology

0 引言

土石混合体边坡作为一种典型的非均质地质体，主要由高强度块石、低强度土体及孔隙构成（Xu et al. 2016；边宏光等，2024），其显著的结构特征表现为材料非均质性、块石空间分布随机性、土石界面不连续性以及尺寸效应（孙华飞等，2014；Wang et al.， 2022；Kang and Wang， 2025）.这种特殊结构特征使得土石混合体边坡在发生大变形时表现出复杂的力学响应，导致其演进范围、冲击特征等与均质边坡有较大差别（Peng et al.， 2023；Abioga et al. 2025）.此类边坡一旦发生大变形破坏，其突发性强、破坏范围广、致灾机理复杂的特点，极易引发重大地质灾害，严重威胁周围居民生命财产安全以及工程建设运维.

数值模拟是研究土石混合体边坡失稳及变形破坏过程的重要手段，然而受限于块石-土体非连续特性，当前在实现土石混合体边坡的高保真建模与大变形破坏过程的跨尺度模拟方面仍存在挑战.一方面，土石混合体的非均质性与破坏行为高度依赖于块石的真实几何形态及其在土体基质中的空间分布，然而传统数值建模方法多采用概化模型（如球形、椭球形或规则多面体）模拟块石（张润涵等， 2023），虽简化了计算，却难以精确反映块石棱角效应、咬合作用及不规则形态对土石界面不连续性、应力集中和尺寸效应的真实影响.另一方面，土石混合体边坡从失稳到剧烈运动堆积过程中土体发生大变形，传统基于网格的方法如有限单元法、有限差分法等，其优势主要体现在边坡失稳前的小变形分析阶段（Liu et al.， 2020；Zhao et al.， 2021），在预测边坡失稳后的大变形破坏过程时，此类方法往往存在网格扭曲、畸变问题，导致计算精度下降甚至计算终止（Zhu et al.， 2022）.尽管存在一些改进措施如网格自适应技术和界面追踪算法等，但这些方法实现的复杂性与高昂的计算成本也限制了其在工程尺度应用的可行性.

近年来，无网格方法或基于点的数值模拟方法，如离散单元法（DEM）、非连续变形分析法（DDA）和光滑粒子流体动力学方法（SPH）等，在岩土工程中大变形分析方面受到越来越多关注（Hu et al.， 2025；李爽等， 2025；彭铭等， 2025）.DEM方法适用于模拟块石等离散介质的运动与接触，而SPH方法作为一种纯拉格朗日、无网格的粒子方法，在处理类流体材料极端大变形、自由表面流动及颗粒破碎方面展现出独特优势（Bui et al.， 2011；Fang et al.， 2024），且便于嵌入各类先进的本构模型.在连续-离散耦合框架下，采用结合DEM与SPH优势的SPH-DEM耦合方法，即采用SPH方法模拟土体大变形过程，采用DEM方法模拟块石的运动和相互碰撞过程，是预测土石混合体边坡大变形过程的有效手段.然而，现有SPH-DEM耦合研究侧重于模拟泥石流或类流体滑坡，多采用流体本构或流变模型描述土体行为，但对于大量以剪切破坏为主、土体仍表现出显著固体力学特性的土石混合体边坡，此类流体本构模型难以准确刻画土体在剪切过程中的塑性屈服、应变软化等关键力学行为，且现有研究较少考虑块石真实形态对土石混合体边坡大变形运动机制与冲击特性的影响.

针对土石混合体边坡研究中存在的精细化建模及块石与土体耦合计算难题，本研究建立了一种三维土石混合体边坡高保真建模技术，并提出基于SPH-DEM耦合的三维土石混合体边坡大变形模拟方法，分析了块石对边坡大变形冲击过程的影响，并以青海省浪加水库滑坡为例，预测了滑坡体再次滑动后冲击大坝附属建筑物的动力过程，以期为类似工程开展边坡失稳大变形预测及灾害风险防控提供技术支撑.

1 土石混合体精细化建模

几何轮廓是块石的关键特征之一，本研究采用三维手持扫描仪，获取真实块石表面的点云数据，并封装为闭合三维曲面，从而实现块石的高分辨率描述和精细化建模.对粒径在30~50 mm范围内的50个真实块石进行扫描，图1给出了扫描不同块石后得到的几何曲面，最终建立的真实块石数据库包括曲面数据文件、质量信息和体积信息.

选取Yang et al.（2022）采用的三维边坡模型进行建模，如图2a所示，距坡脚4 m处设有一混凝土房屋建筑物，如图3b所示.土体密度为2 000 kg/m，摩擦角为21.3°，粘聚力为10 kPa，采用毕肖普法求得土体边坡的稳定系数为1.0，表明边坡处于临界失稳状态，在降雨、地震等因素影响下极易发生滑动，其临界滑移面位置如图3所示.

由土体边坡临界滑移面位置可确定如图4a所示滑坡体，其体积为334.84 m³.以在滑坡体内部生成体积含量为10%的随机块石为例，简要介绍三维土石混合体边坡建模方法.

表1列出的18个块石总体积为33.48 m³，块石粒径区间为［1.5 m， 2.5 m］，填充到滑坡体中满足体积含量为10%的要求.以块石粒径区间中值2.0 m为直径，首先在滑坡体内部生成18个位置随机的球体，如图4b所示，然后将球体全部替换为表1所列的真实块石，即可得到如图4c所示三维土石混合体边坡几何模型.考虑球体位置、块石类型和空间位置变化，可生成多个块石含量为10%的土石混合体边坡模型开展模拟.

2 数值算法

2.1　土体本构模型

考虑边坡失稳大变形运动过程土体屈服后应变软化效应，本研究采用基于Drucker-Prager准则的弹塑性本构模型计算土体应力应变，本构模型的屈服函数y和塑性势函数g分别如下：

y={\alpha }_{\varphi }{I}_{\mathrm{1}}+\sqrt[]{J_{\mathrm{2}}}-k_c,（1）

g={\alpha }_{\psi }{I}_{\mathrm{1}}+\sqrt[]{J_{\mathrm{2}}},（2）

式中：I_{\mathrm{1}}为应力张量第1不变量，J_{\mathrm{2}}为应力张量的第2不变量，k_c、{\alpha }_{\phi }和{\alpha }_{\psi }与粘聚力、内摩擦角和膨胀角相关，定义如下：

{\alpha }_{\phi }=\frac{\mathrm{t}\mathrm{a}\mathrm{n}\phi }{\sqrt[]{\mathrm{9}+\mathrm{12}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\phi }},（3）

k_c=\frac{\mathrm{3}c}{\sqrt[]{\mathrm{9}+\mathrm{12}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\phi }},（4）

{\alpha }_{\psi }=\frac{\mathrm{t}\mathrm{a}\mathrm{n}\psi }{\sqrt[]{\mathrm{9}+\mathrm{12}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\psi }},（5）

式中：\phi为内摩擦角，c为粘聚力，\psi为膨胀角.土体屈服后发生应变软化，采用下式计算土体强度弱化：

c(\kappa )=c_{\mathrm{r}}+\left(c_{\mathrm{p}}-c_{\mathrm{r}}\right)e^{-{\eta }_c\kappa },（6）

\phi (\kappa )={\phi }_{\mathrm{r}}+\left({\phi }_{\mathrm{p}}-{\phi }_{\mathrm{r}}\right)e^{-{\eta }_{\phi }\kappa },（7）

式中：c_{\mathrm{p}}为峰值粘聚力，c_{\mathrm{r}}为残余粘聚力，{\phi }_{\mathrm{p}}为峰值内摩擦角，{\phi }_{\mathrm{r}}为残余内摩擦角，{\eta }_c和{\eta }_{\phi }为控制土体强度衰减速率的系数，\kappa与塑性剪应变张量的增量\mathrm{d}{\mathit{e}}^p有关，\mathrm{d}\kappa =\sqrt[]{\mathrm{2}/\mathrm{3}\mathrm{d}{\mathit{e}}^p:\mathrm{d}{\mathit{e}}^p}.该应变软化本构模型参数如峰值强度、残余强度等可通过三轴压缩试验获得.

2.2　SPH方法基本理论

光滑粒子流体动力学方法（smoothed particle hydrodynamics， SPH）是一种无网格粒子方法，将计算域划分为多个粒子，粒子控制方程如下：

\frac{\mathrm{d}\rho }{\mathrm{d}t}=-\rho \nabla \cdot \mathit{v},（8）

\frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}=\frac{\mathrm{1}}{\rho }\nabla \cdot \mathit{\sigma }+\mathit{f},（9）

式中：\mathit{v}和\rho分别代表速度和密度，\mathit{\sigma }和\mathit{f}分别表示应力张量和体力张量.SPH算法通过采用权重函数得到离散值f来近似函数f(x)及其梯度\nabla f(x)的值，即：

f_i={\displaystyle \sum _{j=\mathrm{1}}^n}\frac{m_j}{{\rho }_j}{f}_j{W}_{ij},（10）

\nabla f_i={\displaystyle \sum _{j=\mathrm{1}}^n}\frac{m_j}{{\rho }_j}{f}_j{\nabla }_i{W}_{ij},（11）

式中：f_i和\nabla f_i分别为x_i位置处的近似函数和函数梯度值，f_j、m_j和{\rho }_j分别为x_j位置处的函数值、质量和密度，W_{ij}是光滑核函数，本研究采用的核函数为Wendland核函数（Wendland， 1995）.SPH算法主要用于模拟土石混合体滑坡中的土体，控制方程如下：

\frac{\mathrm{d}{\rho }_i}{\mathrm{d}t}=\sum _j{m}_j\left({\mathit{v}}_i-{\mathit{v}}_j\right){\nabla }_i{W}_{ij},（12）

\frac{\mathrm{d}{\mathit{v}}_i}{\mathrm{d}t}=\sum _j{m}_j\left(\frac{{\sigma }_i}{{\rho }_i^{\mathrm{2}}}+\frac{{\sigma }_j}{{\rho }_j^{\mathrm{2}}}+{\mathrm{\Pi }}_{ij}\right){\nabla }_i{W}_{ij}+{\mathbf{g}}_i+\frac{{\mathbf{F}}_i^{r\to s}}{m_i},（13）

式中：\mathbf{g}为重力加速度，{\mathrm{\Pi }}_{ij}为确保数值稳定的人工粘性项，{\mathbf{F}}_i^{r\to s}为岩石粒子对土体粒子i的作用力，且{\mathbf{F}}_i^{r\to s}=V_i\sum ({\sigma }_i^r+{\sigma }_j^s){\nabla }_i{W}_{ij}{V}_j，V_j为粒子j的体积.

2.3　SPH-DEM耦合计算理论

在离散单元法框架内，粒子之间通过接触距离传递信息.两个粒子之间的接触力分为法向力和切向力.粒子i和j之间的法向力{\mathit{F}}_{n,ij}和切向力{\mathit{F}}_{t,ij}分别由下式给出（Canelas et al.， 2017；Su et al.， 2024）：

{\mathit{F}}_{n,ij}={\mathit{F}}_n^r+F_n^d=k_{n,ij}{\delta }_{ij}^{\mathrm{3}/\mathrm{2}}{\mathbf{I}}_{ij}-{\gamma }_{n,ij}{\delta }_{ij}^{\mathrm{1}/\mathrm{2}}{\xi }_{ij}{\mathbf{I}}_{ij},(14)

{\mathit{F}}_{t,ij}=\mathrm{m}\mathrm{i}\mathrm{n}\left({\mu }_{f,ij}{\mathit{F}}_{n,ij}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\mathrm{8}{\xi }_{ij}\right){\mathbf{I}}_{ij}^t,{\mathit{F}}_t^r+{\mathit{F}}_t^d\right),(15)

式中：k_{n,ij}=\frac{\mathrm{4}}{\mathrm{3}}{\left((\mathrm{1}-v_i^{\mathrm{2}})/E_i+(\mathrm{1}-v_j^{\mathrm{2}})/E_j\right)}^{-\mathrm{1}}\sqrt[]{{\left(\mathrm{1}/r_i+\mathrm{1}/r_j\right)}^{-\mathrm{1}}}，是与弹性模量、泊松比和粒子半径相关的法向刚度系数，{\gamma }_{n,ij}是法向阻尼常数，k_{t,ij}是切向刚度系数，{\gamma }_{n,ij}是切向阻尼系数，{\delta }_{ij}=\mathrm{m}\mathrm{a}\mathrm{x}\left(\mathrm{0},\left(d_i+d_j\right)/\mathrm{2}-\left|{\mathit{r}}_{\mathbf{i}\mathbf{j}}\right|\right)是粒子间侵入距离，用于表征粒子间碰撞程度.{\mathbf{I}}_{ij}是单位向量，{\xi }_{ij}={\mathit{v}}_{ij}\cdot {\mathbf{I}}_{ij}是法向变形率，{\mu }_{f,ij}为粒子i和粒子j之间的内摩擦系数.土石混合体滑坡中块石和沿途房屋建筑物采用DEM方法模拟，记块石或沿途建筑物为A，则其控制方程为：

M_{\mathrm{A}}\frac{\mathrm{d}{\mathit{V}}_{\mathrm{A}}}{\mathrm{d}t}=\sum _{k\in \mathrm{A}}{m}_k\frac{\mathrm{d}{\mathit{V}}_k}{\mathrm{d}\mathrm{t}}=\sum _{k\in \mathrm{A}}{\mathit{F}}_k^{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}}+\sum _{k\in \mathrm{A}}{\mathbf{F}}_k^{s\to r}+m_k\mathbf{g},(16)

I_{\mathrm{A}}\frac{\mathrm{d}{\mathbf{\Omega }}_{\mathrm{A}}}{\mathrm{d}t}=\sum _{k\in \mathrm{A}}{m}_k\left({\mathit{r}}_k-{\mathit{R}}_{\mathrm{A}}\right)\times \frac{\mathrm{d}{\mathit{V}}_k}{\mathrm{d}t},(17)

式中：M_{\mathrm{A}}是块石或建筑物A的质量，{\mathit{V}}_{\mathrm{A}}其速度，I_{\mathrm{A}}为转动惯量，{\mathrm{\Omega }}_{\mathrm{A}}为角速度，{\mathit{R}}_{\mathrm{A}}为重心，m_k\mathrm{d}{\mathit{V}}_k/\mathrm{d}t为作用在粒子k的作用力，包括重力、DEM块石或建筑物碰撞提供的外力{\mathit{F}}_k^{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}}以及滑坡土粒子的作用力{\mathbf{F}}_k^{s\to r}.DEM表面粒子的应力通过对周围土粒子应力进行插值得到，即{\mathbf{F}}_k^{s\to r}=V_k\sum ({\sigma }_k^r+{\sigma }_j^s){\nabla }_k{W}_{kj}{V}_j.块石或建筑物由DEM粒子组成，通过求解动量方程得到每个粒子的加速度，将其乘以粒子质量得到作用在该粒子上的力，进行矢量求和即为作用在房屋上的冲击力：

\mathit{F}=\sum _{k\in \mathrm{A}}{m}_k\mathrm{d}{\mathit{v}}_k/\mathrm{d}t .（18）

因此，对于三维问题，冲击力大小为F=\sqrt[]{F_x^{\mathrm{2}}+F_y^{\mathrm{2}}+F_z^{\mathrm{2}}}.考虑篇幅限制，对于SPH-DEM耦合算法可靠性验证，读者可参考Su et al.（2026）.

3 三维土石混合体边坡大变形分析

3.1　边坡模型与材料参数

为研究块石分布对边坡大变形冲击过程的影响，针对5%、10%和15%三种块石含量工况，建立土石混合体边坡冲击房屋建筑物数值模型，每种工况随机生成10次土石混合体滑坡模型开展模拟，相关参数如表2所示，其中块石密度为2 731.7kg/m³.在DEM参数设置方面，参考Chau et al.（2002）和Canelas et al.（2017）模型试验及数值分析结果，块石与房屋间的碰撞恢复系数设为0.5，动摩擦系数设为0.35，DEM粒子和SPH粒子半径均设为0.15 m，模型底部基岩和四周为固壁边界（Su et al.， 2026），模拟时长为10 s.

图5给出了含有不同块石含量的土石混合体边坡模型，块石实际含量分别为4.992%、10.000%和15.007%，与目标块石含量5%、10%和15%相比，本研究生成的块石含量误差均低于1%.

3.2　边坡大变形过程模拟

图6和图7分别展示了块石含量为0%（均质边坡）和块石含量为15%的土石混合体边坡在发生大变形过程中的速度云图，整体可知：在自重作用下，边坡沿预设滑移面滑动，短时间内即与房屋发生碰撞.受失稳坡体冲击作用，房屋发生倾倒并向右侧位移；随着运动坡体动能逐渐耗散，房屋最终停止运动.对比图6与图7可知，均质边坡发生大变形时速度场呈现高度对称性，而块石的存在会破坏这种对称性，且块石含量越高，坡体大变形运动时速度场不对称性越显著.此外，相较于均质边坡，土石混合体边坡大变形运动时的局部运动速度更快，具体表现为在相同时刻，土石混合体滑坡速度云图中局部峰值更高.图8给出了边坡t=10 s时累积塑性应变分布情况，可见，均质边坡塑性区主要集中在坡体与基岩接触面以及坡体与房屋接触部位，且剪切带相对较为规则，而土石混合体边坡运动过程中，塑性区扩展受块石阻碍作用，导致剪切带出现绕石和分叉，最终坡体形成一系列不连续且不规则的剪切带.

进一步分析块石含量和位置分布对边坡大变形冲击建筑物动力过程的影响，考虑块石位置分布随机性，对块石含量分别为5%、10%和15%的土石混合体边坡开展10次大变形运动模拟.图9给出了土石混合体边坡大变形冲击作用下房屋建筑物所受冲击力时程曲线，整体可见，土石混合体边坡大变形运动过程及冲击动力响应受块石含量、位置分布等因素综合影响.由冲击力变化趋势可知，冲击力从开始增加（t=2 s）到保持稳定状态（t=5 s），历时约为3 s，表明土石混合体边坡大变形过程发生迅速，能在短时间内对建筑物造成破坏.当块石含量较低时，土石混合体边坡大变形冲击作用下房屋建筑物所受冲击力时程曲线差异较小，如图9a、9b所示.当块石含量较高时，房屋建筑物所受冲击力时程曲线波动较为显著，尤其对于房屋建筑物所受峰值冲击力，如图9c、9d所示，这说明块石含量越高，块石位置分布随机性对冲击力时程曲线波动性的影响增大.

图10给出了房屋建筑物所受峰值冲击力的箱型图，图方框表示在不同块石含量条件下房屋建筑物所受峰值冲击力的平均值，箱子长度则反映了数据的波动范围，整体可见，随着块石含量增加，箱子长度呈现出不规则波动变化，这表明峰值冲击力除了受块石含量影响外，还显著受块石空间分布的影响.在相同块石含量条件下，不同块石分布的边坡模型计算的峰值冲击力仍存在较大差异.例如，在块石含量15%-case8工况下，峰值冲击力（图10突变值）约为0%工况的1.3倍，约为块石含量15%-case1工况的1.6倍.图11给出了块石含量15%-case8工况下三维土石混合体滑坡的运动速度变化，可见块石与房屋发生碰撞（图11），导致冲击力在t=2.5 s到t=3.0 s之间迅速增加，峰值冲击力达到4 280 kN，比不考虑块石时的峰值冲击力3 290 kN高出约30%.

4 土石混合体滑坡冲击过程分析

4.1　滑坡模型与材料参数

本节以青海省浪加水库库区滑坡为例，大坝施工期间，该滑坡体受极端降雨发生滑动，如图12所示，滑坡体最大宽度220 m，最小宽度155 m，滑坡体侵占了大坝部分坝基及坝轴线上游河道，以及原设计导流围堰基坑，严重影响工程建设进度.若遭遇极端暴雨、地震等工况，滑坡可能再次滑动，对沿途的启闭房、施工营地和砼拌合站等大坝附属建筑物造成冲击破坏.为分析浪加滑坡体再次发生大变形的冲击致灾过程，考虑块石含量差异，建立土石混合体滑坡几何模型，如图13所示.为提高计算效率，采用GPU加速技术，模拟滑坡再次滑动冲击大坝附属建筑物动力过程.模型中块石粒径范围为［15 m， 25 m］和［30 m， 50 m］，块石含量误差低于1%，块石相关数据如表3所示.

根据地勘资料，滑坡体在天然状态下的密度为2 100 kg/m³，粘聚力为15 kPa，内摩擦角为29.98°.考虑土体屈服后的应变软化效应，其残余强度参数确定为：残余粘聚力7.5 kPa，残余内摩擦角16.09°.滑坡体材料参数取值详见表4，块石及建筑物（包括启闭房、施工营地和砼拌合站）材料参数见表2.将滑坡体、地基和建筑物离散为粒子，数值模拟中DEM粒子和SPH粒子半径均设为1.75 m，滑坡总粒子数，地形边界粒子数为，模拟计算时间为30 s，模型底部基岩为固壁边界.

4.2　库区土石混合体滑坡模拟

图14与图15分别展示了块石含量为0%和15%时的滑坡运动速度分布.可见，t=5 s时滑坡体已发生显著位移，但受坝址区域平缓地形影响，滑坡体主要沿上下游较为平坦的沟道继续演进.

对于滑坡运动距离，块石含量为0%、5%、15%时，滑坡体前端最大运动距离分别为108 m、102 m、106 m.据此判断，滑坡覆盖范围均未波及117 m外的砼拌合站.对于滑坡运动速度，均质滑坡（块石含量为0%）的局部最大速度值为27 m/s，当块石含量增至5%和15%时，土石混合体滑坡的局部最大速度值分别达到38 m/s和34 m/s，相较于均质滑坡体分别提高了约40%和26%，说明块石一定程度上促进了滑坡体的运动过程.

图16展示了不同块石含量的土石混合体滑坡冲击下，启闭房和施工营地所受冲击力变化过程.t=8 s时两处建筑物开始承受滑坡冲击荷载，施工营地的冲击力随时间首先呈递增趋势，最终趋于稳定值；而启闭房由于距离滑坡体较近，其冲击力响应表现出明显的非线性特征，在t=10 s至t=20 s时段内出现双峰现象.第1个峰值对应于滑坡前缘与建筑物的初始接触碰撞，第2个峰值则源于后续运动土体的持续冲击作用.由模拟结果可知，当块石含量分别为0%、5%和15%时，启闭房所受的峰值冲击力依次为1.89×10⁹ N、1.51×10⁹ N和1.64×10⁹ N；施工营地所受的峰值冲击力则分别为5.37×10⁹ N、4.47×10⁹ N和4.29×10⁹ N.通过对比建筑物自重（启闭房9.42×10⁷ N，施工营地1.28×10⁹ N）可以发现，启闭房承受的冲击力可达其自重的16~20倍，施工营地则为3~4倍.在滑坡冲击荷载作用下，两者均有被滑坡体冲毁的风险.

5 结论

本研究基于连续-离散介质力学框架，采用SPH-DEM耦合计算方法，开展了土石混合体边坡大变形模拟分析，主要结论如下：

（1）基于三维激光扫描技术，开发了三维土石混合体边坡高保真建模技术，采用三维激光扫描设备可构建真实块石数据库，生成指定块石含量的边坡模型，块石含量误差低，实现土石混合体边坡的高保真建模.

（2）建筑物所受冲击力受块石含量、位置分布等因素综合影响，块石含量越高，其位置分布随机性导致冲击力时程曲线的差异性越显著，若块石与建筑物直接碰撞，峰值冲击力较不考虑块石时提高约30%.

（3）青海省浪加水库滑坡大变形计算结果表明，滑坡再次滑动后前缘最大运动距离达108 m，覆盖范围未波及相距117 m的砼拌合站.而启闭房、施工营地所受峰值冲击力分别高达自身重力的20倍、4倍，极易被滑坡冲毁.

（4）本研究所提方法主要考虑土体与岩石相互作用，未来可进一步拓展考虑水、土、岩多相及水、热、力、损伤多场耦合，实现复杂情景下的滑坡灾害精细化模拟预测.

参考文献

[1] Abioga, A., Rachmansyah, A., Zaika, Y., 2025. Rainfall-Induced Slope Stability on Soil-Rock Mixture. Rekayasa Sipil, 19(1): 83-96. https://doi.org/10.21776/ub.rekayasasipil.2025.019.01.11

[2] Bian, H. G., Wang, S., Ma, H. S., et al., 2024. Stability of Soil-Rock Mixture Slopes Based on Random Field Theory. Bulletin of Geological Science and Technology, 43(6): 162-170 (in Chinese with English abstract).

[3] Bui, H. H., Fukagawa, R., Sako, K., et al., 2011. Slope Stability Analysis and Discontinuous Slope Failure Simulation by Elasto-Plastic Smoothed Particle Hydrodynamics (SPH). Géotechnique, 61(7): 565-574. https://doi.org/10.1680/geot.9.p.046

[4] Canelas, R. B., Domínguez, J. M., Crespo, A. J. C., et al., 2017. Resolved Simulation of a Granular-Fluid Flow with a Coupled SPH-DCDEM Model. Journal of Hydraulic Engineering, 143(9): 0601701. https://doi.org/10.1061/(asce)hy.1943-7900.0001331

[5] Chau, K. T., Wong, R. H. C., Wu, J. J., 2002. Coefficient of Restitution and Rotational Motions of Rockfall Impacts. International Journal of Rock Mechanics and Mining Sciences, 39(1): 69-77. https://doi.org/10.1016/s1365-1609(02)00016-3

[6] Fang, H. J., Li, D. Q., Wang, S., et al., 2024. Hypoplastic Modeling Post-Failure Behavior of Sandy Slope with Porosity Spatial Variability. Computers and Geotechnics, 173: 106470. https://doi.org/10.1016/j.compgeo.2024.106470

[7] Hu, Y. Y., Lu, Y., Zheng, Y. W., 2025. Numerical Study on Seepage-Induced Instability of Soil-Rock Mixture Slopes Using CFD-DEM Coupling Method. Computers and Geotechnics, 183: 107206. https://doi.org/10.1016/j.compgeo.2025.107206

[8] Kang, X., Wang, S., et al., 2025. Mechanical Behaviour and Shear Localisation of Gravel Soils: Experimental Study and Numerical Modelling. Géotechnique,76(1):56-71. https://doi.org/10.1680/jgeot.24.01238

[9] Li, S., Peng, M., Shi, Z. M., et al., 2025. Simulation and Analysis of Cascading Hazard Based on Fluid-Soil Coupled SPH Method. Earth Science, 50(10): 3967-3981 (in Chinese with English abstract).

[10] Liu, S. Q., Wang, H. L., Xu, W. Y., et al., 2020. Numerical Investigation of the Influence of Rock Characteristics on the Soil-Rock Mixture (SRM) Slopes Stability. KSCE Journal of Civil Engineering, 24(11): 3247-3256. https://doi.org/10.1007/s12205-020-0034-1

[11] Peng, M., Zhao, Q. X., Li, S., et al., 2025. Two-Phase SPH Simulation of Granular Landslide-Tsunamis Processes Considering Dynamic Seepage. Earth Science, 50(10): 3795-3808 (in Chinese with English abstract).

[12] Peng, X. Y., Yu, P. C., Zhu, H., et al., 2023. Proposal of a Coupled DDA-SPH Method Incorporating a New Contact Algorithm for Soil-Structure Interaction Simulations in Geotechnical Hazards. Computers and Geotechnics, 164: 105849. https://doi.org/10.1016/j.compgeo.2023.105849

[13] Su, Z. Y., Kang, X., Ding, X. C., et al.,2026. SPH-DEM Modeling of Rainfall-Induced Slope Failure in Partially Saturated Soil-Rock Mixture. Computers and Geotechnics, 189: 107635.https://doi.org/10.1016/j.compgeo.2025.107635

[14] Su, Z. Y., Xiang, Y., Li, D. Q., et al., 2024. SPH-DEM Modeling of Cable-Controlled ROVs: Underwater Mobility and Path Planning. Ocean Engineering, 292: 1166. https://doi.org/10.1016/j.oceaneng.20.116623

[15] Sun, H. F., Ju, Y., Wang, X. F., et al., 2014. Review of the Study on Deformation, Failure and the Mesomechanisms of Rock-Soil Mixture(RSM). Scientia Sinica (Technologica), 44(2): 172-181 (in Chinese).

[16] Wang, S., Wu, W., Cui, D. S., 2022. On Mechanical Behaviour of Clastic Soils: Numerical Simulations and Constitutive Modelling. Géotechnique, 72(8): 706-721. https://doi.org/10.1680/jgeot.20.p.184

[17] Wendland, H., 1995. Piecewise Polynomial, Positive Definite and Compactly Supported Radial Functions of Minimal Degree. Advances in Computational Mathematics, 4(1): 389-396. https://doi.org/10.1007/BF02123482

[18] Xu, W. J., Hu, L. M., Gao, W., 2016. Random Generation of the Meso-Structure of a Soil-Rock Mixture and Its Application in the Study of the Mechanical Behavior in a Landslide Dam. International Journal of Rock Mechanics and Mining Sciences, 86: 166-178. https://doi.org/10.1016/j.ijrmms.2016.04.007

[19] Yang, Y. T., Liu, F., Wu, W. A., 2022. Assessing Slope Stability with an Improved 3D Numerical Manifold Method. Rock Mechanics and Rock Engineering, 55(10): 6409-6423. https://doi.org/10.1007/s00603-022-02993-7

[20] Zhang, R. H., Zhang, L. K., Cui, X. C., 2023. Triaxial Numerical Simulation and Macro-Micro Deformation Analysis of Soil-Rock Mixture. Water Resources and Power, 41(7): 171-174, 31 (in Chinese with English abstract).

[21] Zhao, L. H., Huang, D. L., Zhang, S. H., et al., 2021. A New Method for Constructing Finite Difference Model of Soil-Rock Mixture Slope and Its Stability Analysis. International Journal of Rock Mechanics and Mining Sciences, 138: 104605. https://doi.org/10.1016/j.ijrmms.2020.104605

[22] Zhu, C. W., Peng, C., Wu, W., 2022. Lagrangian Meshfree Particle Method (SPH) Based Simulation for Granular Flow in a Rotating Drum with Regularized μ(I) Elastoplastic Model. Powder Technology, 408: 117699. https://doi.org/10.1016/j.powtec.20.117699

基金资助

国家自然科学基金青年学生基础研究项目(523B2091)

国家自然科学基金重点项目(52439007)

国家自然科学基金面上项目(42472355)