考虑动态渗流的散粒体滑坡-涌浪过程两相SPH模拟

Two-Phase SPH Simulation of Granular Landslide-Tsunamis Processes Considering Dynamic Seepage

摘要

滑坡-涌浪是一种典型的多灾种耦合系统，具有跨介质灾种转化的复杂效应.基于黎曼光滑粒子流体动力学（Riemann-SPH），构建了考虑动态渗流的两相SPH滑坡-涌浪分析模型并验证了其准确性：动态渗流作用的引入使得散粒体滑坡-涌浪过程中的动量交换机制更加完整，最大波浪幅值a_m和最大波浪高度H_m的误差分别降低24.72%和41.95%以上.研究发现，滑动面倾角α与滑坡前缘倾角β对涌浪具有协同调控作用：随着α增大，a_m和H_m均呈现先增后减的单峰变化趋势；β的影响则呈现出分段特征：当α+β<90°时，a_m和H_m随角度和增大显著增长，超过该阈值后表现出非单调性变化，表明存在滑坡体积增加与有效作用面积缩减的竞争机制.此外，α的增大强化了渗流、紊流与摩擦等耗散效应，加剧涌浪能量衰减.相关成果可为滑坡-涌浪灾害防治提供科学支撑.

Abstract

The landslide-tsunami is a typical multi-hazard coupled system, characterized by complex effects resulting from the transmedia transformation of hazards. In this paper it proposes a two-phase Riemann-SPH model for landslide-tsunami simulation that incorporates dynamic seepage and is validated against laboratory experiments. The incorporation of dynamic seepage effects enhances the completeness of the momentum exchange mechanism in the granular landslide-tsunami process, reducing the errors in the maximum wave amplitude (a_m) and maximum wave height (H_m) by at least 24.72% and 41.95%, respectively. The results reveal a synergistic regulation of tsunami characteristics by the sliding surface inclination (α) and the landslide leading edge inclination (β): as α increases, the a_m and H_m exhibit a single-peaked, nonlinear increase-then-decrease trend. The influence of β shows a distinct piecewise pattern: when α+β<90°, both a_m and H_m increase significantly with the angle. Beyond this threshold, non-monotonic variations appear, reflecting a competition between the increasing landslide volume and the decreasing effective impact area. Moreover, increasing α enhances seepage, turbulent and frictional dissipation effects, accelerating energy decay. These findings provide scientific support for the mitigation of landslide-tsunami hazards.

Graphical abstract

关键词

滑坡-涌浪灾害链 / 黎曼光滑粒子流体动力学 / 动态渗流 / 滑动面倾角 / 滑坡前缘倾角 / 水土耦合 / 工程地质学.

Key words

landslide-tsunami hazard chain / Riemann-SPH / dynamic seepage / sliding surface inclination / landslide leading edge inclination / fluid-soil coupling / engineering geology

0 引言

滑坡-涌浪是指一定体积的滑坡体在重力作用下高速冲击水体并激发波浪的过程，是近水滑坡灾害链中的重要致灾方式（Evers et al.， 2019； 徐文杰， 2023；刘继芝娴，2023）.随着大规模水电工程的建设与运行，滑坡-涌浪灾害对库区及近水区域的工程稳定性和生命财产安全构成严重威胁，极易引发灾难性后果（李秋旺等， 2024）.如1963年意大利瓦伊昂水库滑坡-涌浪事件，约2.7×10⁸ m³岩土体冲击入水，引发5.0×10⁷ m³水体溢出，造成约2 000人死亡（Barla and Paronuzzi， 2013）；2013年溪洛渡库区黄坪滑坡-涌浪事件导致12人失踪，3人受伤（蒋权， 2019）；2015年三峡库区红岩子滑坡-涌浪事件致使13艘渔船倾覆，2人遇难（Huang et al.， 2016）.相较于单一滑坡灾害，滑坡-涌浪灾害具有显著的时空延展性特征，其灾害演化过程涉及跨尺度灾种转化与跨介质多场耦合的复杂物理机制（李宏伟等，2024）.这种多物理场耦合效应导致灾害影响范围扩大、致灾强度增强、演化过程更加复杂，显著增加了灾害预测与防控的难度（朱艺帆和安超，2024）.

目前关于滑坡-涌浪的研究主要可以分为现场调查、经验公式及解析法、模型试验和数值模拟（Evers and Hager， 2016； Meng et al.， 2023）.在现场调查方面：其信息主要来源于时序监测数据、目击者描述和影像资料等（Fritz，2001； Grilli et al.， 2019）.然而，由于滑坡-涌浪灾害具有突发性强、发生地点难以预测等特点，往往导致灾害现场的监测数据不完整或无法获取.

在经验公式及解析法方面：采用概化的几何模型试验，考虑多个实验变量进行大量实验，分析对比影响因子对滑坡涌浪特征的影响程度和规律，具有概念清楚、适用范围广的特点（Fritz et al.， 2004；殷坤龙等， 2012）.但是由于滑坡-涌浪是一个极其复杂的动力学过程，经验公式法忽略了滑坡体的复杂形态和运动特征，未考虑流固耦合作用，不能重现滑坡-涌浪全过程.

在滑坡-涌浪模型试验方面：主要包括在狭窄水槽中开展的二维试验（Fritz， 2001）和在宽阔水池中开展的三维试验（Heller and Spinneken， 2015）.研究表明，滑坡体质量、几何尺寸（长度、宽度、厚度）、冲击速度以及静水深度等参数对涌浪波形特征具有显著影响，但对滑坡前缘倾角的研究较为有限（Heller and Spinneken， 2013，2015； Tang et al.， 2018）.此外，滑动面倾角（入水冲击角度）对涌浪的影响规律尚存在争议.部分学者认为涌浪首浪峰值和最大波浪幅值与滑动面倾角呈正相关，而另一些学者则认为呈负相关或受复杂因素影响（Cui and Zhu， 2011； Mohammed and Fritz， 2012）.尽管模型试验可通过可靠的数据揭示涌浪的产生机制与传播规律，但受限于现有监测技术，难以实现流场动态演化过程的全域捕捉，从而制约了对水土相互作用细观机理的深入解析.

