泥石流数值模拟方法研究进展

乔成 欧国强 潘华利 王钧 宇岩

摘要：泥石流作为一种多相混合介质，所包含的物理过程和动力学特征非常复杂。针对泥石流问题的数值模拟方法随着数值计算方法和物理计算模型的发展而发展，基于物理过程的数值模拟方法为探究泥石流复杂物理现象背后的机理提供了一种有效手段。回顾了求解泥石流动力问题的数值模拟方法，从连续介质计算方法、离散介质计算方法和混合介质计算方法3个方面分析了不同数值模拟方法的特点和适用情况，介绍了在泥石流分析中常用的数值模拟软件及其特点，展望了求解泥石流动力学问题的数值模拟方法的发展趋势。结果表明：传统的基于网格的计算方法已有长足发展和较长的应用历史，相对比较成熟，但是在处理大变形、快速运移的自由表面流问题时，存在网格容易畸变等问题；基于粒子的计算方法在处理上述问题时无需网格的划分和维护，易于确定自由表面位置和多相间的界面，但存在边界条件处理困难等问题；混合介质计算方法在较小尺度范围内对固体颗粒物质与液相相互作用机理进行探讨时具有重要作用。

关键词：泥石流；数值模拟；动力学模型；多相；连续介质；离散介质；混合介质；深度积分

中图分类号：P642.23文献标志码：A

0引言

泥石流是一种由土、砂、石等固体颗粒物与水组成的，在重力驱动下沿山坡或沟谷运移的混合流体，具有宽级配、高浓度、直进性、大冲大淤、冲击力大、破坏力强等特点，往往由暴雨、融雪、溃坝、滑坡等引发。成因、物源条件等因素的差异导致不同类型泥石流的动力学特征存在很大差异；沿程侵蚀、堆积的发生导致同一场泥石流在不同阶段也存在较大差异。

泥石流动力学特征受多种因素影响，如固体颗粒物体积分数、粒径组成、黏性物质含量等。根据体积分数可以将狭义泥石流分为3类，即黏性泥石流、稀性泥石流和过渡性泥石流，而广义泥石流还包括泥流和水石流。基于内部应力的特点，可将泥石流分为准静力泥石流和动力泥石流[1]。准静力泥石流是指体积分数大于50%（重度为1 830 kg·m-3）的泥石流，相当于黏性泥石流；动力泥石流是指体积分数在20%～50%之间（重度为1 330～1 830 kg·m-3）的泥石流，相当于稀性泥石流。

侵蚀、堆积的发展和变化伴随着泥石流运移的全过程，进而引起泥石流物质构成比例的变化，并进一步导致内部应力构成上的变化。不同类型的应力在泥石流内部具有一种此消彼长的关系。当颗粒碰撞应力占优时，黏性应力就会减小，如水石流，反之亦然。泥石流动力学数学模型应基于泥石流的物质构成和内部应力构成特点来建立。由于模型的复杂性，数学模型一般需要依靠各种数值方法进行求解。

泥石流动力学问题的数值模拟方法总体上可以分为连续介质、离散介质和混合介质3类。本文分析了采用数值模拟方法求解泥石流动力学问题的特殊性和难点，在此基础上进一步分析了上述3类数值方法在求解泥石流动力学问题方面的适用性及相关研究的最新进展，为系统、快速了解泥石流数值模拟研究提供参考。

1连续介质计算方法

连续介质计算方法以被研究对象的质量和变形连续分布为基本特征，建立描述动力学特征的方程组，并将这些方程基于网格或（和）粒子在计算域上进行离散化，结合初始和边界条件来求解。

1.1基于网格的计算方法

传统的基于网格的计算方法有：有限差分法（FDM）、有限元法（FEM）[2]和有限体积法（FVM）[34]。这类方法将计算域划分为很多连续分布的基本单元（网格或体元），在这些离散化的计算域上，控制方程基于某种数值方法被转化为包含节点未知场变量信息的代数方程组。控制方程有微分和积分两种形式，有限差分法求解微分形式的控制方程，在网格离散的基础上采用某种格式的代数差分代替偏微分控制方程中的偏导数，而有限元法和有限体积法则是基于离散控制体求解积分形式的动力学方程。

