基于光滑粒子元的水-沙两相流冲刷数值仿真

Numerical Simulation of Scour by Liquid-sediment Two-phase Flow Based on Smoothed Particle Hydrodynamics

摘要

基础河床冲刷逐渐成为水工建筑物等结构破坏的主要原因，传统分析方法基于欧拉网格法模拟，易出现网格畸变、计算不收敛等情况，降低了计算效率与精度.基于拉格朗日坐标系下无网格光滑粒子元法（Smoothed Particle Hydrodynamics,SPH）对河床冲刷过程进行数值模拟；将泥沙相视为非牛顿流体，分为沉积物、推移质、悬移质三种状态；引入DruckerPrager与Shields应力模型作为泥沙状态转化判断准则，为不同状态泥沙粒子赋予不同流变特性，使计算结果可准确描述水-沙耦合过程.提出了基于SPH多相流模型的河床冲刷数值仿真改进算法，编写了泥沙冲刷计算模块，采用GPU进行计算加速，最后将所得数值水槽模型与Louvain溃坝试验及其他同类仿真结果进行对照.结果表明：本文模型能够更为准确地反映冲刷发展整体趋势；各时刻自由液面、水-沙交界面轮廓均方根误差均在合理范围内；数值模型与试验结果吻合度较高.

Abstract

Scours at the foundations of hydraulic structures have gradually become the major cause of structural damage. The traditional numerical simulation is based on the Euler mesh method， which is not easier to converge or accompanied by large deformation of grids and eventually results in the loss of solution efficiency and accuracy. This paper numerically simulated the scour process of riverbed by utilizing the Smooth Particle Hydrodynamics（SPH）based on the Lagrange coordinate. It treated the sediment phase as a non-Newtonian phase and divided it into three states： sediment， bed load， and suspended load. To accurately describe the effect of liquid-sediment interaction， the Drucker-Prager and Shields stress model was introduced in the numerical model as the criteria for the transformation judgment between three states of sediment， including assigning different rheological properties for sediment particles in different states. In this study， a modification scouring algorithm based on the two-phase flow was proposed， and the sediment scouring calculation module was built and accelerated by GPU. Finally， a numerical flume model was designed for comparison with the Louvain dam-break experiment as well as a similar numerical model. The conclusions are drawn： the current numerical model proposed could more accurately reflect the overall trend of scour development； the RMSE of the free surface profile of the water and sediment-liquid interface were within a reasonable range at the specified moment； the result of numerical model was in good agreement with the experimental data.

Graphical abstract

关键词

桥梁工程 / 数值仿真 / 光滑粒子元 / 河床冲刷 / 两相流 / 侵蚀准则

水流侵蚀导致河床冲刷是水中建筑结构物常见现象，频繁发生于桥梁基础、防波堤、海底管线与堤坝下游等. 2018年，四川省堰塞湖泄洪导致下游河床剧烈冲刷，造成下游竹巴龙金沙江大桥冲毁，多座大桥下部基础被掏空而发生倾斜^［1］；余文畴等^［2］研究发现，大通站长江干流1950—2002年年平均径流量中，每立方米水体平均含沙量达到47.2%，可见河床侵蚀作用十分显著. 因此，有必要针对河床冲刷演进过程进行研究，以指导水工建筑物抗水设计，妥善预防及处置冲刷灾害，减少损失.

床底泥沙受到水流冲击，会在表面产生切应力，当泥沙表面屈服时，细颗粒泥沙将以推移质和悬移质的形式向下游推移、输送，并逐渐形成冲刷坑. 冲刷诱因多种多样：暴雨、泄洪、溃坝等短时间内产生巨大峰值流量皆会造成河床剧烈冲刷. 冲刷本质皆为河床自由表面发生高度非线性变形、破碎以及水和泥沙相互夹带.传统数值模拟方法基于欧拉坐标系建立流域模型，通过数值求解输沙方程和河床变形方程，获得床面各处高程变化，最后采用动网格技术以及体积分数法（VOF）实现河床形态可视化^［3］. 然而，动网格技术对于网格质量具有较高要求，冲刷伴随的河床大变形通常会造成网格畸变、计算不收敛等情况^［4-5］；VOF法通过捕捉水-沙交界面网格中材料体积分数分布情况实现河床形态可视化，因此精度受到交界面网格划分密度限制，对河床剧烈变形的区域，往往不能精准捕捉^［6］. 可见，传统网格法冲刷数值模型仍存在技术缺陷.

基于拉格朗日坐标系的光滑粒子流体动力学（Smoothed Particle Hydrodynamics，SPH）法，克服了传统网格法的固有局限，与欧拉法相比，SPH法粒子本身具有质量，无需额外计算就能保证质量守恒；模拟复杂自由表面流时，无需追踪流体边界及不同流体交界面，对于复杂液面具有优秀的求解能力，可适用于模拟泥沙冲刷、输运等问题^［7-11］. 现阶段，SPH法模拟河床冲刷主要有两种思路：其一，不进行土体粒子建模，使满足条件的边界粒子转化为土粒子进入流域参与计算. 转化后土粒子视为离散元（DEM），不赋予材料本构，其运动方程受牛顿第二定律控制；土体与水体间考虑流固耦合作用力，土体粒子间考虑接触、碰撞作用^{［8， 11］}. 然而，此类模型计算相对复杂，与真实冲刷关联较弱. 其二，基于SPH多相流理论进行水体与土体建模^［9-10］. 土体粒子作为流体相可赋予材料本构，其变形与输运具有真实物理意义.

目前，基于SPH多相流理论的泥沙冲刷研究多集中于河床溃坝冲刷，土体多为非黏性泥沙颗粒（以下简称为泥沙），泥沙粒子流变特性通常采用非牛顿流体进行模拟^［8-14］. Manenti等^［12］在冲淤研究中，基于非牛顿流体理论构建了土体模型，并分别对Mohr-Coulomb屈服理论和Shields理论进行了比较分析，推荐后者作为泥沙模型的材料屈服强度；Fourtakas等^［13］引入了广义Herschel-Bulkley-Papanastasiou（HBP）模型进行土体建模，结合DP强度准则确定材料屈服强度，并应用于水-沙两相流冲刷模拟中，获得了良好效果.

然而，现阶段基于SPH多相流的冲刷数值仿真方法仍处于发展阶段：泥沙本构模型不统一，多数研究采用单一的土体强度准则，并集中于二维模型，对于三维泥沙冲刷模型及组合强度准则研究较少，对泥沙状态转变考虑不够详细^［8-14］. 由于SPH方法对于计算资源要求较高，三维模型需要耗费大量计算时间；同时冲刷过程中土体状态变化急促、剧烈，单一准则难以完整表征泥沙由静止到起动乃至输运的全过程特性变化，最终造成河床冲刷形态与实际仍有差距.