在滑坡-涌浪数值模拟方面：常用的数值模型主要包括深度平均模型、计算流体动力学（CFD）网格模型、多相流模型和粒子法模型四类（Lee et al.， 2022）.深度平均模型以其良好的计算效率得到广泛应用，但其较难考虑流体沿深度方向变化，难以模拟滑坡-涌浪过程中复杂的动力流固耦合作用（Fornaciai et al.， 2019）.CFD模型直接计算深度方向的速度变化并追踪液体表面，可以更加精细地模拟滑坡涌浪过程，但面对涌浪这类复杂的自由液面时，需要在界面处布置高精度网格，影响计算效率（Clous and Abadie，2019）.多相流模型是在CFD基础上改进的欧拉网格多相耦合模型，适用于滑坡-涌浪问题的精细化研究及物理力学机制分析，但其应用仍受限于计算效率（Yu and Lee， 2019； Lee and Huang， 2022）.

粒子法通过一系列粒子描述模拟对象，采用纯拉格朗日形式求解.在处理自由液面时，无需额外方法追踪液面，且液体与固体采用统一形式进行描述与求解，适用于流固耦合计算（Paquier et al.， 2021）.光滑粒子动力学（SPH）是一种典型的无网格拉格朗日方法，广泛用于滑坡-涌浪等大变形过程模拟（Luo et al.， 2021；徐文杰， 2023）.在该方法中，计算域被离散为携带质量、速度等物理属性的粒子单元，粒子间通过核函数进行加权插值，从而近似求解流体动力学控制方程.蒋权（2019）基于SPH两相数值模型，模拟了2013年黄坪滑坡-涌浪事件，获得了较为准确的浪高及完整的涌浪过程.Huang et al.（2023）采用SPH流固耦合框架，模拟了荒田滑坡-涌浪事件，结果与现场数据较为一致.Mao and Guan（2023）发现滑坡潜在能量的减少与最大波高之间存在指数关系.刘继芝娴（2023）基于SPH方法系统研究了巫峡段高陡库岸滑坡-涌浪的形成机理及其动态衰减规律.Wu et al.（2024a）采用MPM-SPH耦合方法模拟了白鹤滩库区王家山滑坡引发的涌浪过程，并开展了相关风险分析.Wu et al.（2024b）基于一个固定前缘倾角的刚体滑块研究了滑坡-涌浪的密度效应，发现滑坡体密度对涌浪特性的影响与滑坡厚度和冲击速度相关.Dai et al.（2025）开发了适用于海底滑坡的多相流SPH模型，研究了流变参数对滑坡运移和涌浪行为的影响.现有的SPH方法通常主要关注滑坡冲击水体过程中的动量与能量交换，而较少考虑水土间的动态渗流效应（Rauter et al.， 2022； Zhu et al.， 2022）.

动态渗流效应是指在滑坡-涌浪过程中，尤其是散粒体滑坡入水阶段，水体以非稳态形式渗入滑坡体内部，并因液相与固相间的相对运动产生渗流拖曳力.该效应显著改变滑坡-涌浪的近场动力学行为，是准确预测涌浪特征参数的关键因素（徐文杰，2023）.此外，已有研究多聚焦于滑动面倾角或滑坡前缘倾角的单一影响（刘继芝娴，2023）.鉴于两者均可直接影响滑坡的有效作用面积，进而影响涌浪生成与演化，其协同作用机制仍有待深入探讨.

为了精确预测散粒体滑坡-涌浪灾害过程中的涌浪参数，本研究基于SPH框架建立了一种考虑动态渗流的水土两相SPH模型，通过与经典试验对比验证了模型的可靠性，并探讨了渗流效应对散粒体滑坡-涌浪数值模拟计算精度的影响.此外，本研究设计了一系列数值试验，系统分析了滑动面倾角（α）与滑坡前缘倾角（β）的协同作用对涌浪特征参数的影响.

1 考虑动态渗流的SPH滑坡-涌浪分析方法

本文采用弱可压缩流体模拟水体流动，并基于Drucker-Prager本构的弹塑性体模拟滑坡体的失稳大变形行为.水体与滑坡体分别由可相互重叠的SPH粒子表示，并结合混合物理论与非线性渗流力改写控制方程，建立了一种考虑动态渗流效应的两相SPH方法，可用于散粒体滑坡-涌浪的分析研究.本文所有开发与计算均基于开源C++代码库SPHinXsys（Zhang et al.， 2021）.

1.1　流体本构方程

在该模型中，水体被视为牛顿弱可压缩流体，状态方程被用来定义压力与密度的对应关系：

p=c_{\mathrm{0}}^{\mathrm{2}}\left(\rho -{\rho }_{\mathrm{0}}\right)，(1)

式中：c₀为数值声速，ρ₀为初始密度.取c₀=10U_max，以保证密度波动小于1%，其中U_max为最大流速.水体剪切应力{\tau }_f的定义为：

{\mathit{\tau }}_f=\mathrm{2}\mu {\dot{\epsilon }}_f=\eta \left[\nabla u_f+{\left(\nabla u_f\right)}^T\right]，(2)

式中：μ为水的动态粘度，{\dot{\mathit{\epsilon }}}_f为应变率张量，η为水的动力粘度.

1.2　土体本构方程

对于土体，应用基于Drucker-Prager屈服准则的弹塑性本构.该模型的屈服条件为：

f\left(I_{\mathrm{1}},J_{\mathrm{2}}\right)=\sqrt[]{J_{\mathrm{2}}}+{\alpha }_{\phi }{I}_{\mathrm{1}}-k_{\mathrm{c}}，(3)

式中：I₁为第一应力不变量，J₂为第二偏应力张量， I 为单位矩阵，α_φ 、k_c为DP准则参数，由粘聚力c和内摩擦角φ决定：

{\alpha }_{\phi }=\frac{\mathrm{t}\mathrm{a}\mathrm{n}\phi }{\sqrt[]{\mathrm{9}+\mathrm{12}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\phi }}，(4)

k_c=\frac{\mathrm{3}c}{\sqrt[]{\mathrm{9}+\mathrm{12}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\phi }}.(5)

非关联流动法则和Jaumann应变率被应用于计算屈服应力，此时应力应变关系为：