基于网格的计算方法在模拟泥石流问题时的主要挑战是：需要处理自由面和锋面移动、侵蚀锋面处物理量的稳定性、激波捕捉、不规则底床和沟道复杂网格的生成和维护等。在泥石流高速运移或岩土体崩滑问题中，急流汇入缓流时会有激波产生。将拉格朗日移动网格格式应用于非守恒形式的动力学方程可以较好地处理自由边界问题，但是无法处理激波。一些改进格式由此诞生，如修正的Godunov算法[5]、PatrovGalerkin法[6]、无振荡算法[7]、近似Riemann求解器[8]、全变差减小格式（TVD）[9]等。上述激波处理方法主要基于有限差分的守恒非振荡形式[10]或者有限体积法对流深守恒的属性[11]。这些算法明显改进了基于网格的计算方法，但这些算法应用于有限体积法时不能保证流深始终为正[12]，为此MangeneyCastelnau基于在粒子微观尺度引入的狄拉克分布函数提出了一种有限体积法的动力格式[12]。这种格式能有效处理不连续问题，再现颗粒崩滑运移中的流动和休止过程，并能保证流深始终为正。任意拉格朗日欧拉法（ALE）使有限差分法和有限元法能够更好地处理运动自由面的位置确定问题，同时该方法与基于固定欧拉网格的处理自由面方法（如流体体积（VOF）法[13]）相比，自由面的离散化精度更高[14]。传统的基于网格的计算方法求解泥石流运移问题的主要困难在于网格的移动和维护以及由网格引起的数值收敛问题。

Savage等研究了有限质量颗粒物质沿粗糙斜面崩滑问题，假设颗粒物质满足库仑摩擦定律，获得了描述崩滑的深度平均方程，即著名的SavageHutter模型[15]。Savage等采用欧拉和拉格朗日两种有限差分形式对模型进行求解，发现拉格朗日型差分格式能更好地处理颗粒物和空气的交界面，求解更简洁有效，结果更稳定可靠[15]。Wieland等采用有限体积法与有限差分法联合的方法求解二维SavageHutter模型[16]。Pitman等在SavageHutter模型中考虑了底床侵蚀的影响，并采用Godunov型有限体积法进行求解，提高了数值稳定性和计算精度[17]。Hungr采用拉格朗日有限差分法求解了滑坡动态模拟（DAN）模型[18]。该模型将运移物质看作等效流体，可以使用不同的流变模型并可近似考虑流变参数受底床物质加入、孔隙压力等因素的影响，近似考虑了运移过程中沟道的侧限作用。Iverson将SavageHutter模型的单相颗粒流模型拓展为二维颗粒流体两相混合模型[19]，又进一步拓展至三维情形，并可以考虑运移过程中孔隙压力的演变[20]，在计算过程中间无需调整模型初始参数便可描述泥石流从起动到堆积的全部动力过程。Denlinger将SavageHutter模型在三维不规则地形上进一步完善，显式地考虑了竖向重力加速度的影响并通过侧应力系数进一步完善了三维空间中库仑应力的表示，采用基于Roe型黎曼求解器的高精度有限体积法求解模型[21]。之后，Iverson等显式地考虑了剪胀和孔隙压力变化之间的相互影响，使模型所包含的方程数增加到5个，分别描述流深、固体体积比、底床孔隙压力和2个速度分量，每个方程包含一个体现颗粒剪胀影响的源项[22]，并采用Godunov型激波捕捉格式的有限体积法进行求解[3]。由于模型假设固液两相具有相同的速度，所以该模型实际是准两相流模型。

为探究泥石流在运移过程中高流动性的成因，He等基于热多孔弹性介质模型原理采用基于TVDMUSCL格式的有限体积法进行分析，发现摩擦生热产生的孔隙压力增加，提高了岩土体的流动性并减小了阻力，摩擦角的改变和剪胀的出现会明显影响运移动力过程和孔隙压力的发展[4]。目前，多数模型对泥石流运移中底床物质挟带的影响采取简化处理，要么不考虑挟带，要么基于经验公式来考虑。Han等基于动力学原理考虑了挟带的影响，建立了基于深度平均原理的动力学模型，并采用显式有限差分法进行求解[23]。

有限元法、有限差分法和有限体积法单独无法或很难处理多相问题，尤其是可变、复杂的相（Phase）间界面的多相问题，为此要与一些新型方法相结合（如流体体积法[13]）。流体体积法对位于相间界面上的体元，通过求解一个额外的偏微分方程来确定每个控制体内液体的填充比例，借此确定界面的位置。但这类方法在处理复杂、急变交界面方面还存在一定困难，而且存在由迁移项引起的数值扩散问题[24]。

基于网格的计算方法采用连续分布的网格对计算域进行空间离散化。对于欧拉网格方法，对复杂的计算域进行网格划分是复杂和耗时的；在固定的欧拉网格上准确定义自由表面、变形边界、移动交界面难度较大。拉格朗日网格可以采用自适应网格和网格再分技术处理交界面的快速移动和捕捉问题，但计算时步在网格细化后需要适当减小，计算耗时会明显增加。

1.2基于粒子的计算方法

基于粒子的计算方法不需要繁琐的网格划分，分布于计算域内的离散粒子并不是物理上的真实粒子而是一种对连续计算域进行离散的方式，这些粒子携带各种场变量（如质量、密度、速度等）以及其他与具体问题相关的派生变量（如涡量）。纯无网格粒子法一般是完全拉格朗日框架下的，如光滑粒子流体动力学法（SPH）和移动粒子半隐式法（MPS）。这两种方法很类似，如粒子的影响都是通过核函数定义的，但不同的是空间导数的近似是基于局部加权平均的，不涉及核函数的梯度，并且动量方程中的压力项是隐式确定的。