针对上述技术难题，本文提出了基于SPH多相流模型的河床冲刷数值仿真改进算法，主要包含三个步骤：①构建泥沙粒子；②根据侵蚀模型和起动准则判断泥沙粒子状态并分类转化；③计算不同状态下泥沙粒子流变特性. 具体来说，基于开源软件 DualSPHysics进行二次开发，编写出泥沙冲刷计算模块，利用开源CUDA框架，进行GPU计算加速；基于SPH多相流理论与HBP非牛顿流体模型，实现河床溃坝冲刷数值建模；引入Drucker-Prager与Shields应力模型作为泥沙强度计算准则，考虑泥沙冲刷全过程中由沉积物到推移质再到悬移质的状态转变；最后，通过溃坝冲刷数值模拟，结合水槽试验验证模型准确性.

1 SPH基本理论

光滑粒子流体动力学（SPH）方法将连续流体离散为具有各种物理量的粒子，通过核函数实现粒子间相互作用，并按Navier-Stokes方程进行运动控制. 本节介绍了SPH法控制方程及其离散形式，结合层流黏性应力与亚粒子（Sub-Particle Scale）紊流应力模型封闭N-S方程.

1.1 SPH法的积分插值理论

SPH法根据相邻粒子物理性质在每个粒子周围进行局部积分，从而更新其物理量. 相邻粒子由核函数检索，用W表示.

核函数W是关于光滑长度h的函数，其作用是将连续微分方程基于插值函数在特定点处进行积分，任意物理量连续形式用F可表示为

