喬成 歐國(guó)強(qiáng) 潘華利 王鈞 宇巖

摘要：泥石流作為一種多相混合介質(zhì)，所包含的物理過(guò)程和動(dòng)力學(xué)特征非常復(fù)雜。針對(duì)泥石流問(wèn)題的數(shù)值模擬方法隨著數(shù)值計(jì)算方法和物理計(jì)算模型的發(fā)展而發(fā)展，基于物理過(guò)程的數(shù)值模擬方法為探究泥石流復(fù)雜物理現(xiàn)象背后的機(jī)理提供了一種有效手段?；仡櫫饲蠼饽嗍鲃?dòng)力問(wèn)題的數(shù)值模擬方法，從連續(xù)介質(zhì)計(jì)算方法、離散介質(zhì)計(jì)算方法和混合介質(zhì)計(jì)算方法3個(gè)方面分析了不同數(shù)值模擬方法的特點(diǎn)和適用情況，介紹了在泥石流分析中常用的數(shù)值模擬軟件及其特點(diǎn)，展望了求解泥石流動(dòng)力學(xué)問(wèn)題的數(shù)值模擬方法的發(fā)展趨勢(shì)。結(jié)果表明：傳統(tǒng)的基于網(wǎng)格的計(jì)算方法已有長(zhǎng)足發(fā)展和較長(zhǎng)的應(yīng)用歷史，相對(duì)比較成熟，但是在處理大變形、快速運(yùn)移的自由表面流問(wèn)題時(shí)，存在網(wǎng)格容易畸變等問(wèn)題；基于粒子的計(jì)算方法在處理上述問(wèn)題時(shí)無(wú)需網(wǎng)格的劃分和維護(hù)，易于確定自由表面位置和多相間的界面，但存在邊界條件處理困難等問(wèn)題；混合介質(zhì)計(jì)算方法在較小尺度范圍內(nèi)對(duì)固體顆粒物質(zhì)與液相相互作用機(jī)理進(jìn)行探討時(shí)具有重要作用。

關(guān)鍵詞：泥石流；數(shù)值模擬；動(dòng)力學(xué)模型；多相；連續(xù)介質(zhì)；離散介質(zhì)；混合介質(zhì)；深度積分

中圖分類(lèi)號(hào)：P642.23文獻(xiàn)標(biāo)志碼：A

0引言

泥石流是一種由土、砂、石等固體顆粒物與水組成的，在重力驅(qū)動(dòng)下沿山坡或溝谷運(yùn)移的混合流體，具有寬級(jí)配、高濃度、直進(jìn)性、大沖大淤、沖擊力大、破壞力強(qiáng)等特點(diǎn)，往往由暴雨、融雪、潰壩、滑坡等引發(fā)。成因、物源條件等因素的差異導(dǎo)致不同類(lèi)型泥石流的動(dòng)力學(xué)特征存在很大差異；沿程侵蝕、堆積的發(fā)生導(dǎo)致同一場(chǎng)泥石流在不同階段也存在較大差異。

泥石流動(dòng)力學(xué)特征受多種因素影響，如固體顆粒物體積分?jǐn)?shù)、粒徑組成、黏性物質(zhì)含量等。根據(jù)體積分?jǐn)?shù)可以將狹義泥石流分為3類(lèi)，即黏性泥石流、稀性泥石流和過(guò)渡性泥石流，而廣義泥石流還包括泥流和水石流?；趦?nèi)部應(yīng)力的特點(diǎn)，可將泥石流分為準(zhǔn)靜力泥石流和動(dòng)力泥石流[1]。準(zhǔn)靜力泥石流是指體積分?jǐn)?shù)大于50%（重度為1 830 kg·m-3）的泥石流，相當(dāng)于黏性泥石流；動(dòng)力泥石流是指體積分?jǐn)?shù)在20%～50%之間（重度為1 330～1 830 kg·m-3）的泥石流，相當(dāng)于稀性泥石流。

侵蝕、堆積的發(fā)展和變化伴隨著泥石流運(yùn)移的全過(guò)程，進(jìn)而引起泥石流物質(zhì)構(gòu)成比例的變化，并進(jìn)一步導(dǎo)致內(nèi)部應(yīng)力構(gòu)成上的變化。不同類(lèi)型的應(yīng)力在泥石流內(nèi)部具有一種此消彼長(zhǎng)的關(guān)系。當(dāng)顆粒碰撞應(yīng)力占優(yōu)時(shí)，黏性應(yīng)力就會(huì)減小，如水石流，反之亦然。泥石流動(dòng)力學(xué)數(shù)學(xué)模型應(yīng)基于泥石流的物質(zhì)構(gòu)成和內(nèi)部應(yīng)力構(gòu)成特點(diǎn)來(lái)建立。由于模型的復(fù)雜性，數(shù)學(xué)模型一般需要依靠各種數(shù)值方法進(jìn)行求解。