光滑粒子流体动力学法是最古老的纯拉格朗日框架下的无网格方法，最初用于求解开放空间天体物理学问题，因为空间中物质点的运动与气体或流体的运动很相似，可以用经典牛顿流体动力学的控制方程求解。也正是由于这种相似性，光滑粒子流体动力学法后来被拓展至计算流体力学领域，如自由表面流动问题[25]。光滑粒子流体动力学法的核心是基于核函数的插值原理，场函数及其导数通过核函数可以转化为连续积分形式，该积分形式又可近似为支持域内所有离散粒子上累加求和的离散形式。在转化为连续积分的过程中，对场函数的微分转化为对光滑核函数的微分，这可以放宽对场函数连续性的要求，与弱形式的算法很相似。与基于欧拉网格的方法相比，光滑粒子流体动力学法严格满足质量守恒，更适于处理大变形、可移动边界、自由表面追踪、移动交界面问题[26]。

基于光滑粒子流体动力学原理的计算模型可以采用不同的流变本构关系。以黏性应力为主的泥石流（通常拜格诺数小于40，且雷诺数小于500）剪胀和湍动特征都不明显，可以基于黏塑性本构模型进行求解（如宾汉模型[27]、Cross模型[28]、HerschelBulkley模型[29]以及Voellmy模型等）；而水石流（通常拜格诺数大于450，且相对流深小于10）主要以固体颗粒间的相互作用（碰撞、摩擦）为主，剪胀效应明显，可以采用剪胀体模型[30]进行分析。基于粒子的计算方法对侵蚀和挟带的影响目前主要是基于经验公式来考虑。McDougall等基于已有的DAN模型将运移物质视作符合摩擦型流变本构关系的等效流体，底床物质按摩擦型或Voellmy型流变特征进行考虑，采用基于光滑粒子流体动力学格式的深度积分模型进行分析，与室内和原位测试结果吻合较好[31]。Pastor等用一个额外的方程来考虑孔隙压力升高导致的增速和孔隙压力消散带来的滑移休止作用，底床摩擦分别按宾汉模型、摩擦流体模型、Voellmy模型和黏性摩擦模型考虑，并按Egashira、Hungr和Blanc等3种侵蚀率经验公式对侵蚀进行分析，并比较不同公式的准确性和特点[32]。

基于粒子的计算方法与基于网格的计算方法相比更适合在复杂的地形上进行计算。对于泥石流体始终与沟床接触的情况，前述方法都适用，但对于泥石流体与沟道某些时刻发生分离的情况（如过坝的泥石流问题），传统有限元法、有限差分法和有限体积法通常很难单独完成计算，而基于粒子的计算方法则能胜任这样的计算。

1.3网格与粒子混合方法

网格与粒子混合方法一般是欧拉和拉格朗日框架相结合的方法，如Particleincell（PIC）法、物质点法（MPM）[33]。该类方法中拉格朗日框架下的粒子用于捕捉变动交界面，而场变量的计算基于欧拉网格进行。倪晋仁等将流团模型引入Particleincell法求解准结构两相流模型，并用于分析阵性泥石流的起动、运移和堆积过程[34]。物质点法是计算流体力学中的Particleincell法在计算固体力学中的拓展，是一种基于有限元法的粒子法。在物质点法中，状态变量以拉格朗日方式在粒子（物质点）上进行处理，物质点在计算过程中随物理量的演进在欧拉背景网格上发生移动，在每一计算循环结束时的迁移步中，背景网格被恢复至初始位置。Mast等通过多尺度、多领域规划算法来处理固液两相，通过双网格实现一般形状的边界处理，使物质点法能分析泥石流问题[35]。

粒子有限元法（PFEM）同样是欧拉与拉格朗日方法的结合，可以考虑自由表面、破碎波、流动分离等。该方法将有限元的节点视为可以自由移动的粒子，粒子的移动基于物理控制方程，粒子移动和网格处理交替进行，每一轮粒子移动后，基于新的计算域构型重新生成网格，因此，可以有效避免大变形中的网格畸变问题。与其他混合方法（如Particleincell法和物质点法）不同的是，粒子有限元法中的粒子是非物质的点，只传递流体强度属性。Zhang等采用粒子有限元法研究流态滑坡问题，分析主要影响因素及流动的演变过程，并与离散元的计算结果进行对比，验证粒子有限元法求解这类问题的有效性[36]。

网格与粒子混合方法的主要缺点是算法较复杂且计算开销较大。在与泥石流动力学问题结合时，如何解决高密度差带来的数值稳定问题以及固体颗粒粒径巨大差异带来的多尺度问题等，都有待进一步深入研究。

2离散介质计算方法

泥石流问题分析中的离散介质计算方法主要是离散元法（DEM）[37]和格子玻尔兹曼法（LBM）。