F\left(\mathit{r}\right)=\int F\left(\mathit{r}\mathrm{\text{'}}\right)W\left(\mathit{r}-\mathit{r}\mathrm{\text{'}},h\right)\mathrm{d}\mathit{r}\mathrm{\text{'}}.（1）

式中：F代表粒子的任意物理量；\mathit{r}和\mathit{r}\mathit{\text{'}}分别表示中心粒子和邻域粒子位置矢量；h为光滑长度. 将F在A点插值并离散化，则其离散形式表达为

F\left({\mathit{r}}_a\right)\approx \sum _{b}F\left({\mathit{r}}_b\right)W\left({\mathit{r}}_a-{\mathit{r}}_b,h\right)\mathrm{\Delta }{v}_b.（2）

式中：下标a表示中心粒子；下标b表示邻域粒子；\mathrm{\Delta }{v}_b表示邻域粒子体积. 核函数选取参照Fourtakas等^［13］的研究采用五次核函数^［15］：

W\left(r,h\right)=\begin{array}{cc}\frac{\mathrm{21}}{\mathrm{16}\mathrm{\pi }{h}^{\mathrm{3}}}{\left(\mathrm{1}-\frac{q}{\mathrm{2}}\right)}^{\mathrm{4}}\left(\mathrm{2}q+\mathrm{1}\right),& \mathrm{0}\le q\le \mathrm{2}\end{array}.（3）

式中：q为粒子间距与平滑长度的比值.

1.2 控制方程的离散形式

N-S方程中质量守恒方程与动量守恒方程分别表示为式（4）与式（5）：

\frac{\mathrm{d}\rho }{\mathrm{d}t}+\rho \frac{\partial \mathit{u}}{\partial x}=\mathrm{0}.（4）

\frac{\mathrm{d}\mathit{u}}{\mathrm{d}t}=-\frac{\mathrm{1}}{\rho }\nabla P+\mathit{g}+\Gamma .（5）

式中：\mathit{u}为流速矢量；\mathit{g}为重力加速度矢量；P为压力项；\rho为流体密度；\Gamma表示流体黏性力耗散项.

将质量守恒方程改写为SPH离散形式：

\frac{\mathrm{d}{\rho }_a}{\mathrm{d}t}={\rho }_a\sum _b\frac{m_b}{{\rho }_b}{\mathit{u}}_{ab}\cdot {\nabla }_a{W}_{ab}.（6）

式中：下标a表示中心粒子；下标b表示邻域粒子；下标ab表示粒子a与b的差；{\nabla }_a为针对粒子a的哈密顿算子；m为粒子质量；{\rho }_a与{\rho }_b分别表示粒子a与粒子b的密度.

通过Favre平均法对密度进行加权平均，可将SPS项引入SPH方法中，动量守恒方程改写为SPH离散形式^［16］：

\begin{array}{ll}\frac{\mathrm{d}{\mathit{u}}_a}{\mathrm{d}t}=& -\sum _b{m}_b\left(\frac{P_b+P_a}{{\rho }_b\cdot {\rho }_a}\right){\nabla }_a{W}_{ab}+\mathit{g}+\\ & \sum _b{m}_b\left[\frac{\mathrm{4}{\mu }_{\mathrm{0}}{\mathit{r}}_{ab}\cdot {\nabla }_a{W}_{ab}}{\left({\rho }_a+{\rho }_b\right)\left(r_{ab}^{\mathrm{2}}+{\eta }^{\mathrm{2}}\right)}\right]{\mathit{u}}_{ab}+\\ & \sum _b{m}_b\left(\frac{{\tau }_{ij}^b}{{\rho }_b^{\mathrm{2}}}+\frac{{\tau }_{ij}^a}{{\rho }_a^{\mathrm{2}}}\right){\nabla }_a{W}_{ab}.\end{array}（7）

式中：{\mu }_{\mathrm{0}}为运动黏度；{\tau }_{ij}为SPS应力矢量.

2 基于SPH的两相流数值模型

基于SPH多相流模型的冲刷模拟，通常将泥沙视为非牛顿流体相，常见模型有Bingham模型^［17］、HB（Herschel-Bulkley）^［18］模型、HBP（Herschel-Bulkley-Papanastasiou）模型^［13］等.此类非牛顿模型在达到屈服强度前后，其表观黏度变化显著：材料所受剪切应力达到屈服强度前，黏度通常较大，流体表现为“难流通”，而在材料屈服后，表观黏度通常急剧下降，逐渐转变为“易流通”，因而适宜模拟冲刷起动前后的泥沙流变行为.然而，Bingham模型及HB模型在低剪切速率时其应力非连续，容易造成数值计算不收敛. 本文采用HBP模型建立水-沙两相流模型，将沉积相分为沉积物、推移质、悬移质三种状态，赋予不同流变特性进行计算.

2.1 泥沙相非牛顿流体模型

HBP模型计算流体表观黏度如下：

{\mu }_{\mathrm{o}\mathrm{r}\mathrm{i}}=\frac{\left|{\tau }_{\mathrm{c}}\right|}{\left|{\mathrm{Ⅱ}}_{\mathrm{D}}\right|}\left(\mathrm{1}-{\mathrm{e}}^{-m\sqrt[]{{\mathrm{Ⅱ}}_{\mathrm{D}}}}\right)+\mathrm{2}\mu {\left|\mathrm{4}{\mathrm{Ⅱ}}_{\mathrm{D}}\right|}^{\frac{n-\mathrm{1}}{\mathrm{2}}}.（8）

式中：{\mu }_{\mathrm{o}\mathrm{r}\mathrm{i}}为泥沙表观黏度；\mu为黏性系数；n为Herschel-Bulkley幂指数参数，与剪切应力相关；m为Papanastasiou参数，控制应力指数增长速率，使应力在未屈服区域保持小剪切速率，在已屈服区域保持线性增长；{\mathrm{Ⅱ}}_{\mathrm{D}}为二阶不变剪切应变率张量；{\tau }_{\mathrm{c}}为材料屈服强度.

Herschel-Bulkley模型是广义Bingham模型，屈服后具有Bingham模型流变特性，而Papanastasiou模型适合描述剪切速率趋于0时材料较大的表观黏度，两模型结合使得HBP模型在模拟冲刷时稳定性更高，可连续表达土体从低剪切速率发展到高剪切速率的全过程^［13］（图1）.此外，当m=0，n=1时，该模型降阶为牛顿流体模型，使得多相流计算中不同流相具有相同广义模型. 因而，在DualSPHysics多相流耦合计算中，对于任意粒子核函数半径内的被检索粒子，仅需基于HBP模型对所有被检索粒子进行黏性力求解，从而规避了对土体粒子的分类检索计算，简化了计算框架，加快了计算效率^［13］.

2.2 泥沙侵蚀模型

为了模拟泥沙侵蚀与输运过程，本文将泥沙划分为三种状态：沉积物、推移质、悬移质^［19］. 泥沙粒子状态转化与邻域粒子相关，受其体积浓度、速度及所受剪切应力等因素影响，根据不同状态，赋予不同黏度参与SPH计算.

2.2.1 泥沙屈服强度模型

Drucker-Prager（DP）屈服准则一般形式为：

\sqrt[]{{\mathit{J}}_{\mathrm{2}}}+(\alpha p-\beta )=\mathrm{0}.（9）

式中：{\mathit{J}}_{\mathrm{2}}为二阶不变剪切应力张量；p为饱和泥沙粒子上作用的静水压力；\alpha与\beta由Mohr-Coulomb屈服准则参数给出：

\alpha =\frac{\mathrm{2}\mathrm{s}\mathrm{i}\mathrm{n}\phi }{\sqrt[]{\mathrm{3}}\left(\mathrm{3}-\mathrm{s}\mathrm{i}\mathrm{n}\phi \right)},\beta =\frac{\mathrm{6}c\mathrm{s}\mathrm{i}\mathrm{n}\phi }{\sqrt[]{\mathrm{3}}\left(\mathrm{3}-\mathrm{s}\mathrm{i}\mathrm{n}\phi \right)}.（10）

式中：\phi为内摩擦角；c为土体黏聚力.

根据Fourtakas等^［13］的研究，{\mathit{J}}_{\mathrm{2}}可由土粒子所受水体剪切应力{\tau }_{\alpha \beta }计算：

{\mathit{J}}_{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{2}}{\tau }_{\alpha \beta }{\tau }_{\alpha \beta }.（11）

剪切应力{\tau }_{\alpha \beta }可由二阶不变剪切应变率张量Ⅱ_D表述：

{\tau }^{\alpha \beta }=\mathrm{2}\sqrt[]{\mathrm{2}}{\mu }_{\mathrm{d}}\sqrt[]{{\mathrm{Ⅱ}}_{\mathrm{D}}}.（12）

联立式（11）、式（12）可以建立{\mathit{J}}_{\mathrm{2}}与{\mathrm{Ⅱ}}_{\mathrm{D}}间联系：

\sqrt[]{{\mathit{J}}_{\mathrm{2}}}=\mathrm{2}{\mu }_{\mathrm{d}}\sqrt[]{{\mathrm{Ⅱ}}_{\mathrm{D}}}.（13）

进而，\left|{\mathit{J}}_{\mathrm{2}}\right|=\left|\alpha p-\beta \right|可作为HBP模型中材料屈服强度限定条件. 计算床面下静止土体屈服强度时，引入Drucker-Prager（DP）屈服准则：

\left|{\tau }_{\mathrm{y}}\right|=\left|\alpha p+\beta \right|.（14）

当沉积物粒子所受切应力大小不超过其屈服强度\left|{\tau }_{\mathrm{y}}\right|时，土体粒子由于高黏度特性几乎不发生剪切变形，其表观黏度{\mu }_{\mathrm{a}\mathrm{p}\mathrm{p}}由\left|{\tau }_{\mathrm{y}}\right|代入式（8）进行计算.

对于河床表面与水接触的土粒子，若满足由沉积物转化为推移质的条件，土体屈服强度采用Shields准则计算并替换，详见2.2.2节.

2.2.2 泥沙起动判定

初始条件下，将泥沙层置于水层下方，当泥沙粒子所受切应力不超过临界值时，泥沙粒子固定不动；达到临界值时，泥沙由沉积物转化为推移质，随水流沿床面输运.

泥沙临界切应力由侵蚀准则计算. Manenti等^［12］在冲淤研究中，对Mohr-Coulomb屈服理论和Shields理论进行了比较分析. 研究表明，Shields理论对于模型参数均具有较高敏感性；同时Shields公式中物理量更易由试验测得，相较Mohr-Coulomb准则中应变率等参数，计算灵敏度更高，本文最终选择Shields理论作为泥沙状态转化判断准则.

转化为推移质的泥沙粒子须位于河床表面，同时考虑到推移质粒子为饱和土的特性，其转化条件设置以下两个：

1）在泥沙粒子的邻域粒子中，至少包括一个水粒子.

2）待转化泥沙粒子实际质量应小于阈值质量.

当泥沙粒子同时满足以上条件时，基于Shields侵蚀准则，计算水平床面临界转化剪切应力：

{\tau }_{\mathrm{b}\mathrm{c}\mathrm{r},\mathrm{0}}={\theta }_{\mathrm{c}\mathrm{r}}({\rho }_{\mathrm{s}}-{\rho }_{\mathrm{w}})gd.（15）

式中：{\tau }_{\mathrm{b}\mathrm{c}\mathrm{r},\mathrm{0}}为水平床面上的临界剪切应力；d为特征粒径，对于非均匀材料取中值粒径d_{\mathrm{50}}；{\rho }_{\mathrm{s}}为泥沙饱和密度；{\theta }_{\mathrm{c}\mathrm{r}}为临界Shields数，是泥沙粒子所受临界起动力与最大抗力之比；{\rho }_{\mathrm{w}}为水的密度；g 为重力加速度.

{\theta }_{\mathrm{c}\mathrm{r}}求解基于雷诺数Re ^［12］：

\left\{\begin{array}{l}{\theta }_{\mathrm{c}\mathrm{r}}=\mathrm{0.010}\mathrm{595}\mathrm{l}\mathrm{n}\left(Re\right)+\frac{\mathrm{0.110}\mathrm{476}}{R_{\mathrm{*}}}+\mathrm{0.002}\mathrm{719}\mathrm{7},Re\le \mathrm{500};\\ {\theta }_{\mathrm{c}\mathrm{r}}=\mathrm{0.068},Re>\mathrm{500}.\end{array}\right.（16）

式中：R_{*}=\frac{u_{*}d}{v}，v为水体运动黏度，u_{*}为摩阻流速，按式（17）计算：

u_{\mathrm{*}}=\mathrm{11.6}\frac{v}{\delta }.（17）

根据Prandtl混合长度理论，假定沉积物表面存在厚度为\delta的层流亚层，且流速u\left(z\right)按线性分布；在其上紊流层中，流速u\left(z\right)按对数分布：

\left\{\begin{array}{ll}u(z)=\frac{u_{\mathrm{*}}^{\mathrm{2}}}{v}z,& z\le \delta ;\\ \frac{u(z)}{u_{\mathrm{*}}}=\frac{\mathrm{1}}{\kappa }\mathrm{l}\mathrm{n}\left(\frac{z}{z_{\mathrm{0}}}\right),& z>\delta .\end{array}\right.（18）

式中：\kappa为冯·卡曼常数，取0.41；z为水粒子到泥沙-水交界处的垂直距离；z_{\mathrm{0}}为床底粗糙度，计算如下：

（19）

式中：k_{\mathrm{s}}为当量砂径糙率，按经验取泥沙中值粒径.

由式（17）、式（19）可知，摩阻流速u_{*}与\delta和z_{\mathrm{0}}有关，因此需迭代求解. 设定u_{*}初始值u_{*\mathrm{0}}从的稳态流开始迭代，直至流速曲线u\left(z\right)在边界层厚度\delta收敛，即u\left({\delta }^{-}\right)=u\left({\delta }^{+}\right)，此时可求得摩阻流速u_{*}，继而求出水平床面的临界剪切应力{\tau }_{\mathrm{b}\mathrm{c}\mathrm{r},\mathrm{0}}，以及斜坡修正临界剪切应力{\tau }_{\mathrm{b}\mathrm{c}\mathrm{r}}.

水流冲刷的起动力{\tau }_{\mathrm{b}}由Einstein对数流速分布公式^［20］计算：

\frac{{\tau }_{\mathrm{b}}}{\rho }={\left(\kappa d\right)}^{\mathrm{2}}{\left(\frac{\mathrm{\Delta }u}{d}\right)}^{\mathrm{2}}.（20）

式中：d为粒子特征粒径；∆u为泥沙粒子与附近水粒子的速度差，取中心泥沙粒子邻域内所有水粒子进行积分插值：

\mathrm{\Delta }u=\sum u_{\mathrm{b}}W\left({\mathit{r}}_{\mathrm{a}}-{\mathit{r}}_{\mathrm{b}},h\right)\mathrm{\Delta }{V}_{\mathrm{b}}.（21）

若泥沙粒子所受剪切应力{\tau }_{\mathrm{b}}超过临界剪切应力{\tau }_{\mathrm{b}\mathrm{c}\mathrm{r}}，则判断泥沙粒子由沉积物转变为推移质. 转化后粒子表观黏度需结合HBP模型进行更新，将{\tau }_{\mathrm{b}\mathrm{c}\mathrm{r}}代入式（8）中代替屈服应力力τ_c，计算得到黏度μ_ori.

2.2.3 泥沙上升与沉降判定

已屈服推移质粒子在满足条件时会转化为悬移质，被水流挟带而发生输运. 当流速减小、大量悬移质粒子聚集时，悬移质再次沉降，转化为推移质.

悬移质周围泥沙粒子体积浓度C_{v,i}不超过临界浓度C_{\mathrm{v}\mathrm{c}\mathrm{r}}：

C_{v,i}=\frac{V_{\mathrm{s}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}}}{V}=\frac{{\sum }_{j_{\mathrm{s}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}}\in \mathrm{2}h}^N\frac{m_j}{{\rho }_j}}{{\sum }_{j\in \mathrm{2}h}^N\frac{m_j}{{\rho }_j}}\le C_{\mathrm{v}\mathrm{c}\mathrm{r}}.（22）

式中：i表示中心泥沙粒子；j表示邻域粒子； 临界体积浓度C_{\mathrm{v}\mathrm{c}\mathrm{r}}取0.3^［12］，当满足式（22）时，泥沙粒子可以作为牛顿流体处理，可以进一步考察流速条件.

引入Mastbergen公式^［21］，计算泥沙临界上升流速与临界沉降流速，如式（23）：

u_{\mathrm{l}\mathrm{i}\mathrm{f}\mathrm{t}}={\alpha }_i{n}_{\mathrm{s}}{d_{\mathrm{*}}}^{\mathrm{0.3}}(\theta -{\theta }_{\mathrm{c}\mathrm{r}}{)}^{\mathrm{1.5}}\sqrt[]{\frac{d({\rho }_{\mathrm{s}}-\rho )g}{\rho }}.（23-a）

u_{\mathrm{s}\mathrm{e}\mathrm{t}}=\frac{v}{d}\left(\sqrt[]{\mathrm{10}.{\mathrm{36}}^{\mathrm{2}}+\mathrm{1.049}{d_{\mathrm{*}}}^{\mathrm{3}}}-\mathrm{10.36}\right).（23-b）

式中：{\alpha }_i为泥沙输运系数；n_{\mathrm{s}}为床面法向量；d_{*}为泥沙粒径系数，d_{\mathrm{*}}=d_{\mathrm{50}}{\left[\rho \left({\rho }_{\mathrm{s}}-\rho \right)g/{\mu }^{\mathrm{2}}\right]}^{\frac{\mathrm{1}}{\mathrm{3}}}；\theta为泥沙粒子当前Shields系数，其余符号同上.

若u\ge u_{\mathrm{l}\mathrm{i}\mathrm{f}\mathrm{t}}，则判断推移质转化为悬移质，作为牛顿流体参与计算.其黏度采用Vand胶体^［22］方程计算：

{\mu }_{\mathrm{s}\mathrm{u}\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}}=v{\mathrm{e}}^{\frac{\mathrm{0.5}{C}_{\mathrm{v}}}{\mathrm{1}-\frac{\mathrm{39}}{\mathrm{64}}{C}_{\mathrm{v}}}}.（24）

若u\le u_{\mathrm{s}\mathrm{e}\mathrm{t}}，则判断悬疑质沉降，转化为推移质，黏度按推移质计算.

2.2.4 临界起动切应力斜坡修正

对于任意坡度下泥沙起动分析，需对水平床面临界剪切应力{\tau }_{\mathrm{b}\mathrm{c}\mathrm{r},\mathrm{0}}进行修正，得到适用于斜坡的临界剪切应力{\tau }_{\mathrm{b}\mathrm{c}\mathrm{r}}. Van Rijn^［23］将修正系数分解为纵向坡和横向坡单独影响；Chen等^［24］根据受力平衡推导出任意三维斜坡临界起动切应力修正系数. 本文在现有工作基础上，考虑水流升力作用影响，进行临界起动切应力斜坡修正.

假设泥沙粒子为一质点，处于斜坡床面上，其受力包括水下重力\mathit{W}（包含重力与浮力）、水流拖曳力{\mathit{F}}_{\mathrm{D}}、水流升力{\mathit{F}}_{\mathrm{L}}、抵抗起动力{\mathit{F}}_{\mathrm{C}}（本文中即土体的库伦摩擦力）. 斜坡、水流方向粒子受力如图2所示.

图2中， i 、 j 、 k 为正交直角坐标系的单位向量； l 、 t 所在平面为斜坡床面，其中 l 处于XOZ平面，与X轴夹角为\alpha， t 处于YOZ平面，与Y轴夹角为\beta； n 为倾斜坡面单位法向量，与Z轴正方向夹角为\gamma，由几何关系可知\mathrm{c}\mathrm{o}\mathrm{s}\gamma =\mathrm{1}/\sqrt[]{\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\alpha +\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\beta +\mathrm{1}}. 水流拖曳力{\mathit{F}}_{\mathrm{D}}方向与流速 u 方向相同. 根据Dey^［25］的研究，假定水流升力{\mathit{F}}_{\mathrm{L}}方向与坡面法向量 n 一致，大小与拖曳力比值η为0.85.

水流拖曳力与水下重力沿坡面的分量合成泥沙起动力\mathit{F}：

\mathit{F}={\mathit{F}}_{\mathrm{D}}+\left[\mathit{W}-\left(\mathit{W}\cdot \mathit{n}\right)\mathit{n}\right].（25）

{\mathit{F}}_{\mathrm{D}}=\left|{\mathit{F}}_{\mathrm{D}}\right|\cdot \mathit{e}=\left(\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_x\mathit{i}+\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_y\mathit{j}+\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_z\mathit{k}\right)\cdot \left|{\mathit{F}}_{\mathrm{D}}\right|.（26）

式中：{\theta }_x、{\theta }_y、{\theta }_z分别为单位向量\mathit{e}与X、Y、Z轴正方向的夹角.

抵抗起动力{\mathit{F}}_{\mathrm{C}}由库伦摩擦力产生，与\mathit{F}方向相反，在极限状态下，

{\mathit{F}}_{\mathrm{C}}=-\left|\left(\mathit{W}\cdot \mathit{n}\right)\mathit{n}+{\mathit{F}}_{\mathrm{L}}\right|\cdot \mathrm{t}\mathrm{a}\mathrm{n}\phi \cdot \frac{\mathit{F}}{\left|\mathit{F}\right|}.（27）

式中：\phi表示泥沙内摩擦角. 当泥沙粒子达到临界起动状态，即满足

\mathit{F}+{\mathit{F}}_{\mathrm{C}}=\dot{\mathrm{0}}或\left|\mathit{F}\right|=\left|{\mathit{F}}_{\mathrm{C}}\right|.（28）

由式（28），推导出\left|{\mathit{F}}_{\mathrm{D}}\right|表达式为：

\left|{\mathit{F}}_{\mathrm{D}}\right|=\sqrt[]{\frac{\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\phi \mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\gamma +{\eta }^{\mathrm{2}}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\phi \mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\gamma -\mathrm{2}\eta \mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\phi \mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_z+\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\theta }_z-\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\gamma -\eta \mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\phi \mathrm{c}\mathrm{o}\mathrm{s}\gamma +\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_z}{\mathrm{1}-{\eta }^{\mathrm{2}}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\phi }}\left|\mathit{W}\right|.（29）

由式（29）得到斜坡床面上的临界起动切应力大小. 对于水平床面，通过受力分析易得：

{\left|\mathit{F}\right|}_{\mathrm{0}}=\left(\left|\mathit{W}\right|-\left|{\mathit{F}}_{\mathrm{L}}\right|\right)\mathrm{t}\mathrm{a}\mathrm{n}\phi .（30）

由式（29）、式（30）可得：

\left|{\mathit{F}}_{\mathrm{D}}\right|/\left|{\mathit{F}}_{\mathrm{0}}\right|=k.（31）

式中：k为临界起动切应力斜坡修正系数：

k=\sqrt[]{\frac{\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\phi \mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\gamma +{\eta }^{\mathrm{2}}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\gamma -\mathrm{2}\eta \mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\phi \mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_z+\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\theta }_z-\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\gamma -\eta \mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\gamma +\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_z}{\mathrm{t}\mathrm{a}\mathrm{n}\phi \left(\mathrm{1}-\eta \mathrm{t}\mathrm{a}\mathrm{n}\phi \right)}}.（32）

水平面的临界起动切应力经过斜坡修正为

{\tau }_{\mathrm{b}\mathrm{c}\mathrm{r}}=k\cdot {\tau }_{\mathrm{b}\mathrm{c}\mathrm{r},\mathrm{0}}.（33）

2.3 泥沙模块运行流程

本文冲刷模型中，初始状态下泥沙粒子皆为沉积物，经由每个时间步进行处理，将依次完成沉积物向推移质转化、推移质向悬移质转化以及悬移质向推移质转化. 转化完成后，分别存储3种状态的泥沙粒子，进入SPH求解步骤，直至计算结束，计算流程示意图如图3所示.

2.3.1 预处理模块

获取上一时间步中泥沙粒子状态信息；迭代计算摩阻流速并应用Shields侵蚀准则得到临界Shields数；按Mastbergen公式计算上升与沉降流速.

2.3.2 推移质模块

沉积物粒子转化为推移质粒子须满足2.2.2节所述起动条件. 因此，需判断其是否处于沙床表面，进而利用Einstein对数流速分布公式计算该粒子所受流体切应力. 若切应力超过屈服应力，则储存为新的推移质粒子，并结合HBP模型利用Shield准则计算屈服强度，获得推移质粒子表观黏度.

2.3.3 悬移质模块

推移质粒子转化为悬移质粒子，须满足2.2.3节所述条件. 由于水流裹挟悬移质输送存在浓度阈值，首先判断粒子邻域内泥沙体积浓度，当浓度小于阈值时，判断其速度大小，若速度大于上升流速，则判断为悬移质粒子并依据Vand模型更新粒子黏度.

对于悬移质粒子，其速度若小于沉降流速，则又转化为推移质粒子，返回推移质模块.

2.3.4 沉积物模块

泥沙粒子不起动，基于DP准则结合HBP模型更新其屈服强度.

3 溃坝冲刷数值仿真

SPH法具有可以模拟流体大变形、复杂液面、不同液相间相互裹挟作用的优势. 在两相流冲刷模型中，水-沙相互作用将根据控制方程自动计算，无需单独求解输运方程. 本文数值仿真分为两步：①基于DualSPHysics开源代码进行二次开发，完成泥沙状态转化功能；②进行溃坝水流冲刷模拟，结合水槽试验结果验证数值模型.

需要注意，本节重点对提出的泥沙模块进行验证，对于DualSPHysics水体计算速度场、压力场和湍流场精确性验证，不在本文的讨论范围内. Sato等^［26］对多种场景下DualSPHysics流域计算精度进行了验证，限于篇幅，不赘述.

3.1 基于DualSPHysics的冲刷仿真计算流程

本文在DualSPHysics开源代码的基础上，实现沉积物粒子起动判断与类型转化. 二次开发模块分为三部分执行：①初始化模块；②判断模块；③黏度计算模块.程序的运行流程如图4所示.

针对不同问题与工程条件，本文二次开发模块具有高度可定制性，可根据研究目标灵活地进行土层参数、泥沙起动理论变换. 精确捕捉两相流体交界面处粒子黏性与屈服特性，需尽可能提高粒子分辨率，即缩小粒子间距、增加粒子数量. 为利用有限计算机资源提高SPH模拟效率，本文基于CUDA框架利用GPU实现泥沙模块加速计算，以解决三维模型对于计算资源要求高的问题.

3.2 数值模型参数设置

Dey^［25］在研究黎曼波理论时引用了Louvain溃坝侵蚀试验，该试验已被多名学者^{［9， 11， 27-28］}用作冲刷模型评估标准，可为本文提供验证结果. 试验中采用长度2.5 m、宽度0.1 m、高度\mathrm{0.35}\mathrm{m}的透明水槽，在底部铺满厚度5~6 cm、等效直径\mathrm{3.5}\mathrm{m}\mathrm{m}、密度\mathrm{1}\mathrm{540}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}的PVC颗粒. 试验开始时，通过高速抬升水闸释放高度为\mathrm{10}\mathrm{c}\mathrm{m}的水柱，并采用高速摄像机记录冲刷过程.

