高山峽谷區山洪溝口泥沙堆積特征數值試驗

聶一品，蘭 玲，王協康

(四川大學山區河流保護與治理全國重點實驗室，四川成都 610065)

中國有近一半的地區被劃分為山洪災害防治區[1]。山洪泥石流作為一種典型的山區地質現象，不僅有著廣泛的分布[2]，且頻繁成災。在以西藏東部為典型代表的高山峽谷地區，常存在高差大、海拔高且狹窄的山洪泥石流溝道[3]。在地質活動和地區氣候環境的影響下，高原冰川發育，局地性暴雨強盛，松散堆積物厚度大[4]，山洪泥石流活動頻繁。受國民經濟發展和可居住條件影響，高山峽谷溝口附近仍有大量的人類活動。山洪泥石流沖出溝道發生堆積會嚴重威脅地區人民生命財產和重大工程安全[5]。因此，眾多學者利用現場調查、物理實驗以及數值模擬的手段對山洪泥石流堆積體的物質組成、發展過程和致災機理進行研究。Cenderelli等[6]通過對歷史山洪泥石流事件的調查，發現堆積體的前緣和表面的泥沙顆粒的粒徑通常是堆積體中最大的。針對2019年汶川“8·20”山洪泥石流事件，Li等[7]通過野外調查并結合影像分析，指出不當的人類活動加劇了山洪泥石流損失。Zheng等[8]基于系列泥沙堆積實驗，表明了山洪泥石流釋放流量的增加會顯著地增大堆積體厚度和長度比值。王協康等[9]以2010年“8·13”文家溝山洪泥石流試驗表明，受主流擠壓作用，堆積體易在溝口形成多元堆積結構。Chen等[10]基于有限元特征的分裂算法建立了山洪泥石流堆積與干流交匯的模型，模擬了文家溝山洪泥石流的堆積過程。Han等[11]引入HBP流體模型，采用SPH方法模擬了2010年日本Yohutagawa山洪泥石流運動，獲得了與現場調查一致的堆積體延伸范圍。

綜上，山洪泥石流堆積數值模擬大都是將多相流動視為單相流動[12]，進而應用庫倫塑性理論或黏塑性理論對流體流動進行描述[13]，但忽略了不同大小顆粒流動過程的相互作用及其運動屬性變化[14]，難以揭示顆粒堆積過程致災機理。離散單元法（DEM）作為一種能準確描述顆粒系統中顆粒之間的復雜的相互作用以及顆粒與周圍環境之間的相互作用的數值模擬方法[15]，能直接得到顆粒在運動過程中的動態信息[16]。計算流體力學方法（CFD）因其具有計算效率高、適用性好等優點，在地球表面的各種流體的計算中發揮著重要的作用[17]。為此，耦合CFD–DEM描述流體–粒子系統運動，已在研究顆粒系統的演進過程中成為了現實[18]，例如滑坡壩破壞[19]以及山區河道泥沙輸移[20]等。Li等[21]采用耦合CFD–DEM方法，研究了山洪泥石流對柔性防護網的影響，指出防護效率與山洪泥石流漿體弗勞德數密切相關。Fang等[22]結合CFD–DEM方法，評估了山洪泥石流對剛性防護結構的沖擊特征，發現山洪泥石流固體體積分數的變化決定了防護結構上靜荷載的分布。然而，耦合CFD–DEM方法對山洪泥石流的研究大都沒有考慮其漿體流變特性的變化對顆粒運動的影響。為此，基于耦合CFD–DEM方法開展數值模擬，擬探究山洪泥石流體積濃度對不同粒徑泥沙顆粒堆積過程中平均堆積距離和分散程度的差異，揭示了顆粒堆積分散程度隨時間的變化規律，為理解山洪泥石流堆積致災機理提供科學依據。

1 計算模型與研究方法

1.1 CFD–DEM耦合模型基本原理

采用CFDEM[23]耦合引擎耦合CFD開源軟件OpenFOAM和DEM開源軟件LIGGGHTS對高山峽谷區山洪溝口泥沙堆積形態及粒度特征進行計算。計算過程中使用的控制方程如下：