2.1离散元法

离散元法以代表真实颗粒物质的理想颗粒体为研究对象，从微观角度出发，通过定义颗粒间的相互作用来反映大量颗粒物质宏观的动力学特征。颗粒间的相互作用在离散元法中有接触作用和非接触作用两类。前者一般从颗粒间接触处的法向和切向来定义力学特征，如线性弹簧黏壶模型、HertzMindlin模型等接触模型；后者体现的是不通过颗粒间的直接接触传递的作用，对于岩土问题主要是孔隙压力，需要在细颗粒含量较高且有孔隙水的情况下考虑。

离散元法适用于对颗粒物质的动力学机理进行研究（如颗粒分选现象[3839]和起动机理[40]），适于分析浆体黏性小的水石流和稀性泥石流。由于离散元法的计算涉及颗粒间复杂接触关系的实时确定、更新和存储等一系列大数据量的计算，因此，计算开销相比连续介质计算方法增加很多，在大尺度的泥石流问题中应用较少。另外，微观粒子间的接触本构模型参数与宏观物质整体的动力学特征参数之间的对应关系复杂。尽管这种微观宏观间的转换最近已经取得了一定进展[4142]，但仍有待进一步深入研究。

2.2格子玻尔兹曼法

格子玻尔兹曼法是一种离散粒子与欧拉网格相结合的计算方法，通过求解带碰撞项（如BhatnagarGrossKrook模型）的离散玻尔兹曼（Boltzmann）方程来模拟流体的流动，而不是求解纳维斯托克斯（NavierStokes）方程。该方法的基本研究对象是粒子团，是一种介观研究方法。它将概率密度函数作为唯一依赖变量，这里的概率定义为在某一时刻一定范围内发现速度满足指定值的粒子的概率。该方法适于处理考虑大颗粒与流体相互作用下的多相、复杂边界、自由表面流问题，易于实现并行计算。王沁等基于格子玻尔兹曼法，结合泥石流流团模型的特点建立特殊格式的格子玻尔兹曼（LB）模型[43]。该模型从平衡和非平衡态分布函数两个层次刻画了泥石流复杂的流变关系和运移特征，尤其是泥石流的堆积过程，并分析了地形对堆积形态的影响。

3混合介质计算方法

混合介质计算方法是连续介质与离散介质的耦合计算，将泥石流中的液态物质用连续介质模型考虑，用离散介质模型模拟泥石流中的固体颗粒。这种方法适合于分析间隙流作用、固体颗粒间的相互作用都比较明显的泥石流。

对包含一定固体颗粒物的流体进行模拟时，离散介质与连续介质耦合的方法是一种很自然的选择。流体中的颗粒物质通过离散元模拟，而携带颗粒物的流体通过连续介质方式模拟。通过混合的计算流体动力学离散元法（CFDDEM）模拟含颗粒流体的流动问题时，流体部分通过局部平均的纳维斯托克斯方程进行描述，之后还出现了大涡模拟（LES）与离散元法的耦合，即LESDEM，连续介质部分采用了过滤的纳维斯托克斯方程。连续介质与离散元耦合的方法还包括有限差分法[44]、有限体积法[45]、光滑粒子流体动力学[46]、格子玻尔兹曼法[47]与离散元法的耦合。

Leonardi等用离散元描述泥石流中固体颗粒的力学行为，并考虑了颗粒的实际粒径分布特征；液相通过格子玻尔兹曼法进行求解，通过一个简化模型考虑颗粒形状对颗粒旋转机制的影响，并考虑了非牛顿流变关系和自由表面，分析了复杂的颗粒物运移、粒径分选和底床侵蚀现象，有利于深入理解和评估泥石流运移机理以及泥石流与防护结构的相互作用[48]。

根据是否考虑固体颗粒对流体的影响，将耦合算法分为单向耦合算法和双向耦合算法。当实际固体颗粒物所占比例较大时，为满足精度要求必须采用双向耦合算法。耦合算法涉及的变量多且算法复杂，当模型较大时，每一时步的离散元计算结果与其他算法间交换的数据量很大，对计算机资源的消耗很大。目前，这类耦合算法还不成熟，处于初期发展阶段，主要结合试验应用于物理现象机理的探讨，较少应用于实际工程问题的分析。

4软件应用

在较大尺度的工程问题分析中，常常将数值计算模型与地理信息系统（GIS）软件结合起来应用，如高桥保模型[49]与GIS的联合应用[50]、FLO2D软件[5155]、EDDA软件[56]等。其他一些则与开源软件GRASS GIS联合应用，如DAN软件[18]、TITAN2D软件[57]、MassMov2D软件[58]以及RAMMS软件等。