泥石流動(dòng)力學(xué)問(wèn)題的數(shù)值模擬方法總體上可以分為連續(xù)介質(zhì)、離散介質(zhì)和混合介質(zhì)3類(lèi)。本文分析了采用數(shù)值模擬方法求解泥石流動(dòng)力學(xué)問(wèn)題的特殊性和難點(diǎn)，在此基礎(chǔ)上進(jìn)一步分析了上述3類(lèi)數(shù)值方法在求解泥石流動(dòng)力學(xué)問(wèn)題方面的適用性及相關(guān)研究的最新進(jìn)展，為系統(tǒng)、快速了解泥石流數(shù)值模擬研究提供參考。

1連續(xù)介質(zhì)計(jì)算方法

連續(xù)介質(zhì)計(jì)算方法以被研究對(duì)象的質(zhì)量和變形連續(xù)分布為基本特征，建立描述動(dòng)力學(xué)特征的方程組，并將這些方程基于網(wǎng)格或（和）粒子在計(jì)算域上進(jìn)行離散化，結(jié)合初始和邊界條件來(lái)求解。

1.1基于網(wǎng)格的計(jì)算方法

傳統(tǒng)的基于網(wǎng)格的計(jì)算方法有：有限差分法（FDM）、有限元法（FEM）[2]和有限體積法（FVM）[34]。這類(lèi)方法將計(jì)算域劃分為很多連續(xù)分布的基本單元（網(wǎng)格或體元），在這些離散化的計(jì)算域上，控制方程基于某種數(shù)值方法被轉(zhuǎn)化為包含節(jié)點(diǎn)未知場(chǎng)變量信息的代數(shù)方程組?？刂品匠逃形⒎趾头e分兩種形式，有限差分法求解微分形式的控制方程，在網(wǎng)格離散的基礎(chǔ)上采用某種格式的代數(shù)差分代替偏微分控制方程中的偏導(dǎo)數(shù)，而有限元法和有限體積法則是基于離散控制體求解積分形式的動(dòng)力學(xué)方程。

基于網(wǎng)格的計(jì)算方法在模擬泥石流問(wèn)題時(shí)的主要挑戰(zhàn)是：需要處理自由面和鋒面移動(dòng)、侵蝕鋒面處物理量的穩(wěn)定性、激波捕捉、不規(guī)則底床和溝道復(fù)雜網(wǎng)格的生成和維護(hù)等。在泥石流高速運(yùn)移或巖土體崩滑問(wèn)題中，急流匯入緩流時(shí)會(huì)有激波產(chǎn)生。將拉格朗日移動(dòng)網(wǎng)格格式應(yīng)用于非守恒形式的動(dòng)力學(xué)方程可以較好地處理自由邊界問(wèn)題，但是無(wú)法處理激波。一些改進(jìn)格式由此誕生，如修正的Godunov算法[5]、PatrovGalerkin法[6]、無(wú)振蕩算法[7]、近似Riemann求解器[8]、全變差減小格式（TVD）[9]等。上述激波處理方法主要基于有限差分的守恒非振蕩形式[10]或者有限體積法對(duì)流深守恒的屬性[11]。這些算法明顯改進(jìn)了基于網(wǎng)格的計(jì)算方法，但這些算法應(yīng)用于有限體積法時(shí)不能保證流深始終為正[12]，為此MangeneyCastelnau基于在粒子微觀尺度引入的狄拉克分布函數(shù)提出了一種有限體積法的動(dòng)力格式[12]。這種格式能有效處理不連續(xù)問(wèn)題，再現(xiàn)顆粒崩滑運(yùn)移中的流動(dòng)和休止過(guò)程，并能保證流深始終為正。任意拉格朗日歐拉法（ALE）使有限差分法和有限元法能夠更好地處理運(yùn)動(dòng)自由面的位置確定問(wèn)題，同時(shí)該方法與基于固定歐拉網(wǎng)格的處理自由面方法（如流體體積（VOF）法[13]）相比，自由面的離散化精度更高[14]。傳統(tǒng)的基于網(wǎng)格的計(jì)算方法求解泥石流運(yùn)移問(wèn)題的主要困難在于網(wǎng)格的移動(dòng)和維護(hù)以及由網(wǎng)格引起的數(shù)值收斂問(wèn)題。

