聶一品，蘭 玲，王協(xié)康

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

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

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

1 計算模型與研究方法

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

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

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

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

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

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

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

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

1.2 模型及計算參數(shù)設(shè)置

1）模型設(shè)置

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

圖1 典型山洪泥石流溝道數(shù)值試驗?zāi)Ｐ虵ig. 1 Numerical test model of a typical landslide gully

2）計算參數(shù)設(shè)置

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

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

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

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

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

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

2 計算結(jié)果與分析

2.1 泥沙堆積形態(tài)的演變

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

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