1）顆粒運動控制方程。顆粒運動主要考慮平動和轉動。運動過程顆?？赡芘c鄰近的顆粒和邊壁發生碰撞或者與周圍的流體產生相互作用進而實現動量和能量交換。因此，單個顆粒的運動可以通過牛頓第二定律進行描述，其控制方程可以寫成：

式（1）～（2）中，mi為顆粒i的質量，v i為顆粒i的平動速度，Ii為顆粒的轉動慣量， ωi為顆粒i的角速度，為其它顆?；蜻叡谧饔糜陬w粒i上的接觸力，為顆粒i通過顆粒–流體相互作用受到的作用力，為顆粒i的重力，M i,j為其它顆?；蜻叡谧饔糜陬w粒i上的力矩。

2）流體運動控制方程。黏性不可壓縮流體運動通過連續性方程和動量方程進行描述，其形式如下：

式中， αf為流體的體積分數，uf為流體的運動速度，ρf為流體的密度，p為流體的壓強， τf為流體的應力張量，g為重力加速度，ffluid為流體平均下的顆粒通過顆粒–流體相互作用受到的阻力。

3）顆粒–流體間相互作用力控制方程。顆粒–流體相互作用力是由顆粒和流體之間相互運動而產生的黏性阻力，其大小與流體的流動狀態和顆粒雷諾數有關。表達式為：

式中：Cd為阻力系數； χ為經驗修正因子，χ=3.7-0.65exp(-(1.5-lgRep)2/2)，其中，Rep為顆粒雷諾數。

1.2 模型及計算參數設置

1）模型設置

采用的概化高山峽谷區山洪泥石流溝道數值模型如圖1所示。通過控制水箱中流體的高度，使不同體積濃度山洪泥石流在擋板下的出口處具有相同的壓強。將粒徑為1至2 mm的泥沙顆粒視為泥石流漿體的組成部分，使用單相流體模擬泥石流漿體的運動并結合流體體積分數（VOF）方法捕捉泥石流與空氣的界面。擋板處于x= 0的平面上。流通區底部鋪有厚度為0.1 m的卵礫石層，其粒徑為6至24 mm。流通區和加速區比降為10%，長度分別為3m和1m。堆積區采用1%的比降，大小為6m×6m。模擬的總時間為20 s。

圖1 典型山洪泥石流溝道數值試驗模型Fig. 1 Numerical test model of a typical landslide gully

2）計算參數設置

已有研究表明，山洪泥石流漿體的流變特性可以使用牛頓流體、賓漢流體、冪律流體以及二項式模式進行表述[24]，而賓漢流體因其形式簡單且具有較高的精度已經在山洪泥石流3維模擬中有著廣泛的應用[11]。此外，已有大量研究對賓漢流體模型的流變參數進行了研究[12]，故在本數值試驗中使用賓漢流體模型對山洪泥石流的流變行為進行表達。試驗中采用4種典型的山洪泥石流體積濃度（Cv）分別為15%、30%、45%和60%。在屈服應力的計算中，當體積濃度較低時采用文獻[24]中的Kang和Zhang公式（0.15

表1 山洪泥石流漿體模擬參數Tab.1 Simulation parameters of debris flow

耦合模型中顆粒計算參數見表2。其中，顆粒之間的計算參數與顆粒和壁面之間的計算參數保持一致。

表2 耦合模型中顆粒計算參數Tab.2 Particle calculation parameters

使用泥漿潰壩實驗來檢驗耦合模型對兩相流體的交界面捕捉能力以驗證模型的適用性。圖2（a）繪制了Komatina[29]實驗中第0.6 s時的泥漿表面和模擬結果。由圖2（a）可見，模型對兩相流體交界面的捕捉是較為準確的。隨后計算顆粒在賓漢流體中沉降時阻力系數–雷諾數的關系來驗證模型的準確性。經過模擬的阻力系數–雷諾數關系如圖2（b）所示，其中理論結果選用Cheng[30]提出的公式。由圖2（b）可見，模型能準確計算顆粒運動過程中顆粒–流體相互作用力。因此，耦合CFD–DEM模型可以準確地計算賓漢流體中顆粒的運動情況。