本文三维数值模型水槽宽度设置为\mathrm{0.1}\mathrm{m}，在厚度为\mathrm{6}\mathrm{c}\mathrm{m}、长度为200 cm的泥沙层上建立深度为\mathrm{10}\mathrm{c}\mathrm{m}、长度为\mathrm{100}\mathrm{c}\mathrm{m}的水体. 当模拟开始时，水体将迅速跌落，水-沙耦合作用将引起泥沙起动与输送，从而发生溃坝侵蚀.

SPH法中粒子流速与其设置间距紧密关联，参照Zubeldia等^［19］的研究，通过设置4组不同间距（0.016 m、0.008 m、0.004 m、0.002 m）模型组别与水槽试验数据分别进行对比，从而确定适宜粒子间距，模型设置如图5所示. 水槽空间布置同溃坝试验相同，长2.5 m，水箱右侧与水槽连接位置（X=0 m）设有0.02 m开口，水槽底部铺设泥沙，为避免泥沙随水流发生输运，其密度、屈服强度及黏度预设值偏大，可视为固定床面. 水粒子动力黏度取1×10^-6 m²·s，密度为\mathrm{1}\mathrm{000}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}，模型求解时长为8 s. 各组别SPH模型与试验水头在t=8 s时前进位置绘制如图6所示. 由图6可知，随着粒子间距不断减小，流体运动模拟结果逐渐逼近真实值. 据此建立数值模型，粒子间距取d_{\mathrm{p}}=\mathrm{0.002}\mathrm{m}，共产生3 087 630个粒子.