FLO2D软件基于体积守恒和水力学原理建立微分形式的质量守恒和动量守恒方程，采用显式中心差分法对守恒方程进行求解。FLO2D软件可以计算流深、流速和影响范围的时间变化，计算受数字高程模型精度以及水力参数的影响明显，无法考虑激波、水跃和沟道侵蚀现象；将泥石流视为洪水或挟沙水流，因此，对泥石流的模拟结果与实际监测结果在某些情况下吻合的并不理想[5960]。Debris2D软件[61]是一款专门分析泥石流的数值软件，控制方程基于浅水假设，主要输入参数是地形图和初始物源分布，与FLO2D软件不同的是，Debris2D软件无需输入曼宁系数和降雨水位图，但需要提供底床物质的屈服应力值，物源分布的输入可以基于现场勘查、航拍或卫星图像。FLO2D软件必须由用户人为终止计算，而Debris2D软件模拟中的泥石流在满足休止条件后自动休止，适于对滑坡引发的泥石流进行分析。泥石流运移模型（DMM）[62]由香港土木工程署拓展土力工程处开发，是对Hungr的滑坡动态模拟模型[18]的改进。泥石流运移模型去除了滑坡动态模拟模型中沟道断面为矩形且侧边界无黏的限制，不需要预先定义滑移体的宽度，泥石流运移所受抗力来自沟道的全部湿周，并可以预测堆积扇的形状。泥石流运移模型和滑坡动态模拟模型共同的缺点是需要预先定义泥石流运移路径，需要在用户输入的运移路径轮廓上预先估计侵蚀和堆积的影响，这两点的完成都有较大难度。McDougall等的DAN3D软件[31]基于光滑粒子流体动力学法，可以在三维地形图上进行计算并且无需事先定义运移路径，侵蚀的影响通过模型内置的经验模型进行考虑。

尽管已有多个基于GIS的软件得到应用，但这些软件目前所采用的泥石流动力学机理相对较简单，在处理复杂地形时都进行了一定的简化。开发基于更加复杂的动力学机理[19，63]以及完全基于GIS来执行并可以考虑复杂地形影响的软件仍然面临诸多挑战。

5讨论

侵蚀和挟带对泥石流的运移动力学特征有重要影响，已有一些经验型侵蚀率计算公式被提出[6466]，结合室内试验进行了验证，并已加入到深度积分模型中

[18，31，6771]。相比经验型公式，基于动力学机理的侵蚀率计算公式[23]物理意义更加明确，计算不受经验参数准确性的影响，公式中的某些参数尽管具有物理意义，但为保证计算结果与实际情况的吻合，参数取值往往也需要在试算后进行调整。上述基于经验或理论的侵蚀率计算公式在引入计算模型时需要在质量守恒方程中加入反映质量演变的边界通量项，需要在动量方程中加入底床剪切应力和底床物质加入引起的动量通量。目前多数数值模拟方法通常只是将侵蚀和挟带引起的物理影响在描述泥石流体的动力学模型中进行了考虑，侵蚀引起的模型几何特征变化并未考虑。当侵蚀深度较大时，由此引起的计算结果与真实情况的偏差将明显影响计算准确性。

目前数值计算方法在与GIS的联合应用中对地形信息都采取了一定简化，如只考虑在运移路径方向上的地形曲率而忽略垂直运移路径方向上的曲率，应用中有要求地形的坡度变化比较缓和等限制。如何突破这些限制，更加准确地考虑地形信息，同时将更加复杂的泥石流动力学模型[22，63] 引入到GIS应用中来，进一步完善激波捕捉等数值计算问题，是今后数值模拟方法与GIS联合应用需进一步解决的问题。目前对泥石流的起动机制研究相对较少，应将考虑降雨入渗的分布式水文模型与泥石流动力学模型在GIS平台内进行整合，在考虑降雨引起孔隙压力变化情况下，研究泥石流从起动到运移直至休止的全过程动力学问题。

6结语

（1）数值模拟方法的选择既要结合具体泥石流的物质组成和动力学特征，还应考虑所研究问题的侧重点。在小尺度范围内配合模型试验对泥石流（尤其混合型）的物理机理进行分析，适于采用离散元与连续介质耦合的方法。基于连续介质的方法计算消耗易于接受，其中有限元法、有限差分法和有限体积法更加成熟。针对基于粒子的计算方法的改进一直积极进行着，如半解析墙边界对光滑粒子流体动力学法边界条件施加效果的改进。拉格朗日框架下的无网格粒子法在数值算法方面还有较大的提升空间，因此，仍将是数值计算理论研究的热点方向之一。

（2）伴随泥石流动力学模型的不断发展完善和计算机技术的进步，以及并行计算技术、GPU加速技术和云计算技术在大尺度泥石流数值模拟计算中的推广，数值模拟技术将在泥石流灾害分析中发挥越来越重要的作用。

参考文献：

References：

[1]TAKAHASHI T.Debris Flow：Mechanics，Prediction and Countermeasures[M].2nd ed.Boca Raton：CRC Press，2014.