圖2 CFD–DEM耦合模型驗證Fig. 2 CFD–DEM coupling model validation results

2 計算結果與分析

2.1 泥沙堆積形態的演變

計算開始后，床面泥沙顆粒在山洪泥石流的沖刷下發生運動。當泥沙顆粒運動到堆積區后，由于山洪泥石流速度降低，泥沙顆粒開始發生沉積。泥沙顆粒在堆積區中不同時刻的堆積的形態如圖3所示，圖3中，體積分數表明該時刻堆積區上山洪泥石流的運動痕跡。從圖3（a）～（c）可以看出：體積濃度較低時（Cv≤45%），山洪泥石流沖刷的顆粒在初始堆積階段（t=5.0 s），由于堆積平面足夠大且沒有阻擋物，泥沙呈現扇形堆積；隨著山洪泥石流將更多的顆粒攜帶到堆積區中（t=7.5 s），堆積體堆積范圍不斷擴大，但依然呈現扇形分布。當Cv≤45%時，山洪泥石流對流通區床面的沖刷隨著體積濃度的增加而增加[31]，因此，在同一時刻下泥沙顆粒的堆積范圍隨著Cv的增加不斷增加。在t=10 s時，由于Cv=30%和45%的模型中的床面已經被侵蝕完畢，因此在堆積區入口處，泥沙不斷受到后續山洪泥石流的推動而產生空隙。在圖3（d）中，溝床顆粒在Cv=60%山洪泥石流沖刷下形成的堆積體形態在堆積區與流通區交界處有一段很短的梯形堆積區域，隨后擴展成半圓狀。此時，由于山洪泥石流體積濃度很高，在運動過程中需要更多地克服黏性阻力，因此，堆積體發展的速度有所減緩。

圖3 泥沙顆粒堆積過程模擬結果Fig. 3 Simulation results of sediment particle accumulation process

2.2 顆粒平均堆積距離的演化

為了探究顆粒的粒徑對其堆積行為的影響，分別計算了4種體積濃度下不同粒徑顆粒與堆積體幾何中心之間的距離（ ?r）隨時間變化情況，計算結果如圖4所示。從圖4中可以看出，各粒徑泥沙顆粒 ?r隨時間變化可分為4個階段，分別為高速增加階段、初次減速階段、增速恢復階段和穩定發展階段。

第1階段為高速增加階段，該階段中山洪泥石流從水箱中流出后在床面產生沿程沖刷，在流通區中產生了“龍頭”，大量泥沙顆粒在聚集在山洪泥石流的前緣。當山洪泥石流到達堆積區時，“龍頭”中的顆粒以很高的速度進入，導致 ?r的高速增加。這一階段中，山洪泥石流的沖刷使粒徑較大泥沙顆粒運動速度顯著快于較小的泥沙顆粒，因而當其進入堆積區后，大粒徑泥沙的 ?r增速大于粒徑較小的泥沙。第2個階段為初次減速階段。在山洪泥石流的持續沖刷下，流通區床面顆粒不斷受到侵蝕。但緊隨“龍頭”的后續漿體中泥沙含量相對較低，導致了 ?r的增速減緩。

當流通區床面泥沙被山洪泥石流完全侵蝕后，處于堆積區中的泥沙顆粒在后續漿體的推動下產生整體性運動， ?r隨時間的變化進入到第3階段即增速恢復階段。這一階段中，更大粒徑泥沙顆粒的 ?r增速依然大于較小的泥沙。此時，山洪泥石流漿體受到外圍顆粒的阻擋，后續的漿體尚未從堆積體中流出。但隨著模擬的進行，在某一方向上顆粒的堆積寬度減小，導致漿體從堆積體中流出，對堆積體的整體推動作用減小，形成了穩定發展的第4個階段。Cv=15%時，穩定發展階段的堆積體中顆粒的 ?r最終不變，但隨著Cv的增加， ?r不再保持不變，其增速顯著增大。在該階段中， ?r的增速隨著粒徑的減小逐漸增加，表明了在溝床的顆粒完全進入到堆積平面中后，由于山洪泥石流的持續作用，各個粒徑的顆粒逐漸地混合在一起。Cv=60%的山洪泥石流的流動速度較低，到模擬結束時，溝床泥沙顆粒沒有被完全侵蝕，因此，模擬中 ?r隨時間的變化僅有高速增加階段和初次減速階段。當Cv≤45%，隨著Cv的增加，相同粒徑的泥沙顆粒在同一時刻下的 ?r逐漸增加，堆積體對堆積區的侵占更加嚴重。