泥沙与水的密度分别设置为{\rho }_{\mathrm{s}}=\mathrm{1}\mathrm{540}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}和{\rho }_{\mathrm{w}}=\mathrm{1}\mathrm{000}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}. 泥沙相作为非牛顿流体，参照Fourtakas等^［13］的研究，HBP模型参数取m=\mathrm{100}，n=\mathrm{1.8}，动力黏度取\mathrm{0.001}{\mathrm{m}}^{\mathrm{2}}/\mathrm{s}，库伦黏聚力c=\mathrm{0}，内摩擦角\phi =\mathrm{38}°，以模拟物理试验中无黏聚力PVC材料特性；水作为牛顿流体，取m=\mathrm{0}，n=\mathrm{1}，动力黏度取\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}{\mathrm{m}}^{\mathrm{2}}/\mathrm{s}. 沉积物转化推移质质量阈值参照文献［29］取1 250 kg/m³.

3.3 溃坝冲刷形态及误差分析

为验证SPH法模拟河床冲刷等大变形及强非线性问题的精度优势，本文采用Flow-3D建立相同空间布置网格计算模型，泥沙及水体材料参数设置与SPH模型相同，顶面为压力边界，其余为壁面，采用LES湍流模型，网格尺寸为0.002 m，模型如图7所示.