Savage等研究了有限質(zhì)量顆粒物質(zhì)沿粗糙斜面崩滑問(wèn)題，假設(shè)顆粒物質(zhì)滿足庫(kù)侖摩擦定律，獲得了描述崩滑的深度平均方程，即著名的SavageHutter模型[15]。Savage等采用歐拉和拉格朗日兩種有限差分形式對(duì)模型進(jìn)行求解，發(fā)現(xiàn)拉格朗日型差分格式能更好地處理顆粒物和空氣的交界面，求解更簡(jiǎn)潔有效，結(jié)果更穩(wěn)定可靠[15]。Wieland等采用有限體積法與有限差分法聯(lián)合的方法求解二維SavageHutter模型[16]。Pitman等在SavageHutter模型中考慮了底床侵蝕的影響，并采用Godunov型有限體積法進(jìn)行求解，提高了數(shù)值穩(wěn)定性和計(jì)算精度[17]。Hungr采用拉格朗日有限差分法求解了滑坡動(dòng)態(tài)模擬（DAN）模型[18]。該模型將運(yùn)移物質(zhì)看作等效流體，可以使用不同的流變模型并可近似考慮流變參數(shù)受底床物質(zhì)加入、孔隙壓力等因素的影響，近似考慮了運(yùn)移過(guò)程中溝道的側(cè)限作用。Iverson將SavageHutter模型的單相顆粒流模型拓展為二維顆粒流體兩相混合模型[19]，又進(jìn)一步拓展至三維情形，并可以考慮運(yùn)移過(guò)程中孔隙壓力的演變[20]，在計(jì)算過(guò)程中間無(wú)需調(diào)整模型初始參數(shù)便可描述泥石流從起動(dòng)到堆積的全部動(dòng)力過(guò)程。Denlinger將SavageHutter模型在三維不規(guī)則地形上進(jìn)一步完善，顯式地考慮了豎向重力加速度的影響并通過(guò)側(cè)應(yīng)力系數(shù)進(jìn)一步完善了三維空間中庫(kù)侖應(yīng)力的表示，采用基于Roe型黎曼求解器的高精度有限體積法求解模型[21]。之后，Iverson等顯式地考慮了剪脹和孔隙壓力變化之間的相互影響，使模型所包含的方程數(shù)增加到5個(gè)，分別描述流深、固體體積比、底床孔隙壓力和2個(gè)速度分量，每個(gè)方程包含一個(gè)體現(xiàn)顆粒剪脹影響的源項(xiàng)[22]，并采用Godunov型激波捕捉格式的有限體積法進(jìn)行求解[3]。由于模型假設(shè)固液兩相具有相同的速度，所以該模型實(shí)際是準(zhǔn)兩相流模型。

為探究泥石流在運(yùn)移過(guò)程中高流動(dòng)性的成因，He等基于熱多孔彈性介質(zhì)模型原理采用基于TVDMUSCL格式的有限體積法進(jìn)行分析，發(fā)現(xiàn)摩擦生熱產(chǎn)生的孔隙壓力增加，提高了巖土體的流動(dòng)性并減小了阻力，摩擦角的改變和剪脹的出現(xiàn)會(huì)明顯影響運(yùn)移動(dòng)力過(guò)程和孔隙壓力的發(fā)展[4]。目前，多數(shù)模型對(duì)泥石流運(yùn)移中底床物質(zhì)挾帶的影響采取簡(jiǎn)化處理，要么不考慮挾帶，要么基于經(jīng)驗(yàn)公式來(lái)考慮。Han等基于動(dòng)力學(xué)原理考慮了挾帶的影響，建立了基于深度平均原理的動(dòng)力學(xué)模型，并采用顯式有限差分法進(jìn)行求解[23]。

有限元法、有限差分法和有限體積法單獨(dú)無(wú)法或很難處理多相問(wèn)題，尤其是可變、復(fù)雜的相（Phase）間界面的多相問(wèn)題，為此要與一些新型方法相結(jié)合（如流體體積法[13]）。流體體積法對(duì)位于相間界面上的體元，通過(guò)求解一個(gè)額外的偏微分方程來(lái)確定每個(gè)控制體內(nèi)液體的填充比例，借此確定界面的位置。但這類(lèi)方法在處理復(fù)雜、急變交界面方面還存在一定困難，而且存在由遷移項(xiàng)引起的數(shù)值擴(kuò)散問(wèn)題[24]。

基于網(wǎng)格的計(jì)算方法采用連續(xù)分布的網(wǎng)格對(duì)計(jì)算域進(jìn)行空間離散化。對(duì)于歐拉網(wǎng)格方法，對(duì)復(fù)雜的計(jì)算域進(jìn)行網(wǎng)格劃分是復(fù)雜和耗時(shí)的；在固定的歐拉網(wǎng)格上準(zhǔn)確定義自由表面、變形邊界、移動(dòng)交界面難度較大。拉格朗日網(wǎng)格可以采用自適應(yīng)網(wǎng)格和網(wǎng)格再分技術(shù)處理交界面的快速移動(dòng)和捕捉問(wèn)題，但計(jì)算時(shí)步在網(wǎng)格細(xì)化后需要適當(dāng)減小，計(jì)算耗時(shí)會(huì)明顯增加。

1.2基于粒子的計(jì)算方法