[2]CROSTA G B，IMPOSIMATO S，RODDEMAN D.Numerical Modelling of Entrainment/Deposition in Rock and Debrisavalanches[J].Engineering Geology，2009，109（1/2）：135145.

[3]GEORGE D L，IVERSON R M.A Depthaveraged Debrisflow Model That Includes the Effects of Evolving Dilatancy：Ⅱ.Numerical Predictions and Experimental Tests[J].Proceedings of the Royal Society A：Mathematical，Physical and Engineering Sciences，2014，470：20130820.

[4]HE S M，LIU W，WANG J.Dynamic Simulation of Landslide Based on Thermoporoelastic Approach[J].Computers and Geosciences，2015，75（1）：2432.

[5]FRACCAROLLO L，CAPART H，ZECH Y.A Godunov Method for the Computation of Erosional Shallow Water Transients[J].International Journal for Numerical Methods in Fluids，2003，41（9）：951976.

[6]AUDUSSE E，BOUCHUT F，BRISTEAU M O，et al.A Fast and Stable Wellbalanced Scheme with Hydrostatic Reconstruction for Shallow Water Flows[J].SIAM Journal on Scientific Computing，2004，25（6）：20502065.

[7]SANDERS B F.Highresolution and Nonoscillatory Solution of the St.Venant Equations in Nonrectangular and Nonprismatic Channels[J].Journal of Hydraulic Research，2001，39（3）：321330.

[8]ARMANINI A，FRACCAROLLO L，ROSATTI G.Twodimensional Simulation of Debris Flows in Erodible Channels[J].Computers and Geosciences，2009，35（5）：9931006.

[9]OUYANG C J，HE S M，XU Q，et al.A MacCormackTVD Finite Difference Method to Simulate the Mass Flow in Mountainous Terrain with Variable Computational Domain[J].Computers and Geosciences，2013，52（1）：110.

[10]TAI Y C，NOELLE S，GRAY J M N T，et al.Shockcapturing and Fronttracking Methods for Granular Avalanches[J].Journal of Computational Physics，2002，175（1）：269301.

[11]LAIGLE D，COUSSOT P.Numerical Modeling of Mudflows[J].Journal of Hydraulic Engineering，1997，123（7）：617623.

[12]MANGENEYCASTELNAU A.Numerical Modeling of Avalanches Based on Saint Venant Equations Using a Kinetic Scheme[J].Journal of Geophysical Research，2003，108（B11）：25272544.

[13]HIRT C W，NICHOLS B D.Volume of Fluid（VOF） Method for the Dynamics of Free Boundaries[J].Journal of Computational Physics，1981，39（1）：201225.

[14]BRAESS H，WRIGGERS P.Arbitrary Lagrangian Eulerian Finite Element Analysis of Free Surface Flow[J].Computer Methods in Applied Mechanics and Engineering，2000，190（1/2）：95109.

[15]SAVAGE S B，HUTTER K.The Motion of a Finite Mass of Granular Material Down a Rough Incline[J].Journal of Fluid Mechanics，1989，199：177215.

[16]WIELAND M，GRAY J，HUTTER K.Channelized Freesurface Flow of Cohesionless Granular Avalanches in a Chute with Shallow Lateral Curvature[J].Journal of Fluid Mechanics，1999，392：73100.

[17]PITMAN E B，NICHITA C C，PATRA A，et al.Computing Granular Avalanches and Landslides[J].Physics of Fluids，2003，15（12）：36383646.

[18]HUNGR O.A Model for the Runout Analysis of Rapid Flow Slides，Debris Flows，and Avalanches[J].Canadian Geotechnical Journal，1995，32（4）：610623.

[19]IVERSON R M.The Physics of Debris Flows[J].Reviews of Geophysics，1997，35（3）：245296.

[20]IVERSON R M，DENLINGER R P.Flow of Variably Fluidized Granular Masses Across Threedimensional Terrain：1.Coulomb Mixture Theory[J].Journal of Geophysical Research，2001，106（B1）：537552.

[21]DENLINGER R P.Granular Avalanches Across Irregular Threedimensional Terrain：1.Theory and Computation[J].Journal of Geophysical Research，2004，109（F1）：337357.

[22]IVERSON R M，GEORGE D L.A Depthaveraged Debrisflow Model that Includes the Effects of Evolving Dilatancy：Ⅰ.Physical Basis[J].Proceedings of the Royal Society A：Mathematical，Physical and Engineering Sciences，2014，470：20130819.

[23]HAN Z，CHEN G Q，LI Y G，et al.Numerical Simulation of Debrisflow Behavior Incorporating a Dynamic Method for Estimating the Entrainment[J].Engineering Geology，2015，190：5264.

[24]SHAKIBAEINIA A，JIN Y C.A Meshfree Particle Model for Simulation of Mobilebed Dam Break[J].Advances in Water Resources，2011，34（6）：794807.

[25]MONAGHAN J J.Simulating Free Surface Flows with SPH[J].Journal of Computational Physics，1994，110（2）：399406.