圖4 不同粒徑顆粒距堆積體幾何中心的平均距離Fig. 4 Average distance between particles of different sizes and geometric center of accumulation

模擬計算了大量（共計30577個）泥沙顆粒的運動數據，因此從統計的角度分析堆積區內顆粒運動的內在規律是可行的。將每一組模擬中顆粒的 ?r按照繪制箱線圖的方法進行統計。如圖5所示，dQ1與dQ3分別代表經過統計得到的 ?r的下和上四分位數。dIQR為四分位距，等于dQ1與dQ3的差值，反映了統計意義上 ?r的集中趨勢，可以表征顆粒在堆積體中分布的分散程度。上邊界使用dQ3+1.5×dIQR表示，代表了非正態分布的 ?r異常值的判別標準，用于表示堆積體外邊界范圍。經過統計分析，模擬過程中任何時刻異常 ?r的數量都不多于8.5%，因此使用統計的方法得到的dQ3與dIQR和dQ3+1.5×dIQR來分別判斷山洪泥石流堆積體中粒度特征和堆積范圍是可以接受的。

圖5 箱線圖組成部分Fig. 5 Components of box plot

圖6為堆積體中泥沙顆粒分散程度（dIQR）與泥沙顆粒的上四分位數（dQ3）的關系隨粒徑變化趨勢，其中，dIQR/dQ3是將模擬過程中統計得到的dIQR與dQ3的數值相除后進行平均。

圖6 顆粒堆積體中d IQR與d Q3之間的關系Fig.6 Relationship between dIQR and dQ3

由圖6可知，當山洪泥石流體積濃度較低（Cv≤45%）時，dIQR/dQ3的值會隨泥沙顆粒粒徑的增加而不斷的減小。這一現象說明，隨著粒徑的增大，泥沙顆粒堆積更加集中。此外，不同的dIQR/dQ3代表了各粒徑顆粒在堆積平面上所處位置的不同。由圖3（a）～（c）可知：在Cv≤45%的山洪泥石流沖刷形成的泥沙堆積體中，較大的泥沙顆粒顯著地位于堆積體的外側且各粒徑顆粒在堆積體中的分選良好；Cv=60%時，堆積體中不同粒徑顆粒的dIQR/dQ3不隨粒徑發生顯著的變化，大致保持在0.45附近。這個現象表明，隨著山洪泥石流體積濃度增加，各種粒徑的顆粒較為均勻的混雜在堆積體當中。結合圖3（d）可發現，即使當t=10 s時，堆積體的中心仍然可以明顯看見大粒徑泥沙的存在，泥沙顆粒之間的分選性較差。除粒徑為6 mm的顆粒外，當Cv≤45%時，dIQR/dQ3隨Cv的增加逐漸減小。這個現象表明隨著山洪泥石流體積濃度的增加，泥沙顆粒的分選性逐漸變差。綜合以上分析，堆積體中dIQR/dQ3的比值可以較為準確地代表各組分顆粒在堆積體內的相對位置，同時不同顆粒dIQR/dQ3數值的變化也表明了堆積體中顆粒之間具有良好的分選性。

2.3 泥沙顆粒分散程度隨時間的變化