基于粒子的計(jì)算方法不需要繁瑣的網(wǎng)格劃分，分布于計(jì)算域內(nèi)的離散粒子并不是物理上的真實(shí)粒子而是一種對(duì)連續(xù)計(jì)算域進(jìn)行離散的方式，這些粒子攜帶各種場(chǎng)變量（如質(zhì)量、密度、速度等）以及其他與具體問(wèn)題相關(guān)的派生變量（如渦量）。純無(wú)網(wǎng)格粒子法一般是完全拉格朗日框架下的，如光滑粒子流體動(dòng)力學(xué)法（SPH）和移動(dòng)粒子半隱式法（MPS）。這兩種方法很類(lèi)似，如粒子的影響都是通過(guò)核函數(shù)定義的，但不同的是空間導(dǎo)數(shù)的近似是基于局部加權(quán)平均的，不涉及核函數(shù)的梯度，并且動(dòng)量方程中的壓力項(xiàng)是隱式確定的。

光滑粒子流體動(dòng)力學(xué)法是最古老的純拉格朗日框架下的無(wú)網(wǎng)格方法，最初用于求解開(kāi)放空間天體物理學(xué)問(wèn)題，因?yàn)榭臻g中物質(zhì)點(diǎn)的運(yùn)動(dòng)與氣體或流體的運(yùn)動(dòng)很相似，可以用經(jīng)典牛頓流體動(dòng)力學(xué)的控制方程求解。也正是由于這種相似性，光滑粒子流體動(dòng)力學(xué)法后來(lái)被拓展至計(jì)算流體力學(xué)領(lǐng)域，如自由表面流動(dòng)問(wèn)題[25]。光滑粒子流體動(dòng)力學(xué)法的核心是基于核函數(shù)的插值原理，場(chǎng)函數(shù)及其導(dǎo)數(shù)通過(guò)核函數(shù)可以轉(zhuǎn)化為連續(xù)積分形式，該積分形式又可近似為支持域內(nèi)所有離散粒子上累加求和的離散形式。在轉(zhuǎn)化為連續(xù)積分的過(guò)程中，對(duì)場(chǎng)函數(shù)的微分轉(zhuǎn)化為對(duì)光滑核函數(shù)的微分，這可以放寬對(duì)場(chǎng)函數(shù)連續(xù)性的要求，與弱形式的算法很相似。與基于歐拉網(wǎng)格的方法相比，光滑粒子流體動(dòng)力學(xué)法嚴(yán)格滿足質(zhì)量守恒，更適于處理大變形、可移動(dòng)邊界、自由表面追蹤、移動(dòng)交界面問(wèn)題[26]。

基于光滑粒子流體動(dòng)力學(xué)原理的計(jì)算模型可以采用不同的流變本構(gòu)關(guān)系。以黏性應(yīng)力為主的泥石流（通常拜格諾數(shù)小于40，且雷諾數(shù)小于500）剪脹和湍動(dòng)特征都不明顯，可以基于黏塑性本構(gòu)模型進(jìn)行求解（如賓漢模型[27]、Cross模型[28]、HerschelBulkley模型[29]以及Voellmy模型等）；而水石流（通常拜格諾數(shù)大于450，且相對(duì)流深小于10）主要以固體顆粒間的相互作用（碰撞、摩擦）為主，剪脹效應(yīng)明顯，可以采用剪脹體模型[30]進(jìn)行分析?；诹Ｗ拥挠?jì)算方法對(duì)侵蝕和挾帶的影響目前主要是基于經(jīng)驗(yàn)公式來(lái)考慮。McDougall等基于已有的DAN模型將運(yùn)移物質(zhì)視作符合摩擦型流變本構(gòu)關(guān)系的等效流體，底床物質(zhì)按摩擦型或Voellmy型流變特征進(jìn)行考慮，采用基于光滑粒子流體動(dòng)力學(xué)格式的深度積分模型進(jìn)行分析，與室內(nèi)和原位測(cè)試結(jié)果吻合較好[31]。Pastor等用一個(gè)額外的方程來(lái)考慮孔隙壓力升高導(dǎo)致的增速和孔隙壓力消散帶來(lái)的滑移休止作用，底床摩擦分別按賓漢模型、摩擦流體模型、Voellmy模型和黏性摩擦模型考慮，并按Egashira、Hungr和Blanc等3種侵蝕率經(jīng)驗(yàn)公式對(duì)侵蝕進(jìn)行分析，并比較不同公式的準(zhǔn)確性和特點(diǎn)[32]。

基于粒子的計(jì)算方法與基于網(wǎng)格的計(jì)算方法相比更適合在復(fù)雜的地形上進(jìn)行計(jì)算。對(duì)于泥石流體始終與溝床接觸的情況，前述方法都適用，但對(duì)于泥石流體與溝道某些時(shí)刻發(fā)生分離的情況（如過(guò)壩的泥石流問(wèn)題），傳統(tǒng)有限元法、有限差分法和有限體積法通常很難單獨(dú)完成計(jì)算，而基于粒子的計(jì)算方法則能勝任這樣的計(jì)算。

1.3網(wǎng)格與粒子混合方法