本次SPH计算平台为Intel i5-12600KF+NIVIDA 3080Ti，模拟时长为1 s，初始状态模型如 图8所示，求解时间为2 h 22 min.图9、图10展示了SPH模型、网格模型与Louvain溃坝试验在0.25 s、0.50 s、0.75 s、1.00 s时剖面比较结果，其中：图9（a）为Fraccarollo物理试验剖面图；图9（b）为本文模拟结果，其中蓝色表示水体，红色表示沉积物；图9 （c）为Flow-3D网格模型模拟结果，其中蓝色为水体，灰色为沉积物.由于篇幅限制，本文仅展示0.25 s时刻试验及各模型结果.图10展示了Louvain试验结果、Fourtakas模拟结果、本文模拟结果以及网格模型模拟结果的水、沙表面高程对比图，通过计算同一时刻试验与模拟结果液面曲线竖直方向高程均方根误差（E_RMS），量化比较两者轮廓线吻合程度.需要注意，在试验与各数值模拟中，水与沉积物均置于三维水槽中，图中展示仅为水槽侧面视图；同时，由于粒子总数量庞大造成其分辨率较高，部分悬疑质粒子紧邻水-沙交界面，影响到真实河床表面高程观测，因此实际河床冲刷后高程应以图10为准.

由图9可知，在t=0.25 s时刻，本文模型和Fourtakas模型均能较好地再现Louvain试验中水头行进位置与形状；水-沙交界面处均较为准确地再现了最大冲刷坑深度与位置.从图9（b）局部放大图中可以看出，交界面处漂移泥沙粒子较多，且受到水流冲击和裹挟作用于近床处发生悬移，河床真实冲刷位置应位于该类粒子下方. 本文模型较好地模拟出沉积物泥沙受水体冲蚀发生状态转化，并在水流裹挟下推移、输送，最终沉降、堆积. 图9（c）网格模型中河床在水流侵蚀下形成冲刷坑，未再现交界面处水-沙裹挟及水头后方位置泥沙堆积行为，其冲刷坑位置更为靠后，水头模拟误差较大.图10（a）反映出本文模型与Fourtakas模型总体均方根误差相近，但本文水头行进位置与试验结果更为相符，而网格模型水面线及冲刷地形均方根误差均为SPH模型数倍以上.