\begin{array}{c}\dot{\mathit{\sigma }}\end{array}=\mathrm{2}G{\dot{\mathit{\epsilon }}}^s+Ktr\left(\dot{\mathit{\epsilon }}\right)\mathit{I}-\left[\dot{\lambda }\mathrm{3}{\alpha }_{\psi }K\mathit{I}+\frac{G}{\sqrt[]{J_{\mathrm{2}}}}{\mathit{\sigma }}^s\right]+\mathit{\sigma }•{\dot{\mathit{\omega }}}^T+\dot{\mathit{\omega }}•\mathit{\sigma }，(6)

式中：K和G分别为体积模量和剪切模量，\dot{\mathit{\epsilon }}是应变张量，{\dot{\mathit{\epsilon }}}^s=\dot{\mathit{\epsilon }}-tr\left(\dot{\mathit{\epsilon }}\right)\mathit{I}/\mathrm{3}是偏应变率张量，\dot{\mathit{\omega }}是旋转张量，λ是塑性应变乘子：

\dot{\lambda }=\frac{\mathrm{3}{\alpha }_{\phi }Ktr(\dot{\mathit{\epsilon }})+(G/\sqrt[]{J_{\mathrm{2}}}){\mathit{\sigma }}^s:\dot{\mathit{\epsilon }}}{\mathrm{9}{\alpha }_{\phi }K{\alpha }_{\psi }+G}，(7)

式中：ψ为土体膨胀角.

在应力更新过程中，应力状态有可能会超过DP准则的屈服面，违背理想弹塑性假设.为此，应力拉回算法被采用，具体为：

，(8)

式中：\tilde{\mathbf{\sigma }}是拉回后的应力张量.

1.3　基于混合物理论的控制方程

基于混合物理论模拟水、土混合，用下标f和s分别代表水体和土体.在混合物理论中，数值空间被内在密度为ρ_f的水及内在密度为ρ_s的土体同时占据，其各自的体积分数为\varphi_f和\varphi_s（Zhu et al.， 2022）.此时，其体积分数满足归一化条件：

{\varphi }_{\mathrm{f}}+{\varphi }_{\mathrm{s}}=\mathrm{1} .(9)

在拉格朗日框架下，两相质量守恒方程为：

\frac{{\mathrm{d}}^{\mathrm{s}}{\rho }_{\mathrm{s}}}{\mathrm{d}t}=-{\rho }_{\mathrm{s}}\nabla \bullet {\mathit{u}}_{\mathrm{s}},(10)

\frac{{\mathrm{d}}^{\mathrm{f}}{\rho }_{\mathrm{f}}}{\mathrm{d}t}=-\frac{{\rho }_{\mathrm{f}}}{{\varphi }_{\mathrm{f}}}\nabla \bullet \left({\varphi }_{\mathrm{f}}{\mathit{u}}_{\mathrm{f}}+{\varphi }_{\mathrm{s}}{\mathit{u}}_{\mathrm{s}}\right),(11)

式中：{\rho }_{\mathrm{s}}和{\rho }_{\mathrm{f}}分别为土体和流体的材料密度.u_s、u_f分别为土体和水体的速度.基于线性动量守恒，水体和土体的动量守恒方程分别为：

{\rho }_{\mathrm{s}}\frac{{\mathrm{d}}^{\mathrm{s}}{\mathit{u}}_{\mathrm{s}}}{\mathrm{d}t}=\nabla \bullet {\mathit{\sigma }}_{\mathrm{s}}+{\rho }_{\mathrm{s}}\mathrm{g}-{\varphi }_{\mathrm{s}}\nabla p_{\mathrm{f}}+\mathit{R}，(12)

{\varphi }_{\mathrm{f}}{\rho }_{\mathrm{f}}\frac{{\mathrm{d}}^{\mathrm{f}}{\mathit{u}}_{\mathrm{f}}}{\mathrm{d}t}=-{\varphi }_{\mathrm{f}}\nabla p_{\mathrm{f}}+\eta {\nabla }^{\mathrm{2}}{\mathit{\mu }}_{\mathrm{f}}+{\varphi }_{\mathrm{f}}{\rho }_{\mathrm{f}}\mathrm{g}-\mathit{R}，(13)

式中：g=9.8 m/s²为重力加速度，{\mathit{\sigma }}_{\mathrm{s}}为土体内部有效应力，p_{\mathrm{f}}为水体压力，\eta为动力粘度， R 为渗流拖曳力，其公式为：

\mathit{R}=\frac{\mu }{k}\left({\mathit{\mu }}_{\mathrm{f}}-{\mathit{\mu }}_{\mathrm{s}}\right)+\frac{\mathrm{1.75}}{\sqrt[]{\mathrm{150}}}\frac{{\rho }_{\mathrm{f}}}{\sqrt[]{k}{\varphi }_{\mathrm{f}}}‖{\mathit{\mu }}_{\mathrm{f}}-{\mathit{\mu }}_{\mathrm{s}}‖\left({\mathit{\mu }}_{\mathrm{f}}-{\mathit{\mu }}_{\mathrm{s}}\right)，(14)

k=\frac{{\varphi }_{\mathrm{f}}^{\mathrm{3}}{D}_{\mathrm{50}}^{\mathrm{2}}}{\mathrm{150}{\left(\mathrm{1}-{\varphi }_{\mathrm{f}}\right)}^{\mathrm{2}}}，(15)

式中：D₅₀为中值粒径.

在SPH模型中，用一系列拉格朗日粒子离散化计算空间.在多相SPH模型中，水和土将会被离散为两组粒子同时计算，相互作用.两种粒子可以重叠，形成两层SPH粒子，每个粒子除携带其物理力学属性如速度、质量等外，还会携带孔隙率属性（图1）.采用SPH核函数W_{ij}=W(r_i-r_j,h_s)可以对物理方程进行离散化，从而对每个SPH粒子进行计算，其中i， j为粒子的序号，h_s为光滑半径.对流体的质量守恒和动量守恒方程进行离散化后，可以得到流体的SPH形式控制方程为：