[26]LIU G R，LIU M B.Smoothed Particle Hydrodynamics：A Meshfree Particle Method[M].Hackensack：World Scientific，2003.

[27]缪吉伦，张文忠，周家俞.基于SPH方法的黏性泥石流堆积形态数值模拟[J].自然灾害学报，2013，22（6）：125130.

MIAO Jilun，ZHANG Wenzhong，ZHOU Jiayu.Numerical Simulation of the Accumulation State of Viscous Debris Flow by Smooth Particle Hydrodynamics Method[J].Journal of Natural Disasters，2013，22（6）：125130.

[28]SHAO S D，LO E Y M.Incompressible SPH Method for Simulating Newtonian and NonNewtonian Flows with a Free Surface[J].Advances in Water Resources，2003，26（7）：787800.

[29]PASCULLI A，SCIARRA N，MINATTI L，et al.SPH Modeling of Fast Muddy Debris Flow：Numerical and Experimental Comparison of Certain Commonly Utilized Approaches[J].Italian Journal of Geosciences，2013，132（3）：350365.

[30]PASTOR M，BLANC T，PASTOR M J.A Depthintegrated Viscoplastic Model for Dilatant Saturated Cohesivefrictional Fluidized Mixtures：Application to Fast Catastrophic Landslides[J].Journal of NonNewtonian Fluid Mechanics，2009，158（1/2/3）：142153.

[31]MCDOUGALL S，HUNGR O.Dynamic Modelling of Entrainment in Rapid Landslides[J].Canadian Geotechnical Journal，2005，42（5）：14371448.

[32]PASTOR M，BLANC T，HADDAD B，et al.Application of a SPH Depthintegrated Model to Landslide Runout Analysis[J].Landslides，2014，11（5）：793812.

[33]SULSKY D，ZHOU S J，SCHREYER H L.Application of a Particleincell Method to Solid Mechanics[J].Computer Physics Communications，1995，87（1/2）：236252.

[34]倪晋仁，廖谦，曲轶众，等.阵性泥石流运动与堆积的欧拉─拉格朗日模型：Ⅰ.理论[J].自然灾害学报，2000，9（3）：814.

NI Jinren，LIAO Qian，QU Yizhong，et al.An EulerLagrangian Model for Noncontinuous Debris Flow：Ⅰ.Theory[J].Journal of Natural Disasters，2000，9（3）：814.

[35]MAST C，MACKENZIEHELNWEIN P，ARDUINO P，et al.Landslide and Debris Flowinduced Static and Dynamic Loads on Protective Structures[M]∥BORJA R I.Multiscale and Multiphysics Processes in Geomechanics.Heidelberg：Springer，2011：169172.

[36]ZHANG X，KRABBENHOFT K，SHENG D C，et al.Numerical Simulation of a Flowlike Landslide Using the Particle Finite Element Method[J].Computational Mechanics，2015，55（1）：167177.

[37]CUNDALL P A，STRACK O D L.A Discrete Numerical Model for Granular Assemblies[J].Geotechnique，1979，29（1）：4765.

[38]ZHOU G D，NG C W W.Numerical Investigation of Reverse Segregation in Debris Flows by DEM[J].Granular Matter，2010，12（5）：507516.

[39]周公旦，孙其诚，崔鹏.泥石流颗粒物质分选机理和效应[J].四川大学学报：工程科学版，2013，45（1）：2836.

ZHOU Gongdan，SUN Qicheng，CUI Peng.Study on the Mechanisms of Solids Segregation in Granular Debris Flows[J].Journal of Sichuan University：Engineering Science Edition，2013，45（1）：2836.

[40]ZHOU J，LI Y X，JIA M C，et al.Numerical Simulation of Failure Behavior of Granular Debris Flows Based on Flume Model Tests[J].The Scientific World Journal，2013（1）：603130.

[41]THORNTON A，WEINHART T，LUDING S，et al.Modeling of Particle Size Segregation：Calibration Using the Discrete Particle Method[J].International Journal of Modern Physics C，2011，23（8）：2934.

[42]WEINHART T，LUDING S，THORNTON A R.From Discrete Particles to Continuum Fields in Mixtures[J].AIP Conference Proceedings，2013，1542（1）：12021205.

[43]王沁，姚令侃.格子Boltzmann方法及其在泥石流堆积研究中的应用[J].灾害学，2007，22（3）：15.

WANG Qin，YAO Lingkan.Lattice Boltzmann Method and Its Application in the Study on Deposition of Debris Flow[J].Journal of Catastrophology，2007，22（3）：15.

[44]BARTZKE G，HUHN K.A Conceptual Model of Porespace Blockage in Mixed Sediments Using a New Numerical Approach，with Implications for Sediment Bed Stabilization[J].Geomarine Letters，2015，35（3）：189202.