網(wǎng)格與粒子混合方法一般是歐拉和拉格朗日框架相結(jié)合的方法，如Particleincell（PIC）法、物質(zhì)點(diǎn)法（MPM）[33]。該類(lèi)方法中拉格朗日框架下的粒子用于捕捉變動(dòng)交界面，而場(chǎng)變量的計(jì)算基于歐拉網(wǎng)格進(jìn)行。倪晉仁等將流團(tuán)模型引入Particleincell法求解準(zhǔn)結(jié)構(gòu)兩相流模型，并用于分析陣性泥石流的起動(dòng)、運(yùn)移和堆積過(guò)程[34]。物質(zhì)點(diǎn)法是計(jì)算流體力學(xué)中的Particleincell法在計(jì)算固體力學(xué)中的拓展，是一種基于有限元法的粒子法。在物質(zhì)點(diǎn)法中，狀態(tài)變量以拉格朗日方式在粒子（物質(zhì)點(diǎn)）上進(jìn)行處理，物質(zhì)點(diǎn)在計(jì)算過(guò)程中隨物理量的演進(jìn)在歐拉背景網(wǎng)格上發(fā)生移動(dòng)，在每一計(jì)算循環(huán)結(jié)束時(shí)的遷移步中，背景網(wǎng)格被恢復(fù)至初始位置。Mast等通過(guò)多尺度、多領(lǐng)域規(guī)劃算法來(lái)處理固液兩相，通過(guò)雙網(wǎng)格實(shí)現(xiàn)一般形狀的邊界處理，使物質(zhì)點(diǎn)法能分析泥石流問(wèn)題[35]。

粒子有限元法（PFEM）同樣是歐拉與拉格朗日方法的結(jié)合，可以考慮自由表面、破碎波、流動(dòng)分離等。該方法將有限元的節(jié)點(diǎn)視為可以自由移動(dòng)的粒子，粒子的移動(dòng)基于物理控制方程，粒子移動(dòng)和網(wǎng)格處理交替進(jìn)行，每一輪粒子移動(dòng)后，基于新的計(jì)算域構(gòu)型重新生成網(wǎng)格，因此，可以有效避免大變形中的網(wǎng)格畸變問(wèn)題。與其他混合方法（如Particleincell法和物質(zhì)點(diǎn)法）不同的是，粒子有限元法中的粒子是非物質(zhì)的點(diǎn)，只傳遞流體強(qiáng)度屬性。Zhang等采用粒子有限元法研究流態(tài)滑坡問(wèn)題，分析主要影響因素及流動(dòng)的演變過(guò)程，并與離散元的計(jì)算結(jié)果進(jìn)行對(duì)比，驗(yàn)證粒子有限元法求解這類(lèi)問(wèn)題的有效性[36]。

網(wǎng)格與粒子混合方法的主要缺點(diǎn)是算法較復(fù)雜且計(jì)算開(kāi)銷(xiāo)較大。在與泥石流動(dòng)力學(xué)問(wèn)題結(jié)合時(shí)，如何解決高密度差帶來(lái)的數(shù)值穩(wěn)定問(wèn)題以及固體顆粒粒徑巨大差異帶來(lái)的多尺度問(wèn)題等，都有待進(jìn)一步深入研究。

2離散介質(zhì)計(jì)算方法

泥石流問(wèn)題分析中的離散介質(zhì)計(jì)算方法主要是離散元法（DEM）[37]和格子玻爾茲曼法（LBM）。

2.1離散元法

離散元法以代表真實(shí)顆粒物質(zhì)的理想顆粒體為研究對(duì)象，從微觀角度出發(fā)，通過(guò)定義顆粒間的相互作用來(lái)反映大量顆粒物質(zhì)宏觀的動(dòng)力學(xué)特征。顆粒間的相互作用在離散元法中有接觸作用和非接觸作用兩類(lèi)。前者一般從顆粒間接觸處的法向和切向來(lái)定義力學(xué)特征，如線性彈簧黏壺模型、HertzMindlin模型等接觸模型；后者體現(xiàn)的是不通過(guò)顆粒間的直接接觸傳遞的作用，對(duì)于巖土問(wèn)題主要是孔隙壓力，需要在細(xì)顆粒含量較高且有孔隙水的情況下考慮。

離散元法適用于對(duì)顆粒物質(zhì)的動(dòng)力學(xué)機(jī)理進(jìn)行研究（如顆粒分選現(xiàn)象[3839]和起動(dòng)機(jī)理[40]），適于分析漿體黏性小的水石流和稀性泥石流。由于離散元法的計(jì)算涉及顆粒間復(fù)雜接觸關(guān)系的實(shí)時(shí)確定、更新和存儲(chǔ)等一系列大數(shù)據(jù)量的計(jì)算，因此，計(jì)算開(kāi)銷(xiāo)相比連續(xù)介質(zhì)計(jì)算方法增加很多，在大尺度的泥石流問(wèn)題中應(yīng)用較少。另外，微觀粒子間的接觸本構(gòu)模型參數(shù)與宏觀物質(zhì)整體的動(dòng)力學(xué)特征參數(shù)之間的對(duì)應(yīng)關(guān)系復(fù)雜。盡管這種微觀宏觀間的轉(zhuǎn)換最近已經(jīng)取得了一定進(jìn)展[4142]，但仍有待進(jìn)一步深入研究。