在t=0.50 s时刻［图10（b）］，本文与Fourtakas模型水头与推移质行进长度均大于试验结果，而本文所得最大冲刷深度相较于Fourtakas模型更为贴近试验结果；网格模型中仍未出现泥沙堆积，仅有冲蚀作用；由高程对比可知，针对自由液面与床沙表面整体形态模拟，本文模型均方根误差均优于Fourtakas模型，网格模型模拟误差最大.

在t=0.75 s时刻［图10（c）］和t=\mathrm{1.00}\mathrm{s}时刻［图 10（d）］，本文模型自由液面均方根误差略大于Fourtakas模型，而水沙交界面显著优于Fourtakas模型，网格模型误差最大.本文模型所得最大冲刷深度略小于试验值. 总体上看，相较于SPH模型，网格模型难以精准捕捉水沙交界面处大变形、强非线性行为，水面及冲刷地形模拟误差均较大.

为进一步对SPH数值仿真结果的局部误差进行分析，选取坝内断面（S1）、最大冲深位置断面（S2）、Louvain试验水头处断面（S3）共3处断面对本文模型与Fourtakas模型模拟误差进行比较，各断面分布汇总如图11所示.分别提取各时刻两种模型所有断面的自由液面高程及泥沙界面高程，并计算局部误差率绘制成图12.由图12（a）可知，两种模型在水头处（S3）液面模拟误差相比于其他断面均偏大，这与水头推移位置模拟略远于试验结果相关. Fourtakas模型液面高程整体上略低于本文模型，此种现象是由其河床冲蚀深度整体大于试验结果造成的，河床冲蚀越深使得水位高程相应地沉降越多.本文模型河床冲蚀趋势与试验结果更为贴近，但液面整体上略高于试验结果，推测原因可能为：Louvain溃坝侵蚀试验采用PVC塑料颗粒模拟泥沙，而颗粒间存在间隙，随着试验的进行，水流会发生渗流，导致水位下降；Fourtakas模型与本文模型均采用SPH方法模拟溃坝冲刷，未能考虑水流向下渗透作用，最终造成本文模型冲刷结果相近而水位偏高的情况，此亦能较好地解释Fourtakas模型冲蚀大于试验结果而水位几乎与试验持平的现象.

图12（b）中S1、S3断面处泥沙高程模拟，本文模型明显优于Fourtakas模型，Fourtakas模型在最大冲深位置前后的模拟误差较大，本文模型改进了这一不足；而各个时刻最大冲深结果，两模型各有优劣. 图13给出了Louvain试验、Fourtakas模型以及本文模型最大冲深处高程随时间变化趋势图. 由图13可知，在0.50 s前，本文模型可以较为精准地预测最大冲刷深度. 本文模型和Fourtakas模型模拟的最大冲刷深度均出现在0.50 s时刻，而实际最大冲深出现在 0.75 s时刻；自0.50 s后，本文模型和Fourtakas模型最大冲刷深度均出现回升，推测原因为冲刷坑周围泥沙随水流回填至最大冲深处，造成高程回升，而此现象亦在原试验0.75 s之后发生. 整体上看，本文模型所得全过程最大冲深略小于试验值，而Fourtakas模型略大于试验值.另外，泥沙回填与其流动性相关，最终由其黏度决定，本文泥沙采用HBP模型，其黏度受m、n值控制，为统一变量从而进行冲刷精度验证，本文参数取值参照了Fourtakas模型；然而，在本文模型中由于引入了Shields准则作为屈服泥沙强度计算依据，原有模型参数取值可能造成黏度计算误差，从而加剧了冲刷坑泥沙回填现象.事实上，Louvain溃坝侵蚀试验对照研究中^{［9，11，13，28］}，即使采用相同泥沙本构，也未出现模型参数取值统一，基于SPH法溃坝冲刷最大深度的精确模拟，仍需在后续研究工作中对泥沙模型及参其数取值深入研究.

综上，相比于Fourtakas模型，本文建立的改进模型能够对试验河床整体冲刷形态、演变及推进过程进行更为精确地模拟；两种模型最大冲刷深度与真实结果均存在一定误差. 总体上看，使用本文提出的两相流SPH改进算法，仿真计算结果较为可信.

4 结论与展望

1）提出一种基于光滑粒子元的水-沙两相流冲刷数值仿真改进方法. 将泥沙视为一种非牛顿流体，引入HBP模型描述其非牛顿流体行为. 根据泥沙运动状态将其划分为沉积物、推移质、悬移质三种类型，引入上升与沉降流速模型作为三种不同状态泥沙转化判断依据. 应用Drucker-Prager与Shields侵蚀准则计算泥沙起动临界切应力，并为不同状态泥沙粒子赋予不同流变特性.

2）考虑水流升力与水流拖曳力，推导了在三维斜坡床面上泥沙起动临界切应力修正系数，作为水平床面泥沙侵蚀准则的补充.

3）基于所提出的SPH多相流模型河床冲刷数值仿真改进算法，采用三维溃坝冲刷算例与试验结果以及其他同类仿真结果对比验证模型精度. 结果表明，相比于其他同类模型，本文建立的改进模型能高精度模拟试验河床整体冲刷形态、演变及推进过程.

4）当前模型最大冲深与试验结果仍存在一定误差，模型中亦未能考虑水体渗透影响，后续工作将针对上述问题开展进一步研究.

参考文献

[1] 唐国汉， 冮大兴， 郑旭峰， 等 ．国道318线金沙江大桥抢通及恢复重建设计［J］．山西建筑， 2021， 47（15）： 143-145．

[2] TANG G H， GANG D X， ZHENG X F，et al ．Design of rehabilitation and reconstruction of Jinsha River Bridge on G318［J］．Shanxi Architecture，2021， 47（15）： 143-145．（in Chinese）

[3] 余文畴， 张志林 ．2002—2018年长江口基本河槽冲刷及形态调整演化趋势［J］．长江科学院院报，2021，38（8）：1-8．

[4] YU W C， ZHANG Z L ．Evolution trend of basic channel scour and morphological adjustment in Yangtze River Estuary from 2002 to 2018［J］．Journal of Yangtze River Scientific Research Institute，2021，38（8）： 1-8．（in Chinese）