圖7為Cv=15%山洪泥石流沖刷下，粒徑為6mm的顆粒在堆積過程中dIQR隨時間的變化過程。由圖7可知，隨著時間的增加，dIQR的增速逐漸減小，且有接近一個最大值的趨勢。由于被沖刷的泥沙顆粒需要一定的時間才能進入堆積區之中，故假設dIQR隨時間的增長滿足形如a×tb+c的冪函數關系。根據模擬數據點的分布情況，a和b都小于0，而c大于0。圖7的預測帶代表了dIQR在95%置信度上可能的范圍，且幾乎所有的模擬數據點都在95%預測帶以內，這代表了冪函數關系對dIQR隨時間的變化過程較為精準。結合圖6的分析結果，可以在任意時刻通過dIQR的擬合結果，推斷出dQ3的數值，從而估算山洪泥石流堆積體影響范圍時間的變化過程。

圖7 泥沙顆粒分散程度（d IQR）隨時間變化過程Fig.7 Variation in the dispersion of sediment particles(d IQR)over time

所有模擬泥沙顆粒分散程度（dIQR）隨時間變化過程的擬合參數與體積濃度（Cv）和粒徑（d）之間的關系如圖8所示。

圖8 擬合指標i fit與碰撞指標i collision的關系Fig.8 Relationship between the ifit and icollision

圖8中箭頭的方向指示了不同模擬中同一粒徑顆粒的擬合指標（ifit，即a×b/c）變化趨勢。Cv的增加代表了山洪泥石流漿體內顆粒碰撞的幾率的增加[32]，而離散的顆粒與漿體中顆粒碰撞幾率會隨著顆粒粒徑d增大而增大，故碰撞指標（icollision，即Cv×d）能代表離散的顆粒通過顆粒碰撞獲得能量的幾率大小。而a×b/c的增加則代表了顆粒的dIQR極值的減小和其增速的增大，也即代表著顆粒達到極限分散程度所需時間也較小。

粒徑越大的泥沙顆?？梢酝ㄟ^碰撞能夠獲得更多的能量，其進入到堆積區后能在更短時間內克服較多摩擦阻力，實現更快的分散。隨著Cv的逐漸減小，粒徑較小的顆粒與漿體發生能量交換的幾率大大減小，因此ifit隨d的增益也越大。而當顆粒完成一定程度的分散之后，堆積體的存在導致了漿體的運動方向發生偏轉，處于堆積體前緣的大顆粒難以被漿體推動，導致其分散狀態難以繼續改變，達到了極限分散狀態。對相同粒徑的泥沙顆粒來說，隨著Cv的增加，漿體在運動過程中因其黏性消耗的能量不斷增加，抵消了一部分傳遞給離散顆粒的能量，抑制了顆粒的分散。因此，如圖8所示，同一Cv下，隨著d的增大，ifit也隨之增加，且Cv越小，a×b/c隨d的增益也越大。而在同一粒徑下，a×b/c會隨著Cv的增加而呈現出先增加后減小的趨勢。

3 結 論

基于高山峽谷區山洪泥石流溝道的特點，借助耦合CFD–DEM方法對溝口泥沙顆粒的運動進行了數值計算，重點研究了堆積體整體的形態及不同粒徑顆粒的分布特征隨時間的變化，主要結論如下：

1）隨著山洪泥石流體積濃度的增加，泥沙堆積體的發展速度呈現先增加后減小的趨勢，當體積濃度較高時堆積體的形態也隨之發生改變。

2）不同粒徑泥沙顆粒與堆積體幾何中心之間的平均堆積距離（ ?r）隨著模擬時間的變化呈現出4個不同的階段。粒徑的增大會導致堆積距離的增速增大，但隨著運動的發展，各粒徑顆粒堆積距離之間的差距會逐漸減小。

3）通過統計得到顆粒堆積距離的四分位距（dIQR）和上邊界（dQ3+1.5×dIQR）可以表示堆積體中顆粒的分散程度以及堆積范圍。不同粒徑顆粒在堆積體中dIQR/dQ3的變化則反映了堆積體中分選性的良好與否。dIQR隨時間變化的趨勢的符合冪函數關系，通過擬合結果與dIQR/dQ3的數值變化規律可以得到任意時刻的泥石流堆積體的堆積范圍。

4）相同體積濃度下，隨著顆粒粒徑的增加，顆粒的極限分散程度減小但能更快地分散。在顆粒粒徑相同時，隨著體積濃度的增加，顆粒的極限分散程度先減小后增加而分散速度先增加后減小。