\begin{array}{c}\begin{array}{c}\frac{{\mathrm{d}}^{\mathrm{f}}{\rho }_i}{\mathrm{d}t}=\end{array}\end{array}-\frac{{\rho }_i}{{\varphi }_i}{\sum }_j^{}\left({\varphi }_j{\mathit{u}}_j-{\varphi }_i{\mathit{u}}_i\right)\bullet {\nabla }_i{W}_{ij}{V}_j-\frac{{\rho }_i}{{\varphi }_i}{\sum }_j^{}\left(\mathrm{1}-{\varphi }_j\right){\mathit{u}}_j^{\mathrm{s}}\bullet {\nabla }_i{W}_{ij}{V}_j+\widehat{\rho }{\nabla }_i{W}_{ij}V，(16)

\frac{{\mathrm{d}}^{\mathrm{f}}{\mathit{u}}_i}{\mathrm{d}t}=-\frac{\mathrm{1}}{{\rho }_i}{\sum }_j^{}\left(p_i+p_j\right){\nabla }_i{W}_{ij}{V}_j+{\mathit{g}}_i-\frac{\mathrm{1}}{{\rho }_i{\varphi }_i}{\sum }_a^{}{f}_{ia}{W}_{ia}{V}_a+\widehat{P}{\nabla }_i{W}_{ij}V，(17)

式中：i，j表示水粒子的序号，a，b表示土粒子的序号，V为粒子体积.基于黎曼相互作用求解（许强和董秀军，2011），\widehat{\rho }和\widehat{P}分别为：

\widehat{\rho }={\rho }_{\mathrm{L}}\frac{P_{\mathrm{L}}-P_{\mathrm{R}}}{{\rho }_{\mathrm{L}}{c}_{\mathrm{L}}+{\rho }_{\mathrm{R}}{c}_{\mathrm{R}}}，(18)

\widehat{P}=\frac{\mathrm{2}{\rho }_{\mathrm{L}}{c}_{\mathrm{L}}{\rho }_{\mathrm{R}}{c}_{\mathrm{R}}\left(U_{\mathrm{L}}-U_{\mathrm{R}}\right)\lambda }{{\rho }_{\mathrm{L}}\left({\rho }_{\mathrm{L}}{c}_{\mathrm{L}}+{\rho }_{\mathrm{R}}{c}_{\mathrm{R}}\right)}，(19)

\lambda =\mathrm{m}\mathrm{i}\mathrm{n}\left(\xi \mathrm{m}\mathrm{a}\mathrm{x}\left(\frac{\left(U_{\mathrm{L}}-U_{\mathrm{R}}\right)\left({\rho }_{\mathrm{L}}+{\rho }_{\mathrm{R}}\right)}{{\rho }_{\mathrm{L}}{c}_{\mathrm{L}}+{\rho }_{\mathrm{R}}{c}_{\mathrm{R}}},\mathrm{0}\right),\mathrm{1}\right)，(20)

式中：\xi=3为耗散系数.下标L，R分别代表黎曼相互作用中的左右项，计算方式为：

