喬成 歐國強 潘華利 王鈞 宇巖

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

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

中圖分類號：P642.23文獻標志碼：A

0引言

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

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

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

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

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

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

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

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

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

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

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

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

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

1.2基于粒子的計算方法

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

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

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

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

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

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

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

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

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

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

2.1離散元法

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

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

2.2格子玻爾茲曼法

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

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

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

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

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

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

4軟件應(yīng)用

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

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

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

5討論

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

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

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

6結(jié)語

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

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

參考文獻：

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].自然災害學報，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]周公旦，孫其誠，崔鵬.泥石流顆粒物質(zhì)分選機理和效應(yīng)[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方法及其在泥石流堆積研究中的應(yīng)用[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，F(xiàn)ULLERTON W T.Twodimensional Water Flood and Mudflow Simulation[J].Journal of Hydraulic Engineering，1993，119（2）：244261.

[52]張鵬，馬金珠，舒和平，等.基于FLO2D模型的泥石流運動沖淤數(shù)值模擬[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與沖量模型的泥石流危險度分區(qū)方法及應(yīng)用[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]黃勛，唐川，周偉.基于數(shù)值模擬的暴雨泥石流暴發(fā)頻率計算模型[J].工程地質(zhì)學報，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ì)災害與環(huán)境保護，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.