2.2格子玻爾茲曼法

格子玻爾茲曼法是一種離散粒子與歐拉網(wǎng)格相結(jié)合的計(jì)算方法，通過(guò)求解帶碰撞項(xiàng)（如BhatnagarGrossKrook模型）的離散玻爾茲曼（Boltzmann）方程來(lái)模擬流體的流動(dòng)，而不是求解納維斯托克斯（NavierStokes）方程。該方法的基本研究對(duì)象是粒子團(tuán)，是一種介觀研究方法。它將概率密度函數(shù)作為唯一依賴變量，這里的概率定義為在某一時(shí)刻一定范圍內(nèi)發(fā)現(xiàn)速度滿足指定值的粒子的概率。該方法適于處理考慮大顆粒與流體相互作用下的多相、復(fù)雜邊界、自由表面流問(wèn)題，易于實(shí)現(xiàn)并行計(jì)算。王沁等基于格子玻爾茲曼法，結(jié)合泥石流流團(tuán)模型的特點(diǎn)建立特殊格式的格子玻爾茲曼（LB）模型[43]。該模型從平衡和非平衡態(tài)分布函數(shù)兩個(gè)層次刻畫(huà)了泥石流復(fù)雜的流變關(guān)系和運(yùn)移特征，尤其是泥石流的堆積過(guò)程，并分析了地形對(duì)堆積形態(tài)的影響。

3混合介質(zhì)計(jì)算方法

混合介質(zhì)計(jì)算方法是連續(xù)介質(zhì)與離散介質(zhì)的耦合計(jì)算，將泥石流中的液態(tài)物質(zhì)用連續(xù)介質(zhì)模型考慮，用離散介質(zhì)模型模擬泥石流中的固體顆粒。這種方法適合于分析間隙流作用、固體顆粒間的相互作用都比較明顯的泥石流。

對(duì)包含一定固體顆粒物的流體進(jìn)行模擬時(shí)，離散介質(zhì)與連續(xù)介質(zhì)耦合的方法是一種很自然的選擇。流體中的顆粒物質(zhì)通過(guò)離散元模擬，而攜帶顆粒物的流體通過(guò)連續(xù)介質(zhì)方式模擬。通過(guò)混合的計(jì)算流體動(dòng)力學(xué)離散元法（CFDDEM）模擬含顆粒流體的流動(dòng)問(wèn)題時(shí)，流體部分通過(guò)局部平均的納維斯托克斯方程進(jìn)行描述，之后還出現(xiàn)了大渦模擬（LES）與離散元法的耦合，即LESDEM，連續(xù)介質(zhì)部分采用了過(guò)濾的納維斯托克斯方程。連續(xù)介質(zhì)與離散元耦合的方法還包括有限差分法[44]、有限體積法[45]、光滑粒子流體動(dòng)力學(xué)[46]、格子玻爾茲曼法[47]與離散元法的耦合。

Leonardi等用離散元描述泥石流中固體顆粒的力學(xué)行為，并考慮了顆粒的實(shí)際粒徑分布特征；液相通過(guò)格子玻爾茲曼法進(jìn)行求解，通過(guò)一個(gè)簡(jiǎn)化模型考慮顆粒形狀對(duì)顆粒旋轉(zhuǎn)機(jī)制的影響，并考慮了非牛頓流變關(guān)系和自由表面，分析了復(fù)雜的顆粒物運(yùn)移、粒徑分選和底床侵蝕現(xiàn)象，有利于深入理解和評(píng)估泥石流運(yùn)移機(jī)理以及泥石流與防護(hù)結(jié)構(gòu)的相互作用[48]。

根據(jù)是否考慮固體顆粒對(duì)流體的影響，將耦合算法分為單向耦合算法和雙向耦合算法。當(dāng)實(shí)際固體顆粒物所占比例較大時(shí)，為滿足精度要求必須采用雙向耦合算法。耦合算法涉及的變量多且算法復(fù)雜，當(dāng)模型較大時(shí)，每一時(shí)步的離散元計(jì)算結(jié)果與其他算法間交換的數(shù)據(jù)量很大，對(duì)計(jì)算機(jī)資源的消耗很大。目前，這類(lèi)耦合算法還不成熟，處于初期發(fā)展階段，主要結(jié)合試驗(yàn)應(yīng)用于物理現(xiàn)象機(jī)理的探討，較少應(yīng)用于實(shí)際工程問(wèn)題的分析。

4軟件應(yīng)用