\left\{\begin{array}{c}\left({\rho }_{\mathrm{L}},U_{\mathrm{L}},P_{\mathrm{L}},c_{\mathrm{L}}\right)=\left({\rho }_i,-{\mathbf{v}}_i\bullet {\mathbf{e}}_{ij},p_i,c_i\right)\\ \left({\rho }_{\mathrm{R}},U_{\mathrm{R}},P_{\mathrm{R}},c_{\mathrm{R}}\right)=\left({\rho }_j,-{\mathbf{v}}_j\bullet {\mathbf{e}}_{ij},p_j,c_j\right)\end{array}\right.，(21)

式中：c为声速.采用同样的方式，可以对土体的质量守恒方程和能量守恒方程进行黎曼SPH离散化，得到：

\frac{{\mathrm{d}}^{\mathrm{s}}{\rho }_{\mathrm{a}}}{\mathrm{d}t}={\rho }_{\mathrm{a}}{\sum }_{\mathrm{b}}^{}\left({\mathit{u}}_{\mathrm{a}}-{\mathit{u}}_{\mathrm{b}}\right){\nabla }_{\mathrm{a}}{W}_{\mathrm{a}\mathrm{b}}{V}_{\mathrm{b}}+\widehat{\rho }{\nabla }_i{W}_{ij}V，(22)

\frac{{\mathrm{d}}^{\mathrm{s}}{\mathit{u}}_{\mathrm{a}}}{\mathrm{d}t}=\frac{\mathrm{1}}{{\rho }_{\mathrm{a}}}{\sum }_{\mathrm{b}}^{}\left({\sigma }_{\mathrm{a}}+{\sigma }_{\mathrm{b}}\right){\nabla }_{\mathrm{a}}{W}_{\mathrm{a}\mathrm{b}}{V}_{\mathrm{b}}+g_{\mathrm{a}}-\frac{{\varphi }_{\mathrm{a}}}{{\rho }_{\mathrm{a}}}{\sum }_i^{}\left(p_i-p_{\mathrm{a}}\right){\nabla }_{\mathrm{a}}{W}_{\mathrm{a}i}{V}_i+\frac{\mathrm{1}}{{\rho }_{\mathrm{a}}}{\sum }_i^{}{\mathbf{f}}_{\mathrm{a}i}{W}_{\mathrm{a}i}{V}_i+\widehat{P}{\nabla }_i{W}_{ij}V，(23)

式中：\widehat{\rho }和\widehat{P}计算方式与水体相同，但在土体中耗散系数ζ取15D， 其中D代表维度，在二维算例下取2，三维算例下取3（Zhang et al.， 2024）.

2 数值验证

2.1　计算模型

基于滑坡涌浪室内模型试验（Viroulet et al.， 2013）对SPH模拟滑坡-涌浪进行了验证，如图2所示.试验中，滑坡体由质量为2 kg的干颗粒组成，初始置于45°的斜坡上，并在挡板限制下保持静止.滑坡体下设置深度为0.15 m的静止水体.涌浪的波面数据通过距离滑坡体位置分别为0.45 m、0.75 m、1.05 m和1.35 m四个浪高仪（WG）进行采集（图2）.SPH模型中散粒体滑坡及水体的详细物理参数列于表1.

为了深入探讨渗流效应在滑坡-涌浪灾害中对涌浪特征参数的影响，本文分别对考虑渗流和未考虑渗流两种情形进行了数值计算.

2.2　计算结果

图3a~3e展示了滑坡涌浪试验的动态演化过程，t=0 s时，挡板被快速撤掉，滑坡体开始失稳并发生变形；t=0.3 s时，滑坡体前缘进入水中，与水体发生剧烈混合和冲击，形成近场涌浪；t=0.9 s时，滑坡体完全落入水中，水体最左端受滑坡抬升效应影响产生二次涌浪.随着滑坡体逐渐趋于稳定，所产生的波浪向远处传播并逐渐衰减，如t=1.5 s时所示.考虑渗流效应的SPH模拟结果（图3f~3j）与试验观测基本一致，成功复现了散粒体滑坡从失稳启动到涌浪传播的全过程动力学行为.相比之下，忽略渗流效应的模拟（图3k~3o）则与试验结果存在显著偏差.具体而言，考虑渗流效应时，滑坡体前缘呈现明显的团聚与堆积特征，运移形态表现为前缘厚后缘薄（图3g），与试验观测一致；而忽略渗流效应时，滑坡前缘形态明显偏离实际.此外，渗流作用下在滑坡体附近形成局部高速流动区，波形特征与试验结果高度一致；而未考虑渗流效应时，波形形态与试验结果存在显著差异.这一对比表明，本文建立的考虑渗流效应的SPH方法能够准确再现滑坡-涌浪事件的动态演化过程，包括滑坡体失稳、入水冲击、波浪生成与传播等关键阶段.

图4进一步对比了数值模拟与模型试验的波形数据，结果表明，考虑渗流效应的SPH计算方法与试验结果具有较高的一致性，而未考虑渗流效应模拟到的浪高显著偏小.具体而言，对于最大波浪幅值（a_m），未考虑渗流效应时，4个监测点的误差分别为48.15%、50.78%、53.14%和60.66%；而在考虑渗流效应后，误差显著降低至23.43%、14.77%、8.53%和1.68%.对于最大波浪高度（H_m），未考虑渗流效应时，4个监测点的误差分别为59.57%、61.99%、67.12%和71.62%；而在考虑渗流效应后，误差分别降至17.62%、10.18%、4.86%和2.69%.这一对比结果充分验证了渗流效应对滑坡-涌浪模拟精度的重要影响.考虑动态渗流的两相SPH方法引入了与水土相对速度相关的渗流拖曳力项，增强了滑坡-涌浪过程中的水土动量耦合，显著提高了涌浪波高的模拟精度，使其更接近真实试验结果.综上结果表明，考虑渗流的两相SPH算法能够更真实地反映土-水两相介质的相互作用机制，可以更好地应用于散粒体滑坡-涌浪的模拟研究.

3 滑动面倾角和滑坡前缘倾角协同作用机制

3.1　工况设置

为研究滑动面倾角（α）和滑坡前缘倾角（β）对涌浪特征参数的协同影响，采用概化模型进行数值计算，模型几何示意图如图5所示（Wu et al.， 2024b）.滑坡体形状设定为直角梯形，底边长为0.42 m，高为0.08 m，初始水位为0.3 m.设定滑动面倾角分别为35°、45°、50°、55°和65°，滑坡前缘倾角分别为20°、30°、40°、50°和60°（表2）.滑坡体由中值粒径为4 mm的散粒体组成，内摩擦角23.3°，密度为2 000 kg/m³.散粒体滑坡在重力作用下沿滑动面滑入水中，模拟过程中在布置了4个浪高仪（WG）以记录涌浪近场波面的变化.通过系统改变α和β的组合，研究α和β协调影响下涌浪特征参数的差异.

3.2　涌浪形成规律

当滑坡前缘倾角β为20°时，滑动面倾角α由35°逐步增至65°，a_m和H_m均表现出先增后减的非线性变化趋势，并在α=55°时达到峰值（图6a）.当β进一步增大至30°、40°、50°和60°时，a_m和H_m依然表现出类似的变化规律，峰值位置集中在α=50°或55°.相比之下，β对涌浪参数的影响相对较弱，但其变化规律更为复杂（图7）.具体而言，当α=35°时，随着β的增大，a_m和H_m呈现先增大后减小的单峰特征；当α=45°时，随着β的增大，a_m和H_m先增大后趋于稳定；而当α=50°、55°或65°时，随着β的增大，a_m和H_m则表现出先增大后减小再增大的多阶段变化特征.

滑坡与水体之间的有效作用面积是传递滑坡动能和动量、进而生成涌浪的关键因素（Li et al.， 2024）.滑坡的厚度（S）本身并不直接影响有效作用面积，关键因素在于其在垂直方向的分量（S_z）.当滑坡前缘倾角为90°时，S_z可用常见的S\times \mathrm{c}\mathrm{o}\mathrm{s}\alpha形式表示，其中\mathrm{c}\mathrm{o}\mathrm{s}\alpha表示滑坡厚度在竖直方向上的投影长度.若前缘倾角不为90°，则滑坡的垂直入水分量同时受到滑动面倾角（α）与前缘倾角（β）的协同影响，S_z则可表示为S\times \mathrm{s}\mathrm{i}\mathrm{n}(\alpha +\beta ).在本研究中，滑坡厚度保持为0.08 m不变，通过调控α与β的组合关系以改变S_z，进而调节滑坡与水体的有效作用面积，对涌浪特征参数产生显著影响.所以a_m和H_m随着α的增加均呈现先增后减的非线性变化趋势；而β的变化除影响有效作用面积外，还在一定程度上改变滑坡体积，使其对涌浪参数的影响机制更加复杂.

3.3　涌浪衰减规律

在每组工况下，4个监测点的a_m和H_m通过拟合关系z=ax^-b （其中x为涌浪传播距离）进行分析，如图8和图9所示.分析表明，a_m和H_m随着涌浪传播距离的增加而逐渐减小，但衰减速率表现出不同的变化趋势，尤其在α和β的不同组合下，涌浪衰减速率呈现出不同的减缓和加速模式.对a_m而言，在α相同条件下，随着α+β趋近于90°，拟合参数b值逐渐减小，表明衰减速率明显减缓；当α+β超过90°后，a_m衰减速率先增强后减弱.进一步分析发现，当β固定为20°或30°时，a_m衰减速率随α的增大大致呈减缓趋势；然而，当β为40°、50°或60°时，a_m衰减速率随α的增大呈现出先减缓、后加速、再减缓的非单调变化趋势.对于H_m，当α=35°时，衰减速率随β增大呈减缓趋势；而当α为45°、50°、55°和65°的情况下，H_m衰减速率则呈现出减缓和加快交替的变化规律.此外，当β固定时，H_m衰减速率随α的增加表现为单调增强.从物理机制的角度看，当α+β接近90°时，滑坡体有效作用面积的推移方向与波浪传播方向趋于一致.导致波浪能量的传播模式与滑坡体的推移作用相互耦合，从而减少了波浪与土体间的能量耗散效应，导致能量耗散速率降低，进而使涌浪的衰减速率相对减缓.

4 讨论

4.1　不同滑坡的涌浪冲击特征对比

与散体滑坡相比，刚体滑坡因其整体性强、不发生形变，且水体无法渗入，具有显著的水土体积置换效应与前缘推动作用，滑坡前缘与水体剧烈碰撞，掀起高幅涌浪，并在前缘形成明显的高速水流区（图10a）.而散体滑坡允许水体渗入滑体内部，削弱了体积置换与推动效应，导致涌浪波高显著降低，涌浪总动能仅为刚体滑坡的21.55%（t=0.3 s）.除水土颗粒间的碰撞导致能量传递外，散体滑坡中还存在与水土相对运动有关的渗流拖曳力.若在数值模拟中忽略拖曳力作用（图10b），则水土间动量传递不充分，涌浪波形偏离实际.考虑渗流拖曳力后，散体滑坡前缘可形成一定范围的高速区（图10c），模拟结果更接近实际，涌浪总动能提升至刚体滑坡的34.01%（t=0.3 s），能量传递机制更为完整.

综上，尽管水体入渗削弱了滑坡体对水体的体积置换与推动效应，降低涌浪强度，但渗流拖曳力在散粒体滑坡-涌浪动力学过程中具有重要作用，其缺失将导致对涌浪强度的显著低估.

4.2　角度协同影响

研究结果表明，当α+β<90°时，随着α+β的增大，a_m和H_m均呈现增加趋势.然而，对于α+β>90°的情况，a_m和H_m的变化规律则表现出更为复杂的非线性特征.例如，当α=50°且β=40°时，a_m和H_m分别为0.081 m和0.126 m；而当β保持不变、α增至55°、65°时，a_m下降了2.59%和7.02%，H_m下降了2.56%和3.02%.进一步分析发现，当α保持不变且β=50°时，β的增大会导致滑坡体积增加约3.79%，但a_m和H_m却分别下降了1.85%和1.84%.然而，当α+β远大于90°时，a_m和H_m的增长规律进一步复杂化.例如，当α=45°且β=60°时，a_m和H_m分别为0.074 m和0.107 m；当α增加至55°时，a_m和H_m分别显著增加了11.42%和19.54%；而当α进一步增加至65°时，a_m和H_m的增幅则分别减缓至6.49%和17.36%.这一结果表明，当α+β>90°时，存在一种竞争机制：一方面，α+β的增大会导致有效作用面积的减少，削弱滑坡对水体的推力，从而抑制a_m和H_m的增长；另一方面，β的增大又会导致滑坡体积的增加，增强动能输入，进而促进a_m和H_m的上升.这种竞争机制揭示了涌浪生成过程中多因素耦合作用的复杂性，其中α和β是导致这种复杂性的两个重要因素.

此外，当β较小时，随着α的增大（即α+β逐渐接近90°），涌浪衰减速率逐渐减缓.然而，当β较大时，初期随着α增大，衰减速率减缓；但随着α进一步增大，滑坡体的滑移速率显著提升，从而输入波浪系统的能量急剧增加（图11）.此时，高速波浪容易引发非线性效应和加剧土体孔隙渗流作用，同时底部摩擦阻力与波浪破碎产生的紊流耗散逐渐占据主导地位，最终导致涌浪衰减速率加速.此外，在α保持不变的情况下，β的变化会影响滑坡体积和有效作用面积，从而对涌浪衰减速率产生影响.定量分析表明，当α固定为50°时，β=40°工况下浪高仪监测域内水体获得的水平动量达到最大值，较β=20°、30°、50°和60°分别提高2.02%、1.20%、0.63%和0.87%，表明此时波浪传播方向与滑坡推移运动方向实现最优耦合（图12）.当β进一步增大至60°时，尽管耦合程度下降，但由于滑坡体积显著增加，系统向水体输入更多、更持续的能量与动量，仍表现出涌浪衰减速率减缓的特征，尽管涌浪传播方向尚未与滑坡推移方向形成最优耦合，仍表现出涌浪衰减速率减缓的特征.上述结果表明，波浪衰减过程中的能量传递与耗散机制高度复杂，且α与β均为其重要影响参数.

4.3　本研究的局限性

本文构建了适用于散粒体滑坡-涌浪过程的动态渗流两相SPH模型，进而初步揭示了滑动面倾角与滑坡前缘倾角对涌浪特征参数的协同影响机制.但需指出的是，当前研究结论尚待物理模型试验进一步验证；同时，本文仅针对近水滑坡自然垮塌入水过程进行分析，尚未涉及高速滑坡情形下滑动面倾角与滑坡前缘倾角协同作用机制的研究.

5 结论与展望

本研究基于动态渗流的两相光滑粒子流体动力学（SPH）模型，通过数值模拟与模型试验的对比验证，系统揭示了渗流效应在滑坡-涌浪过程中的关键作用，并深入分析了滑动面倾角（α）与滑坡前缘倾角（β）对涌浪特征参数的协同影响机制.主要结论如下：

（1）建立了考虑动态渗流的两相SPH模型，该模型基于混合物理论模拟水-土混合体的相互作用，并引入非线性Darcy-Forchheimer模型描述渗流拖曳力.基于室内试验验证了该模型在模拟散粒体滑坡-涌浪过程中的可靠性.与未考虑渗流的SPH模型相比，考虑渗流的SPH模型显著提高了涌浪波高和波形的预测精度，各监测点的最大波浪幅值（a_m）和最大波浪高度（H_m）的误差显著降低，各监测点a_m和H_m的误差分别降低24.72%和41.95%以上.

（2）数值试验结果表明α和β对涌浪特征参数a_m和H_m存在较为显著的协同影响.随着α的增大，a_m和H_m呈现出先递增后递减的非线性趋势；相比之下，β的影响相对较弱，但变化规律更为复杂.当α+β<90°时，a_m和H_m随α+β的增大而增加；而当α+β>90°时，a_m和H_m的变化规律则表现出二者呈现显著的非单调关系，表明存在滑坡体积增加与有效作用面积减少之间的竞争机制.

（3）波浪传播过程中，a_m和H_m的衰减特性受参数α和β的协同调控，呈现显著的非线性特征.研究表明：α增大导致波浪能量增强，与紊流、摩擦耗散效应的耦合作用进一步加剧了衰减速率的复杂性；β主要通过改变有效作用面积影响衰减速率，当α+β趋于90°时，波浪传播方向与滑坡体运动方向趋于一致，能量耗散较弱，致使近场涌浪衰减速率明显降低.

综上所述，本研究基于考虑动态渗流的两相SPH模型，揭示了滑动面倾角与滑坡前缘倾角对涌浪特征参数的协同影响机制.研究成果可为滑坡-涌浪灾害的精确预测与防控提供理论参考，并在工程实践中通过优化滑坡前缘倾角以减弱涌浪危害.上述结论尚未涉及滑坡高速冲击下的波形演化特征及涌浪衰减规律.后续研究应重点突破高速粒子边界穿透问题，发展适用于高速滑坡-涌浪过程的动态渗流两相SPH算法，从而建立更完善的灾害预测理论体系.

参考文献

[1] Barla, G., Paronuzzi, P., 2013. The 1963 Vajont Landslide: 50th Anniversary. Rock Mechanics and Rock Engineering, 46(6): 1267-1270. https://doi.org/10.1007/s00603-013-0483-7

[2] Clous, L., Abadie, S., 2019. Simulation of Energy Transfers in Waves Generated by Granular Slides. Landslides, 16(9): 1663-1679. https://doi.org/10.1007/s10346-019-01180-0

[3] Cui, P., Zhu, X. H., 2011. Surge Generation in Reservoirs by Landslides Triggered by the Wenchuan Earthquake. Journal of Earthquake and Tsunami, 5(5): 461-474. https://doi.org/10.1142/s1793431111001194

[4] Dai, Z. L., Lan, B. S., Jiang, M. T., et al., 2025. Numerical Modeling of Submarine Landslide Motion and Impact Behavior Based on the SPH Method. Journal of Ocean University of China, 24(2): 365-376. https://doi.org/10.1007/s11802-025-5853-8

[5] Evers, F. M., Hager, W. H., 2016. Spatial Impulse Waves: Wave Height Decay Experiments at Laboratory Scale. Landslides, 13(6): 1395-1403. https://doi.org/10.1007/s10346-016-0719-1

[6] Evers, F. M., Hager, W. H., Boes, R. M., 2019. Spatial Impulse Wave Generation and Propagation. Journal of Waterway, Port, Coastal, and Ocean Engineering, 145(3): 04019011. https://doi.org/10.1061/(asce)ww.1943-5460.0000514

[7] Fornaciai, A., Favalli, M., Nannipieri, L., 2019. Numerical Simulation of the Tsunamis Generated by the Sciara Del Fuoco Landslides (Stromboli Island, Italy). Scientific Reports, 9(1): 18542. https://doi.org/10.1038/s41598-019-54949-7

[8] Fritz, H. M., 2001. Lituya Bay Case:Rockslide Impact and Wave Run-up. Science of Tsunami Hazards, 19, 3.

[9] Fritz, H. M., Hager, W. H., Minor, H. E., 2004. Near Field Characteristics of Landslide Generated Impulse Waves. Journal of Waterway, Port, Coastal, and Ocean Engineering, 130(6): 287-302. https://doi.org/10.1061/(asce)0733-950x(2004)130:6(287)

[10] Grilli, S. T., Tappin, D. R., Carey, S., et al., 2019. Modelling of the Tsunami from the December 22, 2018 Lateral Collapse of Anak Krakatau Volcano in the Sunda Straits, Indonesia. Scientific Reports, 9: 11946. https://doi.org/10.1038/s41598-019-48327-6

[11] Heller, V., Spinneken, J., 2013. Improved Landslide-Tsunami Prediction: Effects of Block Model Parameters and Slide Model. Journal of Geophysical Research: Oceans, 118(3): 1489-1507. https://doi.org/10.1002/jgrc.20099

[12] Heller, V., Spinneken, J., 2015. On the Effect of the Water Body Geometry on Landslide–Tsunamis: Physical Insight from Laboratory Tests and 2D to 3D Wave Parameter Transformation. Coastal Engineering, 104: 113-134. https://doi.org/10.1016/j.coastaleng.2015.06.006

[13] Huang, B. L., Yin, Y. P., Du, C. L., 2016. Risk Management Study on Impulse Waves Generated by Hongyanzi Landslide in Three Gorges Reservoir of China on June 24, 2015. Landslides, 13(3): 603-616. https://doi.org/10.1007/s10346-016-0702-x

[14] Huang, C., Hu, C., An, Y., et al., 2023. Numerical Simulation of the Large-Scale Huangtian (China) Landslide-Generated Impulse Waves by a GPU-Accelerated Three-Dimensional Soil‒Water Coupled SPH Model. Water Resources Research, 59(6): e2022WR034157. https://doi.org/10.1029/2022wr034157

[15] Jiang, Q., 2019. Unified Particle Method Research for Simulation of Landslides Generated Waves in Reservoir Bank (Dissertation). Ningbo Institute of Material Technology, Chinese Academy of Sciences, Ningbo(in Chinese with English abstract).

[16] Lee, C. H., Huang, Z. H., 2022. Effects of Grain Size on Subaerial Granular Landslides and Resulting Impulse Waves: Experiment and Multi-Phase Flow Simulation. Landslides, 19(1): 137-153. https://doi.org/10.1007/s10346-021-01760-z

[17] Lee, C. H., Lo, P. H., Shi, H. B., et al., 2022. Numerical Modeling of Generation of Landslide Tsunamis: A Review. Journal of Earthquake and Tsunami, 16(6): 2241001. https://doi.org/10.1142/s1793431122410019

[18] Li, H. W., Xu, Z. G., Shi, J. Y., et al., 2024. Tsunami Potential Threat from the Ryukyu Trench on Chinese Coast Based on Subduction Zone Dynamics Parameters. Earth Science, 49(2): 403-413 (in Chinese with English abstract).

[19] Li, P. F., Jing, H. X., Li, G. D., 2024. Generation and Prediction of Water Waves Induced by Rigid Piston-Like Landslide. Natural Hazards, 120(3): 2683-2704. https://doi.org/10.1007/s11069-023-06300-7

[20] Li, Q. W., Huang, B. L., Zhang, P., et al., 2024. Influence of the Degree of Landslide Fragmentation on the Characteristics of Landslide Impulse Wave. Rock and Soil Mechanics, 45(11): 3345-3354 (in Chinese with English abstract).

[21] Liu, J. X. Z., 2023. Partitioning Prediction Study of Landslide-Tsunamis in the Wu Gorge of the Three Gorges Reservoir Area(Dissertation). China University of Geosciences, Wuhan (in Chinese with English abstract).

[22] Luo, M., Khayyer, A., Lin, P. Z., 2021. Particle Methods in Ocean and Coastal Engineering. Applied Ocean Research, 114: 102734. https://doi.org/10.1016/j.apor.2021.102734

[23] Mao, Y. F., Guan, M. F., 2023. Mesh-Free Simulation of Height and Energy Transfer of Landslide-Induced Tsunami Waves. Ocean Engineering, 284: 115219. https://doi.org/10.1016/j.oceaneng.2023.115219

[24] Meng, Z. Z., Zhang, J. X., Hu, Y. T., et al., 2023. Temporal Prediction of Landslide-Generated Waves Using a Theoretical-Statistical Combined Method. Journal of Marine Science and Engineering, 11(6): 1151. https://doi.org/10.3390/jmse11061151

[25] Mohammed, F., Fritz, H. M., 2012. Physical Modeling of Tsunamis Generated by Three-Dimensional Deformable Granular Landslides. Journal of Geophysical Research (Oceans), 117(C11): C11015. https://doi.org/10.1029/2011JC007850

[26] Paquier, A. E., Oudart, T., Le Bouteiller, C., et al., 2021. 3D Numerical Simulation of Seagrass Movement under Waves and Currents with GPUSPH. International Journal of Sediment Research, 36(6): 711-722. https://doi.org/10.1016/j.ijsrc.2020.08.003

[27] Rauter, M., Viroulet, S., Gylfadóttir, S. S., et al., 2022. Granular Porous Landslide Tsunami Modelling—The 2014 Lake Askja Flank Collapse. Nature Communications, 13(1): 678. https://doi.org/10.1038/s41467-022-28296-7

[28] Tang, G. Q., Lu, L., Teng, Y. F., et al., 2018. Impulse Waves Generated by Subaerial Landslides of Combined Block Mass and Granular Material. Coastal Engineering, 141: 68-85. https://doi.org/10.1016/j.coastaleng.2018.09.003

[29] Viroulet, S., Sauret, A., Kimmoun, O., et al., 2013. Granular Collapse into Water: Toward Tsunami Landslides. Journal of Visualization, 16(3): 189-191. https://doi.org/10.1007/s12650-013-0171-4

[30] Wu, H., Shi, A. C., Ni, W. D., et al., 2024a. Numerical Simulation on Potential Landslide–Induced Wave Hazards by a Novel Hybrid Method. Engineering Geology, 331: 107429. https://doi.org/10.1016/j.enggeo.2024.107429

[31] Wu, H., Zhong, Q. M., Deng, Z., et al., 2024b. Numerical Investigation of the Effect of Landslide Relative Density on the Impulse Wave Amplitude. Ocean Engineering, 309: 118563. https://doi.org/10.1016/j.oceaneng.2024.118563

[32] Xu, Q., Dong, X. J., 2011. Genetic Types of Large-Scale Landslides Induced by Wenchuan Earthquake. Earth Science, 36(6): 1134-1142 (in Chinese with English abstract).

[33] Xu, W. J., 2023. Research Advances in Disaster Dynamics of Landslide Tsunami. Journal of Engineering Geology, 31(6): 1929-1940 (in Chinese with English abstract).

[34] Yin, K. L., Liu, Y. L., Wang, Y., et al., 2012. Physical Model Experiments of Landslide-Induced Surge in Three Gorges Reservoir. Earth Science, 37(5): 1067-1074 (in Chinese with English abstract).

[35] Yu, M. L., Lee, C. H., 2019. Multi-Phase-Flow Modeling of Underwater Landslides on an Inclined Plane and Consequently Generated Waves. Advances in Water Resources, 133: 103421. https://doi.org/10.1016/j.advwatres.2019.103421

[36] Zhang, C., Rezavand, M., Zhu, Y. J., et al., 2021. SPHinXsys: An Open-Source Multi-Physics and Multi-Resolution Library Based on Smoothed Particle Hydrodynamics. Computer Physics Communications, 267: 108066. https://doi.org/10.1016/j.cpc.2021.108066

[37] Zhang, S. H., Zhang, C., Hu, X. Y., et al., 2024. A Riemann-Based SPH Method for Modelling Large Deformation of Granular Materials. Computers and Geotechnics, 167: 106052. https://doi.org/10.1016/j.compgeo.2023.106052

[38] Zhu, C. W., Peng, C., Wu, W., et al., 2022. A Multi-Layer SPH Method for Generic Water–Soil Dynamic Coupling Problems. Part I: Revisit, Theory, and Validation. Computer Methods in Applied Mechanics and Engineering, 396: 115106. https://doi.org/10.1016/j.cma.2022.115106

[39] Zhu, Y. F., An, C., 2024. Application of Uniform Slip Models to Tsunami Early Warning: A Case Study of 2021 Mw 8.2 Alaska Peninsula Earthquake. Earth Science, 49(2): 500-510 (in Chinese with English abstract).

基金资助

国家自然科学基金-联合基金重点项目(U23A2044)

广西重点研发计划项目(No.桂科AB25069121)

国家自然科学基金-青年基金项目(42207238)

福建省自然科学基金项目-面上项目(2022J011253)