[5] 熊文， 蔡春声， 张嵘钊 ．桥梁水毁研究综述［J］．中国公路学报，2021， 34（11）： 10-28．

[6] XIONG W， CAI C S， ZHANG R Z. Review of hydraulic bridge failures［J］．China Journal of Highway and Transport，2021， 34（11）： 10-28．（in Chinese）

[7] XIONG W， TANG P B， KONG B，et al ．Computational simulation of live-bed bridge scour considering suspended sediment loads［J］．Journal of Computing in Civil Engineering，2017，31（5）：4017040．

[8] 熊文，姚浩， CAI C S，等 ．考虑悬移质效应的桥墩动床冲刷精细化分析方法［J］．湖南大学学报（自然科学版），2016， 43（5）：52-60．

[9] XIONG W， YAO H， CAI C S，et al ．Bridge scour simulation in live-bed condition with suspended load［J］．Journal of Hunan University （Natural Sciences），2016，43（5）： 52-60．（in Chinese）

[10] 张曙光， 尹进步， 张根广 ．基于Flow-3D的圆柱形桥墩局部冲刷大涡模拟［J］．泥沙研究，2020，45（1）： 67-73．

[11] ZHANG S G， YIN J B， ZHANG G G ．Large-eddy simulation on local scour of cylindrical piers based on Flow-3D［J］．Journal of Sediment Research， 2020， 45（1）： 67-73．（in Chinese）

[12] 王占彬， 张卫杰， 张健，等 ．基于并行SPH方法的地震滑坡对桥桩的冲击作用［J］．湖南大学学报（自然科学版）， 2022， 49（7）： 54-65．

[13] WANG Z B， ZHANG W J， ZHANG J，et al ．Impact of earthquake-induced landslide on bridge pile based on parallelized SPH method［J］. Journal of Hunan University （Natural Sciences），2022， 49（7）： 54-65．（in Chinese）

[14] KIM J， LEE J H， JANG H，et al ．Numerical investigation of scour by incompressible SPH coupled with coarse-grained DEM［J］. Soil Dynamics and Earthquake Engineering，2021， 151： 106998．

[15] MA X J， ZHANG B W， CHEN J，et al ．The simulation of sediment transport and erosion caused by free-surface flow based on two-phase SPH model with the improved Shields criterion［J］．Ocean Dynamics，2022， 72（2）： 169-186．

[16] NG F C， ZAWAWI M H， AZMAN A，et al ．Smooth particle hydrodynamics modelling of liquid-sediment system and coastal wave breaker［J］．Ocean Dynamics，2022，72（2）：99-114．

[17] WANG D， LI S W， ARIKAWA T， et al ．ISPH simulation of scour behind seawall due to continuous tsunami overflow［J］．Coastal Engineering Journal， 2016， 58（3）： 1650014-1．

[18] MANENTI S， SIBILLA S， GALLATI M， et al ．SPH simulation of sediment flushing induced by a rapid water flow［J］．Journal of Hydraulic Engineering，2012， 138（3）： 272-284．

[19] FOURTAKAS G， ROGERS B D ．Modelling multi-phase liquid-sediment scour and resuspension induced by rapid flows using Smoothed Particle Hydrodynamics （SPH） accelerated with a Graphics Processing Unit （GPU）［J］．Advances in Water Resources，2016， 92： 186-199．

[20] SHAKIBAEINIA A， JIN Y C ．A mesh-free particle model for simulation of mobile-bed dam break［J］．Advances in Water Resources，2011， 34（6）： 794-807．

[21] WENDLAND H ．Piecewise polynomial，positive definite and compactly supported radial functions of minimal degree［J］．Advances in Computational Mathematics， 1995， 4（1）： 389-396．

[22] AWAD B N， TAIT M J ．Macroscopic modelling for screens inside a tuned liquid damper using incompressible smoothed particle hydrodynamics［J］．Ocean Engineering， 2022， 263： 112320．

[23] NIKEGHBALI P， BENJANKAR R ．Bingham-plastic and spring-dashpot models in SPH method：simulation of bed load material beneath the violent flows［C］//Geo-Congress 2022．Charlotte，North Carolina. Reston， VA： American Society of Civil Engineers， 2022： 375-384．

[24] HOSSEINI S， MANZARI M， HANNANI S ．A fully explicit three-step SPH algorithm for simulation of non-Newtonian fluid flow［J］．International Journal of Numerical Methods for Heat & Fluid Flow，2007， 17（7）： 715-735．

[25] ZUBELDIA E H， FOURTAKAS G， ROGERS B D，et al ．Multi-phase SPH model for simulation of erosion and scouring by means of the shields and Drucker-Prager criteria［J］. Advances in Water Resources，2018，117：98-114．

[26] EINSTEIN H A， EL-SAMNI E S A ．Hydrodynamic forces on a rough wall［J］. Reviews of Modern Physics， 1949， 21（3）：520-524．

[27] MASTBERGEN D R， VAN DEN BERG J H ．Breaching in fine sands and the generation of sustained turbidity currents in submarine canyons［J］．Sedimentology，2003， 50（4）： 625-637．

[28] VAND V. Viscosity of solutions and suspensions；theory［J］．The Journal of Physical and Colloid Chemistry，1948，52（2）：277-299．

[29] VAN RIJN L C. Principles of sediment transport in rivers，estuaries，and coastal seas［M］. Amsterdam：Aqua Publications，1993．

[30] CHEN X L， MA J M， DEY S ．Sediment transport on arbitrary slopes：simplified model［J］．Journal of Hydraulic Engineering，2010， 136（5）： 311-317．

[31] DEY S ．Threshold of sediment motion on combined transverse and longitudinal sloping beds［J］．Journal of Hydraulic Research，2003， 41（4）： 405-415．

[32] SATO K， KAWASAKI K， WATANABE K，et al ．Validation of the applicability of the particle-based open-source software DualSPHysics to violent flow fields［J］．Coastal Engineering Journal， 2021， 63（4）： 545-572．

[33] FRACCAROLLO L， CAPART H ．Riemann wave description of erosional dam-break flows［J］．Journal of Fluid Mechanics，2002， 461： 183-228．

[34] MEMARZADEH R， BARANI G， GHAEINI-HESSAROEYEH M ．Numerical modeling of sediment transport based on unsteady and steady flows by incompressible smoothed particle hydrodynamics method［J］．Journal of Hydrodynamics， 2018， 30（5）： 928-942．

[35] 王东. 水流作用下建筑物周围局部冲刷ISPH数值模拟及实验研究［D］. 天津： 天津大学，2017．

[36] WANG D. ISPH simulation and experimental research of local scour around structures under currents［D］. Tianjin：Tianjin University， 2017．（in Chinese）

基金资助

国家自然科学基金资助项目(52022021,51978160)