[45]CHAREYRE B，CATALANO E，BARTHELEMY E.Porescale Modeling of Fluidparticles Interaction and Emerging Poromechanical Effects[J].International Journal for Numerical and Analytical Methods in Geomechanics，2014，38（1）：5171.

[46]REN B，JIN Z，GAO R，et al.SPHDEM Modeling of the Hydraulic Stability of 2D Blocks on a Slope[J].Journal of Waterway，Port，Coastal，and Ocean Engineering，2013，140（6）：04014022.

[47]EL SHAMY U，ABDELHAMID Y.Modeling Granular Soils Liquefaction Using Coupled Lattice Boltzmann Method and Discrete Element Method[J].Soil Dynamics and Earthquake Engineering，2014，67：119132.

[48]LEONARDI A，WITTEL F K，MENDOZA M，et al.Particlefluidstructure Interaction for Debris Flow Impact on Flexible Barriers[Z].arXiv：1409.8034，2014.

[49]TAKAHASHI T，NAKAGAWA H，HARADA T，et al.Routing Debris Flows with Particle Segregation[J].Journal of Hydraulic Engineering，1992，118（11）：14901507.

[50]CHAU K T，LO K H.Hazard Assessment of Debris Flows for Leung King Estate of Hong Kong by Incorporating GIS with Numerical Simulations[J].Natural Hazards and Earth System Sciences，2004，4（1）：103116.

[51]OBRIEN J S，JULIEN P Y，FULLERTON W T.Twodimensional Water Flood and Mudflow Simulation[J].Journal of Hydraulic Engineering，1993，119（2）：244261.

[52]张鹏，马金珠，舒和平，等.基于FLO2D模型的泥石流运动冲淤数值模拟[J].兰州大学学报：自然科学版，2014，50（3）：363368.

ZHANG Peng，MA Jinzhu，SHU Heping，et al.Numerical Simulation of Erosion and Deposition Debris Flow Based on FLO2D Model[J].Journal of Lanzhou University：Natural Sciences，2014，50（3）：363368.

[53]贾涛，唐川，王纳纳.基于FLO2D与冲量模型的泥石流危险度分区方法及应用[J].水电能源科学，2015，33（2）：152155.

JIA Tao，TANG Chuan，WANG Nana.Method and Application of Debris Flow Hazard Zoning Based on FLO2D and Impulse Model[J].Water Resources and Power，2015，33（2）：152155.

[54]黄勋，唐川，周伟.基于数值模拟的暴雨泥石流暴发频率计算模型[J].工程地质学报，2014，22（6）：12711278.

HUANG Xun，TANG Chuan，ZHOU Wei.Numerical Simulation of Occurrence Frequency Estimation Model for Debris Flows[J].Journal of Engineering Geology，2014，22（6）：12711278.

[55]王纳纳，唐川.基于FLO2D的都江堰市龙池镇黄央沟泥石流数值模拟[J].地质灾害与环境保护，2014，25（1）：107112.

WANG Nana，TANG Chuan.Numerical Simulation of Huangyang Gully Debris Flow Longchi Town Dujiangyan City Based on FLO2D[J].Journal of Geological Hazards and Environment Preservation，2014，25（1）：107112.

[56]CHEN H X，ZHANG L M.EDDA 1.0：Integrated Simulation of Debris Flow Erosion，Deposition and Property Changes[J].Geoscientific Model Development，2015，8（3）：829844.

[57]LE L，PITMAN E B.A Model of Granular Flows over an Erodible Surface[J].SIAM Journal on Applied Mathematics，2009，70（5）：14071427.

[58]BEGUERIA S，ASCH T W J V，MALET J P，et al.A GISbased Numerical Model for Simulating the Kinematics of Mud and Debris Flows over Complex Terrain[J].Natural Hazards and Earth System Sciences，2009，9（6）：18971909.

[59]STOLZ A，HUGGEL C.Debris Flows in the Swiss National Park：The Influence of Different Flow Models and Varying DEM Grid Size on Modeling Results[J].Landslides，2008，5（3）：311319.

[60]BONIELLO M A，CALLIGARIS C，LAPASIN R，et al.Rheological Investigation and Simulation of a Debrisflow Event in the Fella Watershed[J].Nature Hazards and Earth System Sciences，2010，10（5）：989997.

[61]LIU K F，HUANG M C.Numerical Simulation of Debris Flow with Application on Hazard Area Mapping[J].Computational Geosciences，2006，10（2）：221240.

[62]KWAN J S H，SUN H W.An Improved Landslide Mobility Model[J].Canadian Geotechnical Journal，2006，43（5）：531539.

[63]PUDASAINI S P.A General Twophase Debris Flow Model[J].Journal of Geophysical Research，2012，117（F3）：F03036.

[64]TAKAHASHI T.Debris Flow[M].Rotterdam：Balkema，1991.

[65]BRUFAU P，GARCIANAVARRO P，GHILARDI P，et al.1D Mathematical Modelling of Debris Flow[J].Journal of Hydraulic Research，2000，38（6）：435446.