在較大尺度的工程問(wèn)題分析中，常常將數(shù)值計(jì)算模型與地理信息系統(tǒng)（GIS）軟件結(jié)合起來(lái)應(yīng)用，如高橋保模型[49]與GIS的聯(lián)合應(yīng)用[50]、FLO2D軟件[5155]、EDDA軟件[56]等。其他一些則與開(kāi)源軟件GRASS GIS聯(lián)合應(yīng)用，如DAN軟件[18]、TITAN2D軟件[57]、MassMov2D軟件[58]以及RAMMS軟件等。

FLO2D軟件基于體積守恒和水力學(xué)原理建立微分形式的質(zhì)量守恒和動(dòng)量守恒方程，采用顯式中心差分法對(duì)守恒方程進(jìn)行求解。FLO2D軟件可以計(jì)算流深、流速和影響范圍的時(shí)間變化，計(jì)算受數(shù)字高程模型精度以及水力參數(shù)的影響明顯，無(wú)法考慮激波、水躍和溝道侵蝕現(xiàn)象；將泥石流視為洪水或挾沙水流，因此，對(duì)泥石流的模擬結(jié)果與實(shí)際監(jiān)測(cè)結(jié)果在某些情況下吻合的并不理想[5960]。Debris2D軟件[61]是一款專門(mén)分析泥石流的數(shù)值軟件，控制方程基于淺水假設(shè)，主要輸入?yún)?shù)是地形圖和初始物源分布，與FLO2D軟件不同的是，Debris2D軟件無(wú)需輸入曼寧系數(shù)和降雨水位圖，但需要提供底床物質(zhì)的屈服應(yīng)力值，物源分布的輸入可以基于現(xiàn)場(chǎng)勘查、航拍或衛(wèi)星圖像。FLO2D軟件必須由用戶人為終止計(jì)算，而Debris2D軟件模擬中的泥石流在滿足休止條件后自動(dòng)休止，適于對(duì)滑坡引發(fā)的泥石流進(jìn)行分析。泥石流運(yùn)移模型（DMM）[62]由香港土木工程署拓展土力工程處開(kāi)發(fā)，是對(duì)Hungr的滑坡動(dòng)態(tài)模擬模型[18]的改進(jìn)。泥石流運(yùn)移模型去除了滑坡動(dòng)態(tài)模擬模型中溝道斷面為矩形且側(cè)邊界無(wú)黏的限制，不需要預(yù)先定義滑移體的寬度，泥石流運(yùn)移所受抗力來(lái)自溝道的全部濕周，并可以預(yù)測(cè)堆積扇的形狀。泥石流運(yùn)移模型和滑坡動(dòng)態(tài)模擬模型共同的缺點(diǎn)是需要預(yù)先定義泥石流運(yùn)移路徑，需要在用戶輸入的運(yùn)移路徑輪廓上預(yù)先估計(jì)侵蝕和堆積的影響，這兩點(diǎn)的完成都有較大難度。McDougall等的DAN3D軟件[31]基于光滑粒子流體動(dòng)力學(xué)法，可以在三維地形圖上進(jìn)行計(jì)算并且無(wú)需事先定義運(yùn)移路徑，侵蝕的影響通過(guò)模型內(nèi)置的經(jīng)驗(yàn)?zāi)Ｐ瓦M(jìn)行考慮。

盡管已有多個(gè)基于GIS的軟件得到應(yīng)用，但這些軟件目前所采用的泥石流動(dòng)力學(xué)機(jī)理相對(duì)較簡(jiǎn)單，在處理復(fù)雜地形時(shí)都進(jìn)行了一定的簡(jiǎn)化。開(kāi)發(fā)基于更加復(fù)雜的動(dòng)力學(xué)機(jī)理[19，63]以及完全基于GIS來(lái)執(zhí)行并可以考慮復(fù)雜地形影響的軟件仍然面臨諸多挑戰(zhàn)。

5討論

侵蝕和挾帶對(duì)泥石流的運(yùn)移動(dòng)力學(xué)特征有重要影響，已有一些經(jīng)驗(yàn)型侵蝕率計(jì)算公式被提出[6466]，結(jié)合室內(nèi)試驗(yàn)進(jìn)行了驗(yàn)證，并已加入到深度積分模型中

[18，31，6771]。相比經(jīng)驗(yàn)型公式，基于動(dòng)力學(xué)機(jī)理的侵蝕率計(jì)算公式[23]物理意義更加明確，計(jì)算不受經(jīng)驗(yàn)參數(shù)準(zhǔn)確性的影響，公式中的某些參數(shù)盡管具有物理意義，但為保證計(jì)算結(jié)果與實(shí)際情況的吻合，參數(shù)取值往往也需要在試算后進(jìn)行調(diào)整。上述基于經(jīng)驗(yàn)或理論的侵蝕率計(jì)算公式在引入計(jì)算模型時(shí)需要在質(zhì)量守恒方程中加入反映質(zhì)量演變的邊界通量項(xiàng)，需要在動(dòng)量方程中加入底床剪切應(yīng)力和底床物質(zhì)加入引起的動(dòng)量通量。目前多數(shù)數(shù)值模擬方法通常只是將侵蝕和挾帶引起的物理影響在描述泥石流體的動(dòng)力學(xué)模型中進(jìn)行了考慮，侵蝕引起的模型幾何特征變化并未考慮。當(dāng)侵蝕深度較大時(shí)，由此引起的計(jì)算結(jié)果與真實(shí)情況的偏差將明顯影響計(jì)算準(zhǔn)確性。

目前數(shù)值計(jì)算方法在與GIS的聯(lián)合應(yīng)用中對(duì)地形信息都采取了一定簡(jiǎn)化，如只考慮在運(yùn)移路徑方向上的地形曲率而忽略垂直運(yùn)移路徑方向上的曲率，應(yīng)用中有要求地形的坡度變化比較緩和等限制。如何突破這些限制，更加準(zhǔn)確地考慮地形信息，同時(shí)將更加復(fù)雜的泥石流動(dòng)力學(xué)模型[22，63] 引入到GIS應(yīng)用中來(lái)，進(jìn)一步完善激波捕捉等數(shù)值計(jì)算問(wèn)題，是今后數(shù)值模擬方法與GIS聯(lián)合應(yīng)用需進(jìn)一步解決的問(wèn)題。目前對(duì)泥石流的起動(dòng)機(jī)制研究相對(duì)較少，應(yīng)將考慮降雨入滲的分布式水文模型與泥石流動(dòng)力學(xué)模型在GIS平臺(tái)內(nèi)進(jìn)行整合，在考慮降雨引起孔隙壓力變化情況下，研究泥石流從起動(dòng)到運(yùn)移直至休止的全過(guò)程動(dòng)力學(xué)問(wèn)題。

6結(jié)語(yǔ)

（1）數(shù)值模擬方法的選擇既要結(jié)合具體泥石流的物質(zhì)組成和動(dòng)力學(xué)特征，還應(yīng)考慮所研究問(wèn)題的側(cè)重點(diǎn)。在小尺度范圍內(nèi)配合模型試驗(yàn)對(duì)泥石流（尤其混合型）的物理機(jī)理進(jìn)行分析，適于采用離散元與連續(xù)介質(zhì)耦合的方法?；谶B續(xù)介質(zhì)的方法計(jì)算消耗易于接受，其中有限元法、有限差分法和有限體積法更加成熟。針對(duì)基于粒子的計(jì)算方法的改進(jìn)一直積極進(jìn)行著，如半解析墻邊界對(duì)光滑粒子流體動(dòng)力學(xué)法邊界條件施加效果的改進(jìn)。拉格朗日框架下的無(wú)網(wǎng)格粒子法在數(shù)值算法方面還有較大的提升空間，因此，仍將是數(shù)值計(jì)算理論研究的熱點(diǎn)方向之一。

（2）伴隨泥石流動(dòng)力學(xué)模型的不斷發(fā)展完善和計(jì)算機(jī)技術(shù)的進(jìn)步，以及并行計(jì)算技術(shù)、GPU加速技術(shù)和云計(jì)算技術(shù)在大尺度泥石流數(shù)值模擬計(jì)算中的推廣，數(shù)值模擬技術(shù)將在泥石流災(zāi)害分析中發(fā)揮越來(lái)越重要的作用。

參考文獻(xiàn)：

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，F(xiàn)RACCAROLLO 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方法的黏性泥石流堆積形態(tài)數(shù)值模擬[J].自然災(zāi)害學(xué)報(bào)，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]倪晉仁，廖謙，曲軼眾，等.陣性泥石流運(yùn)動(dòng)與堆積的歐拉─拉格朗日模型：Ⅰ.理論[J].自然災(zāi)害學(xué)報(bào)，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]周公旦，孫其誠(chéng)，崔鵬.泥石流顆粒物質(zhì)分選機(jī)理和效應(yīng)[J].四川大學(xué)學(xué)報(bào)：工程科學(xué)版，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方法及其在泥石流堆積研究中的應(yīng)用[J].災(zāi)害學(xué)，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，F(xiàn)ULLERTON W T.Twodimensional Water Flood and Mudflow Simulation[J].Journal of Hydraulic Engineering，1993，119（2）：244261.

[52]張鵬，馬金珠，舒和平，等.基于FLO2D模型的泥石流運(yùn)動(dòng)沖淤數(shù)值模擬[J].蘭州大學(xué)學(xué)報(bào)：自然科學(xué)版，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與沖量模型的泥石流危險(xiǎn)度分區(qū)方法及應(yīng)用[J].水電能源科學(xué)，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]黃勛，唐川，周偉.基于數(shù)值模擬的暴雨泥石流暴發(fā)頻率計(jì)算模型[J].工程地質(zhì)學(xué)報(bào)，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的都江堰市龍池鎮(zhèn)黃央溝泥石流數(shù)值模擬[J].地質(zhì)災(zāi)害與環(huán)境保護(hù)，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.